Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNADingfelderChargeIncreaseModel.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/G4DNADingfelderChargeIncreaseModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNADingfelderChargeIncreaseModel.cc (Version 11.2.2)


  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 //                                                 26 //
 27                                                    27 
 28 #include "G4DNADingfelderChargeIncreaseModel.h     28 #include "G4DNADingfelderChargeIncreaseModel.hh"
 29 #include "G4PhysicalConstants.hh"                  29 #include "G4PhysicalConstants.hh"
 30 #include "G4SystemOfUnits.hh"                      30 #include "G4SystemOfUnits.hh"
 31 #include "G4DNAMolecularMaterial.hh"               31 #include "G4DNAMolecularMaterial.hh"
 32 #include "G4Log.hh"                                32 #include "G4Log.hh"
 33 #include "G4Pow.hh"                                33 #include "G4Pow.hh"
 34 #include "G4Alpha.hh"                              34 #include "G4Alpha.hh"
 35                                                    35 
 36 static G4Pow * gpow = G4Pow::GetInstance();        36 static G4Pow * gpow = G4Pow::GetInstance();
 37                                                    37 
 38 //....oooOO0OOooo........oooOO0OOooo........oo     38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 39                                                    39 
 40 using namespace std;                               40 using namespace std;
 41                                                    41 
 42 //....oooOO0OOooo........oooOO0OOooo........oo     42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 43                                                    43 
 44 G4DNADingfelderChargeIncreaseModel::G4DNADingf     44 G4DNADingfelderChargeIncreaseModel::G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition*,
 45                                                    45                                                                        const G4String& nam) :
 46 G4VEmModel(nam)                                    46 G4VEmModel(nam)
 47 {                                                  47 {
 48   if (verboseLevel > 0)                            48   if (verboseLevel > 0)
 49   {                                                49   {
 50     G4cout << "Dingfelder charge increase mode     50     G4cout << "Dingfelder charge increase model is constructed " << G4endl;
 51   }                                                51   }
 52 }                                                  52 }
 53                                                    53 
 54 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 55                                                    55 
 56 void G4DNADingfelderChargeIncreaseModel::Initi     56 void G4DNADingfelderChargeIncreaseModel::Initialise(const G4ParticleDefinition* particle,
 57                                                    57                                                     const G4DataVector& /*cuts*/)
 58 {                                                  58 {
 59                                                    59 
 60   if (verboseLevel > 3)                            60   if (verboseLevel > 3)
 61   {                                                61   {
 62     G4cout << "Calling G4DNADingfelderChargeIn     62     G4cout << "Calling G4DNADingfelderChargeIncreaseModel::Initialise()"
 63         << G4endl;                                 63         << G4endl;
 64   }                                                64   }
 65                                                    65 
 66   // Energy limits                                 66   // Energy limits
 67                                                    67 
 68   G4DNAGenericIonsManager *instance;               68   G4DNAGenericIonsManager *instance;
 69   instance = G4DNAGenericIonsManager::Instance     69   instance = G4DNAGenericIonsManager::Instance();
 70   hydrogenDef = instance->GetIon("hydrogen");      70   hydrogenDef = instance->GetIon("hydrogen");
 71   alphaPlusPlusDef = G4Alpha::Alpha();             71   alphaPlusPlusDef = G4Alpha::Alpha();
 72   alphaPlusDef = instance->GetIon("alpha+");       72   alphaPlusDef = instance->GetIon("alpha+");
 73   heliumDef = instance->GetIon("helium");          73   heliumDef = instance->GetIon("helium");
 74                                                    74 
 75   G4String hydrogen;                               75   G4String hydrogen;
 76   G4String alphaPlus;                              76   G4String alphaPlus;
 77   G4String helium;                                 77   G4String helium;
 78                                                    78 
 79   // Limits                                        79   // Limits
 80                                                    80 
 81   hydrogen = hydrogenDef->GetParticleName();       81   hydrogen = hydrogenDef->GetParticleName();
 82   lowEnergyLimit[hydrogen] = 100. * eV;            82   lowEnergyLimit[hydrogen] = 100. * eV;
 83   highEnergyLimit[hydrogen] = 100. * MeV;          83   highEnergyLimit[hydrogen] = 100. * MeV;
 84                                                    84 
 85   alphaPlus = alphaPlusDef->GetParticleName();     85   alphaPlus = alphaPlusDef->GetParticleName();
 86   lowEnergyLimit[alphaPlus] = 1. * keV;            86   lowEnergyLimit[alphaPlus] = 1. * keV;
 87   highEnergyLimit[alphaPlus] = 400. * MeV;         87   highEnergyLimit[alphaPlus] = 400. * MeV;
 88                                                    88 
 89   helium = heliumDef->GetParticleName();           89   helium = heliumDef->GetParticleName();
 90   lowEnergyLimit[helium] = 1. * keV;               90   lowEnergyLimit[helium] = 1. * keV;
 91   highEnergyLimit[helium] = 400. * MeV;            91   highEnergyLimit[helium] = 400. * MeV;
 92                                                    92 
 93   //                                               93   //
 94                                                    94 
 95   if (particle==hydrogenDef)                       95   if (particle==hydrogenDef)
 96   {                                                96   {
 97     SetLowEnergyLimit(lowEnergyLimit[hydrogen]     97     SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
 98     SetHighEnergyLimit(highEnergyLimit[hydroge     98     SetHighEnergyLimit(highEnergyLimit[hydrogen]);
 99   }                                                99   }
100                                                   100 
101   if (particle==alphaPlusDef)                     101   if (particle==alphaPlusDef)
102   {                                               102   {
103     SetLowEnergyLimit(lowEnergyLimit[alphaPlus    103     SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
104     SetHighEnergyLimit(highEnergyLimit[alphaPl    104     SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
105   }                                               105   }
106                                                   106 
107   if (particle==heliumDef)                        107   if (particle==heliumDef)
108   {                                               108   {
109     SetLowEnergyLimit(lowEnergyLimit[helium]);    109     SetLowEnergyLimit(lowEnergyLimit[helium]);
110     SetHighEnergyLimit(highEnergyLimit[helium]    110     SetHighEnergyLimit(highEnergyLimit[helium]);
111   }                                               111   }
112                                                   112 
113   // Final state                                  113   // Final state
114                                                   114 
115   //ALPHA+                                        115   //ALPHA+
116                                                   116 
117   f0[0][0]=1.;                                    117   f0[0][0]=1.;
118   a0[0][0]=2.25;                                  118   a0[0][0]=2.25;
119   a1[0][0]=-0.75;                                 119   a1[0][0]=-0.75;
120   b0[0][0]=-32.10;                                120   b0[0][0]=-32.10;
121   c0[0][0]=0.600;                                 121   c0[0][0]=0.600;
122   d0[0][0]=2.40;                                  122   d0[0][0]=2.40;
123   x0[0][0]=4.60;                                  123   x0[0][0]=4.60;
124                                                   124 
125   x1[0][0]=-1.;                                   125   x1[0][0]=-1.;
126   b1[0][0]=-1.;                                   126   b1[0][0]=-1.;
127                                                   127 
128   numberOfPartialCrossSections[0]=1;              128   numberOfPartialCrossSections[0]=1;
129                                                   129 
130   //HELIUM                                        130   //HELIUM
131                                                   131 
132   f0[0][1]=1.;                                    132   f0[0][1]=1.;
133   a0[0][1]=2.25;                                  133   a0[0][1]=2.25;
134   a1[0][1]=-0.75;                                 134   a1[0][1]=-0.75;
135   b0[0][1]=-30.93;                                135   b0[0][1]=-30.93;
136   c0[0][1]=0.590;                                 136   c0[0][1]=0.590;
137   d0[0][1]=2.35;                                  137   d0[0][1]=2.35;
138   x0[0][1]=4.29;                                  138   x0[0][1]=4.29;
139                                                   139 
140   f0[1][1]=1.;                                    140   f0[1][1]=1.;
141   a0[1][1]=2.25;                                  141   a0[1][1]=2.25;
142   a1[1][1]=-0.75;                                 142   a1[1][1]=-0.75;
143   b0[1][1]=-32.61;                                143   b0[1][1]=-32.61;
144   c0[1][1]=0.435;                                 144   c0[1][1]=0.435;
145   d0[1][1]=2.70;                                  145   d0[1][1]=2.70;
146   x0[1][1]=4.45;                                  146   x0[1][1]=4.45;
147                                                   147 
148   x1[0][1]=-1.;                                   148   x1[0][1]=-1.;
149   b1[0][1]=-1.;                                   149   b1[0][1]=-1.;
150                                                   150 
151   x1[1][1]=-1.;                                   151   x1[1][1]=-1.;
152   b1[1][1]=-1.;                                   152   b1[1][1]=-1.;
153                                                   153 
154   numberOfPartialCrossSections[1]=2;              154   numberOfPartialCrossSections[1]=2;
155                                                   155 
156   //                                              156   //
157                                                   157 
158   if( verboseLevel>0 )                            158   if( verboseLevel>0 )
159   {                                               159   {
160     G4cout << "Dingfelder charge increase mode    160     G4cout << "Dingfelder charge increase model is initialized " << G4endl
161     << "Energy range: "                           161     << "Energy range: "
162     << LowEnergyLimit() / keV << " keV - "        162     << LowEnergyLimit() / keV << " keV - "
163     << HighEnergyLimit() / MeV << " MeV for "     163     << HighEnergyLimit() / MeV << " MeV for "
164     << particle->GetParticleName()                164     << particle->GetParticleName()
165     << G4endl;                                    165     << G4endl;
166   }                                               166   }
167                                                   167 
168   // Initialize water density pointer             168   // Initialize water density pointer
169   fpMolWaterDensity = G4DNAMolecularMaterial::    169   fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
170                                                   170 
171   if (isInitialised) return;                      171   if (isInitialised) return;
172                                                   172 
173   fParticleChangeForGamma = GetParticleChangeF    173   fParticleChangeForGamma = GetParticleChangeForGamma();
174   isInitialised = true;                           174   isInitialised = true;
175 }                                                 175 }
176                                                   176 
177 //....oooOO0OOooo........oooOO0OOooo........oo    177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178                                                   178 
179 G4double G4DNADingfelderChargeIncreaseModel::C    179 G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material* material,
180                                                   180                                                                    const G4ParticleDefinition* particleDefinition,
181                                                   181                                                                    G4double k,
182                                                   182                                                                    G4double,
183                                                   183                                                                    G4double)
184 {                                                 184 {
185   if (verboseLevel > 3)                           185   if (verboseLevel > 3)
186   {                                               186   {
187     G4cout                                        187     G4cout
188         << "Calling CrossSectionPerVolume() of    188         << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel"
189         << G4endl;                                189         << G4endl;
190   }                                               190   }
191                                                   191 
192   // Calculate total cross section for model      192   // Calculate total cross section for model
193                                                   193 
194   if (                                            194   if (
195       particleDefinition != hydrogenDef           195       particleDefinition != hydrogenDef
196       &&                                          196       &&
197       particleDefinition != alphaPlusDef          197       particleDefinition != alphaPlusDef
198       &&                                          198       &&
199       particleDefinition != heliumDef             199       particleDefinition != heliumDef
200   )                                               200   )
201                                                   201 
202   return 0;                                       202   return 0;
203                                                   203 
204   G4double lowLim = 0;                            204   G4double lowLim = 0;
205   G4double highLim = 0;                           205   G4double highLim = 0;
206   G4double totalCrossSection = 0.;                206   G4double totalCrossSection = 0.;
207                                                   207 
208   G4double waterDensity = (*fpMolWaterDensity)    208   G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
209                                                   209 
210   const G4String& particleName = particleDefin    210   const G4String& particleName = particleDefinition->GetParticleName();
211                                                   211 
212   std::map< G4String,G4double,std::less<G4Stri    212   std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
213   pos1 = lowEnergyLimit.find(particleName);       213   pos1 = lowEnergyLimit.find(particleName);
214                                                   214 
215   if (pos1 != lowEnergyLimit.end())               215   if (pos1 != lowEnergyLimit.end())
216   {                                               216   {
217     lowLim = pos1->second;                        217     lowLim = pos1->second;
218   }                                               218   }
219                                                   219 
220   std::map< G4String,G4double,std::less<G4Stri    220   std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
221   pos2 = highEnergyLimit.find(particleName);      221   pos2 = highEnergyLimit.find(particleName);
222                                                   222 
223   if (pos2 != highEnergyLimit.end())              223   if (pos2 != highEnergyLimit.end())
224   {                                               224   {
225     highLim = pos2->second;                       225     highLim = pos2->second;
226   }                                               226   }
227                                                   227 
228   if (k >= lowLim && k <= highLim)                228   if (k >= lowLim && k <= highLim)
229   {                                               229   {
230     //HYDROGEN                                    230     //HYDROGEN
231     if (particleDefinition == hydrogenDef)        231     if (particleDefinition == hydrogenDef)
232     {                                             232     {
233       const G4double aa = 2.835;                  233       const G4double aa = 2.835;
234       const G4double bb = 0.310;                  234       const G4double bb = 0.310;
235       const G4double cc = 2.100;                  235       const G4double cc = 2.100;
236       const G4double dd = 0.760;                  236       const G4double dd = 0.760;
237       const G4double fac = 1.0e-18;               237       const G4double fac = 1.0e-18;
238       const G4double rr = 13.606 * eV;            238       const G4double rr = 13.606 * eV;
239                                                   239 
240       G4double t = k / (proton_mass_c2/electro    240       G4double t = k / (proton_mass_c2/electron_mass_c2);
241       G4double x = t / rr;                        241       G4double x = t / rr;
242       G4double temp = 4.0 * pi * Bohr_radius/n    242       G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
243       G4double sigmal = temp * cc * (gpow->pow    243       G4double sigmal = temp * cc * (gpow->powA(x,dd));
244       G4double sigmah = temp * (aa * G4Log(1.0    244       G4double sigmah = temp * (aa * G4Log(1.0 + x) + bb) / x;
245       totalCrossSection = 1.0/(1.0/sigmal + 1.    245       totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
246     }                                             246     }
247     else                                          247     else
248     {                                             248     {
249       totalCrossSection = Sum(k,particleDefini    249       totalCrossSection = Sum(k,particleDefinition);
250     }                                             250     }
251   }                                               251   }
252                                                   252 
253   if (verboseLevel > 2)                           253   if (verboseLevel > 2)
254   {                                               254   {
255     G4cout << "_______________________________    255     G4cout << "__________________________________" << G4endl;
256     G4cout << "G4DNADingfelderChargeIncreaseMo    256     G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl;
257     G4cout << "Kinetic energy(eV)=" << k/eV <<    257     G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl;
258     G4cout << "Cross section per water molecul    258     G4cout << "Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl;
259     G4cout << "Cross section per water molecul    259     G4cout << "Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl;
260     //  G4cout << " - Cross section per water     260     //  G4cout << " - Cross section per water molecule (cm^-1)=" 
261     //  << sigma*material->GetAtomicNumDensity    261     //  << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
262     G4cout << "G4DNADingfelderChargeIncreaseMo    262     G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl;
263   }                                               263   }
264                                                   264 
265   return totalCrossSection*waterDensity;          265   return totalCrossSection*waterDensity;
266   // return totalCrossSection*material->GetAto    266   // return totalCrossSection*material->GetAtomicNumDensityVector()[1];
267                                                   267 
268 }                                                 268 }
269                                                   269 
270 //....oooOO0OOooo........oooOO0OOooo........oo    270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
271                                                   271 
272 void G4DNADingfelderChargeIncreaseModel::Sampl    272 void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector<
273                                                   273                                                                G4DynamicParticle*>* fvect,
274                                                   274                                                            const G4MaterialCutsCouple* /*couple*/,
275                                                   275                                                            const G4DynamicParticle* aDynamicParticle,
276                                                   276                                                            G4double,
277                                                   277                                                            G4double)
278 {                                                 278 {
279   if (verboseLevel > 3)                           279   if (verboseLevel > 3)
280   {                                               280   {
281     G4cout                                        281     G4cout
282         << "Calling SampleSecondaries() of G4D    282         << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel"
283         << G4endl;                                283         << G4endl;
284   }                                               284   }
285                                                   285   
286   if (!statCode) fParticleChangeForGamma->Prop    286   if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.);
287                                                   287   
288   G4ParticleDefinition* definition = aDynamicP    288   G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
289                                                   289 
290   G4double particleMass = definition->GetPDGMa    290   G4double particleMass = definition->GetPDGMass();
291                                                   291 
292   G4double inK = aDynamicParticle->GetKineticE    292   G4double inK = aDynamicParticle->GetKineticEnergy();
293                                                   293 
294   G4int finalStateIndex = RandomSelect(inK,def    294   G4int finalStateIndex = RandomSelect(inK,definition);
295                                                   295 
296   G4int n = NumberOfFinalStates(definition,fin    296   G4int n = NumberOfFinalStates(definition,finalStateIndex);
297                                                   297 
298   G4double outK = 0.;                             298   G4double outK = 0.;
299                                                   299   
300   if (!statCode) outK = inK - IncomingParticle    300   if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
301                                                   301 
302   else outK = inK;                                302   else outK = inK;
303                                                   303   
304   if (statCode)                                   304   if (statCode) 
305     fParticleChangeForGamma->                     305     fParticleChangeForGamma->
306       ProposeLocalEnergyDeposit(IncomingPartic    306       ProposeLocalEnergyDeposit(IncomingParticleBindingEnergyConstant(definition,finalStateIndex));
307                                                   307 
308   fParticleChangeForGamma->ProposeTrackStatus(    308   fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
309                                                   309 
310   G4double electronK;                             310   G4double electronK;
311   if (definition == hydrogenDef) electronK = i    311   if (definition == hydrogenDef) electronK = inK*electron_mass_c2/proton_mass_c2;
312   else electronK = inK*electron_mass_c2/(parti    312   else electronK = inK*electron_mass_c2/(particleMass);
313                                                   313 
314   if (outK<0)                                     314   if (outK<0)
315   {                                               315   {
316     G4Exception("G4DNADingfelderChargeIncrease    316     G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004",
317         FatalException,"Final kinetic energy i    317         FatalException,"Final kinetic energy is negative.");
318   }                                               318   }
319                                                   319 
320   auto dp = new G4DynamicParticle(OutgoingPart    320   auto dp = new G4DynamicParticle(OutgoingParticleDefinition(definition,finalStateIndex),
321       aDynamicParticle->GetMomentumDirection()    321       aDynamicParticle->GetMomentumDirection(),
322       outK);                                      322       outK);
323                                                   323 
324   fvect->push_back(dp);                           324   fvect->push_back(dp);
325                                                   325 
326   n = n - 1;                                      326   n = n - 1;
327                                                   327 
328   while (n>0)                                     328   while (n>0)
329   {                                               329   {
330     n--;                                          330     n--;
331     fvect->push_back(new G4DynamicParticle        331     fvect->push_back(new G4DynamicParticle
332         (G4Electron::Electron(), aDynamicParti    332         (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) );
333   }                                               333   }
334 }                                                 334 }
335                                                   335 
336 //....oooOO0OOooo........oooOO0OOooo........oo    336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337                                                   337 
338 G4int G4DNADingfelderChargeIncreaseModel::Numb    338 G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
339                                                   339                                                               G4int finalStateIndex)
340                                                   340 
341 {                                                 341 {
342                                                   342 
343   if (particleDefinition == hydrogenDef)          343   if (particleDefinition == hydrogenDef)
344     return 2;                                     344     return 2;
345                                                   345 
346   if (particleDefinition == alphaPlusDef)         346   if (particleDefinition == alphaPlusDef)
347     return 2;                                     347     return 2;
348                                                   348 
349   if (particleDefinition == heliumDef)            349   if (particleDefinition == heliumDef)
350   {                                               350   {
351     if (finalStateIndex == 0)                     351     if (finalStateIndex == 0)
352       return 2;                                   352       return 2;
353     return 3;                                     353     return 3;
354   }                                               354   }
355                                                   355 
356   return 0;                                       356   return 0;
357 }                                                 357 }
358                                                   358 
359 //....oooOO0OOooo........oooOO0OOooo........oo    359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
360                                                   360 
361 G4ParticleDefinition* G4DNADingfelderChargeInc    361 G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
362                                                   362                                                                                      G4int finalStateIndex)
363 {                                                 363 {
364                                                   364 
365   if (particleDefinition == hydrogenDef)          365   if (particleDefinition == hydrogenDef)
366     return G4Proton::Proton();                    366     return G4Proton::Proton();
367                                                   367 
368   if (particleDefinition == alphaPlusDef)         368   if (particleDefinition == alphaPlusDef)
369     return alphaPlusPlusDef;                      369     return alphaPlusPlusDef;
370                                                   370 
371   if (particleDefinition == heliumDef)            371   if (particleDefinition == heliumDef)
372   {                                               372   {
373     if (finalStateIndex == 0)                     373     if (finalStateIndex == 0)
374       return alphaPlusDef;                        374       return alphaPlusDef;
375     return alphaPlusPlusDef;                      375     return alphaPlusPlusDef;
376   }                                               376   }
377                                                   377 
378   return nullptr;                                 378   return nullptr;
379 }                                                 379 }
380                                                   380 
381 //....oooOO0OOooo........oooOO0OOooo........oo    381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
382                                                   382 
383 G4double G4DNADingfelderChargeIncreaseModel::I    383 G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
384                                                   384                                                                                    G4int finalStateIndex)
385 {                                                 385 {
386                                                   386 
387   if (particleDefinition == hydrogenDef)          387   if (particleDefinition == hydrogenDef)
388     return 13.6 * eV;                             388     return 13.6 * eV;
389                                                   389 
390   if (particleDefinition == alphaPlusDef)         390   if (particleDefinition == alphaPlusDef)
391   {                                               391   {
392     // Binding energy for    He+ -> He++ + e-     392     // Binding energy for    He+ -> He++ + e-    54.509 eV
393     // Binding energy for    He  -> He+  + e-     393     // Binding energy for    He  -> He+  + e-    24.587 eV
394     return 54.509 * eV;                           394     return 54.509 * eV;
395   }                                               395   }
396                                                   396 
397   if (particleDefinition == heliumDef)            397   if (particleDefinition == heliumDef)
398   {                                               398   {
399     // Binding energy for    He+ -> He++ + e-     399     // Binding energy for    He+ -> He++ + e-    54.509 eV
400     // Binding energy for    He  -> He+  + e-     400     // Binding energy for    He  -> He+  + e-    24.587 eV
401                                                   401 
402     if (finalStateIndex == 0)                     402     if (finalStateIndex == 0)
403       return 24.587 * eV;                         403       return 24.587 * eV;
404     return (54.509 + 24.587) * eV;                404     return (54.509 + 24.587) * eV;
405   }                                               405   }
406                                                   406 
407   return 0;                                       407   return 0;
408 }                                                 408 }
409                                                   409 
410 //....oooOO0OOooo........oooOO0OOooo........oo    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411                                                   411 
412 G4double G4DNADingfelderChargeIncreaseModel::P    412 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(const G4double& k,
413                                                   413                                                                  const G4int& index,
414                                                   414                                                                  const G4ParticleDefinition* particleDefinition)
415 {                                                 415 {
416   G4int particleTypeIndex = 0;                    416   G4int particleTypeIndex = 0;
417                                                   417 
418   if (particleDefinition == alphaPlusDef)         418   if (particleDefinition == alphaPlusDef)
419     particleTypeIndex = 0;                        419     particleTypeIndex = 0;
420                                                   420 
421   if (particleDefinition == heliumDef)            421   if (particleDefinition == heliumDef)
422     particleTypeIndex = 1;                        422     particleTypeIndex = 1;
423                                                   423 
424   //                                              424   //
425   // sigma(T) = f0 10 ^ y(log10(T/eV))            425   // sigma(T) = f0 10 ^ y(log10(T/eV))
426   //                                              426   //
427   //         /  a0 x + b0                    x    427   //         /  a0 x + b0                    x < x0
428   //         |                                    428   //         |
429   // y(x) = <   a0 x + b0 - c0 (x - x0)^d0   x    429   // y(x) = <   a0 x + b0 - c0 (x - x0)^d0   x0 <= x < x1
430   //         |                                    430   //         |
431   //         \  a1 x + b1                    x    431   //         \  a1 x + b1                    x >= x1
432   //                                              432   //
433   //                                              433   //
434   // f0, a0, a1, b0, b1, c0, d0, x0, x1 are pa    434   // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
435   //                                              435   //
436   // f0 has been added to the code in order to    436   // f0 has been added to the code in order to manage partial (shell-dependent) cross sections 
437   // (if no shell dependence is present. f0=1.    437   // (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
438   //                                              438   //
439   // From Rad. Phys. and Chem. 59 (2000) 255-2    439   // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
440   // Inelastic-collision cross sections of liq    440   // Inelastic-collision cross sections of liquid water for interactions of energetic proton
441   //                                              441   //
442                                                   442 
443   if (x1[index][particleTypeIndex] < x0[index]    443   if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
444   {                                               444   {
445     //                                            445     //
446     // if x1 < x0 means that x1 and b1 will be    446     // if x1 < x0 means that x1 and b1 will be calculated with the following formula 
447     // (this piece of code is run on all alpha    447     // (this piece of code is run on all alphas and not on protons)
448     //                                            448     //
449     // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 /     449     // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
450     //                                            450     //
451     // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x    451     // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
452     //                                            452     //
453                                                   453 
454     x1[index][particleTypeIndex] = x0[index][p    454     x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
455         + gpow->powA((a0[index][particleTypeIn    455         + gpow->powA((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
456                        / (c0[index][particleTy    456                        / (c0[index][particleTypeIndex]
457                            * d0[index][particl    457                            * d0[index][particleTypeIndex]),
458                    1. / (d0[index][particleTyp    458                    1. / (d0[index][particleTypeIndex] - 1.));
459     b1[index][particleTypeIndex] = (a0[index][    459     b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
460         - a1[index][particleTypeIndex]) * x1[i    460         - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
461         + b0[index][particleTypeIndex]            461         + b0[index][particleTypeIndex]
462         - c0[index][particleTypeIndex]            462         - c0[index][particleTypeIndex]
463             * gpow->powA(x1[index][particleTyp    463             * gpow->powA(x1[index][particleTypeIndex]
464                            - x0[index][particl    464                            - x0[index][particleTypeIndex],
465                        d0[index][particleTypeI    465                        d0[index][particleTypeIndex]);
466   }                                               466   }
467                                                   467 
468   G4double x(G4Log(k / eV)/gpow->logZ(10));       468   G4double x(G4Log(k / eV)/gpow->logZ(10));
469   G4double y;                                     469   G4double y;
470                                                   470 
471   if (x < x0[index][particleTypeIndex])           471   if (x < x0[index][particleTypeIndex])
472     y = a0[index][particleTypeIndex] * x + b0[    472     y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
473   else if (x < x1[index][particleTypeIndex])      473   else if (x < x1[index][particleTypeIndex])
474     y = a0[index][particleTypeIndex] * x + b0[    474     y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
475         - c0[index][particleTypeIndex]            475         - c0[index][particleTypeIndex]
476             * gpow->powA(x - x0[index][particl    476             * gpow->powA(x - x0[index][particleTypeIndex],
477                        d0[index][particleTypeI    477                        d0[index][particleTypeIndex]);
478   else                                            478   else
479     y = a1[index][particleTypeIndex] * x + b1[    479     y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
480                                                   480 
481   return f0[index][particleTypeIndex] * gpow->    481   return f0[index][particleTypeIndex] * gpow->powA(10., y) * m * m;
482                                                   482 
483 }                                                 483 }
484                                                   484 
485 //....oooOO0OOooo........oooOO0OOooo........oo    485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
486                                                   486 
487 G4int G4DNADingfelderChargeIncreaseModel::Rand    487 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(const G4double& k,
488                                                   488                                                        const G4ParticleDefinition* particleDefinition)
489 {                                                 489 {
490   G4int particleTypeIndex = 0;                    490   G4int particleTypeIndex = 0;
491                                                   491 
492   if (particleDefinition == hydrogenDef)          492   if (particleDefinition == hydrogenDef)
493     return 0;                                     493     return 0;
494                                                   494 
495   if (particleDefinition == alphaPlusDef)         495   if (particleDefinition == alphaPlusDef)
496     particleTypeIndex = 0;                        496     particleTypeIndex = 0;
497                                                   497 
498   if (particleDefinition == heliumDef)            498   if (particleDefinition == heliumDef)
499     particleTypeIndex = 1;                        499     particleTypeIndex = 1;
500                                                   500 
501   const G4int n = numberOfPartialCrossSections    501   const G4int n = numberOfPartialCrossSections[particleTypeIndex];
502   auto values(new G4double[n]);                   502   auto values(new G4double[n]);
503   G4double value = 0;                             503   G4double value = 0;
504   G4int i = n;                                    504   G4int i = n;
505                                                   505 
506   while (i > 0)                                   506   while (i > 0)
507   {                                               507   {
508     i--;                                          508     i--;
509     values[i] = PartialCrossSection(k, i, part    509     values[i] = PartialCrossSection(k, i, particleDefinition);
510     value += values[i];                           510     value += values[i];
511   }                                               511   }
512                                                   512 
513   value *= G4UniformRand();                       513   value *= G4UniformRand();
514                                                   514 
515   i = n;                                          515   i = n;
516   while (i > 0)                                   516   while (i > 0)
517   {                                               517   {
518     i--;                                          518     i--;
519                                                   519 
520     if (values[i] > value)                        520     if (values[i] > value)
521       break;                                      521       break;
522                                                   522 
523     value -= values[i];                           523     value -= values[i];
524   }                                               524   }
525                                                   525 
526   delete[] values;                                526   delete[] values;
527                                                   527 
528   return i;                                       528   return i;
529 }                                                 529 }
530                                                   530 
531 //....oooOO0OOooo........oooOO0OOooo........oo    531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
532                                                   532 
533 G4double G4DNADingfelderChargeIncreaseModel::S    533 G4double G4DNADingfelderChargeIncreaseModel::Sum(const G4double& k,
534                                                   534                                                  const G4ParticleDefinition* particleDefinition)
535 {                                                 535 {
536   G4int particleTypeIndex = 0;                    536   G4int particleTypeIndex = 0;
537                                                   537 
538   if (particleDefinition == alphaPlusDef)         538   if (particleDefinition == alphaPlusDef)
539     particleTypeIndex = 0;                        539     particleTypeIndex = 0;
540                                                   540 
541   if (particleDefinition == heliumDef)            541   if (particleDefinition == heliumDef)
542     particleTypeIndex = 1;                        542     particleTypeIndex = 1;
543                                                   543 
544   G4double totalCrossSection = 0.;                544   G4double totalCrossSection = 0.;
545                                                   545 
546   for (G4int i = 0; i < numberOfPartialCrossSe    546   for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
547   {                                               547   {
548     totalCrossSection += PartialCrossSection(k    548     totalCrossSection += PartialCrossSection(k, i, particleDefinition);
549   }                                               549   }
550                                                   550 
551   return totalCrossSection;                       551   return totalCrossSection;
552 }                                                 552 }
553                                                   553