Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNAMillerGreenExcitationModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 #include "G4Pow.hh"
 34 #include "G4Alpha.hh"
 35 
 36 static G4Pow * gpow = G4Pow::GetInstance();
 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 38 
 39 using namespace std;
 40 
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 42 
 43 G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(const G4ParticleDefinition*,
 44                                                                  const G4String& nam)
 45 :G4VEmModel(nam)
 46 {
 47   fpMolWaterDensity = nullptr;
 48 
 49   nLevels=0;
 50   kineticEnergyCorrection[0]=0.;
 51   kineticEnergyCorrection[1]=0.;
 52   kineticEnergyCorrection[2]=0.;
 53   kineticEnergyCorrection[3]=0.;
 54 
 55   verboseLevel= 0;
 56   // Verbosity scale:
 57   // 0 = nothing
 58   // 1 = warning for energy non-conservation
 59   // 2 = details of energy budget
 60   // 3 = calculation of cross sections, file openings, sampling of atoms
 61   // 4 = entering in methods
 62 
 63   if( verboseLevel>0 )
 64   {
 65     G4cout << "Miller & Green excitation model is constructed " << G4endl;
 66   }
 67   fParticleChangeForGamma = nullptr;
 68 
 69   // Selection of stationary mode
 70 
 71   statCode = false;
 72 }
 73 
 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75 
 76 G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel()
 77 = default;
 78 
 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 80 
 81 void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle,
 82                                                  const G4DataVector& /*cuts*/)
 83 {
 84 
 85   if (verboseLevel > 3)
 86     G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl;
 87 
 88   // Energy limits
 89 
 90   G4DNAGenericIonsManager *instance;
 91   instance = G4DNAGenericIonsManager::Instance();
 92   protonDef = G4Proton::ProtonDefinition();
 93   hydrogenDef = instance->GetIon("hydrogen");
 94   alphaPlusPlusDef = G4Alpha::Alpha();
 95   alphaPlusDef = instance->GetIon("alpha+");
 96   heliumDef = instance->GetIon("helium");
 97 
 98   G4String proton;
 99   G4String hydrogen;
100   G4String alphaPlusPlus;
101   G4String alphaPlus;
102   G4String helium;
103 
104   // LIMITS AND CONSTANTS
105 
106   proton = protonDef->GetParticleName();
107   lowEnergyLimit[proton] = 10. * eV;
108   highEnergyLimit[proton] = 500. * keV;
109   
110   kineticEnergyCorrection[0] = 1.;
111   slaterEffectiveCharge[0][0] = 0.;
112   slaterEffectiveCharge[1][0] = 0.;
113   slaterEffectiveCharge[2][0] = 0.;
114   sCoefficient[0][0] = 0.;
115   sCoefficient[1][0] = 0.;
116   sCoefficient[2][0] = 0.;
117 
118   hydrogen = hydrogenDef->GetParticleName();
119   lowEnergyLimit[hydrogen] = 10. * eV;
120   highEnergyLimit[hydrogen] = 500. * keV;
121   
122   kineticEnergyCorrection[0] = 1.;
123   slaterEffectiveCharge[0][0] = 0.;
124   slaterEffectiveCharge[1][0] = 0.;
125   slaterEffectiveCharge[2][0] = 0.;
126   sCoefficient[0][0] = 0.;
127   sCoefficient[1][0] = 0.;
128   sCoefficient[2][0] = 0.;
129 
130   alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
131   lowEnergyLimit[alphaPlusPlus] = 1. * keV;
132   highEnergyLimit[alphaPlusPlus] = 400. * MeV;
133 
134   kineticEnergyCorrection[1] = 0.9382723/3.727417;
135   slaterEffectiveCharge[0][1]=0.;
136   slaterEffectiveCharge[1][1]=0.;
137   slaterEffectiveCharge[2][1]=0.;
138   sCoefficient[0][1]=0.;
139   sCoefficient[1][1]=0.;
140   sCoefficient[2][1]=0.;
141 
142   alphaPlus = alphaPlusDef->GetParticleName();
143   lowEnergyLimit[alphaPlus] = 1. * keV;
144   highEnergyLimit[alphaPlus] = 400. * MeV;
145 
146   kineticEnergyCorrection[2] = 0.9382723/3.727417;
147   slaterEffectiveCharge[0][2]=2.0;
148 
149   // Following values provided by M. Dingfelder
150   slaterEffectiveCharge[1][2]=2.00;
151   slaterEffectiveCharge[2][2]=2.00;
152   //
153   sCoefficient[0][2]=0.7;
154   sCoefficient[1][2]=0.15;
155   sCoefficient[2][2]=0.15;
156 
157   helium = heliumDef->GetParticleName();
158   lowEnergyLimit[helium] = 1. * keV;
159   highEnergyLimit[helium] = 400. * MeV;
160   
161   kineticEnergyCorrection[3] = 0.9382723/3.727417;
162   slaterEffectiveCharge[0][3]=1.7;
163   slaterEffectiveCharge[1][3]=1.15;
164   slaterEffectiveCharge[2][3]=1.15;
165   sCoefficient[0][3]=0.5;
166   sCoefficient[1][3]=0.25;
167   sCoefficient[2][3]=0.25;
168 
169   //
170 
171   if (particle==protonDef)
172   {
173     SetLowEnergyLimit(lowEnergyLimit[proton]);
174     SetHighEnergyLimit(highEnergyLimit[proton]);
175   }
176 
177   if (particle==hydrogenDef)
178   {
179     SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
180     SetHighEnergyLimit(highEnergyLimit[hydrogen]);
181   }
182 
183   if (particle==alphaPlusPlusDef)
184   {
185     SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
186     SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
187   }
188 
189   if (particle==alphaPlusDef)
190   {
191     SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
192     SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
193   }
194 
195   if (particle==heliumDef)
196   {
197     SetLowEnergyLimit(lowEnergyLimit[helium]);
198     SetHighEnergyLimit(highEnergyLimit[helium]);
199   }
200 
201   //
202 
203   nLevels = waterExcitation.NumberOfLevels();
204 
205   //
206   if( verboseLevel>0 )
207   {
208     G4cout << "Miller & Green excitation model is initialized " << G4endl
209            << "Energy range: "
210            << LowEnergyLimit() / eV << " eV - "
211            << HighEnergyLimit() / keV << " keV for "
212            << particle->GetParticleName()
213            << G4endl;
214   }
215 
216   // Initialize water density pointer
217   fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
218 
219   if (isInitialised) { return; }
220   fParticleChangeForGamma = GetParticleChangeForGamma();
221   isInitialised = true;
222 
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 
227 G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material,
228                                                                 const G4ParticleDefinition* particleDefinition,
229                                                                 G4double k,
230                                                                 G4double,
231                                                                 G4double)
232 {
233   if (verboseLevel > 3)
234     G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl;
235 
236   // Calculate total cross section for model
237 
238   if (
239        particleDefinition != protonDef
240        &&
241        particleDefinition != hydrogenDef
242        &&
243        particleDefinition != alphaPlusPlusDef
244        &&
245        particleDefinition != alphaPlusDef
246        &&
247        particleDefinition != heliumDef
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 == alphaPlusDef
283            ||
284            particleDefinition == heliumDef
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 == alphaPlusDef )
299         crossSection = crossSection +  sigmaExcitation ;
300 
301       if ( particleDefinition == heliumDef )
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 == alphaPlusDef )
319         crossSection = crossSection +  sigmaExcitation ;
320 
321       if ( particleDefinition == heliumDef )
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  
425   if (particleDefinition == protonDef) particleTypeIndex=0;
426   if (particleDefinition == hydrogenDef) particleTypeIndex=0;
427   if (particleDefinition == alphaPlusPlusDef) particleTypeIndex=1;
428   if (particleDefinition == alphaPlusDef) particleTypeIndex=2;
429   if (particleDefinition == heliumDef) particleTypeIndex=3;
430 
431   G4double tCorrected;
432   tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
433 
434   // SI - added protection
435   if (tCorrected < Eliq[excitationLevel]) return 0;
436   //
437 
438   G4int z = 10;
439 
440   G4double numerator;
441   numerator = gpow->powA(z * aj[excitationLevel], omegaj[excitationLevel]) *
442               gpow->powA(tCorrected - Eliq[excitationLevel], nu);
443 
444   // H case : see S. Uehara et al. IJRB 77, 2, 139-154 (2001) - section 3.3
445 
446   if (particleDefinition == hydrogenDef)
447       numerator = gpow->powA(z * 0.75*aj[excitationLevel], omegaj[excitationLevel]) *
448                   gpow->powA(tCorrected - Eliq[excitationLevel], nu);
449 
450 
451   G4double power;
452   power = omegaj[excitationLevel] + nu;
453 
454   G4double denominator;
455   denominator = gpow->powA(jj[excitationLevel], power) + gpow->powA(tCorrected, power);
456 
457   G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
458 
459   zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, Eliq[excitationLevel], slaterEffectiveCharge[0][particleTypeIndex], 1.) +
460           sCoefficient[1][particleTypeIndex] * S_2s(k, Eliq[excitationLevel], slaterEffectiveCharge[1][particleTypeIndex], 2.) +
461           sCoefficient[2][particleTypeIndex] * S_2p(k, Eliq[excitationLevel], slaterEffectiveCharge[2][particleTypeIndex], 2.) );
462 
463   if (particleDefinition == hydrogenDef) zEff = 1.;
464 
465   G4double cross = sigma0 * zEff * zEff * numerator / denominator;
466 
467 
468   return cross;
469 }
470 
471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
472 
473 G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle)
474 {
475   G4int i = nLevels;
476   G4double value = 0.;
477   std::deque<G4double> values;
478 
479   if ( particle == alphaPlusPlusDef ||
480        particle == protonDef||
481        particle == hydrogenDef  ||
482        particle == alphaPlusDef  ||
483        particle == heliumDef
484        )
485   {
486     while (i > 0)
487     {
488       i--;
489       G4double partial = PartialCrossSection(k,i,particle);
490       values.push_front(partial);
491       value += partial;
492     }
493 
494     value *= G4UniformRand();
495 
496     i = nLevels;
497 
498     while (i > 0)
499     {
500       i--;
501       if (values[i] > value) return i;
502       value -= values[i];
503     }
504   }
505 
506   /*
507   // add ONE or TWO electron-water excitation for alpha+ and helium
508 
509   if ( particle == alphaPlusDef
510        ||
511        particle == heliumDef
512      )
513   {
514     while (i>0)
515     {
516       i--;
517 
518           G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
519           excitationXS->Initialise(G4Electron::ElectronDefinition());
520 
521       G4double sigmaExcitation=0;
522 
523       if (k*0.511/3728 > 8.23*eV && k*0.511/3728 < 10*MeV ) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i);
524 
525       G4double partial = PartialCrossSection(k,i,particle);
526 
527       if (particle == alphaPlusDef) partial = PartialCrossSection(k,i,particle) + sigmaExcitation;
528       if (particle == heliumDef) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation;
529 
530       values.push_front(partial);
531       value += partial;
532       delete excitationXS;
533     }
534 
535     value*=G4UniformRand();
536 
537     i=5;
538     while (i>0)
539     {
540       i--;
541 
542       if (values[i]>value) return i;
543 
544       value-=values[i];
545     }
546   }
547   */
548 
549   return 0;
550 }
551 
552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 
554 G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle)
555 {
556   G4double totalCrossSection = 0.;
557 
558   for (G4int i=0; i<nLevels; i++)
559   {
560     totalCrossSection += PartialCrossSection(k,i,particle);
561   }
562   return totalCrossSection;
563 }
564 
565 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
566 
567 G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t,
568                                                G4double energyTransferred,
569                                                G4double _slaterEffectiveCharge,
570                                                G4double shellNumber)
571 {
572   // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
573   // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
574 
575   G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
576   G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
577 
578   return value;
579 }
580 
581 
582 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
583 
584 G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t,
585                                                G4double energyTransferred,
586                                                G4double _slaterEffectiveCharge,
587                                                G4double shellNumber)
588 {
589   // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
590   // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
591 
592   G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
593   G4double value =  1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
594 
595   return value;
596 
597 }
598 
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
600 
601 G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t,
602                                                G4double energyTransferred,
603                                                G4double _slaterEffectiveCharge,
604                                                G4double shellNumber)
605 {
606   // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
607   // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
608 
609   G4double r = R(t, energyTransferred, _slaterEffectiveCharge, shellNumber);
610   G4double value =  1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
611 
612   return value;
613 }
614 
615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
616 
617 G4double G4DNAMillerGreenExcitationModel::R(G4double t,
618                                             G4double energyTransferred,
619                                             G4double _slaterEffectiveCharge,
620                                             G4double shellNumber)
621 {
622   // tElectron = m_electron / m_alpha * t
623   // Dingfelder, in Chattanooga 2005 proceedings, p 4
624 
625   G4double tElectron = 0.511/3728. * t;
626 
627   // The following is provided by M. Dingfelder
628   G4double H = 2.*13.60569172 * eV;
629   G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (_slaterEffectiveCharge/shellNumber);
630 
631   return value;
632 }
633 
634