ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNARuddIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNARuddIonisationModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include "G4DNARuddIonisationModel.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNARuddAngle.hh"
36 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh"
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
41 using namespace std;
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 
45 G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
46  const G4String& nam) :
47 G4VEmModel(nam), isInitialised(false)
48 {
49  fpWaterDensity = 0;
50 
51  slaterEffectiveCharge[0] = 0.;
52  slaterEffectiveCharge[1] = 0.;
53  slaterEffectiveCharge[2] = 0.;
54  sCoefficient[0] = 0.;
55  sCoefficient[1] = 0.;
56  sCoefficient[2] = 0.;
57 
58  lowEnergyLimitForZ1 = 0 * eV;
59  lowEnergyLimitForZ2 = 0 * eV;
60  lowEnergyLimitOfModelForZ1 = 100 * eV;
61  lowEnergyLimitOfModelForZ2 = 1 * keV;
62  killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
63  killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
64 
65  verboseLevel = 0;
66  // Verbosity scale:
67  // 0 = nothing
68  // 1 = warning for energy non-conservation
69  // 2 = details of energy budget
70  // 3 = calculation of cross sections, file openings, sampling of atoms
71  // 4 = entering in methods
72 
73  if (verboseLevel > 0)
74  {
75  G4cout << "Rudd ionisation model is constructed " << G4endl;
76  }
77 
78  // Define default angular generator
79  SetAngularDistribution(new G4DNARuddAngle());
80 
81  // Mark this model as "applicable" for atomic deexcitation
82  SetDeexcitationFlag(true);
83  fAtomDeexcitation = 0;
84  fParticleChangeForGamma = 0;
85 
86  // Selection of stationary mode
87 
88  statCode = false;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
94 {
95  // Cross section
96 
97  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
98  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99  {
100  G4DNACrossSectionDataSet* table = pos->second;
101  delete table;
102  }
103 
104  // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
105  // Coverity however will signal this as an error
106  // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107 
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
112 void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
113  const G4DataVector& /*cuts*/)
114 {
115 
116  if (verboseLevel > 3)
117  {
118  G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
119  }
120 
121  // Energy limits
122 
123  G4String fileProton("dna/sigma_ionisation_p_rudd");
124  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
125  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
126  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
127  G4String fileHelium("dna/sigma_ionisation_he_rudd");
128 
129  G4DNAGenericIonsManager *instance;
130  instance = G4DNAGenericIonsManager::Instance();
131  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
132  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
133  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
134  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
135  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
136 
137  G4String proton;
138  G4String hydrogen;
139  G4String alphaPlusPlus;
140  G4String alphaPlus;
141  G4String helium;
142 
143  G4double scaleFactor = 1 * m*m;
144 
145  // LIMITS AND DATA
146 
147  // ********************************************************
148 
149  proton = protonDef->GetParticleName();
150  tableFile[proton] = fileProton;
151 
152  lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153  highEnergyLimit[proton] = 500. * keV;
154 
155  // Cross section
156 
157  G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
158  eV,
159  scaleFactor );
160  tableProton->LoadData(fileProton);
161  tableData[proton] = tableProton;
162 
163  // ********************************************************
164 
165  hydrogen = hydrogenDef->GetParticleName();
166  tableFile[hydrogen] = fileHydrogen;
167 
168  lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
169  highEnergyLimit[hydrogen] = 100. * MeV;
170 
171  // Cross section
172 
173  G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
174  eV,
175  scaleFactor );
176  tableHydrogen->LoadData(fileHydrogen);
177 
178  tableData[hydrogen] = tableHydrogen;
179 
180  // ********************************************************
181 
182  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
183  tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
184 
185  lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
186  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
187 
188  // Cross section
189 
190  G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
191  eV,
192  scaleFactor );
193  tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
194 
195  tableData[alphaPlusPlus] = tableAlphaPlusPlus;
196 
197  // ********************************************************
198 
199  alphaPlus = alphaPlusDef->GetParticleName();
200  tableFile[alphaPlus] = fileAlphaPlus;
201 
202  lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
203  highEnergyLimit[alphaPlus] = 400. * MeV;
204 
205  // Cross section
206 
207  G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
208  eV,
209  scaleFactor );
210  tableAlphaPlus->LoadData(fileAlphaPlus);
211  tableData[alphaPlus] = tableAlphaPlus;
212 
213  // ********************************************************
214 
215  helium = heliumDef->GetParticleName();
216  tableFile[helium] = fileHelium;
217 
218  lowEnergyLimit[helium] = lowEnergyLimitForZ2;
219  highEnergyLimit[helium] = 400. * MeV;
220 
221  // Cross section
222 
223  G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
224  eV,
225  scaleFactor );
226  tableHelium->LoadData(fileHelium);
227  tableData[helium] = tableHelium;
228 
229  //
230 
231  if (particle==protonDef)
232  {
233  SetLowEnergyLimit(lowEnergyLimit[proton]);
234  SetHighEnergyLimit(highEnergyLimit[proton]);
235  }
236 
237  if (particle==hydrogenDef)
238  {
239  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
240  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
241  }
242 
243  if (particle==heliumDef)
244  {
245  SetLowEnergyLimit(lowEnergyLimit[helium]);
246  SetHighEnergyLimit(highEnergyLimit[helium]);
247  }
248 
249  if (particle==alphaPlusDef)
250  {
251  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
252  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
253  }
254 
255  if (particle==alphaPlusPlusDef)
256  {
257  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
258  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
259  }
260 
261  if( verboseLevel>0 )
262  {
263  G4cout << "Rudd ionisation model is initialized " << G4endl
264  << "Energy range: "
265  << LowEnergyLimit() / eV << " eV - "
266  << HighEnergyLimit() / keV << " keV for "
267  << particle->GetParticleName()
268  << G4endl;
269  }
270 
271  // Initialize water density pointer
272  fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
273 
274  //
275 
276  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
277 
278  if (isInitialised)
279  { return;}
280  fParticleChangeForGamma = GetParticleChangeForGamma();
281  isInitialised = true;
282 
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 
287 G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
288  const G4ParticleDefinition* particleDefinition,
289  G4double k,
290  G4double,
291  G4double)
292 {
293  if (verboseLevel > 3)
294  {
295  G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
296  << G4endl;
297  }
298 
299  // Calculate total cross section for model
300 
301  G4DNAGenericIonsManager *instance;
302  instance = G4DNAGenericIonsManager::Instance();
303 
304  if (
305  particleDefinition != G4Proton::ProtonDefinition()
306  &&
307  particleDefinition != instance->GetIon("hydrogen")
308  &&
309  particleDefinition != instance->GetIon("alpha++")
310  &&
311  particleDefinition != instance->GetIon("alpha+")
312  &&
313  particleDefinition != instance->GetIon("helium")
314  )
315 
316  return 0;
317 
318  G4double lowLim = 0;
319 
320  if ( particleDefinition == G4Proton::ProtonDefinition()
321  || particleDefinition == instance->GetIon("hydrogen")
322  )
323 
324  lowLim = lowEnergyLimitOfModelForZ1;
325 
326  if ( particleDefinition == instance->GetIon("alpha++")
327  || particleDefinition == instance->GetIon("alpha+")
328  || particleDefinition == instance->GetIon("helium")
329  )
330 
331  lowLim = lowEnergyLimitOfModelForZ2;
332 
333  G4double highLim = 0;
334  G4double sigma=0;
335 
336  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
337 
338  const G4String& particleName = particleDefinition->GetParticleName();
339 
340  // SI - the following is useless since lowLim is already defined
341  /*
342  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
343  pos1 = lowEnergyLimit.find(particleName);
344 
345  if (pos1 != lowEnergyLimit.end())
346  {
347  lowLim = pos1->second;
348  }
349  */
350 
351  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
352  pos2 = highEnergyLimit.find(particleName);
353 
354  if (pos2 != highEnergyLimit.end())
355  {
356  highLim = pos2->second;
357  }
358 
359  if (k <= highLim)
360  {
361  //SI : XS must not be zero otherwise sampling of secondaries method ignored
362 
363  if (k < lowLim) k = lowLim;
364 
365  //
366 
367  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
368  pos = tableData.find(particleName);
369 
370  if (pos != tableData.end())
371  {
372  G4DNACrossSectionDataSet* table = pos->second;
373  if (table != 0)
374  {
375  sigma = table->FindValue(k);
376  }
377  }
378  else
379  {
380  G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
381  FatalException,"Model not applicable to particle type.");
382  }
383 
384  }
385 
386  if (verboseLevel > 2)
387  {
388  G4cout << "__________________________________" << G4endl;
389  G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
390  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
391  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
392  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
393  //G4cout << " - Cross section per water molecule (cm^-1)="
394  //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
395  G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
396  }
397 
398  return sigma*waterDensity;
399 
400 }
401 
402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
403 
404 void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
405  const G4MaterialCutsCouple* couple,
406  const G4DynamicParticle* particle,
407  G4double,
408  G4double)
409 {
410  if (verboseLevel > 3)
411  {
412  G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
413  << G4endl;
414  }
415 
416  G4double lowLim = 0;
417  G4double highLim = 0;
418 
419  G4DNAGenericIonsManager *instance;
420  instance = G4DNAGenericIonsManager::Instance();
421 
422  if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
423  || particle->GetDefinition() == instance->GetIon("hydrogen")
424  )
425 
426  lowLim = killBelowEnergyForZ1;
427 
428  if ( particle->GetDefinition() == instance->GetIon("alpha++")
429  || particle->GetDefinition() == instance->GetIon("alpha+")
430  || particle->GetDefinition() == instance->GetIon("helium")
431  )
432 
433  lowLim = killBelowEnergyForZ2;
434 
435  G4double k = particle->GetKineticEnergy();
436 
437  const G4String& particleName = particle->GetDefinition()->GetParticleName();
438 
439  // SI - the following is useless since lowLim is already defined
440  /*
441  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
442  pos1 = lowEnergyLimit.find(particleName);
443 
444  if (pos1 != lowEnergyLimit.end())
445  {
446  lowLim = pos1->second;
447  }
448  */
449 
450  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
451  pos2 = highEnergyLimit.find(particleName);
452 
453  if (pos2 != highEnergyLimit.end())
454  {
455  highLim = pos2->second;
456  }
457 
458  if (k >= lowLim && k <= highLim)
459  {
460  G4ParticleDefinition* definition = particle->GetDefinition();
461  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
462  /*
463  G4double particleMass = definition->GetPDGMass();
464  G4double totalEnergy = k + particleMass;
465  G4double pSquare = k*(totalEnergy+particleMass);
466  G4double totalMomentum = std::sqrt(pSquare);
467  */
468 
469  G4int ionizationShell = RandomSelect(k,particleName);
470 
471  G4double bindingEnergy = 0;
472  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
473 
474  //SI: additional protection if tcs interpolation method is modified
475  if (k<bindingEnergy) return;
476  //
477 
478  // SI - For atom. deexc. tagging - 23/05/2017
479  G4int Z = 8;
480  //
481 
482  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
483 
484  G4ThreeVector deltaDirection =
485  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
486  Z, ionizationShell,
487  couple->GetMaterial());
488 
489  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
490  fvect->push_back(dp);
491 
492  // Ignored for ions on electrons
493  /*
494  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
495 
496  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
497  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
498  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
499  G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
500  finalPx /= finalMomentum;
501  finalPy /= finalMomentum;
502  finalPz /= finalMomentum;
503 
504  G4ThreeVector direction;
505  direction.set(finalPx,finalPy,finalPz);
506 
507  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
508  */
509 
510  fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
511 
512  // sample deexcitation
513  // here we assume that H_{2}O electronic levels are the same of Oxigen.
514  // this can be considered true with a rough 10% error in energy on K-shell,
515 
516  size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
517  size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
518 
519  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
520 
521  // SI: only atomic deexcitation from K shell is considered
522  if(fAtomDeexcitation && ionizationShell == 4)
523  {
524  const G4AtomicShell* shell
525  = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
526  secNumberInit = fvect->size();
527  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
528  secNumberFinal = fvect->size();
529 
530  if(secNumberFinal > secNumberInit)
531  {
532  for (size_t i=secNumberInit; i<secNumberFinal; ++i)
533  {
534  //Check if there is enough residual energy
535  if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
536  {
537  //Ok, this is a valid secondary: keep it
538  bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
539  }
540  else
541  {
542  //Invalid secondary: not enough energy to create it!
543  //Keep its energy in the local deposit
544  delete (*fvect)[i];
545  (*fvect)[i]=0;
546  }
547  }
548  }
549 
550  }
551 
552  //This should never happen
553  if(bindingEnergy < 0.0)
554  G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
555  "em2050",FatalException,"Negative local energy deposit");
556 
557  //bindingEnergy has been decreased
558  //by the amount of energy taken away by deexc. products
559  if (!statCode)
560  {
561  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
562  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
563  }
564  else
565  {
566  fParticleChangeForGamma->SetProposedKineticEnergy(k);
567  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
568  }
569 
570  // debug
571  // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
572  // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
573  // = bindingEnergy-deexSecEnergy
574  // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
575 
576  // TEST //////////////////////////
577  // if (secondaryKinetic<0) abort();
578  // if (scatteredEnergy<0) abort();
579  // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
580  // if (k-scatteredEnergy<0) abort();
582 
583  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
584  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
585  ionizationShell,
586  theIncomingTrack);
587  }
588 
589  // SI - not useful since low energy of model is 0 eV
590 
591  if (k < lowLim)
592  {
593  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
594  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
595  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
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 
608  G4DNAGenericIonsManager *instance;
609  instance = G4DNAGenericIonsManager::Instance();
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 
763  G4DNAGenericIonsManager* instance;
764  instance = G4DNAGenericIonsManager::Instance();
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 
847  // if ( particleDefinition == G4Proton::ProtonDefinition()
848  // || particleDefinition == instance->GetIon("hydrogen")
849  // )
850 
851  if (isProtonOrHydrogen)
852  {
853  return (sigma);
854  }
855 
856  if (particleDefinition == instance->GetIon("alpha++"))
857  {
858  slaterEffectiveCharge[0] = 0.;
859  slaterEffectiveCharge[1] = 0.;
860  slaterEffectiveCharge[2] = 0.;
861  sCoefficient[0] = 0.;
862  sCoefficient[1] = 0.;
863  sCoefficient[2] = 0.;
864  }
865 
866  else if (particleDefinition == instance->GetIon("alpha+"))
867  {
868  slaterEffectiveCharge[0] = 2.0;
869  // The following values are provided by M. Dingfelder (priv. comm)
870  slaterEffectiveCharge[1] = 2.0;
871  slaterEffectiveCharge[2] = 2.0;
872  //
873  sCoefficient[0] = 0.7;
874  sCoefficient[1] = 0.15;
875  sCoefficient[2] = 0.15;
876  }
877 
878  else if (particleDefinition == instance->GetIon("helium"))
879  {
880  slaterEffectiveCharge[0] = 1.7;
881  slaterEffectiveCharge[1] = 1.15;
882  slaterEffectiveCharge[2] = 1.15;
883  sCoefficient[0] = 0.5;
884  sCoefficient[1] = 0.25;
885  sCoefficient[2] = 0.25;
886  }
887 
888  // if ( particleDefinition == instance->GetIon("helium")
889  // || particleDefinition == instance->GetIon("alpha+")
890  // || particleDefinition == instance->GetIon("alpha++")
891  // )
892  if (isHelium)
893  {
894  sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
895  * ((F1 + w * F2)
896  / (std::pow((1. + w), 3)
897  * (1. + G4Exp(alphaConst * (w - wc) / v))));
898 
899  if (j == 4)
900  sigma = Gj[j]
901  * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
902  * ((F1 + w * F2)
903  / (std::pow((1. + w), 3)
904  * (1. + G4Exp(alphaConst * (w - wc) / v))));
905 
906  G4double zEff = particleDefinition->GetPDGCharge() / eplus
907  + particleDefinition->GetLeptonNumber();
908 
909  zEff -= (sCoefficient[0]
910  * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
911  + sCoefficient[1]
912  * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
913  + sCoefficient[2]
914  * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
915 
916  return zEff * zEff * sigma;
917  }
918 
919  return 0;
920 }
921 
922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
923 
924 G4double G4DNARuddIonisationModel::S_1s(G4double t,
925  G4double energyTransferred,
926  G4double slaterEffectiveChg,
927  G4double shellNumber)
928 {
929  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
930  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
931 
932  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
933  G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
934 
935  return value;
936 }
937 
938 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
939 
940 G4double G4DNARuddIonisationModel::S_2s(G4double t,
941  G4double energyTransferred,
942  G4double slaterEffectiveChg,
943  G4double shellNumber)
944 {
945  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
946  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
947 
948  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
949  G4double value = 1.
950  - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
951 
952  return value;
953 
954 }
955 
956 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
957 
958 G4double G4DNARuddIonisationModel::S_2p(G4double t,
959  G4double energyTransferred,
960  G4double slaterEffectiveChg,
961  G4double shellNumber)
962 {
963  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
964  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
965 
966  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
967  G4double value = 1.
968  - G4Exp(-2 * r)
969  * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
970 
971  return value;
972 }
973 
974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
975 
976 G4double G4DNARuddIonisationModel::R(G4double t,
977  G4double energyTransferred,
978  G4double slaterEffectiveChg,
979  G4double shellNumber)
980 {
981  // tElectron = m_electron / m_alpha * t
982  // Dingfelder, in Chattanooga 2005 proceedings, p 4
983 
984  G4double tElectron = 0.511 / 3728. * t;
985  // The following values are provided by M. Dingfelder (priv. comm)
986  G4double H = 2. * 13.60569172 * eV;
987  G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
988  * (slaterEffectiveChg / shellNumber);
989 
990  return value;
991 }
992 
993 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
994 
995 G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
996  G4double k)
997 {
998  G4DNAGenericIonsManager *instance;
999  instance = G4DNAGenericIonsManager::Instance();
1000 
1001  if (particleDefinition == G4Proton::Proton())
1002  {
1003  return (1.);
1004  } else if (particleDefinition == instance->GetIon("hydrogen"))
1005  {
1006  G4double value = (std::log10(k / eV) - 4.2) / 0.5;
1007  // The following values are provided by M. Dingfelder (priv. comm)
1008  return ((0.6 / (1 + G4Exp(value))) + 0.9);
1009  } else
1010  {
1011  return (1.);
1012  }
1013 }
1014 
1015 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1016 
1017 G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1018  const G4String& particle)
1019 {
1020 
1021  // BEGIN PART 1/2 OF ELECTRON CORRECTION
1022 
1023  // add ONE or TWO electron-water ionisation for alpha+ and helium
1024 
1025  G4int level = 0;
1026 
1027  // Retrieve data table corresponding to the current particle type
1028 
1029  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1030  pos = tableData.find(particle);
1031 
1032  if (pos != tableData.end())
1033  {
1034  G4DNACrossSectionDataSet* table = pos->second;
1035 
1036  if (table != 0)
1037  {
1038  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1039 
1040  const size_t n(table->NumberOfComponents());
1041  size_t i(n);
1042  G4double value = 0.;
1043 
1044  while (i > 0)
1045  {
1046  i--;
1047  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1048  value += valuesBuffer[i];
1049  }
1050 
1051  value *= G4UniformRand();
1052 
1053  i = n;
1054 
1055  while (i > 0)
1056  {
1057  i--;
1058 
1059  if (valuesBuffer[i] > value)
1060  {
1061  delete[] valuesBuffer;
1062  return i;
1063  }
1064  value -= valuesBuffer[i];
1065  }
1066 
1067  if (valuesBuffer)
1068  delete[] valuesBuffer;
1069 
1070  }
1071  } else
1072  {
1073  G4Exception("G4DNARuddIonisationModel::RandomSelect",
1074  "em0002",
1075  FatalException,
1076  "Model not applicable to particle type.");
1077  }
1078 
1079  return level;
1080 }
1081 
1082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1083 
1084 G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1085 {
1086  G4double sigma = 0.;
1087 
1088  const G4DynamicParticle* particle = track.GetDynamicParticle();
1089  G4double k = particle->GetKineticEnergy();
1090 
1091  G4double lowLim = 0;
1092  G4double highLim = 0;
1093 
1094  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1095 
1096  std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1097  pos1 = lowEnergyLimit.find(particleName);
1098 
1099  if (pos1 != lowEnergyLimit.end())
1100  {
1101  lowLim = pos1->second;
1102  }
1103 
1104  std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1105  pos2 = highEnergyLimit.find(particleName);
1106 
1107  if (pos2 != highEnergyLimit.end())
1108  {
1109  highLim = pos2->second;
1110  }
1111 
1112  if (k >= lowLim && k <= highLim)
1113  {
1114  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1115  pos = tableData.find(particleName);
1116 
1117  if (pos != tableData.end())
1118  {
1119  G4DNACrossSectionDataSet* table = pos->second;
1120  if (table != 0)
1121  {
1122  sigma = table->FindValue(k);
1123  }
1124  } else
1125  {
1126  G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1127  "em0002",
1128  FatalException,
1129  "Model not applicable to particle type.");
1130  }
1131  }
1132 
1133  return sigma;
1134 }
1135 
1136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1137 
1138 G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1139  const G4String& /* particle */)
1140 {
1141  return 0;
1142 }
1143