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