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