ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNARuddIonisationExtendedModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNARuddIonisationExtendedModel.cc
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 
45 using namespace std;
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 G4DNARuddIonisationExtendedModel::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 
101 G4DNARuddIonisationExtendedModel::~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 
120 void 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 
535 G4double 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 
672 void 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();
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 
887 G4double 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 /*
921 void 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 
960 G4double 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 
1095 G4double 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;
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 
1213 G4double 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 
1229 G4double 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 
1246 G4double 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 
1262 G4double 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 
1280 G4double 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 
1300 G4int 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 
1361 G4double 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 
1414 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1415 {
1416  return 0;
1417 }
1418