Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNARuddIonisationModel.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: G4DNARuddIonisationModel.cc 96060 2016-03-11 12:58:04Z gcosmo $
27 // GEANT4 tag $Name: $
28 //
29 
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4UAtomicDeexcitation.hh"
34 #include "G4LossTableManager.hh"
35 #include "G4DNAChemistryManager.hh"
37 #include "G4DNARuddAngle.hh"
38 #include "G4DeltaAngle.hh"
39 #include "G4Exp.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 using namespace std;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
48  const G4String& nam) :
49  G4VEmModel(nam), isInitialised(false)
50 {
51  fpWaterDensity = 0;
52 
53  slaterEffectiveCharge[0] = 0.;
54  slaterEffectiveCharge[1] = 0.;
55  slaterEffectiveCharge[2] = 0.;
56  sCoefficient[0] = 0.;
57  sCoefficient[1] = 0.;
58  sCoefficient[2] = 0.;
59 
60  lowEnergyLimitForZ1 = 0 * eV;
61  lowEnergyLimitForZ2 = 0 * eV;
62  lowEnergyLimitOfModelForZ1 = 100 * eV;
63  lowEnergyLimitOfModelForZ2 = 1 * keV;
64  killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
65  killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
66 
67  verboseLevel = 0;
68  // Verbosity scale:
69  // 0 = nothing
70  // 1 = warning for energy non-conservation
71  // 2 = details of energy budget
72  // 3 = calculation of cross sections, file openings, sampling of atoms
73  // 4 = entering in methods
74 
75  if (verboseLevel > 0)
76  {
77  G4cout << "Rudd ionisation model is constructed " << G4endl;
78  }
79 
80  // define default angular generator
82 
83  //Mark this model as "applicable" for atomic deexcitation
84  SetDeexcitationFlag(true);
85  fAtomDeexcitation = 0;
87 
88  // Selection of stationary mode
89 
90  statCode = false;
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
96 {
97  // Cross section
98 
99  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
100  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
101  {
102  G4DNACrossSectionDataSet* table = pos->second;
103  delete table;
104  }
105 
106  // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
107  // Coverity however will signal this as an error
108  //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
109 
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
115  const G4DataVector& /*cuts*/)
116 {
117 
118  if (verboseLevel > 3)
119  {
120  G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
121  }
122 
123  // Energy limits
124 
125  G4String fileProton("dna/sigma_ionisation_p_rudd");
126  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
127  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
128  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
129  G4String fileHelium("dna/sigma_ionisation_he_rudd");
130 
134  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
135  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
136  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
137  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
138 
140  G4String hydrogen;
141  G4String alphaPlusPlus;
142  G4String alphaPlus;
143  G4String helium;
144 
145  G4double scaleFactor = 1 * m*m;
146 
147  // LIMITS AND DATA
148 
149  // ********************************************************
150 
151  proton = protonDef->GetParticleName();
152  tableFile[proton] = fileProton;
153 
154  lowEnergyLimit[proton] = lowEnergyLimitForZ1;
155  highEnergyLimit[proton] = 500. * keV;
156 
157  // Cross section
158 
160  eV,
161  scaleFactor );
162  tableProton->LoadData(fileProton);
163  tableData[proton] = tableProton;
164 
165  // ********************************************************
166 
167  hydrogen = hydrogenDef->GetParticleName();
168  tableFile[hydrogen] = fileHydrogen;
169 
170  lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
171  highEnergyLimit[hydrogen] = 100. * MeV;
172 
173  // Cross section
174 
176  eV,
177  scaleFactor );
178  tableHydrogen->LoadData(fileHydrogen);
179 
180  tableData[hydrogen] = tableHydrogen;
181 
182  // ********************************************************
183 
184  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
185  tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
186 
187  lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
188  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
189 
190  // Cross section
191 
193  eV,
194  scaleFactor );
195  tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
196 
197  tableData[alphaPlusPlus] = tableAlphaPlusPlus;
198 
199  // ********************************************************
200 
201  alphaPlus = alphaPlusDef->GetParticleName();
202  tableFile[alphaPlus] = fileAlphaPlus;
203 
204  lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
205  highEnergyLimit[alphaPlus] = 400. * MeV;
206 
207  // Cross section
208 
210  eV,
211  scaleFactor );
212  tableAlphaPlus->LoadData(fileAlphaPlus);
213  tableData[alphaPlus] = tableAlphaPlus;
214 
215  // ********************************************************
216 
217  helium = heliumDef->GetParticleName();
218  tableFile[helium] = fileHelium;
219 
220  lowEnergyLimit[helium] = lowEnergyLimitForZ2;
221  highEnergyLimit[helium] = 400. * MeV;
222 
223  // Cross section
224 
226  eV,
227  scaleFactor );
228  tableHelium->LoadData(fileHelium);
229  tableData[helium] = tableHelium;
230 
231  //
232 
233  if (particle==protonDef)
234  {
235  SetLowEnergyLimit(lowEnergyLimit[proton]);
236  SetHighEnergyLimit(highEnergyLimit[proton]);
237  }
238 
239  if (particle==hydrogenDef)
240  {
241  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
242  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
243  }
244 
245  if (particle==heliumDef)
246  {
247  SetLowEnergyLimit(lowEnergyLimit[helium]);
248  SetHighEnergyLimit(highEnergyLimit[helium]);
249  }
250 
251  if (particle==alphaPlusDef)
252  {
253  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
254  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
255  }
256 
257  if (particle==alphaPlusPlusDef)
258  {
259  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
260  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
261  }
262 
263  if( verboseLevel>0 )
264  {
265  G4cout << "Rudd ionisation model is initialized " << G4endl
266  << "Energy range: "
267  << LowEnergyLimit() / eV << " eV - "
268  << HighEnergyLimit() / keV << " keV for "
269  << particle->GetParticleName()
270  << G4endl;
271  }
272 
273  // Initialize water density pointer
275 
276  //
277 
278  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
279 
280  if (isInitialised)
281  { return;}
283  isInitialised = true;
284 
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
288 
290  const G4ParticleDefinition* particleDefinition,
291  G4double k,
292  G4double,
293  G4double)
294 {
295  if (verboseLevel > 3)
296  {
297  G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
298  << G4endl;
299  }
300 
301  // Calculate total cross section for model
302 
305 
306  if (
307  particleDefinition != G4Proton::ProtonDefinition()
308  &&
309  particleDefinition != instance->GetIon("hydrogen")
310  &&
311  particleDefinition != instance->GetIon("alpha++")
312  &&
313  particleDefinition != instance->GetIon("alpha+")
314  &&
315  particleDefinition != instance->GetIon("helium")
316  )
317 
318  return 0;
319 
320  G4double lowLim = 0;
321 
322  if ( particleDefinition == G4Proton::ProtonDefinition()
323  || particleDefinition == instance->GetIon("hydrogen")
324  )
325 
326  lowLim = lowEnergyLimitOfModelForZ1;
327 
328  if ( particleDefinition == instance->GetIon("alpha++")
329  || particleDefinition == instance->GetIon("alpha+")
330  || particleDefinition == instance->GetIon("helium")
331  )
332 
333  lowLim = lowEnergyLimitOfModelForZ2;
334 
335  G4double highLim = 0;
336  G4double sigma=0;
337 
338  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
339 
340  if(waterDensity!= 0.0)
341  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
342  {
343  const G4String& particleName = particleDefinition->GetParticleName();
344 
345  // SI - the following is useless since lowLim is already defined
346  /*
347  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
348  pos1 = lowEnergyLimit.find(particleName);
349 
350  if (pos1 != lowEnergyLimit.end())
351  {
352  lowLim = pos1->second;
353  }
354  */
355 
356  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
357  pos2 = highEnergyLimit.find(particleName);
358 
359  if (pos2 != highEnergyLimit.end())
360  {
361  highLim = pos2->second;
362  }
363 
364  if (k <= highLim)
365  {
366  //SI : XS must not be zero otherwise sampling of secondaries method ignored
367 
368  if (k < lowLim) k = lowLim;
369 
370  //
371 
372  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
373  pos = tableData.find(particleName);
374 
375  if (pos != tableData.end())
376  {
377  G4DNACrossSectionDataSet* table = pos->second;
378  if (table != 0)
379  {
380  sigma = table->FindValue(k);
381  }
382  }
383  else
384  {
385  G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
386  FatalException,"Model not applicable to particle type.");
387  }
388 
389  } // if (k >= lowLim && k < highLim)
390 
391  if (verboseLevel > 2)
392  {
393  G4cout << "__________________________________" << G4endl;
394  G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
395  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
396  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
397  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
398  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
399  G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
400  }
401 
402  } // if (waterMaterial)
403  else
404  {
405  if (verboseLevel > 2)
406  {
407  G4cout << "Warning : RuddIonisationModel: WATER DENSITY IS NULL" << G4endl;
408  }
409  }
410 
411  return sigma*waterDensity;
412  // return sigma*material->GetAtomicNumDensityVector()[1];
413 
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417 
418 void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
419  const G4MaterialCutsCouple* couple,
420  const G4DynamicParticle* particle,
421  G4double,
422  G4double)
423 {
424  if (verboseLevel > 3)
425  {
426  G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
427  << G4endl;
428  }
429 
430  G4double lowLim = 0;
431  G4double highLim = 0;
432 
435 
436  if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
437  || particle->GetDefinition() == instance->GetIon("hydrogen")
438  )
439 
440  lowLim = killBelowEnergyForZ1;
441 
442  if ( particle->GetDefinition() == instance->GetIon("alpha++")
443  || particle->GetDefinition() == instance->GetIon("alpha+")
444  || particle->GetDefinition() == instance->GetIon("helium")
445  )
446 
447  lowLim = killBelowEnergyForZ2;
448 
449  G4double k = particle->GetKineticEnergy();
450 
451  const G4String& particleName = particle->GetDefinition()->GetParticleName();
452 
453  // SI - the following is useless since lowLim is already defined
454  /*
455  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
456  pos1 = lowEnergyLimit.find(particleName);
457 
458  if (pos1 != lowEnergyLimit.end())
459  {
460  lowLim = pos1->second;
461  }
462  */
463 
464  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
465  pos2 = highEnergyLimit.find(particleName);
466 
467  if (pos2 != highEnergyLimit.end())
468  {
469  highLim = pos2->second;
470  }
471 
472  if (k >= lowLim && k <= highLim)
473  {
474  G4ParticleDefinition* definition = particle->GetDefinition();
475  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
476  /*
477  G4double particleMass = definition->GetPDGMass();
478  G4double totalEnergy = k + particleMass;
479  G4double pSquare = k*(totalEnergy+particleMass);
480  G4double totalMomentum = std::sqrt(pSquare);
481  */
482 
483  G4int ionizationShell = RandomSelect(k,particleName);
484 
485  // sample deexcitation
486  // here we assume that H_{2}O electronic levels are the same of Oxigen.
487  // this can be considered true with a rough 10% error in energy on K-shell,
488 
489  G4int secNumberInit = 0;// need to know at a certain point the enrgy of secondaries
490  G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
491 
493  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
494 
495  //SI: additional protection if tcs interpolation method is modified
496  if (k<bindingEnergy) return;
497  //
498 
499  G4int Z = 8;
500  if(fAtomDeexcitation)
501  {
503 
504  if (ionizationShell <5 && ionizationShell >1)
505  {
506  as = G4AtomicShellEnumerator(4-ionizationShell);
507  }
508  else if (ionizationShell <2)
509  {
510  as = G4AtomicShellEnumerator(3);
511  }
512 
513  // DEBUG
514  // if (ionizationShell == 4) {
515  //
516  // G4cout << "Z: " << Z << " as: " << as
517  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
518  // G4cout << "Press <Enter> key to continue..." << G4endl;
519  // G4cin.ignore();
520  // }
521 
522  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
523  secNumberInit = fvect->size();
524  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
525  secNumberFinal = fvect->size();
526  }
527 
528  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
529 
530  G4ThreeVector deltaDirection =
531  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
532  Z, ionizationShell,
533  couple->GetMaterial());
534 
535  // Ignored for ions on electrons
536  /*
537  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
538 
539  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
540  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
541  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
542  G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
543  finalPx /= finalMomentum;
544  finalPy /= finalMomentum;
545  finalPz /= finalMomentum;
546 
547  G4ThreeVector direction;
548  direction.set(finalPx,finalPy,finalPz);
549 
550  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
551  */
553 
554  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
555  G4double deexSecEnergy = 0;
556  for (G4int j=secNumberInit; j < secNumberFinal; j++)
557  {
558 
559  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
560 
561  }
562 
563  if (!statCode)
564  {
566  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
567  }
568  else
569  {
572  }
573 
574  // debug
575  // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
576  // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
577  // = bindingEnergy-deexSecEnergy
578  // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
579 
580  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
581  fvect->push_back(dp);
582 
583  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
585  ionizationShell,
586  theIncomingTrack);
587  }
588 
589  // SI - not useful since low energy of model is 0 eV
590 
591  if (k < lowLim)
592  {
596  }
597 
598 }
599 
600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
601 
602 G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
603  G4double k,
604  G4int shell)
605 {
606  G4double maximumKineticEnergyTransfer = 0.;
607 
610 
611  if (particleDefinition == G4Proton::ProtonDefinition()
612  || particleDefinition == instance->GetIon("hydrogen"))
613  {
614  maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
615  }
616 
617  else if (particleDefinition == instance->GetIon("helium")
618  || particleDefinition == instance->GetIon("alpha+")
619  || particleDefinition == instance->GetIon("alpha++"))
620  {
621  maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
622  }
623 
624  G4double crossSectionMaximum = 0.;
625 
626  for (G4double value = waterStructure.IonisationEnergy(shell);
627  value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
628  value += 0.1 * eV)
629  {
630  G4double differentialCrossSection =
631  DifferentialCrossSection(particleDefinition, k, value, shell);
632  if (differentialCrossSection >= crossSectionMaximum)
633  crossSectionMaximum = differentialCrossSection;
634  }
635 
636  G4double secElecKinetic = 0.;
637 
638  do
639  {
640  secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
641  }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
642  k,
643  secElecKinetic+waterStructure.IonisationEnergy(shell),
644  shell));
645 
646  return (secElecKinetic);
647 }
648 
649 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650 
651 // The following section is not used anymore but is kept for memory
652 // GetAngularDistribution()->SampleDirectionForShell is used instead
653 
654 /*
655  void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
656  G4double k,
657  G4double secKinetic,
658  G4double & cosTheta,
659  G4double & phi )
660  {
661  G4DNAGenericIonsManager *instance;
662  instance = G4DNAGenericIonsManager::Instance();
663 
664  G4double maxSecKinetic = 0.;
665 
666  if (particleDefinition == G4Proton::ProtonDefinition()
667  || particleDefinition == instance->GetIon("hydrogen"))
668  {
669  maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
670  }
671 
672  else if (particleDefinition == instance->GetIon("helium")
673  || particleDefinition == instance->GetIon("alpha+")
674  || particleDefinition == instance->GetIon("alpha++"))
675  {
676  maxSecKinetic = 4.* (0.511 / 3728) * k;
677  }
678 
679  phi = twopi * G4UniformRand();
680 
681  // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
682 
683  // Restriction below 100 eV from Emfietzoglou (2000)
684 
685  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
686  else cosTheta = (2.*G4UniformRand())-1.;
687 
688  }
689  */
690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
691 G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
692  G4double k,
693  G4double energyTransfer,
694  G4int ionizationLevelIndex)
695 {
696  // Shells ids are 0 1 2 3 4 (4 is k shell)
697  // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
698  // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
699  //
700  // ds S F1(nu) + w * F2(nu)
701  // ---- = G(k) * ---- -------------------------------------------
702  // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
703  //
704  // w is the secondary electron kinetic Energy in eV
705  //
706  // All the other parameters can be found in Rudd's Papers
707  //
708  // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
709  // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
710  //
711 
712  const G4int j = ionizationLevelIndex;
713 
714  G4double A1;
715  G4double B1;
716  G4double C1;
717  G4double D1;
718  G4double E1;
719  G4double A2;
720  G4double B2;
721  G4double C2;
722  G4double D2;
723  G4double alphaConst;
724 
725  // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
726  // The following values are provided by M. dingfelder (priv. comm)
727  const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
728  * eV };
729 
730  if (j == 4)
731  {
732  //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
733  A1 = 1.25;
734  B1 = 0.5;
735  C1 = 1.00;
736  D1 = 1.00;
737  E1 = 3.00;
738  A2 = 1.10;
739  B2 = 1.30;
740  C2 = 1.00;
741  D2 = 0.00;
742  alphaConst = 0.66;
743  } else
744  {
745  //Data For Liquid Water from Dingfelder (Protons in Water)
746  A1 = 1.02;
747  B1 = 82.0;
748  C1 = 0.45;
749  D1 = -0.80;
750  E1 = 0.38;
751  A2 = 1.07;
752  // Value provided by M. Dingfelder (priv. comm)
753  B2 = 11.6;
754  //
755  C2 = 0.60;
756  D2 = 0.04;
757  alphaConst = 0.64;
758  }
759 
760  const G4double n = 2.;
761  const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
762 
765 
766  G4double wBig = (energyTransfer
767  - waterStructure.IonisationEnergy(ionizationLevelIndex));
768  if (wBig < 0)
769  return 0.;
770 
771  G4double w = wBig / Bj[ionizationLevelIndex];
772  // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
773  if (j == 4)
774  w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
775 
776  G4double Ry = 13.6 * eV;
777 
778  G4double tau = 0.;
779 
780  G4bool isProtonOrHydrogen = false;
781  G4bool isHelium = false;
782 
783  if (particleDefinition == G4Proton::ProtonDefinition()
784  || particleDefinition == instance->GetIon("hydrogen"))
785  {
786  isProtonOrHydrogen = true;
787  tau = (electron_mass_c2 / proton_mass_c2) * k;
788  }
789 
790  else if (particleDefinition == instance->GetIon("helium")
791  || particleDefinition == instance->GetIon("alpha+")
792  || particleDefinition == instance->GetIon("alpha++"))
793  {
794  isHelium = true;
795  tau = (0.511 / 3728.) * k;
796  }
797 
798  G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
799  * std::pow((Ry / Bj[ionizationLevelIndex]), 2);
800  if (j == 4)
801  S = 4. * pi * Bohr_radius * Bohr_radius * n
802  * std::pow((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
803  2);
804 
805  G4double v2 = tau / Bj[ionizationLevelIndex];
806  if (j == 4)
807  v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
808 
809  G4double v = std::sqrt(v2);
810  G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
811  if (j == 4)
812  wc = 4. * v2 - 2. * v
813  - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
814 
815  G4double L1 = (C1 * std::pow(v, (D1))) / (1. + E1 * std::pow(v, (D1 + 4.)));
816  G4double L2 = C2 * std::pow(v, (D2));
817  G4double H1 = (A1 * std::log(1. + v2)) / (v2 + (B1 / v2));
818  G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
819 
820  G4double F1 = L1 + H1;
821  G4double F2 = (L2 * H2) / (L2 + H2);
822 
823  G4double sigma =
824  CorrectionFactor(particleDefinition, k) * Gj[j]
825  * (S / Bj[ionizationLevelIndex])
826  * ((F1 + w * F2)
827  / (std::pow((1. + w), 3)
828  * (1. + G4Exp(alphaConst * (w - wc) / v))));
829 
830  if (j == 4)
831  sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
832  * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
833  * ((F1 + w * F2)
834  / (std::pow((1. + w), 3)
835  * (1. + G4Exp(alphaConst * (w - wc) / v))));
836 
837  if ((particleDefinition == instance->GetIon("hydrogen"))
838  && (ionizationLevelIndex == 4))
839  {
840  // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
841  sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
842  * ((F1 + w * F2)
843  / (std::pow((1. + w), 3)
844  * (1. + G4Exp(alphaConst * (w - wc) / v))));
845  }
846 // if ( particleDefinition == G4Proton::ProtonDefinition()
847 // || particleDefinition == instance->GetIon("hydrogen")
848 // )
849  if (isProtonOrHydrogen)
850  {
851  return (sigma);
852  }
853 
854  if (particleDefinition == instance->GetIon("alpha++"))
855  {
856  slaterEffectiveCharge[0] = 0.;
857  slaterEffectiveCharge[1] = 0.;
858  slaterEffectiveCharge[2] = 0.;
859  sCoefficient[0] = 0.;
860  sCoefficient[1] = 0.;
861  sCoefficient[2] = 0.;
862  }
863 
864  else if (particleDefinition == instance->GetIon("alpha+"))
865  {
866  slaterEffectiveCharge[0] = 2.0;
867  // The following values are provided by M. Dingfelder (priv. comm)
868  slaterEffectiveCharge[1] = 2.0;
869  slaterEffectiveCharge[2] = 2.0;
870  //
871  sCoefficient[0] = 0.7;
872  sCoefficient[1] = 0.15;
873  sCoefficient[2] = 0.15;
874  }
875 
876  else if (particleDefinition == instance->GetIon("helium"))
877  {
878  slaterEffectiveCharge[0] = 1.7;
879  slaterEffectiveCharge[1] = 1.15;
880  slaterEffectiveCharge[2] = 1.15;
881  sCoefficient[0] = 0.5;
882  sCoefficient[1] = 0.25;
883  sCoefficient[2] = 0.25;
884  }
885 
886 // if ( particleDefinition == instance->GetIon("helium")
887 // || particleDefinition == instance->GetIon("alpha+")
888 // || particleDefinition == instance->GetIon("alpha++")
889 // )
890  if (isHelium)
891  {
892  sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
893  * ((F1 + w * F2)
894  / (std::pow((1. + w), 3)
895  * (1. + G4Exp(alphaConst * (w - wc) / v))));
896 
897  if (j == 4)
898  sigma = Gj[j]
899  * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
900  * ((F1 + w * F2)
901  / (std::pow((1. + w), 3)
902  * (1. + G4Exp(alphaConst * (w - wc) / v))));
903 
904  G4double zEff = particleDefinition->GetPDGCharge() / eplus
905  + particleDefinition->GetLeptonNumber();
906 
907  zEff -= (sCoefficient[0]
908  * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
909  + sCoefficient[1]
910  * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
911  + sCoefficient[2]
912  * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
913 
914  return zEff * zEff * sigma;
915  }
916 
917  return 0;
918 }
919 
920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
921 
922 G4double G4DNARuddIonisationModel::S_1s(G4double t,
923  G4double energyTransferred,
924  G4double slaterEffectiveChg,
925  G4double shellNumber)
926 {
927  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
928  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
929 
930  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
931  G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
932 
933  return value;
934 }
935 
936 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
937 
938 G4double G4DNARuddIonisationModel::S_2s(G4double t,
939  G4double energyTransferred,
940  G4double slaterEffectiveChg,
941  G4double shellNumber)
942 {
943  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
944  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
945 
946  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
947  G4double value = 1.
948  - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
949 
950  return value;
951 
952 }
953 
954 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
955 
956 G4double G4DNARuddIonisationModel::S_2p(G4double t,
957  G4double energyTransferred,
958  G4double slaterEffectiveChg,
959  G4double shellNumber)
960 {
961  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
962  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
963 
964  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
965  G4double value = 1.
966  - G4Exp(-2 * r)
967  * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
968 
969  return value;
970 }
971 
972 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
973 
974 G4double G4DNARuddIonisationModel::R(G4double t,
975  G4double energyTransferred,
976  G4double slaterEffectiveChg,
977  G4double shellNumber)
978 {
979  // tElectron = m_electron / m_alpha * t
980  // Dingfelder, in Chattanooga 2005 proceedings, p 4
981 
982  G4double tElectron = 0.511 / 3728. * t;
983  // The following values are provided by M. Dingfelder (priv. comm)
984  G4double H = 2. * 13.60569172 * eV;
985  G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
986  * (slaterEffectiveChg / shellNumber);
987 
988  return value;
989 }
990 
991 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
992 
993 G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
994  G4double k)
995 {
998 
999  if (particleDefinition == G4Proton::Proton())
1000  {
1001  return (1.);
1002  } else if (particleDefinition == instance->GetIon("hydrogen"))
1003  {
1004  G4double value = (std::log10(k / eV) - 4.2) / 0.5;
1005  // The following values are provided by M. Dingfelder (priv. comm)
1006  return ((0.6 / (1 + G4Exp(value))) + 0.9);
1007  } else
1008  {
1009  return (1.);
1010  }
1011 }
1012 
1013 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1014 
1015 G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1016  const G4String& particle)
1017 {
1018 
1019  // BEGIN PART 1/2 OF ELECTRON CORRECTION
1020 
1021  // add ONE or TWO electron-water ionisation for alpha+ and helium
1022 
1023  G4int level = 0;
1024 
1025  // Retrieve data table corresponding to the current particle type
1026 
1027  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1028  pos = tableData.find(particle);
1029 
1030  if (pos != tableData.end())
1031  {
1032  G4DNACrossSectionDataSet* table = pos->second;
1033 
1034  if (table != 0)
1035  {
1036  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1037 
1038  const size_t n(table->NumberOfComponents());
1039  size_t i(n);
1040  G4double value = 0.;
1041 
1042  while (i > 0)
1043  {
1044  i--;
1045  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1046  value += valuesBuffer[i];
1047  }
1048 
1049  value *= G4UniformRand();
1050 
1051  i = n;
1052 
1053  while (i > 0)
1054  {
1055  i--;
1056 
1057  if (valuesBuffer[i] > value)
1058  {
1059  delete[] valuesBuffer;
1060  return i;
1061  }
1062  value -= valuesBuffer[i];
1063  }
1064 
1065  if (valuesBuffer)
1066  delete[] valuesBuffer;
1067 
1068  }
1069  } else
1070  {
1071  G4Exception("G4DNARuddIonisationModel::RandomSelect",
1072  "em0002",
1074  "Model not applicable to particle type.");
1075  }
1076 
1077  return level;
1078 }
1079 
1080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1081 
1082 G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1083 {
1084  G4double sigma = 0.;
1085 
1086  const G4DynamicParticle* particle = track.GetDynamicParticle();
1087  G4double k = particle->GetKineticEnergy();
1088 
1089  G4double lowLim = 0;
1090  G4double highLim = 0;
1091 
1092  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1093 
1094  std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1095  pos1 = lowEnergyLimit.find(particleName);
1096 
1097  if (pos1 != lowEnergyLimit.end())
1098  {
1099  lowLim = pos1->second;
1100  }
1101 
1102  std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1103  pos2 = highEnergyLimit.find(particleName);
1104 
1105  if (pos2 != highEnergyLimit.end())
1106  {
1107  highLim = pos2->second;
1108  }
1109 
1110  if (k >= lowLim && k <= highLim)
1111  {
1112  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1113  pos = tableData.find(particleName);
1114 
1115  if (pos != tableData.end())
1116  {
1117  G4DNACrossSectionDataSet* table = pos->second;
1118  if (table != 0)
1119  {
1120  sigma = table->FindValue(k);
1121  }
1122  } else
1123  {
1124  G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1125  "em0002",
1127  "Model not applicable to particle type.");
1128  }
1129  }
1130 
1131  return sigma;
1132 }
1133 
1134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1135 
1136 G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1137  const G4String& /* particle */)
1138 {
1139  return 0;
1140 }
1141 
const double C2
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
const double C1
double S(double temp)
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
G4ParticleChangeForGamma * fParticleChangeForGamma
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
size_t GetIndex() const
Definition: G4Material.hh:262
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
G4DNARuddIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationModel")
string material
Definition: eplot.py:19
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char int const XML_Char * value
Definition: expat.h:331
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double eV
Definition: G4SIunits.hh:215
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4DNAGenericIonsManager * Instance(void)
const G4int n
virtual G4double FindValue(G4double e, G4int componentId=0) const
tuple v
Definition: test.py:18
virtual size_t NumberOfComponents(void) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
const G4Track * GetCurrentTrack() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
static MCTruthManager * instance
static constexpr double pi
Definition: G4SIunits.hh:75
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
G4double GetPDGCharge() const
G4double bindingEnergy(G4int A, G4int Z)
static constexpr double keV
Definition: G4SIunits.hh:216
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
static const G4double pos
G4AtomicShellEnumerator
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const
G4ParticleDefinition * GetIon(const G4String &name)