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 10.7.p2)


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