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 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNAMillerGreenExcitationModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNAMillerGreenExcitationModel.cc (Version 9.5)


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