Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNARuddIonisationExtendedModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// Modified by Z. Francis, S. Incerti to handle HZE
28// && inverse rudd function sampling 26-10-2010
29
30#include "G4DNARuddIonisationExtendedModel.hh"
31#include "G4PhysicalConstants.hh"
32#include "G4SystemOfUnits.hh"
33#include "G4UAtomicDeexcitation.hh"
34#include "G4LossTableManager.hh"
35#include "G4DNAChemistryManager.hh"
36#include "G4DNAMolecularMaterial.hh"
37
38#include "G4IonTable.hh"
39#include "G4DNARuddAngle.hh"
40#include "G4DeltaAngle.hh"
41#include "G4Exp.hh"
42
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45using namespace std;
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
50 const G4String& nam)
51:G4VEmModel(nam),isInitialised(false)
52{
53 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
54 fpWaterDensity = 0;
55
56 slaterEffectiveCharge[0]=0.;
57 slaterEffectiveCharge[1]=0.;
58 slaterEffectiveCharge[2]=0.;
59 sCoefficient[0]=0.;
60 sCoefficient[1]=0.;
61 sCoefficient[2]=0.;
62
63 lowEnergyLimitForA[1] = 0 * eV;
64 lowEnergyLimitForA[2] = 0 * eV;
65 lowEnergyLimitForA[3] = 0 * eV;
66 lowEnergyLimitOfModelForA[1] = 100 * eV;
67 lowEnergyLimitOfModelForA[4] = 1 * keV;
68 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
69 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
70 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
71 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
72
73 verboseLevel= 0;
74 // Verbosity scale:
75 // 0 = nothing
76 // 1 = warning for energy non-conservation
77 // 2 = details of energy budget
78 // 3 = calculation of cross sections, file openings, sampling of atoms
79 // 4 = entering in methods
80
81 if( verboseLevel>0 )
82 {
83 G4cout << "Rudd ionisation model is constructed " << G4endl;
84 }
85
86 // Define default angular generator
87 SetAngularDistribution(new G4DNARuddAngle());
88
89 // Mark this model as "applicable" for atomic deexcitation
90 SetDeexcitationFlag(true);
91 fAtomDeexcitation = 0;
92 fParticleChangeForGamma = 0;
93
94 // Selection of stationary mode
95
96 statCode = false;
97}
98
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
102{
103 // Cross section
104
105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
107 {
108 G4DNACrossSectionDataSet* table = pos->second;
109 delete table;
110 }
111
112 // The following removal is forbidden G4VEnergyLossModel takes care of deletion
113 // however coverity will signal this as an error
114 // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
115
116}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119
120void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
121 const G4DataVector& /*cuts*/)
122{
123 if (verboseLevel > 3)
124 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
125
126 // Energy limits
127
128 G4String fileProton("dna/sigma_ionisation_p_rudd");
129 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
130 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
131 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
132 G4String fileHelium("dna/sigma_ionisation_he_rudd");
133 G4String fileLithium("dna/sigma_ionisation_li_rudd");
134 G4String fileBeryllium("dna/sigma_ionisation_be_rudd");
135 G4String fileBoron("dna/sigma_ionisation_b_rudd");
136 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
137 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
138 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
139 G4String fileSilicon("dna/sigma_ionisation_si_rudd");
140 G4String fileIron("dna/sigma_ionisation_fe_rudd");
141
142 G4DNAGenericIonsManager *instance;
143 instance = G4DNAGenericIonsManager::Instance();
144 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
145 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
146 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
147 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
148 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
149
150 //G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
151 //G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
152 //G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
153 //G4ParticleDefinition* siliconDef = instance->GetIon("silicon");
154 //G4ParticleDefinition* ironDef = instance->GetIon("iron");
155 G4ParticleDefinition* lithiumDef = G4IonTable::GetIonTable()->GetIon(3,7);
156 G4ParticleDefinition* berylliumDef = G4IonTable::GetIonTable()->GetIon(4,9);
157 G4ParticleDefinition* boronDef = G4IonTable::GetIonTable()->GetIon(5,11);
158 G4ParticleDefinition* carbonDef = G4IonTable::GetIonTable()->GetIon(6,12);
159 G4ParticleDefinition* nitrogenDef = G4IonTable::GetIonTable()->GetIon(7,14);
160 G4ParticleDefinition* oxygenDef = G4IonTable::GetIonTable()->GetIon(8,16);
161 G4ParticleDefinition* siliconDef = G4IonTable::GetIonTable()->GetIon(14,28);
162 G4ParticleDefinition* ironDef = G4IonTable::GetIonTable()->GetIon(26,56);
163 //
164
165 G4String proton;
166 G4String hydrogen;
167 G4String alphaPlusPlus;
168 G4String alphaPlus;
169 G4String helium;
170 G4String lithium;
171 G4String beryllium;
172 G4String boron;
173 G4String carbon;
174 G4String nitrogen;
175 G4String oxygen;
176 G4String silicon;
177 G4String iron;
178
179 G4double scaleFactor = 1 * m*m;
180
181 // LIMITS AND DATA
182
183 // **********************************************************************************************
184
185 proton = protonDef->GetParticleName();
186 tableFile[proton] = fileProton;
187 lowEnergyLimit[proton] = lowEnergyLimitForA[1];
188 highEnergyLimit[proton] = 500. * keV;
189
190 // Cross section
191
192 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
193 eV,
194 scaleFactor );
195 tableProton->LoadData(fileProton);
196 tableData[proton] = tableProton;
197
198 // **********************************************************************************************
199
200 hydrogen = hydrogenDef->GetParticleName();
201 tableFile[hydrogen] = fileHydrogen;
202
203 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
204 highEnergyLimit[hydrogen] = 100. * MeV;
205
206 // Cross section
207
208 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
209 eV,
210 scaleFactor );
211 tableHydrogen->LoadData(fileHydrogen);
212
213 tableData[hydrogen] = tableHydrogen;
214
215 // **********************************************************************************************
216
217 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
218 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
219
220 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
221 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
222
223 // Cross section
224
225 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
226 eV,
227 scaleFactor );
228 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
229
230 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
231
232 // **********************************************************************************************
233
234 alphaPlus = alphaPlusDef->GetParticleName();
235 tableFile[alphaPlus] = fileAlphaPlus;
236
237 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
238 highEnergyLimit[alphaPlus] = 400. * MeV;
239
240 // Cross section
241
242 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
243 eV,
244 scaleFactor );
245 tableAlphaPlus->LoadData(fileAlphaPlus);
246 tableData[alphaPlus] = tableAlphaPlus;
247
248 // **********************************************************************************************
249
250 helium = heliumDef->GetParticleName();
251 tableFile[helium] = fileHelium;
252
253 lowEnergyLimit[helium] = lowEnergyLimitForA[4];
254 highEnergyLimit[helium] = 400. * MeV;
255
256 // Cross section
257
258 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
259 eV,
260 scaleFactor );
261 tableHelium->LoadData(fileHelium);
262 tableData[helium] = tableHelium;
263
264 // **********************************************************************************************
265
266 lithium = lithiumDef->GetParticleName();
267 tableFile[lithium] = fileLithium;
268
269 //SI
270 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
271 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
272 lowEnergyLimit[lithium] = 0.5*7*MeV;
273 highEnergyLimit[lithium] = 1e6*7*MeV;
274 //
275
276 // Cross section
277
278 G4DNACrossSectionDataSet* tableLithium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
279 eV,
280 scaleFactor );
281 tableLithium->LoadData(fileLithium);
282 tableData[lithium] = tableLithium;
283
284 // **********************************************************************************************
285
286 beryllium = berylliumDef->GetParticleName();
287 tableFile[beryllium] = fileBeryllium;
288
289 //SI
290 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
291 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
292 lowEnergyLimit[beryllium] = 0.5*9*MeV;
293 highEnergyLimit[beryllium] = 1e6*9*MeV;
294 //
295
296 // Cross section
297
298 G4DNACrossSectionDataSet* tableBeryllium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
299 eV,
300 scaleFactor );
301 tableBeryllium->LoadData(fileBeryllium);
302 tableData[beryllium] = tableBeryllium;
303
304 // **********************************************************************************************
305
306 boron = boronDef->GetParticleName();
307 tableFile[boron] = fileBoron;
308
309 //SI
310 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
311 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
312 lowEnergyLimit[boron] = 0.5*11*MeV;
313 highEnergyLimit[boron] = 1e6*11*MeV;
314 //
315
316 // Cross section
317
318 G4DNACrossSectionDataSet* tableBoron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
319 eV,
320 scaleFactor );
321 tableBoron->LoadData(fileBoron);
322 tableData[boron] = tableBoron;
323
324 // **********************************************************************************************
325
326 carbon = carbonDef->GetParticleName();
327 tableFile[carbon] = fileCarbon;
328
329 //SI
330 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
331 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
332 lowEnergyLimit[carbon] = 0.5*12*MeV;
333 highEnergyLimit[carbon] = 1e6*12*MeV;
334 //
335
336 // Cross section
337
338 G4DNACrossSectionDataSet* tableCarbon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
339 eV,
340 scaleFactor );
341 tableCarbon->LoadData(fileCarbon);
342 tableData[carbon] = tableCarbon;
343
344 // **********************************************************************************************
345
346 oxygen = oxygenDef->GetParticleName();
347 tableFile[oxygen] = fileOxygen;
348
349 //SI
350 //lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
351 //highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
352 lowEnergyLimit[oxygen] = 0.5*16*MeV;
353 highEnergyLimit[oxygen] = 1e6*16*MeV;
354 //
355
356 // Cross section
357
358 G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
359 eV,
360 scaleFactor );
361 tableOxygen->LoadData(fileOxygen);
362 tableData[oxygen] = tableOxygen;
363
364 // **********************************************************************************************
365
366 nitrogen = nitrogenDef->GetParticleName();
367 tableFile[nitrogen] = fileNitrogen;
368
369 //SI
370 //lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
371 //highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
372 lowEnergyLimit[nitrogen] = 0.5*14*MeV;
373 highEnergyLimit[nitrogen] = 1e6*14*MeV;
374 //
375
376 // Cross section
377
378 G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
379 eV,
380 scaleFactor );
381 tableNitrogen->LoadData(fileNitrogen);
382 tableData[nitrogen] = tableNitrogen;
383
384 // **********************************************************************************************
385
386 silicon = siliconDef->GetParticleName();
387 tableFile[silicon] = fileSilicon;
388
389 //lowEnergyLimit[silicon] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
390 //highEnergyLimit[silicon] = 1e6* particle->GetAtomicMass()* MeV;
391 lowEnergyLimit[silicon] = 0.5*28*MeV;
392 highEnergyLimit[silicon] = 1e6*28*MeV;
393 //
394
395 // Cross section
396
397 G4DNACrossSectionDataSet* tableSilicon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
398 eV,
399 scaleFactor );
400 tableSilicon->LoadData(fileSilicon);
401 tableData[silicon] = tableSilicon;
402
403 // **********************************************************************************************
404
405 iron = ironDef->GetParticleName();
406 tableFile[iron] = fileIron;
407
408 //SI
409 //lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
410 //highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
411 lowEnergyLimit[iron] = 0.5*56*MeV;
412 highEnergyLimit[iron] = 1e6*56*MeV;
413 //
414
415 // Cross section
416
417 G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
418 eV,
419 scaleFactor );
420 tableIron->LoadData(fileIron);
421 tableData[iron] = tableIron;
422
423 // **********************************************************************************************
424
425 // SI: not anymore
426 // ZF Following lines can be replaced by:
427 // SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
428 // SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
429 // at least for HZE
430
431 if (particle==protonDef)
432 {
433 SetLowEnergyLimit(lowEnergyLimit[proton]);
434 SetHighEnergyLimit(highEnergyLimit[proton]);
435 }
436
437 if (particle==hydrogenDef)
438 {
439 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
440 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
441 }
442
443 if (particle==heliumDef)
444 {
445 SetLowEnergyLimit(lowEnergyLimit[helium]);
446 SetHighEnergyLimit(highEnergyLimit[helium]);
447 }
448
449 if (particle==alphaPlusDef)
450 {
451 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
452 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
453 }
454
455 if (particle==alphaPlusPlusDef)
456 {
457 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
458 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
459 }
460
461 if (particle==lithiumDef)
462 {
463 SetLowEnergyLimit(lowEnergyLimit[lithium]);
464 SetHighEnergyLimit(highEnergyLimit[lithium]);
465 }
466
467 if (particle==berylliumDef)
468 {
469 SetLowEnergyLimit(lowEnergyLimit[beryllium]);
470 SetHighEnergyLimit(highEnergyLimit[beryllium]);
471 }
472
473 if (particle==boronDef)
474 {
475 SetLowEnergyLimit(lowEnergyLimit[boron]);
476 SetHighEnergyLimit(highEnergyLimit[boron]);
477 }
478
479 if (particle==carbonDef)
480 {
481 SetLowEnergyLimit(lowEnergyLimit[carbon]);
482 SetHighEnergyLimit(highEnergyLimit[carbon]);
483 }
484
485 if (particle==nitrogenDef)
486 {
487 SetLowEnergyLimit(lowEnergyLimit[nitrogen]);
488 SetHighEnergyLimit(highEnergyLimit[nitrogen]);
489 }
490
491 if (particle==oxygenDef)
492 {
493 SetLowEnergyLimit(lowEnergyLimit[oxygen]);
494 SetHighEnergyLimit(highEnergyLimit[oxygen]);
495 }
496
497 if (particle==siliconDef)
498 {
499 SetLowEnergyLimit(lowEnergyLimit[silicon]);
500 SetHighEnergyLimit(highEnergyLimit[silicon]);
501 }
502
503 if (particle==ironDef)
504 {
505 SetLowEnergyLimit(lowEnergyLimit[iron]);
506 SetHighEnergyLimit(highEnergyLimit[iron]);
507 }
508
509 //----------------------------------------------------------------------
510
511 if( verboseLevel>0 )
512 {
513 G4cout << "Rudd ionisation model is initialized " << G4endl
514 << "Energy range: "
515 << LowEnergyLimit() / eV << " eV - "
516 << HighEnergyLimit() / keV << " keV for "
517 << particle->GetParticleName()
518 << G4endl;
519 }
520
521 // Initialize water density pointer
522 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
523
524 //
525
526 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
527
528 if (isInitialised) { return; }
529 fParticleChangeForGamma = GetParticleChangeForGamma();
530 isInitialised = true;
531}
532
533//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
534
535G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
536 const G4ParticleDefinition* particleDefinition,
537 G4double k,
538 G4double,
539 G4double)
540{
541 //SI: particleDefinition->GetParticleName() is for eg. Fe56
542 // particleDefinition->GetPDGMass() is correct
543 // particleDefinition->GetAtomicNumber() is correct
544
545 if (verboseLevel > 3)
546 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
547
548 // Calculate total cross section for model
549
550 G4DNAGenericIonsManager *instance;
551 instance = G4DNAGenericIonsManager::Instance();
552
553 if (
554 particleDefinition != G4Proton::ProtonDefinition()
555 &&
556 particleDefinition != instance->GetIon("hydrogen")
557 &&
558 particleDefinition != instance->GetIon("alpha++")
559 &&
560 particleDefinition != instance->GetIon("alpha+")
561 &&
562 particleDefinition != instance->GetIon("helium")
563 &&
564 // SI
565 //particleDefinition != instance->GetIon("carbon")
566 //&&
567 //particleDefinition != instance->GetIon("nitrogen")
568 //&&
569 //particleDefinition != instance->GetIon("oxygen")
570 //&&
571 //particleDefinition != instance->GetIon("iron")
572 particleDefinition != G4IonTable::GetIonTable()->GetIon(3,7)
573 &&
574 particleDefinition != G4IonTable::GetIonTable()->GetIon(4,9)
575 &&
576 particleDefinition != G4IonTable::GetIonTable()->GetIon(5,11)
577 &&
578 particleDefinition != G4IonTable::GetIonTable()->GetIon(6,12)
579 &&
580 particleDefinition != G4IonTable::GetIonTable()->GetIon(7,14)
581 &&
582 particleDefinition != G4IonTable::GetIonTable()->GetIon(8,16)
583 &&
584 particleDefinition != G4IonTable::GetIonTable()->GetIon(14,28)
585 &&
586 particleDefinition != G4IonTable::GetIonTable()->GetIon(26,56)
587 //
588 )
589
590 return 0;
591
592 G4double lowLim = 0;
593
594 if ( particleDefinition == G4Proton::ProtonDefinition()
595 || particleDefinition == instance->GetIon("hydrogen")
596 )
597
598 lowLim = lowEnergyLimitOfModelForA[1];
599
600 else if ( particleDefinition == instance->GetIon("alpha++")
601 || particleDefinition == instance->GetIon("alpha+")
602 || particleDefinition == instance->GetIon("helium")
603 )
604
605 lowLim = lowEnergyLimitOfModelForA[4];
606
607 else lowLim = lowEnergyLimitOfModelForA[5];
608
609 G4double highLim = 0;
610 G4double sigma=0;
611
612
613 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
614
615 const G4String& particleName = particleDefinition->GetParticleName();
616
617 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
618 pos2 = highEnergyLimit.find(particleName);
619
620 if (pos2 != highEnergyLimit.end())
621 {
622 highLim = pos2->second;
623 }
624
625 if (k <= highLim)
626 {
627
628 //SI : XS must not be zero otherwise sampling of secondaries method ignored
629
630 if (k < lowLim) k = lowLim;
631
632 //
633
634 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
635 pos = tableData.find(particleName);
636
637 if (pos != tableData.end())
638 {
639 G4DNACrossSectionDataSet* table = pos->second;
640 if (table != 0)
641 {
642 sigma = table->FindValue(k);
643 }
644 }
645 else
646 {
647 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
648 FatalException,"Model not applicable to particle type.");
649 }
650
651 } // if (k >= lowLim && k < highLim)
652
653 if (verboseLevel > 2)
654 {
655 G4cout << "__________________________________" << G4endl;
656 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
657 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
658 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
659 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
660 //G4cout << " - Cross section per water molecule (cm^-1)="
661 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
662 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
663
664 }
665
666 return sigma*waterDensity;
667
668}
669
670//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
671
672void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
673 const G4MaterialCutsCouple* couple,
674 const G4DynamicParticle* particle,
675 G4double,
676 G4double)
677{
678 //SI: particle->GetDefinition()->GetParticleName() is for eg. Fe56
679 // particle->GetDefinition()->GetPDGMass() is correct
680 // particle->GetDefinition()->GetAtomicNumber() is correct
681 // particle->GetDefinition()->GetAtomicMass() is correct
682
683 if (verboseLevel > 3)
684 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
685
686 G4double lowLim = 0;
687 G4double highLim = 0;
688
689 // ZF: the following line summarizes the commented part
690
691 if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
692
693 else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
694
695 /*
696
697 if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
698
699 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
700 || particle->GetDefinition() == instance->GetIon("hydrogen")
701 )
702
703 lowLim = killBelowEnergyForA[1];
704
705 if ( particle->GetDefinition() == instance->GetIon("alpha++")
706 || particle->GetDefinition() == instance->GetIon("alpha+")
707 || particle->GetDefinition() == instance->GetIon("helium")
708 )
709
710 lowLim = killBelowEnergyForA[4];
711
712 */
713 //
714
715 G4double k = particle->GetKineticEnergy();
716
717 const G4String& particleName = particle->GetDefinition()->GetParticleName();
718
719 // SI - the following is useless since lowLim is already defined
720 /*
721 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
722 pos1 = lowEnergyLimit.find(particleName);
723
724 if (pos1 != lowEnergyLimit.end())
725 {
726 lowLim = pos1->second;
727 }
728 */
729
730 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
731 pos2 = highEnergyLimit.find(particleName);
732
733 if (pos2 != highEnergyLimit.end()) highLim = pos2->second;
734
735 if (k >= lowLim && k <= highLim)
736
737 // SI: no strict limits, like in the non extended version of the model
738 {
739 G4ParticleDefinition* definition = particle->GetDefinition();
740 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
741 /*
742 G4double particleMass = definition->GetPDGMass();
743 G4double totalEnergy = k + particleMass;
744 G4double pSquare = k*(totalEnergy+particleMass);
745 G4double totalMomentum = std::sqrt(pSquare);
746 */
747
748 G4int ionizationShell = RandomSelect(k,particleName);
749
750 // sample deexcitation
751 // here we assume that H_{2}O electronic levels are the same as Oxygen.
752 // this can be considered true with a rough 10% error in energy on K-shell,
753
754 G4double bindingEnergy = 0;
755 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
756
757 //SI: additional protection if tcs interpolation method is modified
758 if (k<bindingEnergy) return;
759 //
760
761 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
762
763 G4int Z = 8;
764
765 G4ThreeVector deltaDirection =
766 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
767 Z, ionizationShell,
768 couple->GetMaterial());
769
770 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
771 fvect->push_back(dp);
772
773 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
774
775 // SI: the following lines are not needed anymore
776 /*
777 G4double cosTheta = 0.;
778 G4double phi = 0.;
779 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
780
781 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
782 G4double dirX = sinTheta*std::cos(phi);
783 G4double dirY = sinTheta*std::sin(phi);
784 G4double dirZ = cosTheta;
785 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
786 deltaDirection.rotateUz(primaryDirection);
787 */
788
789 // Ignored for ions on electrons
790 /*
791 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
792
793 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
794 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
795 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
796 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
797 finalPx /= finalMomentum;
798 finalPy /= finalMomentum;
799 finalPz /= finalMomentum;
800
801 G4ThreeVector direction;
802 direction.set(finalPx,finalPy,finalPz);
803
804 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
805 */
806
807 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
808 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
809
810 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
811
812 // SI: only atomic deexcitation from K shell is considered
813 if(fAtomDeexcitation && ionizationShell == 4)
814 {
815 const G4AtomicShell* shell
816 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
817 secNumberInit = fvect->size();
818 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
819 secNumberFinal = fvect->size();
820
821 if(secNumberFinal > secNumberInit)
822 {
823 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
824 {
825 //Check if there is enough residual energy
826 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
827 {
828 //Ok, this is a valid secondary: keep it
829 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
830 }
831 else
832 {
833 //Invalid secondary: not enough energy to create it!
834 //Keep its energy in the local deposit
835 delete (*fvect)[i];
836 (*fvect)[i]=0;
837 }
838 }
839 }
840
841 }
842
843 //This should never happen
844 if(bindingEnergy < 0.0)
845 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
846 "em2050",FatalException,"Negative local energy deposit");
847
848 //bindingEnergy has been decreased
849 //by the amount of energy taken away by deexc. products
850 if (!statCode)
851 {
852 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
853 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
854 }
855 else
856 {
857 fParticleChangeForGamma->SetProposedKineticEnergy(k);
858 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
859 }
860
861 // TEST //////////////////////////
862 // if (secondaryKinetic<0) abort();
863 // if (scatteredEnergy<0) abort();
864 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
865 // if (k-scatteredEnergy<0) abort();
866 /////////////////////////////////
867
868 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
869 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
870 ionizationShell,
871 theIncomingTrack);
872 }
873
874 // SI - not useful since low energy of model is 0 eV
875
876 if (k < lowLim)
877 {
878 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
879 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
880 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
881 }
882
883}
884
885//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
886
887G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
888 G4double k,
889 G4int shell)
890{
891 //-- Fast sampling method -----
892 G4double proposed_energy;
893 G4double random1;
894 G4double value_sampling;
895 G4double max1;
896
897 do
898 {
899 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
900
901 max1=0.;
902
903 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
904 max1=RejectionFunction(particleDefinition, k, en, shell);
905
906 random1 = G4UniformRand()*max1;
907
908 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
909
910 } while(random1 > value_sampling);
911
912 return(proposed_energy);
913}
914
915//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
916
917// The following section is not used anymore but is kept for memory
918// GetAngularDistribution()->SampleDirectionForShell is used instead
919
920/*
921void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
922 G4double k,
923 G4double secKinetic,
924 G4double & cosTheta,
925 G4double & phi,
926 G4int shell )
927{
928 G4double maxSecKinetic = 0.;
929 G4double maximumEnergyTransfer = 0.;
930
931 // ZF. generalized & relativistic version
932
933 if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
934 {
935 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
936 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
937 }
938 else
939 {
940 G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
941 G4double en_per_nucleon = k/approx_nuc_number;
942 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
943 G4double gamma = 1./sqrt(1.-beta2);
944 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
945 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
946 }
947
948 maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
949
950 phi = twopi * G4UniformRand();
951
952 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
953 else cosTheta = (2.*G4UniformRand())-1.;
954
955}
956*/
957
958//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
959
960G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
961 G4double k,
962 G4double proposed_ws,
963 G4int ionizationLevelIndex)
964{
965 const G4int j=ionizationLevelIndex;
966 G4double Bj_energy, alphaConst;
967 G4double Ry = 13.6*eV;
968 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
969
970 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
971
972 // Following values provided by M. Dingfelder (priv. comm)
973 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
974
975 if (j == 4)
976 {
977 alphaConst = 0.66;
978 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
979 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
980 //---
981 }
982 else
983 {
984 alphaConst = 0.64;
985 Bj_energy = Bj[ionizationLevelIndex];
986 }
987
988 G4double energyTransfer = proposed_ws + Bj_energy;
989 proposed_ws/=Bj_energy;
990 G4DNAGenericIonsManager *instance;
991 instance = G4DNAGenericIonsManager::Instance();
992 G4double tau = 0.;
993 G4double A_ion = 0.;
994 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
995 A_ion = particleDefinition->GetAtomicMass();
996
997 G4double v2;
998 G4double beta2;
999
1000 if((tau/MeV)<5.447761194e-2)
1001 {
1002 v2 = tau / Bj_energy;
1003 beta2 = 2.*tau / electron_mass_c2;
1004 }
1005 // Relativistic
1006 else
1007 {
1008 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
1009 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
1010 }
1011
1012 G4double v = std::sqrt(v2);
1013 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
1014 G4double rejection_term = 1.+G4Exp(alphaConst*(proposed_ws - wc) / v);
1015 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
1016 //* (S/Bj_energy) ; Not needed anymore
1017
1018 G4bool isHelium = false;
1019
1020 if ( particleDefinition == G4Proton::ProtonDefinition()
1021 || particleDefinition == instance->GetIon("hydrogen")
1022 )
1023 {
1024 return(rejection_term);
1025 }
1026
1027 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
1028 {
1029 G4double Z = particleDefinition->GetAtomicNumber();
1030
1031 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
1032 G4double Zeffion = Z*(1.-G4Exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
1033 rejection_term*=Zeffion*Zeffion;
1034 }
1035
1036 else if (particleDefinition == instance->GetIon("alpha++") )
1037 {
1038 isHelium = true;
1039 slaterEffectiveCharge[0]=0.;
1040 slaterEffectiveCharge[1]=0.;
1041 slaterEffectiveCharge[2]=0.;
1042 sCoefficient[0]=0.;
1043 sCoefficient[1]=0.;
1044 sCoefficient[2]=0.;
1045 }
1046
1047 else if (particleDefinition == instance->GetIon("alpha+") )
1048 {
1049 isHelium = true;
1050 slaterEffectiveCharge[0]=2.0;
1051 // The following values are provided by M. Dingfelder (priv. comm)
1052 slaterEffectiveCharge[1]=2.0;
1053 slaterEffectiveCharge[2]=2.0;
1054 //
1055 sCoefficient[0]=0.7;
1056 sCoefficient[1]=0.15;
1057 sCoefficient[2]=0.15;
1058 }
1059
1060 else if (particleDefinition == instance->GetIon("helium") )
1061 {
1062 isHelium = true;
1063 slaterEffectiveCharge[0]=1.7;
1064 slaterEffectiveCharge[1]=1.15;
1065 slaterEffectiveCharge[2]=1.15;
1066 sCoefficient[0]=0.5;
1067 sCoefficient[1]=0.25;
1068 sCoefficient[2]=0.25;
1069 }
1070
1071 // if ( particleDefinition == instance->GetIon("helium")
1072 // || particleDefinition == instance->GetIon("alpha+")
1073 // || particleDefinition == instance->GetIon("alpha++")
1074 // )
1075
1076 if (isHelium)
1077 {
1078
1079 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
1080
1081 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
1082 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
1083 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
1084
1085 rejection_term*= zEff * zEff;
1086 }
1087
1088 return (rejection_term);
1089}
1090
1091
1092//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1093
1094
1095G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
1096 G4double k,
1097 G4int ionizationLevelIndex)
1098{
1099
1100 const G4int j=ionizationLevelIndex;
1101
1102 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
1103 //G4double alphaConst ;
1104 G4double Bj_energy;
1105
1106 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
1107 // Following values provided by M. Dingfelder (priv. comm)
1108
1109 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
1110
1111 if (j == 4)
1112 {
1113 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
1114 A1 = 1.25;
1115 B1 = 0.5;
1116 C1 = 1.00;
1117 D1 = 1.00;
1118 E1 = 3.00;
1119 A2 = 1.10;
1120 B2 = 1.30;
1121 C2 = 1.00;
1122 D2 = 0.00;
1123 //alphaConst = 0.66;
1124 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
1125 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
1126 //---
1127 }
1128 else
1129 {
1130 //Data For Liquid Water from Dingfelder (Protons in Water)
1131 A1 = 1.02;
1132 B1 = 82.0;
1133 C1 = 0.45;
1134 D1 = -0.80;
1135 E1 = 0.38;
1136 A2 = 1.07;
1137 //B2 = 14.6; From Ding Paper
1138 // Value provided by M. Dingfelder (priv. comm)
1139 B2 = 11.6;
1140 //
1141 C2 = 0.60;
1142 D2 = 0.04;
1143 //alphaConst = 0.64;
1144
1145 Bj_energy = Bj[ionizationLevelIndex];
1146 }
1147
1148 G4double tau = 0.;
1149 G4double A_ion = 0.;
1150 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
1151
1152 A_ion = particle->GetAtomicMass();
1153
1154 G4double v2;
1155 G4double beta2;
1156 if((tau/MeV)<5.447761194e-2)
1157 {
1158 v2 = tau / Bj_energy;
1159 beta2 = 2.*tau / electron_mass_c2;
1160 }
1161 // Relativistic
1162 else
1163 {
1164 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
1165 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
1166 }
1167
1168 G4double v = std::sqrt(v2);
1169 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
1170 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
1171 G4double L2 = C2*std::pow(v,(D2));
1172 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
1173 G4double H2 = (A2/v2) + (B2/(v2*v2));
1174 G4double F1 = L1+H1;
1175 G4double F2 = (L2*H2)/(L2+H2);
1176
1177 // ZF. generalized & relativistic version
1178 G4double maximumEnergy;
1179
1180 //---- maximum kinetic energy , non relativistic ------
1181 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
1182 {
1183 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
1184 }
1185 //---- relativistic -----------------------------------
1186 else
1187 {
1188 G4double gamma = 1./sqrt(1.-beta2);
1189 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
1190 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
1191 }
1192
1193 //either it is transfered energy or secondary electron energy ...
1194 //maximumEnergy-=Bj_energy;
1195
1196 //-----------------------------------------------------
1197 G4double wmax = maximumEnergy/Bj_energy;
1198 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
1199 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
1200 G4double randVal = G4UniformRand();
1201 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
1202 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
1203 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
1204 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1205 proposed_ws*=Bj_energy;
1206
1207 return(proposed_ws);
1208
1209}
1210
1211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1212
1213G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
1214 G4double energyTransferred,
1215 G4double slaterEffectiveChg,
1216 G4double shellNumber)
1217{
1218 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
1219 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
1220
1221 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1222 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1223
1224 return value;
1225}
1226
1227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1228
1229G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
1230 G4double energyTransferred,
1231 G4double slaterEffectiveChg,
1232 G4double shellNumber)
1233{
1234 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
1235 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
1236
1237 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1238 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1239
1240 return value;
1241
1242}
1243
1244//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1245
1246G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
1247 G4double energyTransferred,
1248 G4double slaterEffectiveChg,
1249 G4double shellNumber)
1250{
1251 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
1252 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
1253
1254 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1255 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1256
1257 return value;
1258}
1259
1260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1261
1262G4double G4DNARuddIonisationExtendedModel::R(G4double t,
1263 G4double energyTransferred,
1264 G4double slaterEffectiveChg,
1265 G4double shellNumber)
1266{
1267 // tElectron = m_electron / m_alpha * t
1268 // Dingfelder, in Chattanooga 2005 proceedings, p 4
1269
1270 G4double tElectron = 0.511/3728. * t;
1271 // The following values are provided by M. Dingfelder (priv. comm)
1272 G4double H = 2.*13.60569172 * eV;
1273 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1274
1275 return value;
1276}
1277
1278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1279
1280G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
1281{
1282 // ZF Shortened
1283 G4DNAGenericIonsManager *instance;
1284 instance = G4DNAGenericIonsManager::Instance();
1285
1286 if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
1287 {
1288 G4double value = (std::log10(k/eV)-4.2)/0.5;
1289 // The following values are provided by M. Dingfelder (priv. comm)
1290 return((0.6/(1+G4Exp(value))) + 0.9);
1291 }
1292 else
1293 {
1294 return(1.);
1295 }
1296}
1297
1298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1299
1300G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
1301{
1302
1303 G4int level = 0;
1304
1305 // Retrieve data table corresponding to the current particle type
1306
1307 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1308 pos = tableData.find(particle);
1309
1310 if (pos != tableData.end())
1311 {
1312 G4DNACrossSectionDataSet* table = pos->second;
1313
1314 if (table != 0)
1315 {
1316 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1317
1318 const size_t n(table->NumberOfComponents());
1319 size_t i(n);
1320 G4double value = 0.;
1321
1322 while (i>0)
1323 {
1324 i--;
1325 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1326
1327 value += valuesBuffer[i];
1328 }
1329
1330 value *= G4UniformRand();
1331
1332 i = n;
1333
1334 while (i > 0)
1335 {
1336 i--;
1337
1338 if (valuesBuffer[i] > value)
1339 {
1340 delete[] valuesBuffer;
1341 return i;
1342 }
1343 value -= valuesBuffer[i];
1344 }
1345
1346 if (valuesBuffer) delete[] valuesBuffer;
1347
1348 }
1349 }
1350 else
1351 {
1352 G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
1353 FatalException,"Model not applicable to particle type.");
1354 }
1355
1356 return level;
1357}
1358
1359//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1360
1361G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
1362{
1363 G4double sigma = 0.;
1364
1365 const G4DynamicParticle* particle = track.GetDynamicParticle();
1366 G4double k = particle->GetKineticEnergy();
1367
1368 G4double lowLim = 0;
1369 G4double highLim = 0;
1370
1371 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1372
1373 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1374 pos1 = lowEnergyLimit.find(particleName);
1375
1376 if (pos1 != lowEnergyLimit.end())
1377 {
1378 lowLim = pos1->second;
1379 }
1380
1381 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1382 pos2 = highEnergyLimit.find(particleName);
1383
1384 if (pos2 != highEnergyLimit.end())
1385 {
1386 highLim = pos2->second;
1387 }
1388
1389 if (k >= lowLim && k <= highLim)
1390 {
1391 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1392 pos = tableData.find(particleName);
1393
1394 if (pos != tableData.end())
1395 {
1396 G4DNACrossSectionDataSet* table = pos->second;
1397 if (table != 0)
1398 {
1399 sigma = table->FindValue(k);
1400 }
1401 }
1402 else
1403 {
1404 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
1405 FatalException,"Model not applicable to particle type.");
1406 }
1407 }
1408
1409 return sigma;
1410}
1411
1412//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1413
1414G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1415{
1416 return 0;
1417}
1418
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
eIonizedMolecule
@ eIonizedMolecule
Definition: G4DNAChemistryManager.hh:70
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:46
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
C2
const double C2
Definition: MCGIDI_KalbachMann.cc:15
C1
#define C1
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:132
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:606
G4DNACrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4DNACrossSectionDataSet.cc:483
G4DNACrossSectionDataSet::NumberOfComponents
virtual size_t NumberOfComponents(void) const
Definition: G4DNACrossSectionDataSet.hh:74
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:68
G4DNACrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4DNACrossSectionDataSet.cc:83
G4DNAGenericIonsManager::Instance
static G4DNAGenericIonsManager * Instance(void)
Definition: G4DNAGenericIonsManager.cc:35
G4DNAGenericIonsManager::GetIon
G4ParticleDefinition * GetIon(const G4String &name)
Definition: G4DNAGenericIonsManager.cc:45
G4DNAMolecularMaterial::GetNumMolPerVolTableFor
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
Definition: G4DNAMolecularMaterial.cc:418
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:80
G4DNARuddIonisationExtendedModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4DNARuddIonisationExtendedModel.cc:672
G4DNARuddIonisationExtendedModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4DNARuddIonisationExtendedModel.cc:535
G4DNARuddIonisationExtendedModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4DNARuddIonisationExtendedModel.cc:120
G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
Definition: G4DNARuddIonisationExtendedModel.cc:49
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
G4IonTable::GetIonTable
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:404
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:258
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:651
G4ParticleChangeForGamma::GetCurrentTrack
const G4Track * GetCurrentTrack() const
Definition: G4ParticleChangeForGamma.hh:149
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:119
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:139
G4ParticleDefinition::GetAtomicNumber
G4int GetAtomicNumber() const
G4ParticleDefinition::GetAtomicMass
G4int GetAtomicMass() const
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:112
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.cc:254
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:69
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:611
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:764
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:813
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)