ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hICRU49He.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hICRU49He.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: G4hICRU49He
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/09/2000 V.Ivanchenko clean up - all variable are the same as in ICRU
41 // 03/10/2000 V.Ivanchenko clean up accoding to CodeWizard
42 // 10/05/2001 V.Ivanchenko Clean up againist Linux compilation with -Wall
43 // 26/08/2004 V.Ivanchenko Fix a problem of effective charge
44 //
45 // Class Description:
46 //
47 // Electronic stopping power parametrised according to
48 // ICRU Report N49, 1993. J.F. Ziegler model for He ion.
49 //
50 // Class Description: End
51 //
52 // -------------------------------------------------------------------
53 //
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 
56 #include "G4hICRU49He.hh"
57 
58 #include "globals.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4UnitsTable.hh"
62 #include "G4Material.hh"
63 #include "G4Exp.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 G4hICRU49He::G4hICRU49He():G4VhElectronicStoppingPower(),
68  rateMass(4.0026/1.007276),
69  iMolecula(0)
70 {;}
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 G4hICRU49He::~G4hICRU49He()
75 {;}
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
79 G4bool G4hICRU49He::HasMaterial(const G4Material* material)
80 {
81  G4String chFormula = material->GetChemicalFormula() ;
82  G4String myFormula = G4String(" ");
83 
84  if (myFormula == chFormula ) {
85  if(1 == (material->GetNumberOfElements())) {return true;}
86  return false ;
87  }
88 
89  // ICRU Report N49, 1993. Power's model for He.
90  const size_t numberOfMolecula = 30 ;
91  static const G4String name[numberOfMolecula] = {
92  "H_2", "Be-Solid", "C-Solid", "Graphite", "N_2",
93  "O_2", "Al-Solid", "Si-Solid", "Ar-Solid", "Cu-Solid",
94  "Ge", "W-Solid", "Au-Solid", "Pb-Solid", "C_2H_2",
95  "CO_2", "Cellulose-Nitrat", "C_2H_4", "LiF",
96  "CH_4", "Nylon", "Polycarbonate", "(CH_2)_N-Polyetilene", "PMMA",
97  "(C_8H_8)_N", "SiO_2", "CsI", "H_2O", "H_2O-Gas"};
98 
99  // Special treatment for water in gas state
100 
101  myFormula = G4String("H_2O") ;
102  const G4State theState = material->GetState() ;
103  if( theState == kStateGas && myFormula == chFormula) {
104  chFormula = G4String("H_2O-Gas");
105  }
106 
107  // Search for the material in the table
108  for (size_t i=0; i<numberOfMolecula; i++) {
109  if (chFormula == name[i]) {
110  SetMoleculaNumber(i) ;
111  return true ;
112  }
113  }
114  return false ;
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
119 G4double G4hICRU49He::StoppingPower(const G4Material* material,
120  G4double kineticEnergy)
121 {
122  G4double ionloss = 0.0 ;
123 
124  // pure material (normally not the case for this function)
125  if(1 == (material->GetNumberOfElements())) {
126  G4double z = material->GetZ() ;
127  ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
128 
129  // The data and the fit from:
130  // ICRU Report N49, 1993. Power's model for He.
131  } else if ( iMolecula < 30 ) {
132 
133  // Reduced kinetic energy
134  // in internal units of parametrisation formula (MeV)
135  G4double T = kineticEnergy*rateMass/MeV ;
136 
137  static const G4double c[30][7] = {
138  {8.0080, 3.6287, 23.0700, 14.9900, 0.8507, 0.60, 2.0
139  },{ 13.3100, 3.7432, 39.4130, 12.1990, 1.0950, 0.38, 1.4
140  },{ 22.7240, 3.6040, 47.1810, 17.5490, 0.9040, 0.40, 1.4
141  },{ 24.4040, 2.4032, 48.9440, 27.9730, 1.2933, 0.40, 1.6
142  },{ 58.4719, 1.5115, 77.6421, 102.490, 1.5811, 0.50, 2.0
143  },{ 60.5408, 1.6297, 91.7601, 94.1260, 1.3662, 0.50, 2.0
144  },{ 48.4480, 6.4323, 59.2890, 18.3810, 0.4937, 0.48, 1.6
145  },{ 59.0346, 5.1305, 47.0866, 30.0857, 0.3500, 0.60, 2.0
146  },{ 71.8691, 2.8250, 51.1658, 57.1235, 0.4477, 0.60, 2.0
147  },{ 78.3520, 4.0961, 136.731, 28.4470, 1.0621, 0.52, 1.2
148  },{ 120.553, 1.5374, 49.8740, 82.2980, 0.8733, 0.45, 1.6
149  },{ 249.896, 0.6996, -37.274, 248.592, 1.1052, 0.50, 1.5
150  },{ 246.698, 0.6219, -58.391, 292.921, 0.8186, 0.56, 1.8
151  },{ 248.563, 0.6235, -36.8968, 306.960, 1.3214, 0.50, 2.0
152  },{ 25.5860, 1.7125, 154.723, 118.620, 2.2580, 0.50, 2.0
153  },{ 138.294, 25.6413, 231.873, 17.3780, 0.3218, 0.58, 1.3
154  },{ 83.2091, 1.1294, 135.7457, 190.865, 2.3461, 0.50, 2.0
155  },{ 263.542, 1.4754, 1541.446, 781.898, 1.9209, 0.40, 2.0
156  },{ 59.5545, 1.5354, 132.1523, 153.3537, 2.0262, 0.50, 2.0
157  },{ 31.7380, 19.820, 125.2100, 6.8910, 0.7242, 0.50, 1.1
158  },{ 31.7549, 1.5682, 97.4777, 106.0774, 2.3204, 0.50, 2.0
159  },{ 230.465, 4.8967, 1845.320, 358.641, 1.0774, 0.46, 1.2
160  },{ 423.444, 5.3761, 1189.114, 319.030, 0.7652, 0.48, 1.5
161  },{ 86.3410, 3.3322, 91.0433, 73.1091, 0.4650, 0.50, 2.0
162  },{ 146.105, 9.4344, 515.1500, 82.8860, 0.6239, 0.55, 1.5
163  },{ 238.050, 5.6901, 372.3575, 146.1835, 0.3992, 0.50, 2.0
164  },{ 124.2338, 2.6730, 133.8175, 99.4109, 0.7776, 0.50, 2.0
165  },{ 221.723, 1.5415, 87.7315, 192.5266, 1.0742, 0.50, 2.0
166  },{ 26.7537, 1.3717, 90.8007, 77.1587, 2.3264, 0.50, 2.0
167  },{ 37.6121, 1.8052, 73.0250, 66.2070, 1.4038, 0.50, 2.0} };
168 
169  G4double a1,a2 ;
170 
171  // Free electron gas model
172  if ( T < 0.001 ) {
173  G4double T0 = 0.001 ;
174  a1 = 1.0 - G4Exp(-c[iMolecula][1]*std::pow(T0,-2.0+c[iMolecula][5])) ;
175  a2 = (c[iMolecula][0]*std::log(T0)/T0 + c[iMolecula][2]/T0) *
176  G4Exp(-c[iMolecula][4]*std::pow(T0,-c[iMolecula][6])) +
177  c[iMolecula][3]/(T0*T0) ;
178 
179  ionloss = a1*a2*std::sqrt(T/T0) ;
180 
181  // Main parametrisation
182  } else {
183  a1 = 1.0 - G4Exp(-c[iMolecula][1]*std::pow(T,-2.0+c[iMolecula][5])) ;
184  a2 = (c[iMolecula][0]*std::log(T)/T + c[iMolecula][2]/T) *
185  G4Exp(-c[iMolecula][4]*std::pow(T,-c[iMolecula][6])) +
186  c[iMolecula][3]/(T*T) ;
187  ionloss = a1*a2;
188  }
189 
190  // He effective charge
191  G4double z = (material->GetTotNbOfElectPerVolume()) /
192  (material->GetTotNbOfAtomsPerVolume()) ;
193 
194  ionloss /= HeEffChargeSquare(z, kineticEnergy*rateMass) ;
195 
196  if ( ionloss < 0.0) ionloss = 0.0 ;
197  }
198 
199  return ionloss ;
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 
204 G4double G4hICRU49He::ElectronicStoppingPower(G4double z,
205  G4double kineticEnergy) const
206 {
207  G4double ionloss ;
208  G4int i = G4int(z)-1 ; // index of atom
209  if(i < 0) i = 0 ;
210  if(i > 91) i = 91 ;
211 
212  // The data and the fit from:
213  // ICRU Report 49, 1993. Ziegler's type of parametrisations
214  // Reduced kinetic energy
215 
216  // He energy in internal units of parametrisation formula (MeV)
217  G4double T = kineticEnergy*rateMass/MeV ;
218 
219  static const G4double a[92][5] = {
220  {0.35485, 0.6456, 6.01525, 20.8933, 4.3515
221  },{ 0.58, 0.59, 6.3, 130.0, 44.07
222  },{ 1.42, 0.49, 12.25, 32.0, 9.161
223  },{ 2.1895, 0.47183,7.2362, 134.30, 197.96
224  },{ 3.691, 0.4128, 18.48, 50.72, 9.0
225  },{ 3.83523, 0.42993,12.6125, 227.41, 188.97
226  },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768
227  },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382
228  },{ 1.533, 0.531, 40.44, 18.41, 2.718
229  },{ 2.303, 0.4861, 37.01, 37.96, 5.092
230  },{ 9.894, 0.3081, 23.65, 0.384, 92.93
231  },{ 4.3, 0.47, 34.3, 3.3, 12.74
232  },{ 2.5, 0.625, 45.7, 0.1, 4.359
233  },{ 2.1, 0.65, 49.34, 1.788, 4.133
234  },{ 1.729, 0.6562, 53.41, 2.405, 3.845
235  },{ 1.402, 0.6791, 58.98, 3.528, 3.211
236  },{ 1.117, 0.7044, 69.69, 3.705, 2.156
237  },{ 2.291, 0.6284, 73.88, 4.478, 2.066
238  },{ 8.554, 0.3817, 83.61, 11.84, 1.875
239  },{ 6.297, 0.4622, 65.39, 10.14, 5.036
240  },{ 5.307, 0.4918, 61.74, 12.4, 6.665
241  },{ 4.71, 0.5087, 65.28, 8.806, 5.948
242  },{ 6.151, 0.4524, 83.0, 18.31, 2.71
243  },{ 6.57, 0.4322, 84.76, 15.53, 2.779
244  },{ 5.738, 0.4492, 84.6, 14.18, 3.101
245  },{ 5.013, 0.4707, 85.8, 16.55, 3.211
246  },{ 4.32, 0.4947, 76.14, 10.85, 5.441
247  },{ 4.652, 0.4571, 80.73, 22.0, 4.952
248  },{ 3.114, 0.5236, 76.67, 7.62, 6.385
249  },{ 3.114, 0.5236, 76.67, 7.62, 7.502
250  },{ 3.114, 0.5236, 76.67, 7.62, 8.514
251  },{ 5.746, 0.4662, 79.24, 1.185, 7.993
252  },{ 2.792, 0.6346, 106.1, 0.2986, 2.331
253  },{ 4.667, 0.5095, 124.3, 2.102, 1.667
254  },{ 2.44, 0.6346, 105.0, 0.83, 2.851
255  },{ 1.413, 0.7377, 147.9, 1.466, 1.016
256  },{ 11.72, 0.3826, 102.8, 9.231, 4.371
257  },{ 7.126, 0.4804, 119.3, 5.784, 2.454
258  },{ 11.61, 0.3955, 146.7, 7.031, 1.423
259  },{ 10.99, 0.41, 163.9, 7.1, 1.052
260  },{ 9.241, 0.4275, 163.1, 7.954, 1.102
261  },{ 9.276, 0.418, 157.1, 8.038, 1.29
262  },{ 3.999, 0.6152, 97.6, 1.297, 5.792
263  },{ 4.306, 0.5658, 97.99, 5.514, 5.754
264  },{ 3.615, 0.6197, 86.26, 0.333, 8.689
265  },{ 5.8, 0.49, 147.2, 6.903, 1.289
266  },{ 5.6, 0.49, 130.0, 10.0, 2.844
267  },{ 3.55, 0.6068, 124.7, 1.112, 3.119
268  },{ 3.6, 0.62, 105.8, 0.1692, 6.026
269  },{ 5.4, 0.53, 103.1, 3.931, 7.767
270  },{ 3.97, 0.6459, 131.8, 0.2233, 2.723
271  },{ 3.65, 0.64, 126.8, 0.6834, 3.411
272  },{ 3.118, 0.6519, 164.9, 1.208, 1.51
273  },{ 3.949, 0.6209, 200.5, 1.878, 0.9126
274  },{ 14.4, 0.3923, 152.5, 8.354, 2.597
275  },{ 10.99, 0.4599, 138.4, 4.811, 3.726
276  },{ 16.6, 0.3773, 224.1, 6.28, 0.9121
277  },{ 10.54, 0.4533, 159.3, 4.832, 2.529
278  },{ 10.33, 0.4502, 162.0, 5.132, 2.444
279  },{ 10.15, 0.4471, 165.6, 5.378, 2.328
280  },{ 9.976, 0.4439, 168.0, 5.721, 2.258
281  },{ 9.804, 0.4408, 176.2, 5.675, 1.997
282  },{ 14.22, 0.363, 228.4, 7.024, 1.016
283  },{ 9.952, 0.4318, 233.5, 5.065, 0.9244
284  },{ 9.272, 0.4345, 210.0, 4.911, 1.258
285  },{ 10.13, 0.4146, 225.7, 5.525, 1.055
286  },{ 8.949, 0.4304, 213.3, 5.071, 1.221
287  },{ 11.94, 0.3783, 247.2, 6.655, 0.849
288  },{ 8.472, 0.4405, 195.5, 4.051, 1.604
289  },{ 8.301, 0.4399, 203.7, 3.667, 1.459
290  },{ 6.567, 0.4858, 193.0, 2.65, 1.66
291  },{ 5.951, 0.5016, 196.1, 2.662, 1.589
292  },{ 7.495, 0.4523, 251.4, 3.433, 0.8619
293  },{ 6.335, 0.4825, 255.1, 2.834, 0.8228
294  },{ 4.314, 0.5558, 214.8, 2.354, 1.263
295  },{ 4.02, 0.5681, 219.9, 2.402, 1.191
296  },{ 3.836, 0.5765, 210.2, 2.742, 1.305
297  },{ 4.68, 0.5247, 244.7, 2.749, 0.8962
298  },{ 3.223, 0.5883, 232.7, 2.954, 1.05
299  },{ 2.892, 0.6204, 208.6, 2.415, 1.416
300  },{ 4.728, 0.5522, 217.0, 3.091, 1.386
301  },{ 6.18, 0.52, 170.0, 4.0, 3.224
302  },{ 9.0, 0.47, 198.0, 3.8, 2.032
303  },{ 2.324, 0.6997, 216.0, 1.599, 1.399
304  },{ 1.961, 0.7286, 223.0, 1.621, 1.296
305  },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507
306  },{ 10.31, 0.4613, 261.2, 4.738, 0.9899
307  },{ 7.962, 0.519, 235.7, 4.347, 1.313
308  },{ 6.227, 0.5645, 231.9, 3.961, 1.379
309  },{ 5.246, 0.5947, 228.6, 4.027, 1.432
310  },{ 5.408, 0.5811, 235.7, 3.961, 1.358
311  },{ 5.218, 0.5828, 245.0, 3.838, 1.25}
312  };
313 
314  // Free electron gas model
315  if ( T < 0.001 ) {
316  G4double slow = a[i][0] ;
317  G4double shigh = std::log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
318  * a[i][2]*1000.0 ;
319  ionloss = slow*shigh / (slow + shigh) ;
320  ionloss *= std::sqrt(T*1000.0) ;
321 
322  // Main parametrisation
323  } else {
324  G4double slow = a[i][0] * std::pow((T*1000.0), a[i][1]) ;
325  G4double shigh = std::log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
326  ionloss = slow*shigh / (slow + shigh) ;
327 
328  }
329  if ( ionloss < 0.0) ionloss = 0.0 ;
330 
331  // He effective charge
332  ionloss /= HeEffChargeSquare(z, kineticEnergy*rateMass) ;
333 
334  return ionloss;
335 }