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