Geant4  10.01.p02
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 90769 2015-06-09 10:33:41Z 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 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
100 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
101 {
102  // Cross section
103 
104  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
105  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
106  {
107  G4DNACrossSectionDataSet* table = pos->second;
108  delete table;
109  }
110 
111  // The following removal is forbidden G4VEnergyLossModel takes care of deletion
112  // however coverity will signal this as an error
113  // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
114 
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
119 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
120  const G4DataVector& /*cuts*/)
121 {
122  if (verboseLevel > 3)
123  G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
124 
125  // Energy limits
126 
127  G4String fileProton("dna/sigma_ionisation_p_rudd");
128  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
129  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
130  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
131  G4String fileHelium("dna/sigma_ionisation_he_rudd");
132  G4String fileLithium("dna/sigma_ionisation_li_rudd");
133  G4String fileBeryllium("dna/sigma_ionisation_be_rudd");
134  G4String fileBoron("dna/sigma_ionisation_b_rudd");
135  G4String fileCarbon("dna/sigma_ionisation_c_rudd");
136  G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
137  G4String fileOxygen("dna/sigma_ionisation_o_rudd");
138  G4String fileSilicon("dna/sigma_ionisation_si_rudd");
139  G4String fileIron("dna/sigma_ionisation_fe_rudd");
140 
141  G4DNAGenericIonsManager *instance;
142  instance = G4DNAGenericIonsManager::Instance();
143  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
144  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
145  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
146  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
147  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
148 
149  //SEB
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,6);
156  G4ParticleDefinition* berylliumDef = G4IonTable::GetIonTable()->GetIon(4,8);
157  G4ParticleDefinition* boronDef = G4IonTable::GetIonTable()->GetIon(5,10);
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  //SEB
270  //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
271  //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
272  lowEnergyLimit[lithium] = 0.5*6*MeV;
273  highEnergyLimit[lithium] = 1e6*6*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  //SEB
290  //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
291  //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
292  lowEnergyLimit[beryllium] = 0.5*8*MeV;
293  highEnergyLimit[beryllium] = 1e6*8*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  //SEB
310  //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
311  //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
312  lowEnergyLimit[boron] = 0.5*10*MeV;
313  highEnergyLimit[boron] = 1e6*10*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  //SEB
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  //SEB
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  //SEB
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  //SEB
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  //SEB: 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  //SEB: 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  //SEB
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,6)
573  &&
574  particleDefinition != G4IonTable::GetIonTable()->GetIon(4,8)
575  &&
576  particleDefinition != G4IonTable::GetIonTable()->GetIon(5,10)
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  if(waterDensity!= 0.0)
616 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
617  {
618  const G4String& particleName = particleDefinition->GetParticleName();
619 
620  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
621  pos2 = highEnergyLimit.find(particleName);
622 
623  if (pos2 != highEnergyLimit.end())
624  {
625  highLim = pos2->second;
626  }
627 
628  if (k <= highLim)
629  {
630 
631  //SI : XS must not be zero otherwise sampling of secondaries method ignored
632 
633  if (k < lowLim) k = lowLim;
634 
635  //
636 
637  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
638  pos = tableData.find(particleName);
639 
640  if (pos != tableData.end())
641  {
642  G4DNACrossSectionDataSet* table = pos->second;
643  if (table != 0)
644  {
645  sigma = table->FindValue(k);
646  }
647  }
648  else
649  {
650  G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
651  FatalException,"Model not applicable to particle type.");
652  }
653 
654  } // if (k >= lowLim && k < highLim)
655 
656  if (verboseLevel > 2)
657  {
658  G4cout << "__________________________________" << G4endl;
659  G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
660  G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
661  G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
662  G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
663  //G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
664  G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
665 
666  }
667 
668  } // if (waterMaterial)
669 
670  return sigma*waterDensity;
671 // return sigma*material->GetAtomicNumDensityVector()[1];
672 
673 }
674 
675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
676 
677 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
678  const G4MaterialCutsCouple* couple,
679  const G4DynamicParticle* particle,
680  G4double,
681  G4double)
682 {
683  //SEB: particle->GetDefinition()->GetParticleName() is for eg. Fe56
684  // particle->GetDefinition()->GetPDGMass() is correct
685  // particle->GetDefinition()->GetAtomicNumber() is correct
686  // particle->GetDefinition()->GetAtomicMass() is correct
687 
688  if (verboseLevel > 3)
689  G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
690 
691  G4double lowLim = 0;
692  G4double highLim = 0;
693 
694  // ZF: the following line summarizes the commented part
695 
696  if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
697 
698  else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
699 
700  /*
701 
702  if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
703 
704  if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
705  || particle->GetDefinition() == instance->GetIon("hydrogen")
706  )
707 
708  lowLim = killBelowEnergyForA[1];
709 
710  if ( particle->GetDefinition() == instance->GetIon("alpha++")
711  || particle->GetDefinition() == instance->GetIon("alpha+")
712  || particle->GetDefinition() == instance->GetIon("helium")
713  )
714 
715  lowLim = killBelowEnergyForA[4];
716 
717  */
718  //
719 
720  G4double k = particle->GetKineticEnergy();
721 
722  const G4String& particleName = particle->GetDefinition()->GetParticleName();
723 
724  // SI - the following is useless since lowLim is already defined
725  /*
726  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
727  pos1 = lowEnergyLimit.find(particleName);
728 
729  if (pos1 != lowEnergyLimit.end())
730  {
731  lowLim = pos1->second;
732  }
733  */
734 
735  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
736  pos2 = highEnergyLimit.find(particleName);
737 
738  if (pos2 != highEnergyLimit.end())highLim = pos2->second;
739 
740  if (k >= lowLim && k < highLim)
741  {
742  G4ParticleDefinition* definition = particle->GetDefinition();
743  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
744  /*
745  G4double particleMass = definition->GetPDGMass();
746  G4double totalEnergy = k + particleMass;
747  G4double pSquare = k*(totalEnergy+particleMass);
748  G4double totalMomentum = std::sqrt(pSquare);
749  */
750 
751  G4int ionizationShell = RandomSelect(k,particleName);
752 
753  // sample deexcitation
754  // here we assume that H_{2}O electronic levels are the same of Oxigen.
755  // this can be considered true with a rough 10% error in energy on K-shell,
756 
757  G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
758  G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
759  G4double bindingEnergy = 0;
760  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
761 
762  G4int Z = 8;
763 
764  if(fAtomDeexcitation) {
765  G4AtomicShellEnumerator as = fKShell;
766 
767  if (ionizationShell <5 && ionizationShell >1)
768  {
769  as = G4AtomicShellEnumerator(4-ionizationShell);
770  }
771  else if (ionizationShell <2)
772  {
773  as = G4AtomicShellEnumerator(3);
774  }
775 
776  // DEBUG
777  // if (ionizationShell == 4) {
778  //
779  // G4cout << "Z: " << Z << " as: " << as
780  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
781  // G4cout << "Press <Enter> key to continue..." << G4endl;
782  // G4cin.ignore();
783  // }
784 
785 
786  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
787  secNumberInit = fvect->size();
788  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
789  secNumberFinal = fvect->size();
790  }
791 
792  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
793 
794  G4ThreeVector deltaDirection =
795  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
796  Z, ionizationShell,
797  couple->GetMaterial());
798 
799  // SI: the following lines are not needed anymore
800  /*
801  G4double cosTheta = 0.;
802  G4double phi = 0.;
803  RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
804 
805  G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
806  G4double dirX = sinTheta*std::cos(phi);
807  G4double dirY = sinTheta*std::sin(phi);
808  G4double dirZ = cosTheta;
809  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
810  deltaDirection.rotateUz(primaryDirection);
811  */
812 
813  // Ignored for ions on electrons
814  /*
815  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
816 
817  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
818  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
819  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
820  G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
821  finalPx /= finalMomentum;
822  finalPy /= finalMomentum;
823  finalPz /= finalMomentum;
824 
825  G4ThreeVector direction;
826  direction.set(finalPx,finalPy,finalPz);
827 
828  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
829  */
830 
831  fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
832  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
833  G4double deexSecEnergy = 0;
834  for (G4int j=secNumberInit; j < secNumberFinal; j++) {
835 
836  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
837 
838  }
839 
840  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
841  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
842 
843  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
844  fvect->push_back(dp);
845 
846  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
847  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
848  ionizationShell,
849  theIncomingTrack);
850  }
851 
852  // SI - not useful since low energy of model is 0 eV
853 
854  if (k < lowLim)
855  {
856  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
857  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
858  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
859  }
860 
861 }
862 
863 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
864 
865 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
866  G4double k,
867  G4int shell)
868 {
869  //-- Fast sampling method -----
870  G4double proposed_energy;
871  G4double random1;
872  G4double value_sampling;
873  G4double max1;
874 
875  do
876  {
877  proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
878 
879  max1=0.;
880 
881  for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
882  max1=RejectionFunction(particleDefinition, k, en, shell);
883 
884  random1 = G4UniformRand()*max1;
885 
886  value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
887 
888  } while(random1 > value_sampling);
889 
890  return(proposed_energy);
891 }
892 
893 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
894 
895 // The following section is not used anymore but is kept for memory
896 // GetAngularDistribution()->SampleDirectionForShell is used instead
897 
898 /*
899 void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
900  G4double k,
901  G4double secKinetic,
902  G4double & cosTheta,
903  G4double & phi,
904  G4int shell )
905 {
906  G4double maxSecKinetic = 0.;
907  G4double maximumEnergyTransfer = 0.;
908 
909  // ZF. generalized & relativistic version
910 
911  if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
912  {
913  maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
914  maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
915  }
916  else
917  {
918  G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
919  G4double en_per_nucleon = k/approx_nuc_number;
920  G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
921  G4double gamma = 1./sqrt(1.-beta2);
922  maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
923  maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
924  }
925 
926  maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
927 
928  phi = twopi * G4UniformRand();
929 
930  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
931  else cosTheta = (2.*G4UniformRand())-1.;
932 
933 }
934 */
935 
936 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
937 
938 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
939  G4double k,
940  G4double proposed_ws,
941  G4int ionizationLevelIndex)
942 {
943  const G4int j=ionizationLevelIndex;
944  G4double Bj_energy, alphaConst;
945  G4double Ry = 13.6*eV;
946  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
947 
948  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
949 
950  // Following values provided by M. Dingfelder (priv. comm)
951  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
952 
953  if (j == 4)
954  {
955  alphaConst = 0.66;
956  //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
957  Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
958  //---
959  }
960  else
961  {
962  alphaConst = 0.64;
963  Bj_energy = Bj[ionizationLevelIndex];
964  }
965 
966  G4double energyTransfer = proposed_ws + Bj_energy;
967  proposed_ws/=Bj_energy;
968  G4DNAGenericIonsManager *instance;
969  instance = G4DNAGenericIonsManager::Instance();
970  G4double tau = 0.;
971  G4double A_ion = 0.;
972  tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
973  A_ion = particleDefinition->GetAtomicMass();
974 
975  G4double v2;
976  G4double beta2;
977 
978  if((tau/MeV)<5.447761194e-2)
979  {
980  v2 = tau / Bj_energy;
981  beta2 = 2.*tau / electron_mass_c2;
982  }
983  // Relativistic
984  else
985  {
986  v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
987  beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
988  }
989 
990  G4double v = std::sqrt(v2);
991  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
992  G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
993  rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
994  //* (S/Bj_energy) ; Not needed anymore
995 
996  G4bool isHelium = false;
997 
998  if ( particleDefinition == G4Proton::ProtonDefinition()
999  || particleDefinition == instance->GetIon("hydrogen")
1000  )
1001  {
1002  return(rejection_term);
1003  }
1004 
1005  else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
1006  {
1007  G4double Z = particleDefinition->GetAtomicNumber();
1008 
1009  G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
1010  G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
1011  rejection_term*=Zeffion*Zeffion;
1012  }
1013 
1014  else if (particleDefinition == instance->GetIon("alpha++") )
1015  {
1016  isHelium = true;
1017  slaterEffectiveCharge[0]=0.;
1018  slaterEffectiveCharge[1]=0.;
1019  slaterEffectiveCharge[2]=0.;
1020  sCoefficient[0]=0.;
1021  sCoefficient[1]=0.;
1022  sCoefficient[2]=0.;
1023  }
1024 
1025  else if (particleDefinition == instance->GetIon("alpha+") )
1026  {
1027  isHelium = true;
1028  slaterEffectiveCharge[0]=2.0;
1029  // The following values are provided by M. Dingfelder (priv. comm)
1030  slaterEffectiveCharge[1]=2.0;
1031  slaterEffectiveCharge[2]=2.0;
1032  //
1033  sCoefficient[0]=0.7;
1034  sCoefficient[1]=0.15;
1035  sCoefficient[2]=0.15;
1036  }
1037 
1038  else if (particleDefinition == instance->GetIon("helium") )
1039  {
1040  isHelium = true;
1041  slaterEffectiveCharge[0]=1.7;
1042  slaterEffectiveCharge[1]=1.15;
1043  slaterEffectiveCharge[2]=1.15;
1044  sCoefficient[0]=0.5;
1045  sCoefficient[1]=0.25;
1046  sCoefficient[2]=0.25;
1047  }
1048 
1049  // if ( particleDefinition == instance->GetIon("helium")
1050  // || particleDefinition == instance->GetIon("alpha+")
1051  // || particleDefinition == instance->GetIon("alpha++")
1052  // )
1053 
1054  if (isHelium)
1055  {
1056 
1057  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
1058 
1059  zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
1060  sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
1061  sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
1062 
1063  rejection_term*= zEff * zEff;
1064  }
1065 
1066  return (rejection_term);
1067 }
1068 
1069 
1070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1071 
1072 
1073 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
1074  G4double k,
1075  G4int ionizationLevelIndex)
1076 {
1077 
1078  const G4int j=ionizationLevelIndex;
1079 
1080  G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
1081  //G4double alphaConst ;
1082  G4double Bj_energy;
1083 
1084  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
1085  // Following values provided by M. Dingfelder (priv. comm)
1086 
1087  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
1088 
1089  if (j == 4)
1090  {
1091  //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
1092  A1 = 1.25;
1093  B1 = 0.5;
1094  C1 = 1.00;
1095  D1 = 1.00;
1096  E1 = 3.00;
1097  A2 = 1.10;
1098  B2 = 1.30;
1099  C2 = 1.00;
1100  D2 = 0.00;
1101  //alphaConst = 0.66;
1102  //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
1103  Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
1104  //---
1105  }
1106  else
1107  {
1108  //Data For Liquid Water from Dingfelder (Protons in Water)
1109  A1 = 1.02;
1110  B1 = 82.0;
1111  C1 = 0.45;
1112  D1 = -0.80;
1113  E1 = 0.38;
1114  A2 = 1.07;
1115  //B2 = 14.6; From Ding Paper
1116  // Value provided by M. Dingfelder (priv. comm)
1117  B2 = 11.6;
1118  //
1119  C2 = 0.60;
1120  D2 = 0.04;
1121  //alphaConst = 0.64;
1122 
1123  Bj_energy = Bj[ionizationLevelIndex];
1124  }
1125 
1126  G4double tau = 0.;
1127  G4double A_ion = 0.;
1128  tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
1129 
1130  A_ion = particle->GetAtomicMass();
1131 
1132  G4double v2;
1133  G4double beta2;
1134  if((tau/MeV)<5.447761194e-2)
1135  {
1136  v2 = tau / Bj_energy;
1137  beta2 = 2.*tau / electron_mass_c2;
1138  }
1139  // Relativistic
1140  else
1141  {
1142  v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
1143  beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
1144  }
1145 
1146  G4double v = std::sqrt(v2);
1147  //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
1148  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
1149  G4double L2 = C2*std::pow(v,(D2));
1150  G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
1151  G4double H2 = (A2/v2) + (B2/(v2*v2));
1152  G4double F1 = L1+H1;
1153  G4double F2 = (L2*H2)/(L2+H2);
1154 
1155  // ZF. generalized & relativistic version
1156  G4double maximumEnergy;
1157 
1158  //---- maximum kinetic energy , non relativistic ------
1159  if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
1160  {
1161  maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
1162  }
1163  //---- relativistic -----------------------------------
1164  else
1165  {
1166  G4double gamma = 1./sqrt(1.-beta2);
1167  maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
1168  (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
1169  }
1170 
1171  //either it is transfered energy or secondary electron energy ...
1172  //maximumEnergy-=Bj_energy;
1173 
1174  //-----------------------------------------------------
1175  G4double wmax = maximumEnergy/Bj_energy;
1176  G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
1177  c=1./c;
1178  G4double randVal = G4UniformRand();
1179  G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
1180  proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
1181  // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
1182  proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1183  proposed_ws*=Bj_energy;
1184 
1185  return(proposed_ws);
1186 
1187 }
1188 
1189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1190 
1191 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
1192  G4double energyTransferred,
1193  G4double slaterEffectiveChg,
1194  G4double shellNumber)
1195 {
1196  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
1197  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
1198 
1199  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1200  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1201 
1202  return value;
1203 }
1204 
1205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1206 
1207 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
1208  G4double energyTransferred,
1209  G4double slaterEffectiveChg,
1210  G4double shellNumber)
1211 {
1212  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
1213  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
1214 
1215  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1216  G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1217 
1218  return value;
1219 
1220 }
1221 
1222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1223 
1224 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
1225  G4double energyTransferred,
1226  G4double slaterEffectiveChg,
1227  G4double shellNumber)
1228 {
1229  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
1230  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
1231 
1232  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1233  G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1234 
1235  return value;
1236 }
1237 
1238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1239 
1240 G4double G4DNARuddIonisationExtendedModel::R(G4double t,
1241  G4double energyTransferred,
1242  G4double slaterEffectiveChg,
1243  G4double shellNumber)
1244 {
1245  // tElectron = m_electron / m_alpha * t
1246  // Dingfelder, in Chattanooga 2005 proceedings, p 4
1247 
1248  G4double tElectron = 0.511/3728. * t;
1249  // The following values are provided by M. Dingfelder (priv. comm)
1250  G4double H = 2.*13.60569172 * eV;
1251  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1252 
1253  return value;
1254 }
1255 
1256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1257 
1258 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
1259 {
1260  // ZF Shortened
1261  G4DNAGenericIonsManager *instance;
1262  instance = G4DNAGenericIonsManager::Instance();
1263 
1264  if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
1265  {
1266  G4double value = (std::log10(k/eV)-4.2)/0.5;
1267  // The following values are provided by M. Dingfelder (priv. comm)
1268  return((0.6/(1+std::exp(value))) + 0.9);
1269  }
1270  else
1271  {
1272  return(1.);
1273  }
1274 }
1275 
1276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1277 
1278 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
1279 {
1280 
1281  G4int level = 0;
1282 
1283  // Retrieve data table corresponding to the current particle type
1284 
1285  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1286  pos = tableData.find(particle);
1287 
1288  if (pos != tableData.end())
1289  {
1290  G4DNACrossSectionDataSet* table = pos->second;
1291 
1292  if (table != 0)
1293  {
1294  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1295 
1296  const size_t n(table->NumberOfComponents());
1297  size_t i(n);
1298  G4double value = 0.;
1299 
1300  while (i>0)
1301  {
1302  i--;
1303  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1304 
1305  value += valuesBuffer[i];
1306  }
1307 
1308  value *= G4UniformRand();
1309 
1310  i = n;
1311 
1312  while (i > 0)
1313  {
1314  i--;
1315 
1316  if (valuesBuffer[i] > value)
1317  {
1318  delete[] valuesBuffer;
1319  return i;
1320  }
1321  value -= valuesBuffer[i];
1322  }
1323 
1324  if (valuesBuffer) delete[] valuesBuffer;
1325 
1326  }
1327  }
1328  else
1329  {
1330  G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
1331  FatalException,"Model not applicable to particle type.");
1332  }
1333 
1334  return level;
1335 }
1336 
1337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1338 
1339 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
1340 {
1341  G4double sigma = 0.;
1342 
1343  const G4DynamicParticle* particle = track.GetDynamicParticle();
1344  G4double k = particle->GetKineticEnergy();
1345 
1346  G4double lowLim = 0;
1347  G4double highLim = 0;
1348 
1349  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1350 
1351  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1352  pos1 = lowEnergyLimit.find(particleName);
1353 
1354  if (pos1 != lowEnergyLimit.end())
1355  {
1356  lowLim = pos1->second;
1357  }
1358 
1359  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1360  pos2 = highEnergyLimit.find(particleName);
1361 
1362  if (pos2 != highEnergyLimit.end())
1363  {
1364  highLim = pos2->second;
1365  }
1366 
1367  if (k >= lowLim && k <= highLim)
1368  {
1369  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1370  pos = tableData.find(particleName);
1371 
1372  if (pos != tableData.end())
1373  {
1374  G4DNACrossSectionDataSet* table = pos->second;
1375  if (table != 0)
1376  {
1377  sigma = table->FindValue(k);
1378  }
1379  }
1380  else
1381  {
1382  G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
1383  FatalException,"Model not applicable to particle type.");
1384  }
1385  }
1386 
1387  return sigma;
1388 }
1389 
1390 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1391 
1392 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1393 {
1394  return 0;
1395 }
1396 
cm
static const double cm
Definition: G4SIunits.hh:106
G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
Definition: G4DNARuddIonisationExtendedModel.cc:52
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:623
MeV
static const double MeV
Definition: G4SIunits.hh:193
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:114
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:616
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:603
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:463
G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *particleDefinition, G4double incomingParticleEnergy, G4int shell)
Definition: G4DNARuddIonisationExtendedModel.cc:865
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:599
G4DNARuddIonisationExtendedModel::lowEnergyLimitOfModelForA
std::map< G4double, G4double > lowEnergyLimitOfModelForA
Definition: G4DNARuddIonisationExtendedModel.hh:85
G4DNARuddIonisationExtendedModel::fpWaterDensity
const std::vector< G4double > * fpWaterDensity
Definition: G4DNARuddIonisationExtendedModel.hh:76
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:677
G4DNARuddIonisationExtendedModel::S_1s
G4double S_1s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
Definition: G4DNARuddIonisationExtendedModel.cc:1191
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:707
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:93
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:938
G4DNARuddIonisationExtendedModel::Sum
G4double Sum(G4double energy, const G4String &particle)
Definition: G4DNARuddIonisationExtendedModel.cc:1392
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:1224
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
C1
#define C1
G4DNARuddIonisationExtendedModel::R
G4double R(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
Definition: G4DNARuddIonisationExtendedModel.cc:1240
G4DNAChemistryManager.hh
G4DNARuddAngle.hh
G4DNARuddIonisationExtendedModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4DNARuddIonisationExtendedModel.cc:119
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:106
eV
static const double eV
Definition: G4SIunits.hh:194
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:69
eIonizedMolecule
Definition: G4DNAChemistryManager.hh:69
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4DNARuddIonisationExtendedModel::PartialCrossSection
G4double PartialCrossSection(const G4Track &track)
Definition: G4DNARuddIonisationExtendedModel.cc:1339
G4DNARuddIonisationExtendedModel::ProposedSampledEnergy
G4double ProposedSampledEnergy(G4ParticleDefinition *particle, G4double k, G4int ionizationLevelIndex)
Definition: G4DNARuddIonisationExtendedModel.cc:1073
G4DeltaAngle.hh
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Method used by DNA physics model to create a water molecule.
Definition: G4DNAChemistryManager.cc:518
G4DNARuddIonisationExtendedModel::highEnergyLimit
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
Definition: G4DNARuddIonisationExtendedModel.hh:82
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:606
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:110
G4DNARuddIonisationExtendedModel.hh
instance
static MCTruthManager * instance
Definition: MCTruthManager.cc:43
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.cc:1198
keV
static const double keV
Definition: G4SIunits.hh:195
G4DNARuddIonisationExtendedModel::S_2s
G4double S_2s(G4double t, G4double energyTransferred, G4double slaterEffectiveChg, G4double shellNumber)
Definition: G4DNARuddIonisationExtendedModel.cc:1207
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4SystemOfUnits.hh
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4DNARuddIonisationExtendedModel::killBelowEnergyForA
std::map< G4double, G4double > killBelowEnergyForA
Definition: G4DNARuddIonisationExtendedModel.hh:85
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:309
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:763
eplus
static const double eplus
Definition: G4SIunits.hh:178
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4DNARuddIonisationExtendedModel::CorrectionFactor
G4double CorrectionFactor(G4ParticleDefinition *particleDefinition, G4double k, G4int shell)
Definition: G4DNARuddIonisationExtendedModel.cc:1258
G4DNARuddIonisationExtendedModel::lowEnergyLimitForA
std::map< G4double, G4double > lowEnergyLimitForA
Definition: G4DNARuddIonisationExtendedModel.hh:85
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:535
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:1278
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
G4DNARuddIonisationExtendedModel::lowEnergyLimit
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
Definition: G4DNARuddIonisationExtendedModel.hh:81
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4DNAGenericIonsManager::GetIon
G4ParticleDefinition * GetIon(const G4String &name)
Definition: G4DNAGenericIonsManager.cc:46