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.1.1)


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