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