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