ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAMillerGreenExcitationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAMillerGreenExcitationModel.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 "G4DNAMillerGreenExcitationModel.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4DNAChemistryManager.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
40 G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(const G4ParticleDefinition*,
41  const G4String& nam)
42 :G4VEmModel(nam),isInitialised(false)
43 {
44  fpMolWaterDensity = 0;
45 
46  nLevels=0;
47  kineticEnergyCorrection[0]=0.;
48  kineticEnergyCorrection[1]=0.;
49  kineticEnergyCorrection[2]=0.;
50  kineticEnergyCorrection[3]=0.;
51 
52  verboseLevel= 0;
53  // Verbosity scale:
54  // 0 = nothing
55  // 1 = warning for energy non-conservation
56  // 2 = details of energy budget
57  // 3 = calculation of cross sections, file openings, sampling of atoms
58  // 4 = entering in methods
59 
60  if( verboseLevel>0 )
61  {
62  G4cout << "Miller & Green excitation model is constructed " << G4endl;
63  }
64  fParticleChangeForGamma = 0;
65 
66  // Selection of stationary mode
67 
68  statCode = false;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
73 G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel()
74 {}
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
78 void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle,
79  const G4DataVector& /*cuts*/)
80 {
81 
82  if (verboseLevel > 3)
83  G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl;
84 
85  // Energy limits
86 
87  G4DNAGenericIonsManager *instance;
88  instance = G4DNAGenericIonsManager::Instance();
89  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
90  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
91  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
92  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
93  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
94 
95  G4String proton;
96  G4String hydrogen;
97  G4String alphaPlusPlus;
98  G4String alphaPlus;
99  G4String helium;
100 
101  // LIMITS AND CONSTANTS
102 
103  proton = protonDef->GetParticleName();
104  lowEnergyLimit[proton] = 10. * eV;
105  highEnergyLimit[proton] = 500. * keV;
106 
107  kineticEnergyCorrection[0] = 1.;
108  slaterEffectiveCharge[0][0] = 0.;
109  slaterEffectiveCharge[1][0] = 0.;
110  slaterEffectiveCharge[2][0] = 0.;
111  sCoefficient[0][0] = 0.;
112  sCoefficient[1][0] = 0.;
113  sCoefficient[2][0] = 0.;
114 
115  hydrogen = hydrogenDef->GetParticleName();
116  lowEnergyLimit[hydrogen] = 10. * eV;
117  highEnergyLimit[hydrogen] = 500. * keV;
118 
119  kineticEnergyCorrection[0] = 1.;
120  slaterEffectiveCharge[0][0] = 0.;
121  slaterEffectiveCharge[1][0] = 0.;
122  slaterEffectiveCharge[2][0] = 0.;
123  sCoefficient[0][0] = 0.;
124  sCoefficient[1][0] = 0.;
125  sCoefficient[2][0] = 0.;
126 
127  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
128  lowEnergyLimit[alphaPlusPlus] = 1. * keV;
129  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
130 
131  kineticEnergyCorrection[1] = 0.9382723/3.727417;
132  slaterEffectiveCharge[0][1]=0.;
133  slaterEffectiveCharge[1][1]=0.;
134  slaterEffectiveCharge[2][1]=0.;
135  sCoefficient[0][1]=0.;
136  sCoefficient[1][1]=0.;
137  sCoefficient[2][1]=0.;
138 
139  alphaPlus = alphaPlusDef->GetParticleName();
140  lowEnergyLimit[alphaPlus] = 1. * keV;
141  highEnergyLimit[alphaPlus] = 400. * MeV;
142 
143  kineticEnergyCorrection[2] = 0.9382723/3.727417;
144  slaterEffectiveCharge[0][2]=2.0;
145 
146  // Following values provided by M. Dingfelder
147  slaterEffectiveCharge[1][2]=2.00;
148  slaterEffectiveCharge[2][2]=2.00;
149  //
150  sCoefficient[0][2]=0.7;
151  sCoefficient[1][2]=0.15;
152  sCoefficient[2][2]=0.15;
153 
154  helium = heliumDef->GetParticleName();
155  lowEnergyLimit[helium] = 1. * keV;
156  highEnergyLimit[helium] = 400. * MeV;
157 
158  kineticEnergyCorrection[3] = 0.9382723/3.727417;
159  slaterEffectiveCharge[0][3]=1.7;
160  slaterEffectiveCharge[1][3]=1.15;
161  slaterEffectiveCharge[2][3]=1.15;
162  sCoefficient[0][3]=0.5;
163  sCoefficient[1][3]=0.25;
164  sCoefficient[2][3]=0.25;
165 
166  //
167 
168  if (particle==protonDef)
169  {
170  SetLowEnergyLimit(lowEnergyLimit[proton]);
171  SetHighEnergyLimit(highEnergyLimit[proton]);
172  }
173 
174  if (particle==hydrogenDef)
175  {
176  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
177  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
178  }
179 
180  if (particle==alphaPlusPlusDef)
181  {
182  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
183  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
184  }
185 
186  if (particle==alphaPlusDef)
187  {
188  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
189  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
190  }
191 
192  if (particle==heliumDef)
193  {
194  SetLowEnergyLimit(lowEnergyLimit[helium]);
195  SetHighEnergyLimit(highEnergyLimit[helium]);
196  }
197 
198  //
199 
200  nLevels = waterExcitation.NumberOfLevels();
201 
202  //
203  if( verboseLevel>0 )
204  {
205  G4cout << "Miller & Green excitation model is initialized " << G4endl
206  << "Energy range: "
207  << LowEnergyLimit() / eV << " eV - "
208  << HighEnergyLimit() / keV << " keV for "
209  << particle->GetParticleName()
210  << G4endl;
211  }
212 
213  // Initialize water density pointer
214  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
215 
216  if (isInitialised) { return; }
217  fParticleChangeForGamma = GetParticleChangeForGamma();
218  isInitialised = true;
219 
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
223 
224 G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material,
225  const G4ParticleDefinition* particleDefinition,
226  G4double k,
227  G4double,
228  G4double)
229 {
230  if (verboseLevel > 3)
231  G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl;
232 
233  // Calculate total cross section for model
234 
235  G4DNAGenericIonsManager *instance;
236  instance = G4DNAGenericIonsManager::Instance();
237 
238  if (
239  particleDefinition != G4Proton::ProtonDefinition()
240  &&
241  particleDefinition != instance->GetIon("hydrogen")
242  &&
243  particleDefinition != instance->GetIon("alpha++")
244  &&
245  particleDefinition != instance->GetIon("alpha+")
246  &&
247  particleDefinition != instance->GetIon("helium")
248  )
249 
250  return 0;
251 
252  G4double lowLim = 0;
253  G4double highLim = 0;
254  G4double crossSection = 0.;
255 
256  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
257 
258  const G4String& particleName = particleDefinition->GetParticleName();
259 
260  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
261  pos1 = lowEnergyLimit.find(particleName);
262 
263  if (pos1 != lowEnergyLimit.end())
264  {
265  lowLim = pos1->second;
266  }
267 
268  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
269  pos2 = highEnergyLimit.find(particleName);
270 
271  if (pos2 != highEnergyLimit.end())
272  {
273  highLim = pos2->second;
274  }
275 
276  if (k >= lowLim && k <= highLim)
277  {
278  crossSection = Sum(k,particleDefinition);
279 
280  // add ONE or TWO electron-water excitation for alpha+ and helium
281  /*
282  if ( particleDefinition == instance->GetIon("alpha+")
283  ||
284  particleDefinition == instance->GetIon("helium")
285  )
286  {
287 
288  G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
289  excitationXS->Initialise(G4Electron::ElectronDefinition());
290 
291  G4double sigmaExcitation=0;
292  G4double tmp =0.;
293 
294  if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation =
295  excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp)
296  /material->GetAtomicNumDensityVector()[1];
297 
298  if ( particleDefinition == instance->GetIon("alpha+") )
299  crossSection = crossSection + sigmaExcitation ;
300 
301  if ( particleDefinition == instance->GetIon("helium") )
302  crossSection = crossSection + 2*sigmaExcitation ;
303 
304  delete excitationXS;
305 
306  // Alternative excitation model
307 
308  G4DNABornExcitationModel * excitationXS = new G4DNABornExcitationModel();
309  excitationXS->Initialise(G4Electron::ElectronDefinition());
310 
311  G4double sigmaExcitation=0;
312  G4double tmp=0;
313 
314  if (k*0.511/3728 > 9*eV && k*0.511/3728 < 1*MeV ) sigmaExcitation =
315  excitationXS->CrossSectionPerVolume(material,G4Electron::ElectronDefinition(),k*0.511/3728,tmp,tmp)
316  /material->GetAtomicNumDensityVector()[1];
317 
318  if ( particleDefinition == instance->GetIon("alpha+") )
319  crossSection = crossSection + sigmaExcitation ;
320 
321  if ( particleDefinition == instance->GetIon("helium") )
322  crossSection = crossSection + 2*sigmaExcitation ;
323 
324  delete excitationXS;
325 
326  }
327  */
328 
329  }
330 
331  if (verboseLevel > 2)
332  {
333  G4cout << "__________________________________" << G4endl;
334  G4cout << "G4DNAMillerGreenExcitationModel - XS INFO START" << G4endl;
335  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
336  G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
337  G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*waterDensity/(1./cm) << G4endl;
338  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
339  G4cout << "G4DNAMillerGreenExcitationModel - XS INFO END" << G4endl;
340  }
341 
342  return crossSection*waterDensity;
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346 
347 void G4DNAMillerGreenExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
348  const G4MaterialCutsCouple* /*couple*/,
349  const G4DynamicParticle* aDynamicParticle,
350  G4double,
351  G4double)
352 {
353 
354  if (verboseLevel > 3)
355  G4cout << "Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" << G4endl;
356 
357  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
358 
359  G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition());
360 
361  // Dingfelder's excitation levels
362  const G4double excitation[]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
363  G4double excitationEnergy = excitation[level];
364 
365  G4double newEnergy = 0.;
366 
367  if (!statCode) newEnergy = particleEnergy0 - excitationEnergy;
368 
369  else newEnergy = particleEnergy0;
370 
371  if (newEnergy>0)
372  {
373  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
374  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
375  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
376 
377  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
378  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
379  level, theIncomingTrack);
380 
381  }
382 
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386 
387 G4double G4DNAMillerGreenExcitationModel::GetPartialCrossSection(const G4Material*,
388  G4int level,
389  const G4ParticleDefinition* particleDefinition,
390  G4double kineticEnergy)
391 {
392  return PartialCrossSection(kineticEnergy, level, particleDefinition);
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396 
397 G4double G4DNAMillerGreenExcitationModel::PartialCrossSection(G4double k, G4int excitationLevel,
398  const G4ParticleDefinition* particleDefinition)
399 {
400  // ( ( z * aj ) ^ omegaj ) * ( t - ej ) ^ nu
401  // sigma(t) = zEff^2 * sigma0 * --------------------------------------------
402  // jj ^ ( omegaj + nu ) + t ^ ( omegaj + nu )
403  //
404  // where t is the kinetic energy corrected by Helium mass over proton mass for Helium ions
405  //
406  // zEff is:
407  // 1 for protons
408  // 2 for alpha++
409  // and 2 - c1 S_1s - c2 S_2s - c3 S_2p for alpha+ and He
410  //
411  // Dingfelder et al., RPC 59, 255-275, 2000 from Miller and Green (1973)
412  // Formula (34) and Table 2
413 
414  const G4double sigma0(1.E+8 * barn);
415  const G4double nu(1.);
416  const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
417  const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
418  const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
419 
420  // Dingfelder's excitation levels
421  const G4double Eliq[5]={ 8.17*eV, 10.13*eV, 11.31*eV, 12.91*eV, 14.50*eV};
422 
423  G4int particleTypeIndex = 0;
424  G4DNAGenericIonsManager* instance;
425  instance = G4DNAGenericIonsManager::Instance();
426 
427  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
428  if (particleDefinition == instance->GetIon("hydrogen")) particleTypeIndex=0;
429  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
430  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
431  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3;
432 
433  G4double tCorrected;
434  tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
435 
436  // SI - added protection
437  if (tCorrected < Eliq[excitationLevel]) return 0;
438  //
439 
440  G4int z = 10;
441 
442  G4double numerator;
443  numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) *
444  std::pow(tCorrected - Eliq[excitationLevel], nu);
445 
446  // H case : see S. Uehara et al. IJRB 77, 2, 139-154 (2001) - section 3.3
447 
448  if (particleDefinition == instance->GetIon("hydrogen"))
449  numerator = std::pow(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
450  std::pow(tCorrected - Eliq[excitationLevel], nu);
451 
452 
453  G4double power;
454  power = omegaj[excitationLevel] + nu;
455 
456  G4double denominator;
457  denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
458 
459  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
460 
461  zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
462  sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
463  sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
464 
465  if (particleDefinition == instance->GetIon("hydrogen")) zEff = 1.;
466 
467  G4double cross = sigma0 * zEff * zEff * numerator / denominator;
468 
469 
470  return cross;
471 }
472 
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 
475 G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle)
476 {
477  G4int i = nLevels;
478  G4double value = 0.;
479  std::deque<G4double> values;
480 
481  G4DNAGenericIonsManager *instance;
482  instance = G4DNAGenericIonsManager::Instance();
483 
484  if ( particle == instance->GetIon("alpha++") ||
485  particle == G4Proton::ProtonDefinition()||
486  particle == instance->GetIon("hydrogen") ||
487  particle == instance->GetIon("alpha+") ||
488  particle == instance->GetIon("helium")
489  )
490  {
491  while (i > 0)
492  {
493  i--;
494  G4double partial = PartialCrossSection(k,i,particle);
495  values.push_front(partial);
496  value += partial;
497  }
498 
499  value *= G4UniformRand();
500 
501  i = nLevels;
502 
503  while (i > 0)
504  {
505  i--;
506  if (values[i] > value) return i;
507  value -= values[i];
508  }
509  }
510 
511  /*
512  // add ONE or TWO electron-water excitation for alpha+ and helium
513 
514  if ( particle == instance->GetIon("alpha+")
515  ||
516  particle == instance->GetIon("helium")
517  )
518  {
519  while (i>0)
520  {
521  i--;
522 
523  G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
524  excitationXS->Initialise(G4Electron::ElectronDefinition());
525 
526  G4double sigmaExcitation=0;
527 
528  if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i);
529 
530  G4double partial = PartialCrossSection(k,i,particle);
531 
532  if (particle == instance->GetIon("alpha+")) partial = PartialCrossSection(k,i,particle) + sigmaExcitation;
533  if (particle == instance->GetIon("helium")) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation;
534 
535  values.push_front(partial);
536  value += partial;
537  delete excitationXS;
538  }
539 
540  value*=G4UniformRand();
541 
542  i=5;
543  while (i>0)
544  {
545  i--;
546 
547  if (values[i]>value) return i;
548 
549  value-=values[i];
550  }
551  }
552  */
553 
554  return 0;
555 }
556 
557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
558 
559 G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle)
560 {
561  G4double totalCrossSection = 0.;
562 
563  for (G4int i=0; i<nLevels; i++)
564  {
565  totalCrossSection += PartialCrossSection(k,i,particle);
566  }
567  return totalCrossSection;
568 }
569 
570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
571 
572 G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t,
573  G4double energyTransferred,
574  G4double _slaterEffectiveCharge,
575  G4double shellNumber)
576 {
577  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
578  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
579 
580  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
581  G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
582 
583  return value;
584 }
585 
586 
587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
588 
589 G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t,
590  G4double energyTransferred,
591  G4double _slaterEffectiveCharge,
592  G4double shellNumber)
593 {
594  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
595  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
596 
597  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
598  G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
599 
600  return value;
601 
602 }
603 
604 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
605 
606 G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t,
607  G4double energyTransferred,
608  G4double _slaterEffectiveCharge,
609  G4double shellNumber)
610 {
611  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
612  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
613 
614  G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
615  G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
616 
617  return value;
618 }
619 
620 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
621 
622 G4double G4DNAMillerGreenExcitationModel::R(G4double t,
623  G4double energyTransferred,
624  G4double _slaterEffectiveCharge,
625  G4double shellNumber)
626 {
627  // tElectron = m_electron / m_alpha * t
628  // Dingfelder, in Chattanooga 2005 proceedings, p 4
629 
630  G4double tElectron = 0.511/3728. * t;
631 
632  // The following is provided by M. Dingfelder
633  G4double H = 2.*13.60569172 * eV;
634  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (_slaterEffectiveCharge/shellNumber);
635 
636  return value;
637 }
638