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