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.6.p1)


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