Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4hParametrisedLossModel.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/lowenergy/src/G4hParametrisedLossModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4hParametrisedLossModel.cc (Version 10.6.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 // -------------------------------------------     27 // -------------------------------------------------------------------
 28 //                                                 28 //
 29 // GEANT4 Class file                               29 // GEANT4 Class file
 30 //                                                 30 //
 31 //                                                 31 //
 32 // File name:     G4hParametrisedLossModel         32 // File name:     G4hParametrisedLossModel
 33 //                                                 33 //
 34 // Author:        V.Ivanchenko (Vladimir.Ivanc     34 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
 35 //                                                 35 //
 36 // Creation date: 20 July 2000                     36 // Creation date: 20 July 2000
 37 //                                                 37 //
 38 // Modifications:                                  38 // Modifications:
 39 // 20/07/2000  V.Ivanchenko First implementati     39 // 20/07/2000  V.Ivanchenko First implementation
 40 // 18/08/2000  V.Ivanchenko TRIM85 model is ad     40 // 18/08/2000  V.Ivanchenko TRIM85 model is added
 41 // 03/10/2000  V.Ivanchenko CodeWizard clean u     41 // 03/10/2000  V.Ivanchenko CodeWizard clean up
 42 // 10/05/2001  V.Ivanchenko Clean up againist      42 // 10/05/2001  V.Ivanchenko Clean up againist Linux compilation with -Wall
 43 // 30/12/2003  V.Ivanchenko SRIM2003 model is      43 // 30/12/2003  V.Ivanchenko SRIM2003 model is added
 44 // 07/05/2004  V.Ivanchenko Fix Graphite probl     44 // 07/05/2004  V.Ivanchenko Fix Graphite problem, add QAO model
 45 //                                                 45 //
 46 // Class Description:                              46 // Class Description:
 47 //                                                 47 //
 48 // Low energy protons/ions electronic stopping     48 // Low energy protons/ions electronic stopping power parametrisation
 49 //                                                 49 //
 50 // Class Description: End                          50 // Class Description: End
 51 //                                                 51 //
 52 // -------------------------------------------     52 // -------------------------------------------------------------------
 53 //                                                 53 //
 54 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 55 //....oooOO0OOooo........oooOO0OOooo........oo     55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 56                                                    56 
 57 #include "G4hParametrisedLossModel.hh"             57 #include "G4hParametrisedLossModel.hh"
 58                                                    58 
 59 #include "globals.hh"                              59 #include "globals.hh"
 60 #include "G4PhysicalConstants.hh"                  60 #include "G4PhysicalConstants.hh"
 61 #include "G4SystemOfUnits.hh"                      61 #include "G4SystemOfUnits.hh"
 62 #include "G4UnitsTable.hh"                         62 #include "G4UnitsTable.hh"
 63 #include "G4hZiegler1985p.hh"                      63 #include "G4hZiegler1985p.hh"
 64 #include "G4hICRU49p.hh"                           64 #include "G4hICRU49p.hh"
 65 #include "G4hICRU49He.hh"                          65 #include "G4hICRU49He.hh"
 66 #include "G4DynamicParticle.hh"                    66 #include "G4DynamicParticle.hh"
 67 #include "G4ParticleDefinition.hh"                 67 #include "G4ParticleDefinition.hh"
 68 #include "G4ElementVector.hh"                      68 #include "G4ElementVector.hh"
 69 #include "G4Material.hh"                           69 #include "G4Material.hh"
 70 #include "G4Exp.hh"                                70 #include "G4Exp.hh"
 71                                                    71 
 72 //....oooOO0OOooo........oooOO0OOooo........oo     72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 73                                                    73 
 74 G4hParametrisedLossModel::G4hParametrisedLossM     74 G4hParametrisedLossModel::G4hParametrisedLossModel(const G4String& name)
 75   :G4VLowEnergyModel(name), modelName(name)        75   :G4VLowEnergyModel(name), modelName(name)
 76 {                                                  76 {
 77   InitializeMe();                                  77   InitializeMe();
 78 }                                                  78 }
 79                                                    79 
 80 //....oooOO0OOooo........oooOO0OOooo........oo     80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81                                                    81 
 82 void G4hParametrisedLossModel::InitializeMe()      82 void G4hParametrisedLossModel::InitializeMe()
 83 {                                                  83 {
 84   expStopPower125 = 0.;                            84   expStopPower125 = 0.;
 85                                                    85 
 86   theZieglerFactor = eV*cm2*1.0e-15 ;              86   theZieglerFactor = eV*cm2*1.0e-15 ;
 87                                                    87 
 88   // Registration of parametrisation models        88   // Registration of parametrisation models
 89   const G4String& blank(" ");                  <<  89   G4String blank  = G4String(" ") ;
 90   const G4String& ir49p("ICRU_R49p");          <<  90   G4String ir49p  = G4String("ICRU_R49p") ;
 91   const G4String& ir49He("ICRU_R49He");        <<  91   G4String ir49He = G4String("ICRU_R49He") ;
 92   const G4String& zi85p("Ziegler1985p");       <<  92   G4String zi85p  = G4String("Ziegler1985p") ;
 93   if(zi85p == modelName) {                         93   if(zi85p == modelName) {
 94       eStopingPowerTable = new G4hZiegler1985p     94       eStopingPowerTable = new G4hZiegler1985p();
 95       highEnergyLimit = 100.0*MeV;                 95       highEnergyLimit = 100.0*MeV;
 96       lowEnergyLimit  = 1.0*keV;                   96       lowEnergyLimit  = 1.0*keV;
 97                                                    97 
 98   } else if(ir49p == modelName || blank == mod     98   } else if(ir49p == modelName || blank == modelName) {
 99       eStopingPowerTable = new G4hICRU49p();       99       eStopingPowerTable = new G4hICRU49p();
100       highEnergyLimit = 2.0*MeV;                  100       highEnergyLimit = 2.0*MeV;
101       lowEnergyLimit  = 1.0*keV;                  101       lowEnergyLimit  = 1.0*keV;
102                                                   102 
103   } else if(ir49He == modelName) {                103   } else if(ir49He == modelName) {
104       eStopingPowerTable = new G4hICRU49He();     104       eStopingPowerTable = new G4hICRU49He();
105       highEnergyLimit = 10.0*MeV/4.0;             105       highEnergyLimit = 10.0*MeV/4.0;
106       lowEnergyLimit  = 1.0*keV/4.0;              106       lowEnergyLimit  = 1.0*keV/4.0;
107                                                   107 
108   } else {                                        108   } else {
109       eStopingPowerTable = new G4hICRU49p();      109       eStopingPowerTable = new G4hICRU49p();
110       highEnergyLimit = 2.0*MeV;                  110       highEnergyLimit = 2.0*MeV;
111       lowEnergyLimit  = 1.0*keV;                  111       lowEnergyLimit  = 1.0*keV;
112       G4cout << "G4hParametrisedLossModel Warn    112       G4cout << "G4hParametrisedLossModel Warning: <" << modelName
113              << "> is unknown - default <"        113              << "> is unknown - default <"
114              << ir49p << ">" << " is used for     114              << ir49p << ">" << " is used for Electronic Stopping"
115              << G4endl;                           115              << G4endl;
116       modelName = ir49p;                          116       modelName = ir49p;
117   }                                               117   }
118       /*                                          118       /*
119       G4cout << "G4hParametrisedLossModel: the    119       G4cout << "G4hParametrisedLossModel: the model <"
120              << modelName << ">" << " is used     120              << modelName << ">" << " is used for Electronic Stopping"
121              << G4endl;                           121              << G4endl;
122 */                                                122 */
123 }                                                 123 }
124                                                   124 
125 //....oooOO0OOooo........oooOO0OOooo........oo    125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126                                                   126 
127 G4hParametrisedLossModel::~G4hParametrisedLoss    127 G4hParametrisedLossModel::~G4hParametrisedLossModel()
128 {                                                 128 {
129   delete eStopingPowerTable;                      129   delete eStopingPowerTable;
130 }                                                 130 }
131                                                   131 
132 //....oooOO0OOooo........oooOO0OOooo........oo    132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133                                                   133 
134 G4double G4hParametrisedLossModel::TheValue(co    134 G4double G4hParametrisedLossModel::TheValue(const G4DynamicParticle* particle,
135               const G4Material* material)         135               const G4Material* material)
136 {                                                 136 {
137   G4double scaledEnergy = (particle->GetKineti    137   G4double scaledEnergy = (particle->GetKineticEnergy())
138                         * proton_mass_c2/(part    138                         * proton_mass_c2/(particle->GetMass());
139   G4double factor = theZieglerFactor;             139   G4double factor = theZieglerFactor;
140   if (scaledEnergy < lowEnergyLimit) {            140   if (scaledEnergy < lowEnergyLimit) {
141     if (modelName != "QAO") factor *= std::sqr    141     if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
142     scaledEnergy = lowEnergyLimit;                142     scaledEnergy = lowEnergyLimit;
143   }                                               143   }
144   G4double eloss = StoppingPower(material,scal    144   G4double eloss = StoppingPower(material,scaledEnergy) * factor;
145                                                   145 
146   return eloss;                                   146   return eloss;
147 }                                                 147 }
148                                                   148 
149 //....oooOO0OOooo........oooOO0OOooo........oo    149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150                                                   150 
151 G4double G4hParametrisedLossModel::TheValue(co    151 G4double G4hParametrisedLossModel::TheValue(const G4ParticleDefinition* aParticle,
152               const G4Material* material,         152               const G4Material* material,
153                     G4double kineticEnergy)       153                     G4double kineticEnergy)
154 {                                                 154 {
155   G4double scaledEnergy = kineticEnergy           155   G4double scaledEnergy = kineticEnergy
156                         * proton_mass_c2/(aPar    156                         * proton_mass_c2/(aParticle->GetPDGMass());
157                                                   157 
158   G4double factor = theZieglerFactor;             158   G4double factor = theZieglerFactor;
159   if (scaledEnergy < lowEnergyLimit) {            159   if (scaledEnergy < lowEnergyLimit) {
160     if (modelName != "QAO") factor *= std::sqr    160     if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
161     scaledEnergy = lowEnergyLimit;                161     scaledEnergy = lowEnergyLimit;
162   }                                               162   }
163   G4double eloss = StoppingPower(material,scal    163   G4double eloss = StoppingPower(material,scaledEnergy) * factor;
164                                                   164 
165   return eloss;                                   165   return eloss;
166 }                                                 166 }
167                                                   167 
168 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169                                                   169 
170 G4double G4hParametrisedLossModel::LowEnergyLi    170 G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ,
171               const G4Material*) const            171               const G4Material*) const
172 {                                                 172 {
173   return lowEnergyLimit;                          173   return lowEnergyLimit;
174 }                                                 174 }
175                                                   175 
176 //....oooOO0OOooo........oooOO0OOooo........oo    176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177                                                   177 
178 G4double G4hParametrisedLossModel::HighEnergyL    178 G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ,
179                const G4Material*) const           179                const G4Material*) const
180 {                                                 180 {
181   return highEnergyLimit;                         181   return highEnergyLimit;
182 }                                                 182 }
183 //....oooOO0OOooo........oooOO0OOooo........oo    183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184                                                   184 
185 G4double G4hParametrisedLossModel::LowEnergyLi    185 G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ) const
186 {                                                 186 {
187   return lowEnergyLimit;                          187   return lowEnergyLimit;
188 }                                                 188 }
189                                                   189 
190 //....oooOO0OOooo........oooOO0OOooo........oo    190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191                                                   191 
192 G4double G4hParametrisedLossModel::HighEnergyL    192 G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ) const
193 {                                                 193 {
194   return highEnergyLimit;                         194   return highEnergyLimit;
195 }                                                 195 }
196                                                   196 
197 //....oooOO0OOooo........oooOO0OOooo........oo    197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198                                                   198 
199 G4bool G4hParametrisedLossModel::IsInCharge(co    199 G4bool G4hParametrisedLossModel::IsInCharge(const G4DynamicParticle* ,
200               const G4Material*) const            200               const G4Material*) const
201 {                                                 201 {
202   return true;                                    202   return true;
203 }                                                 203 }
204                                                   204 
205 //....oooOO0OOooo........oooOO0OOooo........oo    205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206                                                   206 
207 G4bool G4hParametrisedLossModel::IsInCharge(co    207 G4bool G4hParametrisedLossModel::IsInCharge(const G4ParticleDefinition* ,
208               const G4Material*) const            208               const G4Material*) const
209 {                                                 209 {
210   return true;                                    210   return true;
211 }                                                 211 }
212                                                   212 
213 //....oooOO0OOooo........oooOO0OOooo........oo    213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214                                                   214 
215 G4double G4hParametrisedLossModel::StoppingPow    215 G4double G4hParametrisedLossModel::StoppingPower(const G4Material* material,
216                    G4double kineticEnergy)        216                    G4double kineticEnergy)
217 {                                                 217 {
218   G4double eloss = 0.0;                           218   G4double eloss = 0.0;
219                                                   219 
220   const std::size_t numberOfElements = materia << 220   const G4int numberOfElements = material->GetNumberOfElements() ;
221   const G4double* theAtomicNumDensityVector =     221   const G4double* theAtomicNumDensityVector =
222     material->GetAtomicNumDensityVector() ;       222     material->GetAtomicNumDensityVector() ;
223                                                   223 
224                                                   224 
225   // compound material with parametrisation       225   // compound material with parametrisation
226   if( (eStopingPowerTable->HasMaterial(materia    226   if( (eStopingPowerTable->HasMaterial(material)) ) {
227                                                   227 
228     eloss = eStopingPowerTable->StoppingPower(    228     eloss = eStopingPowerTable->StoppingPower(material, kineticEnergy);
229     if ("QAO" != modelName) {                     229     if ("QAO" != modelName) {
230       eloss *=  material->GetTotNbOfAtomsPerVo    230       eloss *=  material->GetTotNbOfAtomsPerVolume();
231       if(1 < numberOfElements) {                  231       if(1 < numberOfElements) {
232         G4int nAtoms = 0;                         232         G4int nAtoms = 0;
233                                                   233 
234         const G4int* theAtomsVector = material    234         const G4int* theAtomsVector = material->GetAtomsVector();
235         for (std::size_t iel=0; iel<numberOfEl << 235         for (G4int iel=0; iel<numberOfElements; iel++) {
236           nAtoms += theAtomsVector[iel];          236           nAtoms += theAtomsVector[iel];
237         }                                         237         }
238         eloss /= nAtoms;                          238         eloss /= nAtoms;
239       }                                           239       }
240     }                                             240     }
241                                                   241 
242   // pure material                                242   // pure material
243   } else if(1 == numberOfElements) {              243   } else if(1 == numberOfElements) {
244                                                   244 
245     G4double z = material->GetZ();                245     G4double z = material->GetZ();
246     eloss = (eStopingPowerTable->ElectronicSto    246     eloss = (eStopingPowerTable->ElectronicStoppingPower(z, kineticEnergy))
247                                * (material->Ge    247                                * (material->GetTotNbOfAtomsPerVolume()) ;
248                                                   248 
249   // Experimental data exist only for kinetic     249   // Experimental data exist only for kinetic energy 125 keV
250   } else if( MolecIsInZiegler1988(material)) {    250   } else if( MolecIsInZiegler1988(material)) {
251                                                   251 
252   // Cycle over elements - calculation based o    252   // Cycle over elements - calculation based on Bragg's rule
253     G4double eloss125 = 0.0 ;                     253     G4double eloss125 = 0.0 ;
254     const G4ElementVector* theElementVector =     254     const G4ElementVector* theElementVector =
255                            material->GetElemen    255                            material->GetElementVector() ;
256                                                   256 
257                                                   257 
258     //  loop for the elements in the material     258     //  loop for the elements in the material
259     for (std::size_t i=0; i<numberOfElements;  << 259     for (G4int i=0; i<numberOfElements; i++) {
260       const G4Element* element = (*theElementV    260       const G4Element* element = (*theElementVector)[i] ;
261       G4double z = element->GetZ() ;              261       G4double z = element->GetZ() ;
262       eloss    +=(eStopingPowerTable->Electron    262       eloss    +=(eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
263                                     * theAtomi    263                                     * theAtomicNumDensityVector[i] ;
264       eloss125 +=(eStopingPowerTable->Electron    264       eloss125 +=(eStopingPowerTable->ElectronicStoppingPower(z,125.0*keV))
265                                     * theAtomi    265                                     * theAtomicNumDensityVector[i] ;
266     }                                             266     }
267                                                   267 
268     // Chemical factor is taken into account      268     // Chemical factor is taken into account
269     if (eloss125 > 0.0) {                      << 269     eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
270       eloss *= ChemicalFactor(kineticEnergy, e << 
271     }                                          << 
272                                                   270 
273   // Brugg's rule calculation                     271   // Brugg's rule calculation
274   } else {                                        272   } else {
275     const G4ElementVector* theElementVector =     273     const G4ElementVector* theElementVector =
276                            material->GetElemen    274                            material->GetElementVector() ;
277                                                   275 
278     //  loop for the elements in the material     276     //  loop for the elements in the material
279     for (std::size_t i=0; i<numberOfElements;  << 277     for (G4int i=0; i<numberOfElements; i++)
280     {                                             278     {
281       const G4Element* element = (*theElementV    279       const G4Element* element = (*theElementVector)[i] ;
282       G4double z = element->GetZ() ;              280       G4double z = element->GetZ() ;
283       eloss   += (eStopingPowerTable->Electron    281       eloss   += (eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
284                                    * theAtomic    282                                    * theAtomicNumDensityVector[i];
285     }                                             283     }
286   }                                               284   }
287   return eloss;                                   285   return eloss;
288 }                                                 286 }
289                                                   287 
290 //....oooOO0OOooo........oooOO0OOooo........oo    288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291                                                   289 
292 G4bool G4hParametrisedLossModel::MolecIsInZieg    290 G4bool G4hParametrisedLossModel::MolecIsInZiegler1988(
293                                  const G4Mater    291                                  const G4Material* material)
294 {                                                 292 {
295   // The list of molecules from                   293   // The list of molecules from
296   // J.F.Ziegler and J.M.Manoyan, The stopping    294   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
297   // Nucl. Inst. & Meth. in Phys. Res. B35 (19    295   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
298                                                   296 
299   G4String myFormula = G4String(" ") ;            297   G4String myFormula = G4String(" ") ;
300   const G4String chFormula = material->GetChem    298   const G4String chFormula = material->GetChemicalFormula() ;
301   if (myFormula == chFormula ) return false ;     299   if (myFormula == chFormula ) return false ;
302                                                   300 
303   //  There are no evidence for difference of     301   //  There are no evidence for difference of stopping power depended on
304   //  phase of the compound except for water.     302   //  phase of the compound except for water. The stopping power of the
305   //  water in gas phase can be predicted usin    303   //  water in gas phase can be predicted using Bragg's rule.
306   //                                              304   //
307   //  No chemical factor for water-gas            305   //  No chemical factor for water-gas
308                                                   306 
309   myFormula = G4String("H_2O") ;                  307   myFormula = G4String("H_2O") ;
310   const G4State theState = material->GetState(    308   const G4State theState = material->GetState() ;
311   if( theState == kStateGas && myFormula == ch    309   if( theState == kStateGas && myFormula == chFormula) return false ;
312                                                   310 
313   const std::size_t numberOfMolecula = 53 ;    << 311   const size_t numberOfMolecula = 53 ;
314                                                   312 
315   // The coffecient from Table.4 of Ziegler &     313   // The coffecient from Table.4 of Ziegler & Manoyan
316   static const G4double HeEff = 2.8735 ;          314   static const G4double HeEff = 2.8735 ;
317                                                   315 
318   static const G4String name[numberOfMolecula]    316   static const G4String name[numberOfMolecula] = {
319     "H_2O",     "C_2H_4O",    "C_3H_6O",  "C_2    317     "H_2O",     "C_2H_4O",    "C_3H_6O",  "C_2H_2",             "C_H_3OH",
320     "C_2H_5OH",  "C_3H_7OH",   "C_3H_4",   "NH    318     "C_2H_5OH",  "C_3H_7OH",   "C_3H_4",   "NH_3",               "C_14H_10",
321     "C_6H_6",    "C_4H_10",    "C_4H_6",   "C_    319     "C_6H_6",    "C_4H_10",    "C_4H_6",   "C_4H_8O",            "CCl_4",
322     "CF_4",      "C_6H_8",     "C_6H_12",  "C_    320     "CF_4",      "C_6H_8",     "C_6H_12",  "C_6H_10O",           "C_6H_10",
323     "C_8H_16",   "C_5H_10",    "C_5H_8",   "C_    321     "C_8H_16",   "C_5H_10",    "C_5H_8",   "C_3H_6-Cyclopropane","C_2H_4F_2",
324     "C_2H_2F_2", "C_4H_8O_2",  "C_2H_6",   "C_    322     "C_2H_2F_2", "C_4H_8O_2",  "C_2H_6",   "C_2F_6",             "C_2H_6O",
325     "C_3H_6O",   "C_4H_10O",   "C_2H_4",   "C_    323     "C_3H_6O",   "C_4H_10O",   "C_2H_4",   "C_2H_4O",            "C_2H_4S",
326     "SH_2",      "CH_4",       "CCLF_3",   "CC    324     "SH_2",      "CH_4",       "CCLF_3",   "CCl_2F_2",           "CHCl_2F",
327     "(CH_3)_2S", "N_2O",       "C_5H_10O", "C_    325     "(CH_3)_2S", "N_2O",       "C_5H_10O", "C_8H_6",          "(CH_2)_N",
328     "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8",   "C_    326     "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8",   "C_3H_6-Propylene",   "C_3H_6O",
329     "C_3H_6S",   "C_4H_4S",    "C_7H_8"           327     "C_3H_6S",   "C_4H_4S",    "C_7H_8"
330   };                                              328   };
331                                                   329 
332   static const G4double expStopping[numberOfMo    330   static const G4double expStopping[numberOfMolecula] = {
333      66.1,  190.4, 258.7,  42.2, 141.5,           331      66.1,  190.4, 258.7,  42.2, 141.5,
334     210.9,  279.6, 198.8,  31.0, 267.5,           332     210.9,  279.6, 198.8,  31.0, 267.5,
335     122.8,  311.4, 260.3, 328.9, 391.3,           333     122.8,  311.4, 260.3, 328.9, 391.3,
336     206.6,  374.0, 422.0, 432.0, 398.0,           334     206.6,  374.0, 422.0, 432.0, 398.0,
337     554.0,  353.0, 326.0,  74.6, 220.5,           335     554.0,  353.0, 326.0,  74.6, 220.5,
338     197.4,  362.0, 170.0, 330.5, 211.3,           336     197.4,  362.0, 170.0, 330.5, 211.3,
339     262.3,  349.6,  51.3, 187.0, 236.9,           337     262.3,  349.6,  51.3, 187.0, 236.9,
340     121.9,   35.8, 247.0, 292.6, 268.0,           338     121.9,   35.8, 247.0, 292.6, 268.0,
341     262.3,   49.0, 398.9, 444.0,  22.91,          339     262.3,   49.0, 398.9, 444.0,  22.91,
342      68.0,  155.0,  84.0,  74.2, 254.7,           340      68.0,  155.0,  84.0,  74.2, 254.7,
343     306.8,  324.4, 420.0                          341     306.8,  324.4, 420.0
344   } ;                                             342   } ;
345                                                   343 
346   static const G4double expCharge[numberOfMole    344   static const G4double expCharge[numberOfMolecula] = {
347     HeEff, HeEff, HeEff,   1.0, HeEff,            345     HeEff, HeEff, HeEff,   1.0, HeEff,
348     HeEff, HeEff, HeEff,   1.0,   1.0,            346     HeEff, HeEff, HeEff,   1.0,   1.0,
349       1.0, HeEff, HeEff, HeEff, HeEff,            347       1.0, HeEff, HeEff, HeEff, HeEff,
350     HeEff, HeEff, HeEff, HeEff, HeEff,            348     HeEff, HeEff, HeEff, HeEff, HeEff,
351     HeEff, HeEff, HeEff,   1.0, HeEff,            349     HeEff, HeEff, HeEff,   1.0, HeEff,
352     HeEff, HeEff, HeEff, HeEff, HeEff,            350     HeEff, HeEff, HeEff, HeEff, HeEff,
353     HeEff, HeEff,   1.0, HeEff, HeEff,            351     HeEff, HeEff,   1.0, HeEff, HeEff,
354     HeEff,   1.0, HeEff, HeEff, HeEff,            352     HeEff,   1.0, HeEff, HeEff, HeEff,
355     HeEff,   1.0, HeEff, HeEff,   1.0,            353     HeEff,   1.0, HeEff, HeEff,   1.0,
356       1.0,   1.0,   1.0,   1.0, HeEff,            354       1.0,   1.0,   1.0,   1.0, HeEff,
357     HeEff, HeEff, HeEff                           355     HeEff, HeEff, HeEff
358   } ;                                             356   } ;
359                                                   357 
360   static const G4double numberOfAtomsPerMolecu    358   static const G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
361     3.0,  7.0, 10.0,  4.0,  6.0,                  359     3.0,  7.0, 10.0,  4.0,  6.0,
362     9.0, 12.0,  7.0,  4.0, 24.0,                  360     9.0, 12.0,  7.0,  4.0, 24.0,
363     12.0, 14.0, 10.0, 13.0,  5.0,                 361     12.0, 14.0, 10.0, 13.0,  5.0,
364     5.0, 14.0, 18.0, 17.0, 17.0,                  362     5.0, 14.0, 18.0, 17.0, 17.0,
365     24.0, 15.0, 13.0,  9.0,  8.0,                 363     24.0, 15.0, 13.0,  9.0,  8.0,
366     6.0, 14.0,  8.0,  8.0,  9.0,                  364     6.0, 14.0,  8.0,  8.0,  9.0,
367     10.0, 15.0,  6.0,  7.0,  7.0,                 365     10.0, 15.0,  6.0,  7.0,  7.0,
368     3.0,  5.0,  5.0,  5.0,  5.0,                  366     3.0,  5.0,  5.0,  5.0,  5.0,
369     9.0,  3.0, 16.0, 14.0,  3.0,                  367     9.0,  3.0, 16.0, 14.0,  3.0,
370     9.0, 16.0, 11.0,  9.0, 10.0,                  368     9.0, 16.0, 11.0,  9.0, 10.0,
371     10.0,  9.0, 15.0                              369     10.0,  9.0, 15.0
372   } ;                                             370   } ;
373                                                   371 
374   // Search for the compaund in the table         372   // Search for the compaund in the table
375   for (std::size_t i=0; i<numberOfMolecula; ++ << 373   for (size_t i=0; i<numberOfMolecula; i++)
376     {                                             374     {
377       if(chFormula == name[i]) {                  375       if(chFormula == name[i]) {
378         G4double exp125 = expStopping[i] *        376         G4double exp125 = expStopping[i] *
379                     (material->GetTotNbOfAtoms    377                     (material->GetTotNbOfAtomsPerVolume()) /
380                     (expCharge[i] * numberOfAt    378                     (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
381         SetExpStopPower125(exp125) ;              379         SetExpStopPower125(exp125) ;
382         return true ;                             380         return true ;
383       }                                           381       }
384     }                                             382     }
385                                                   383 
386   return false ;                                  384   return false ;
387 }                                                 385 }
388                                                   386 
389 //....oooOO0OOooo........oooOO0OOooo........oo    387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390                                                   388 
391 G4double G4hParametrisedLossModel::ChemicalFac    389 G4double G4hParametrisedLossModel::ChemicalFactor(
392                             G4double kineticEn    390                             G4double kineticEnergy, G4double eloss125) const
393 {                                                 391 {
394   // Approximation of Chemical Factor accordin    392   // Approximation of Chemical Factor according to
395   // J.F.Ziegler and J.M.Manoyan, The stopping    393   // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
396   // Nucl. Inst. & Meth. in Phys. Res. B35 (19    394   // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
397                                                   395 
398   G4double gamma    = 1.0 + kineticEnergy/prot    396   G4double gamma    = 1.0 + kineticEnergy/proton_mass_c2 ;
399   G4double gamma25  = 1.0 + 25.0*keV /proton_m    397   G4double gamma25  = 1.0 + 25.0*keV /proton_mass_c2 ;
400   G4double gamma125 = 1.0 + 125.0*keV/proton_m    398   G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
401   G4double beta     = std::sqrt(1.0 - 1.0/(gam    399   G4double beta     = std::sqrt(1.0 - 1.0/(gamma*gamma)) ;
402   G4double beta25   = std::sqrt(1.0 - 1.0/(gam    400   G4double beta25   = std::sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
403   G4double beta125  = std::sqrt(1.0 - 1.0/(gam    401   G4double beta125  = std::sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
404                                                   402 
405   G4double factor = 1.0 + (expStopPower125/elo    403   G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
406                    (1.0 + G4Exp( 1.48 * ( beta    404                    (1.0 + G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
407                    (1.0 + G4Exp( 1.48 * ( beta    405                    (1.0 + G4Exp( 1.48 * ( beta/beta25    - 7.0 ) ) ) ;
408                                                   406 
409   return factor ;                                 407   return factor ;
410 }                                                 408 }
411                                                   409 
412 //....oooOO0OOooo........oooOO0OOooo........oo    410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413                                                   411