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