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