| Geant4 Cross Reference |
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" <<
71 70
72 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 72
74 G4hParametrisedLossModel::G4hParametrisedLossM 73 G4hParametrisedLossModel::G4hParametrisedLossModel(const G4String& name)
75 :G4VLowEnergyModel(name), modelName(name) 74 :G4VLowEnergyModel(name), modelName(name)
76 { 75 {
77 InitializeMe(); 76 InitializeMe();
78 } 77 }
79 78
80 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 80
82 void G4hParametrisedLossModel::InitializeMe() 81 void G4hParametrisedLossModel::InitializeMe()
83 { 82 {
84 expStopPower125 = 0.; 83 expStopPower125 = 0.;
85 84
86 theZieglerFactor = eV*cm2*1.0e-15 ; 85 theZieglerFactor = eV*cm2*1.0e-15 ;
87 86
88 // Registration of parametrisation models 87 // Registration of parametrisation models
89 const G4String& blank(" "); << 88 G4String blank = G4String(" ") ;
90 const G4String& ir49p("ICRU_R49p"); << 89 G4String ir49p = G4String("ICRU_R49p") ;
91 const G4String& ir49He("ICRU_R49He"); << 90 G4String ir49He = G4String("ICRU_R49He") ;
92 const G4String& zi85p("Ziegler1985p"); << 91 G4String zi85p = G4String("Ziegler1985p") ;
93 if(zi85p == modelName) { 92 if(zi85p == modelName) {
94 eStopingPowerTable = new G4hZiegler1985p 93 eStopingPowerTable = new G4hZiegler1985p();
95 highEnergyLimit = 100.0*MeV; 94 highEnergyLimit = 100.0*MeV;
96 lowEnergyLimit = 1.0*keV; 95 lowEnergyLimit = 1.0*keV;
97 96
98 } else if(ir49p == modelName || blank == mod 97 } else if(ir49p == modelName || blank == modelName) {
99 eStopingPowerTable = new G4hICRU49p(); 98 eStopingPowerTable = new G4hICRU49p();
100 highEnergyLimit = 2.0*MeV; 99 highEnergyLimit = 2.0*MeV;
101 lowEnergyLimit = 1.0*keV; 100 lowEnergyLimit = 1.0*keV;
102 101
103 } else if(ir49He == modelName) { 102 } else if(ir49He == modelName) {
104 eStopingPowerTable = new G4hICRU49He(); 103 eStopingPowerTable = new G4hICRU49He();
105 highEnergyLimit = 10.0*MeV/4.0; 104 highEnergyLimit = 10.0*MeV/4.0;
106 lowEnergyLimit = 1.0*keV/4.0; 105 lowEnergyLimit = 1.0*keV/4.0;
107 106
108 } else { 107 } else {
109 eStopingPowerTable = new G4hICRU49p(); 108 eStopingPowerTable = new G4hICRU49p();
110 highEnergyLimit = 2.0*MeV; 109 highEnergyLimit = 2.0*MeV;
111 lowEnergyLimit = 1.0*keV; 110 lowEnergyLimit = 1.0*keV;
112 G4cout << "G4hParametrisedLossModel Warn << 111 G4cout << "G4hParametrisedLossModel Warning: <" << modelName
113 << "> is unknown - default <" 112 << "> is unknown - default <"
114 << ir49p << ">" << " is used for 113 << ir49p << ">" << " is used for Electronic Stopping"
115 << G4endl; 114 << G4endl;
116 modelName = ir49p; 115 modelName = ir49p;
117 } 116 }
118 /* 117 /*
119 G4cout << "G4hParametrisedLossModel: the 118 G4cout << "G4hParametrisedLossModel: the model <"
120 << modelName << ">" << " is used 119 << modelName << ">" << " is used for Electronic Stopping"
121 << G4endl; 120 << G4endl;
122 */ 121 */
123 } 122 }
124 123
125 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126 125
127 G4hParametrisedLossModel::~G4hParametrisedLoss 126 G4hParametrisedLossModel::~G4hParametrisedLossModel()
128 { 127 {
129 delete eStopingPowerTable; 128 delete eStopingPowerTable;
130 } 129 }
131 130
132 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 132
134 G4double G4hParametrisedLossModel::TheValue(co 133 G4double G4hParametrisedLossModel::TheValue(const G4DynamicParticle* particle,
135 const G4Material* material) 134 const G4Material* material)
136 { 135 {
137 G4double scaledEnergy = (particle->GetKineti 136 G4double scaledEnergy = (particle->GetKineticEnergy())
138 * proton_mass_c2/(part 137 * proton_mass_c2/(particle->GetMass());
139 G4double factor = theZieglerFactor; 138 G4double factor = theZieglerFactor;
140 if (scaledEnergy < lowEnergyLimit) { 139 if (scaledEnergy < lowEnergyLimit) {
141 if (modelName != "QAO") factor *= std::sqr 140 if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
142 scaledEnergy = lowEnergyLimit; 141 scaledEnergy = lowEnergyLimit;
143 } 142 }
144 G4double eloss = StoppingPower(material,scal 143 G4double eloss = StoppingPower(material,scaledEnergy) * factor;
145 144
146 return eloss; 145 return eloss;
147 } 146 }
148 147
149 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 149
151 G4double G4hParametrisedLossModel::TheValue(co 150 G4double G4hParametrisedLossModel::TheValue(const G4ParticleDefinition* aParticle,
152 const G4Material* material, 151 const G4Material* material,
153 G4double kineticEnergy) 152 G4double kineticEnergy)
154 { 153 {
155 G4double scaledEnergy = kineticEnergy 154 G4double scaledEnergy = kineticEnergy
156 * proton_mass_c2/(aPar 155 * proton_mass_c2/(aParticle->GetPDGMass());
157 156
158 G4double factor = theZieglerFactor; 157 G4double factor = theZieglerFactor;
159 if (scaledEnergy < lowEnergyLimit) { 158 if (scaledEnergy < lowEnergyLimit) {
160 if (modelName != "QAO") factor *= std::sqr 159 if (modelName != "QAO") factor *= std::sqrt(scaledEnergy/lowEnergyLimit);
161 scaledEnergy = lowEnergyLimit; 160 scaledEnergy = lowEnergyLimit;
162 } 161 }
163 G4double eloss = StoppingPower(material,scal 162 G4double eloss = StoppingPower(material,scaledEnergy) * factor;
164 163
165 return eloss; 164 return eloss;
166 } 165 }
167 166
168 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 168
170 G4double G4hParametrisedLossModel::LowEnergyLi 169 G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ,
171 const G4Material*) const 170 const G4Material*) const
172 { 171 {
173 return lowEnergyLimit; 172 return lowEnergyLimit;
174 } 173 }
175 174
176 //....oooOO0OOooo........oooOO0OOooo........oo 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 176
178 G4double G4hParametrisedLossModel::HighEnergyL 177 G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ,
179 const G4Material*) const 178 const G4Material*) const
180 { 179 {
181 return highEnergyLimit; 180 return highEnergyLimit;
182 } 181 }
183 //....oooOO0OOooo........oooOO0OOooo........oo 182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184 183
185 G4double G4hParametrisedLossModel::LowEnergyLi 184 G4double G4hParametrisedLossModel::LowEnergyLimit(const G4ParticleDefinition* ) const
186 { 185 {
187 return lowEnergyLimit; 186 return lowEnergyLimit;
188 } 187 }
189 188
190 //....oooOO0OOooo........oooOO0OOooo........oo 189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191 190
192 G4double G4hParametrisedLossModel::HighEnergyL 191 G4double G4hParametrisedLossModel::HighEnergyLimit(const G4ParticleDefinition* ) const
193 { 192 {
194 return highEnergyLimit; 193 return highEnergyLimit;
195 } 194 }
196 195
197 //....oooOO0OOooo........oooOO0OOooo........oo 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 197
199 G4bool G4hParametrisedLossModel::IsInCharge(co 198 G4bool G4hParametrisedLossModel::IsInCharge(const G4DynamicParticle* ,
200 const G4Material*) const 199 const G4Material*) const
201 { 200 {
202 return true; 201 return true;
203 } 202 }
204 203
205 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 205
207 G4bool G4hParametrisedLossModel::IsInCharge(co 206 G4bool G4hParametrisedLossModel::IsInCharge(const G4ParticleDefinition* ,
208 const G4Material*) const 207 const G4Material*) const
209 { 208 {
210 return true; 209 return true;
211 } 210 }
212 211
213 //....oooOO0OOooo........oooOO0OOooo........oo 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
214 213
215 G4double G4hParametrisedLossModel::StoppingPow 214 G4double G4hParametrisedLossModel::StoppingPower(const G4Material* material,
216 G4double kineticEnergy) 215 G4double kineticEnergy)
217 { 216 {
218 G4double eloss = 0.0; 217 G4double eloss = 0.0;
219 218
220 const std::size_t numberOfElements = materia << 219 const G4int numberOfElements = material->GetNumberOfElements() ;
221 const G4double* theAtomicNumDensityVector = 220 const G4double* theAtomicNumDensityVector =
222 material->GetAtomicNumDensityVector() ; 221 material->GetAtomicNumDensityVector() ;
223 222
224 223
225 // compound material with parametrisation 224 // compound material with parametrisation
226 if( (eStopingPowerTable->HasMaterial(materia 225 if( (eStopingPowerTable->HasMaterial(material)) ) {
227 226
228 eloss = eStopingPowerTable->StoppingPower( 227 eloss = eStopingPowerTable->StoppingPower(material, kineticEnergy);
229 if ("QAO" != modelName) { 228 if ("QAO" != modelName) {
230 eloss *= material->GetTotNbOfAtomsPerVo 229 eloss *= material->GetTotNbOfAtomsPerVolume();
231 if(1 < numberOfElements) { 230 if(1 < numberOfElements) {
232 G4int nAtoms = 0; 231 G4int nAtoms = 0;
233 232
234 const G4int* theAtomsVector = material 233 const G4int* theAtomsVector = material->GetAtomsVector();
235 for (std::size_t iel=0; iel<numberOfEl << 234 for (G4int iel=0; iel<numberOfElements; iel++) {
236 nAtoms += theAtomsVector[iel]; 235 nAtoms += theAtomsVector[iel];
237 } 236 }
238 eloss /= nAtoms; 237 eloss /= nAtoms;
239 } 238 }
240 } 239 }
241 240
242 // pure material 241 // pure material
243 } else if(1 == numberOfElements) { 242 } else if(1 == numberOfElements) {
244 243
245 G4double z = material->GetZ(); 244 G4double z = material->GetZ();
246 eloss = (eStopingPowerTable->ElectronicSto 245 eloss = (eStopingPowerTable->ElectronicStoppingPower(z, kineticEnergy))
247 * (material->Ge 246 * (material->GetTotNbOfAtomsPerVolume()) ;
248 247
249 // Experimental data exist only for kinetic 248 // Experimental data exist only for kinetic energy 125 keV
250 } else if( MolecIsInZiegler1988(material)) { 249 } else if( MolecIsInZiegler1988(material)) {
251 250
252 // Cycle over elements - calculation based o 251 // Cycle over elements - calculation based on Bragg's rule
253 G4double eloss125 = 0.0 ; 252 G4double eloss125 = 0.0 ;
254 const G4ElementVector* theElementVector = 253 const G4ElementVector* theElementVector =
255 material->GetElemen 254 material->GetElementVector() ;
256 255
257 256
258 // loop for the elements in the material 257 // loop for the elements in the material
259 for (std::size_t i=0; i<numberOfElements; << 258 for (G4int i=0; i<numberOfElements; i++) {
260 const G4Element* element = (*theElementV 259 const G4Element* element = (*theElementVector)[i] ;
261 G4double z = element->GetZ() ; 260 G4double z = element->GetZ() ;
262 eloss +=(eStopingPowerTable->Electron 261 eloss +=(eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
263 * theAtomi 262 * theAtomicNumDensityVector[i] ;
264 eloss125 +=(eStopingPowerTable->Electron 263 eloss125 +=(eStopingPowerTable->ElectronicStoppingPower(z,125.0*keV))
265 * theAtomi 264 * theAtomicNumDensityVector[i] ;
266 } 265 }
267 266
268 // Chemical factor is taken into account 267 // Chemical factor is taken into account
269 if (eloss125 > 0.0) { << 268 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
270 eloss *= ChemicalFactor(kineticEnergy, e <<
271 } <<
272 269
273 // Brugg's rule calculation 270 // Brugg's rule calculation
274 } else { 271 } else {
275 const G4ElementVector* theElementVector = 272 const G4ElementVector* theElementVector =
276 material->GetElemen 273 material->GetElementVector() ;
277 274
278 // loop for the elements in the material 275 // loop for the elements in the material
279 for (std::size_t i=0; i<numberOfElements; << 276 for (G4int i=0; i<numberOfElements; i++)
280 { 277 {
281 const G4Element* element = (*theElementV 278 const G4Element* element = (*theElementVector)[i] ;
282 G4double z = element->GetZ() ; 279 G4double z = element->GetZ() ;
283 eloss += (eStopingPowerTable->Electron 280 eloss += (eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
284 * theAtomic 281 * theAtomicNumDensityVector[i];
285 } 282 }
286 } 283 }
287 return eloss; 284 return eloss;
288 } 285 }
289 286
290 //....oooOO0OOooo........oooOO0OOooo........oo 287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 288
292 G4bool G4hParametrisedLossModel::MolecIsInZieg 289 G4bool G4hParametrisedLossModel::MolecIsInZiegler1988(
293 const G4Mater 290 const G4Material* material)
294 { 291 {
295 // The list of molecules from 292 // The list of molecules from
296 // J.F.Ziegler and J.M.Manoyan, The stopping 293 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
297 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 294 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
298 295
299 G4String myFormula = G4String(" ") ; 296 G4String myFormula = G4String(" ") ;
300 const G4String chFormula = material->GetChem 297 const G4String chFormula = material->GetChemicalFormula() ;
301 if (myFormula == chFormula ) return false ; 298 if (myFormula == chFormula ) return false ;
302 << 299
303 // There are no evidence for difference of 300 // There are no evidence for difference of stopping power depended on
304 // phase of the compound except for water. << 301 // phase of the compound except for water. The stopping power of the
305 // water in gas phase can be predicted usin 302 // water in gas phase can be predicted using Bragg's rule.
306 // << 303 //
307 // No chemical factor for water-gas << 304 // No chemical factor for water-gas
308 << 305
309 myFormula = G4String("H_2O") ; 306 myFormula = G4String("H_2O") ;
310 const G4State theState = material->GetState( 307 const G4State theState = material->GetState() ;
311 if( theState == kStateGas && myFormula == ch 308 if( theState == kStateGas && myFormula == chFormula) return false ;
312 << 309
313 const std::size_t numberOfMolecula = 53 ; << 310 const size_t numberOfMolecula = 53 ;
314 311
315 // The coffecient from Table.4 of Ziegler & 312 // The coffecient from Table.4 of Ziegler & Manoyan
316 static const G4double HeEff = 2.8735 ; << 313 static const G4double HeEff = 2.8735 ;
317 << 314
318 static const G4String name[numberOfMolecula] 315 static const G4String name[numberOfMolecula] = {
319 "H_2O", "C_2H_4O", "C_3H_6O", "C_2 316 "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 317 "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_ 318 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
322 "CF_4", "C_6H_8", "C_6H_12", "C_ 319 "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_ 320 "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_ 321 "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_ 322 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
326 "SH_2", "CH_4", "CCLF_3", "CC 323 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
327 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_ << 324 "(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_ 325 "(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" 326 "C_3H_6S", "C_4H_4S", "C_7H_8"
330 }; 327 };
331 << 328
332 static const G4double expStopping[numberOfMo 329 static const G4double expStopping[numberOfMolecula] = {
333 66.1, 190.4, 258.7, 42.2, 141.5, << 330 66.1, 190.4, 258.7, 42.2, 141.5,
334 210.9, 279.6, 198.8, 31.0, 267.5, 331 210.9, 279.6, 198.8, 31.0, 267.5,
335 122.8, 311.4, 260.3, 328.9, 391.3, 332 122.8, 311.4, 260.3, 328.9, 391.3,
336 206.6, 374.0, 422.0, 432.0, 398.0, 333 206.6, 374.0, 422.0, 432.0, 398.0,
337 554.0, 353.0, 326.0, 74.6, 220.5, 334 554.0, 353.0, 326.0, 74.6, 220.5,
338 197.4, 362.0, 170.0, 330.5, 211.3, 335 197.4, 362.0, 170.0, 330.5, 211.3,
339 262.3, 349.6, 51.3, 187.0, 236.9, 336 262.3, 349.6, 51.3, 187.0, 236.9,
340 121.9, 35.8, 247.0, 292.6, 268.0, 337 121.9, 35.8, 247.0, 292.6, 268.0,
341 262.3, 49.0, 398.9, 444.0, 22.91, 338 262.3, 49.0, 398.9, 444.0, 22.91,
342 68.0, 155.0, 84.0, 74.2, 254.7, 339 68.0, 155.0, 84.0, 74.2, 254.7,
343 306.8, 324.4, 420.0 340 306.8, 324.4, 420.0
344 } ; 341 } ;
345 342
346 static const G4double expCharge[numberOfMole 343 static const G4double expCharge[numberOfMolecula] = {
347 HeEff, HeEff, HeEff, 1.0, HeEff, << 344 HeEff, HeEff, HeEff, 1.0, HeEff,
348 HeEff, HeEff, HeEff, 1.0, 1.0, 345 HeEff, HeEff, HeEff, 1.0, 1.0,
349 1.0, HeEff, HeEff, HeEff, HeEff, 346 1.0, HeEff, HeEff, HeEff, HeEff,
350 HeEff, HeEff, HeEff, HeEff, HeEff, 347 HeEff, HeEff, HeEff, HeEff, HeEff,
351 HeEff, HeEff, HeEff, 1.0, HeEff, 348 HeEff, HeEff, HeEff, 1.0, HeEff,
352 HeEff, HeEff, HeEff, HeEff, HeEff, 349 HeEff, HeEff, HeEff, HeEff, HeEff,
353 HeEff, HeEff, 1.0, HeEff, HeEff, 350 HeEff, HeEff, 1.0, HeEff, HeEff,
354 HeEff, 1.0, HeEff, HeEff, HeEff, 351 HeEff, 1.0, HeEff, HeEff, HeEff,
355 HeEff, 1.0, HeEff, HeEff, 1.0, 352 HeEff, 1.0, HeEff, HeEff, 1.0,
356 1.0, 1.0, 1.0, 1.0, HeEff, 353 1.0, 1.0, 1.0, 1.0, HeEff,
357 HeEff, HeEff, HeEff 354 HeEff, HeEff, HeEff
358 } ; 355 } ;
359 356
360 static const G4double numberOfAtomsPerMolecu 357 static const G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
361 3.0, 7.0, 10.0, 4.0, 6.0, << 358 3.0, 7.0, 10.0, 4.0, 6.0,
362 9.0, 12.0, 7.0, 4.0, 24.0, 359 9.0, 12.0, 7.0, 4.0, 24.0,
363 12.0, 14.0, 10.0, 13.0, 5.0, 360 12.0, 14.0, 10.0, 13.0, 5.0,
364 5.0, 14.0, 18.0, 17.0, 17.0, 361 5.0, 14.0, 18.0, 17.0, 17.0,
365 24.0, 15.0, 13.0, 9.0, 8.0, 362 24.0, 15.0, 13.0, 9.0, 8.0,
366 6.0, 14.0, 8.0, 8.0, 9.0, 363 6.0, 14.0, 8.0, 8.0, 9.0,
367 10.0, 15.0, 6.0, 7.0, 7.0, 364 10.0, 15.0, 6.0, 7.0, 7.0,
368 3.0, 5.0, 5.0, 5.0, 5.0, 365 3.0, 5.0, 5.0, 5.0, 5.0,
369 9.0, 3.0, 16.0, 14.0, 3.0, 366 9.0, 3.0, 16.0, 14.0, 3.0,
370 9.0, 16.0, 11.0, 9.0, 10.0, 367 9.0, 16.0, 11.0, 9.0, 10.0,
371 10.0, 9.0, 15.0 368 10.0, 9.0, 15.0
372 } ; 369 } ;
373 370
374 // Search for the compaund in the table 371 // Search for the compaund in the table
375 for (std::size_t i=0; i<numberOfMolecula; ++ << 372 for (size_t i=0; i<numberOfMolecula; i++)
376 { << 373 {
377 if(chFormula == name[i]) { << 374 if(chFormula == name[i]) {
378 G4double exp125 = expStopping[i] * << 375 G4double exp125 = expStopping[i] *
379 (material->GetTotNbOfAtoms 376 (material->GetTotNbOfAtomsPerVolume()) /
380 (expCharge[i] * numberOfAt 377 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
381 SetExpStopPower125(exp125) ; 378 SetExpStopPower125(exp125) ;
382 return true ; 379 return true ;
383 } 380 }
384 } 381 }
385 << 382
386 return false ; 383 return false ;
387 } 384 }
388 385
389 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 387
391 G4double G4hParametrisedLossModel::ChemicalFac 388 G4double G4hParametrisedLossModel::ChemicalFactor(
392 G4double kineticEn 389 G4double kineticEnergy, G4double eloss125) const
393 { 390 {
394 // Approximation of Chemical Factor accordin 391 // Approximation of Chemical Factor according to
395 // J.F.Ziegler and J.M.Manoyan, The stopping 392 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
396 // Nucl. Inst. & Meth. in Phys. Res. B35 (19 393 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
397 << 394
398 G4double gamma = 1.0 + kineticEnergy/prot << 395 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
399 G4double gamma25 = 1.0 + 25.0*keV /proton_m 396 G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
400 G4double gamma125 = 1.0 + 125.0*keV/proton_m 397 G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
401 G4double beta = std::sqrt(1.0 - 1.0/(gam 398 G4double beta = std::sqrt(1.0 - 1.0/(gamma*gamma)) ;
402 G4double beta25 = std::sqrt(1.0 - 1.0/(gam 399 G4double beta25 = std::sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
403 G4double beta125 = std::sqrt(1.0 - 1.0/(gam 400 G4double beta125 = std::sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
404 << 401
405 G4double factor = 1.0 + (expStopPower125/elo 402 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
406 (1.0 + G4Exp( 1.48 * ( beta << 403 (1.0 + std::exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
407 (1.0 + G4Exp( 1.48 * ( beta << 404 (1.0 + std::exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
408 << 405
409 return factor ; 406 return factor ;
410 } 407 }
411 408
412 //....oooOO0OOooo........oooOO0OOooo........oo 409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413 410