ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hParametrisedLossModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hParametrisedLossModel.cc
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  G4String blank = G4String(" ") ;
90  G4String ir49p = G4String("ICRU_R49p") ;
91  G4String ir49He = G4String("ICRU_R49He") ;
92  G4String zi85p = G4String("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 G4int 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 (G4int 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 (G4int 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  eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
270 
271  // Brugg's rule calculation
272  } else {
273  const G4ElementVector* theElementVector =
274  material->GetElementVector() ;
275 
276  // loop for the elements in the material
277  for (G4int i=0; i<numberOfElements; i++)
278  {
279  const G4Element* element = (*theElementVector)[i] ;
280  G4double z = element->GetZ() ;
281  eloss += (eStopingPowerTable->ElectronicStoppingPower(z,kineticEnergy))
282  * theAtomicNumDensityVector[i];
283  }
284  }
285  return eloss;
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 
290 G4bool G4hParametrisedLossModel::MolecIsInZiegler1988(
291  const G4Material* material)
292 {
293  // The list of molecules from
294  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
295  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
296 
297  G4String myFormula = G4String(" ") ;
298  const G4String chFormula = material->GetChemicalFormula() ;
299  if (myFormula == chFormula ) return false ;
300 
301  // There are no evidence for difference of stopping power depended on
302  // phase of the compound except for water. The stopping power of the
303  // water in gas phase can be predicted using Bragg's rule.
304  //
305  // No chemical factor for water-gas
306 
307  myFormula = G4String("H_2O") ;
308  const G4State theState = material->GetState() ;
309  if( theState == kStateGas && myFormula == chFormula) return false ;
310 
311  const size_t numberOfMolecula = 53 ;
312 
313  // The coffecient from Table.4 of Ziegler & Manoyan
314  static const G4double HeEff = 2.8735 ;
315 
316  static const G4String name[numberOfMolecula] = {
317  "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
318  "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
319  "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
320  "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
321  "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
322  "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
323  "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
324  "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
325  "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
326  "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
327  "C_3H_6S", "C_4H_4S", "C_7H_8"
328  };
329 
330  static const G4double expStopping[numberOfMolecula] = {
331  66.1, 190.4, 258.7, 42.2, 141.5,
332  210.9, 279.6, 198.8, 31.0, 267.5,
333  122.8, 311.4, 260.3, 328.9, 391.3,
334  206.6, 374.0, 422.0, 432.0, 398.0,
335  554.0, 353.0, 326.0, 74.6, 220.5,
336  197.4, 362.0, 170.0, 330.5, 211.3,
337  262.3, 349.6, 51.3, 187.0, 236.9,
338  121.9, 35.8, 247.0, 292.6, 268.0,
339  262.3, 49.0, 398.9, 444.0, 22.91,
340  68.0, 155.0, 84.0, 74.2, 254.7,
341  306.8, 324.4, 420.0
342  } ;
343 
344  static const G4double expCharge[numberOfMolecula] = {
345  HeEff, HeEff, HeEff, 1.0, HeEff,
346  HeEff, HeEff, HeEff, 1.0, 1.0,
347  1.0, HeEff, HeEff, HeEff, HeEff,
348  HeEff, HeEff, HeEff, HeEff, HeEff,
349  HeEff, HeEff, HeEff, 1.0, HeEff,
350  HeEff, HeEff, HeEff, HeEff, HeEff,
351  HeEff, HeEff, 1.0, HeEff, HeEff,
352  HeEff, 1.0, HeEff, HeEff, HeEff,
353  HeEff, 1.0, HeEff, HeEff, 1.0,
354  1.0, 1.0, 1.0, 1.0, HeEff,
355  HeEff, HeEff, HeEff
356  } ;
357 
358  static const G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
359  3.0, 7.0, 10.0, 4.0, 6.0,
360  9.0, 12.0, 7.0, 4.0, 24.0,
361  12.0, 14.0, 10.0, 13.0, 5.0,
362  5.0, 14.0, 18.0, 17.0, 17.0,
363  24.0, 15.0, 13.0, 9.0, 8.0,
364  6.0, 14.0, 8.0, 8.0, 9.0,
365  10.0, 15.0, 6.0, 7.0, 7.0,
366  3.0, 5.0, 5.0, 5.0, 5.0,
367  9.0, 3.0, 16.0, 14.0, 3.0,
368  9.0, 16.0, 11.0, 9.0, 10.0,
369  10.0, 9.0, 15.0
370  } ;
371 
372  // Search for the compaund in the table
373  for (size_t i=0; i<numberOfMolecula; i++)
374  {
375  if(chFormula == name[i]) {
376  G4double exp125 = expStopping[i] *
377  (material->GetTotNbOfAtomsPerVolume()) /
378  (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
379  SetExpStopPower125(exp125) ;
380  return true ;
381  }
382  }
383 
384  return false ;
385 }
386 
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388 
389 G4double G4hParametrisedLossModel::ChemicalFactor(
390  G4double kineticEnergy, G4double eloss125) const
391 {
392  // Approximation of Chemical Factor according to
393  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
394  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
395 
396  G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
397  G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
398  G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
399  G4double beta = std::sqrt(1.0 - 1.0/(gamma*gamma)) ;
400  G4double beta25 = std::sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
401  G4double beta125 = std::sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
402 
403  G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
404  (1.0 + G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
405  (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
406 
407  return factor ;
408 }
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....