| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4hIonEffChargeSquare 32 // File name: G4hIonEffChargeSquare
33 // 33 //
34 // Author: V.Ivanchenko (Vladimir.Ivanc 34 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35 // << 35 //
36 // Creation date: 20 July 2000 36 // Creation date: 20 July 2000
37 // 37 //
38 // Modifications: << 38 // Modifications:
39 // 20/07/2000 V.Ivanchenko First implementati 39 // 20/07/2000 V.Ivanchenko First implementation
40 // 18/06/2001 V.Ivanchenko Continuation for e 40 // 18/06/2001 V.Ivanchenko Continuation for eff.charge (small change of y)
41 // 08/10/2002 V.Ivanchenko The charge of the << 41 // 08/10/2002 V.Ivanchenko The charge of the nucleus is used not charge of
42 // DynamicParticle 42 // DynamicParticle
43 // 43 //
44 // Class Description: << 44 // Class Description:
45 // 45 //
46 // Ion effective charge model 46 // Ion effective charge model
47 // J.F.Ziegler and J.M.Manoyan, The stopping o 47 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988 48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
49 // 49 //
50 // Class Description: End << 50 // Class Description: End
51 // 51 //
52 // ------------------------------------------- 52 // -------------------------------------------------------------------
53 // 53 //
54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 55
56 #include "G4hIonEffChargeSquare.hh" 56 #include "G4hIonEffChargeSquare.hh"
57 #include "G4PhysicalConstants.hh" <<
58 #include "G4SystemOfUnits.hh" <<
59 #include "G4DynamicParticle.hh" 57 #include "G4DynamicParticle.hh"
60 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
61 #include "G4Material.hh" 59 #include "G4Material.hh"
62 #include "G4Element.hh" 60 #include "G4Element.hh"
63 #include "G4Exp.hh" <<
64 61
65 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 63
67 G4hIonEffChargeSquare::G4hIonEffChargeSquare(c 64 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name)
68 : G4VLowEnergyModel(name), << 65 : G4VLowEnergyModel(name),
69 theHeMassAMU(4.0026) 66 theHeMassAMU(4.0026)
70 {;} 67 {;}
71 68
72 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 70
74 G4hIonEffChargeSquare::~G4hIonEffChargeSquare( << 71 G4hIonEffChargeSquare::~G4hIonEffChargeSquare()
75 {;} 72 {;}
76 73
77 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 75
79 G4double G4hIonEffChargeSquare::TheValue(const 76 G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle,
80 const G4Mat << 77 const G4Material* material)
81 { 78 {
82 G4double energy = particle->GetKineticEnergy 79 G4double energy = particle->GetKineticEnergy() ;
83 G4double particleMass = particle->GetMass() 80 G4double particleMass = particle->GetMass() ;
84 G4double charge = (particle->GetDefinition() 81 G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ;
85 82
86 G4double q = IonEffChargeSquare(material,ene 83 G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ;
87 84
88 return q ; 85 return q ;
89 } 86 }
90 87
91 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 89
93 G4double G4hIonEffChargeSquare::TheValue(const 90 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle,
94 const G4Mat 91 const G4Material* material,
95 G4double kineticEnergy) << 92 G4double kineticEnergy)
96 { 93 {
97 // SetRateMass(aParticle) ; 94 // SetRateMass(aParticle) ;
98 G4double particleMass = aParticle->GetPDGMas 95 G4double particleMass = aParticle->GetPDGMass() ;
99 G4double charge = (aParticle->GetPDGCharge() 96 G4double charge = (aParticle->GetPDGCharge())/eplus ;
100 97
101 G4double q = IonEffChargeSquare(material,kin 98 G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ;
102 99
103 return q ; 100 return q ;
104 } 101 }
105 102
106 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 104
108 G4double G4hIonEffChargeSquare::HighEnergyLimi 105 G4double G4hIonEffChargeSquare::HighEnergyLimit(
109 const G4ParticleDefinition* , 106 const G4ParticleDefinition* ,
110 const G4Material* ) const 107 const G4Material* ) const
111 { 108 {
112 return 1.0*TeV ; 109 return 1.0*TeV ;
113 } 110 }
114 111
115 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 113
117 G4double G4hIonEffChargeSquare::LowEnergyLimit 114 G4double G4hIonEffChargeSquare::LowEnergyLimit(
118 const G4ParticleDefinition* , 115 const G4ParticleDefinition* ,
119 const G4Material* ) const 116 const G4Material* ) const
120 { 117 {
121 return 0.0 ; 118 return 0.0 ;
122 } 119 }
123 120
124 //....oooOO0OOooo........oooOO0OOooo........oo 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125 122
126 G4double G4hIonEffChargeSquare::HighEnergyLimi 123 G4double G4hIonEffChargeSquare::HighEnergyLimit(
127 const G4ParticleDefinition* ) cons 124 const G4ParticleDefinition* ) const
128 { 125 {
129 return 1.0*TeV ; 126 return 1.0*TeV ;
130 } 127 }
131 128
132 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 130
134 G4double G4hIonEffChargeSquare::LowEnergyLimit 131 G4double G4hIonEffChargeSquare::LowEnergyLimit(
135 const G4ParticleDefinition* ) 132 const G4ParticleDefinition* ) const
136 { 133 {
137 return 0.0 ; 134 return 0.0 ;
138 } 135 }
139 136
140 //....oooOO0OOooo........oooOO0OOooo........oo 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141 << 138
142 G4bool G4hIonEffChargeSquare::IsInCharge(const 139 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* ,
143 const G4Material* 140 const G4Material* ) const
144 { 141 {
145 return true ; 142 return true ;
146 } 143 }
147 144
148 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 << 146
150 G4bool G4hIonEffChargeSquare::IsInCharge(const 147 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* ,
151 const G4Mat 148 const G4Material* ) const
152 { 149 {
153 return true ; 150 return true ;
154 } 151 }
155 152
156 //....oooOO0OOooo........oooOO0OOooo........oo 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 154
158 G4double G4hIonEffChargeSquare::IonEffChargeSq 155 G4double G4hIonEffChargeSquare::IonEffChargeSquare(
159 const G4Material* material, 156 const G4Material* material,
160 G4double kineticEnergy, 157 G4double kineticEnergy,
161 G4double particleMass, 158 G4double particleMass,
162 G4double ionCharge) const 159 G4double ionCharge) const
163 { 160 {
164 // The aproximation of ion effective charge << 161 // The aproximation of ion effective charge from:
165 // J.F.Ziegler, J.P. Biersack, U. Littmark 162 // J.F.Ziegler, J.P. Biersack, U. Littmark
166 // The Stopping and Range of Ions in Matter, 163 // The Stopping and Range of Ions in Matter,
167 // Vol.1, Pergamon Press, 1985 164 // Vol.1, Pergamon Press, 1985
168 165
169 // Fast ions or hadrons 166 // Fast ions or hadrons
170 G4double reducedEnergy = kineticEnergy * pro 167 G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ;
171 if(reducedEnergy < 1.0*keV) reducedEnergy = 168 if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV;
172 if( (reducedEnergy > ionCharge * 10.0 * MeV) << 169 if( (reducedEnergy > ionCharge * 10.0 * MeV) ||
173 (ionCharge < 1.5) ) return ionCharge*ion 170 (ionCharge < 1.5) ) return ionCharge*ionCharge ;
174 171
175 static const G4double vFermi[92] = { << 172 static G4double vFermi[92] = {
176 1.0309, 0.15976, 0.59782, 1.0781, 1.0486 173 1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424,
177 0.45259, 0.71074, 0.90519, 0.97411, 0.9718 174 0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
178 0.81707, 0.9943, 1.1423, 1.2381, 1.1222 175 0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411,
179 0.84912, 0.95, 1.0903, 1.0429, 0.4971 176 0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
180 1.029, 1.2542, 1.122, 1.1241, 1.0882 177 1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054,
181 0.93155, 1.0047, 0.55379, 0.43289, 0.3263 178 0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
182 0.71418, 0.71453, 0.5911, 0.70263, 0.6804 179 0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
183 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 180 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
184 0.66401, 0.84912, 0.88433, 0.80746, 0.4335 181 0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
185 1.9578, 1.0257} ; 182 1.9578, 1.0257} ;
186 183
187 static const G4double c[6] = {0.2865, 0.126 << 184 static G4double lFactor[92] = {
>> 185 1.0, 1.0, 1.1, 1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
>> 186 0.82, 0.81, 0.83, 0.88, 1.0, 0.95, 0.97, 0.99, 0.98, 0.97,
>> 187 0.98, 0.97, 0.96, 0.93, 0.91, 0.9, 0.88, 0.9, 0.9, 0.9,
>> 188 0.9, 0.85, 0.9, 0.9, 0.91, 0.92, 0.9, 0.9, 0.9, 0.9,
>> 189 0.9, 0.88, 0.9, 0.88, 0.88, 0.9, 0.9, 0.88, 0.9, 0.9,
>> 190 0.9, 0.9, 0.96, 1.2, 0.9, 0.88, 0.88, 0.85, 0.9, 0.9,
>> 191 0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1, 1.08, 1.08,
>> 192 1.08, 1.08, 1.09, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15,
>> 193 1.17, 1.2, 1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
>> 194 1.16, 1.16} ;
>> 195
>> 196 static G4double c[6] = {0.2865, 0.1266, -0.001429,
188 0.02402,-0.01135, 0. 197 0.02402,-0.01135, 0.001475} ;
189 198
190 // get elements in the actual material, 199 // get elements in the actual material,
191 const G4ElementVector* theElementVector = ma 200 const G4ElementVector* theElementVector = material->GetElementVector() ;
192 const G4double* theAtomicNumDensityVector = << 201 const G4double* theAtomicNumDensityVector =
193 material->GetAtomicNu 202 material->GetAtomicNumDensityVector() ;
194 const G4int NumberOfElements = (G4int)materi << 203 const G4int NumberOfElements = material->GetNumberOfElements() ;
195 << 204
196 // loop for the elements in the material 205 // loop for the elements in the material
197 // to find out average values Z, vF, lF 206 // to find out average values Z, vF, lF
198 G4double z = 0.0, vF = 0.0, norm = 0.0 ; << 207 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
199 208
200 if( 1 == NumberOfElements ) { 209 if( 1 == NumberOfElements ) {
201 z = material->GetZ() ; 210 z = material->GetZ() ;
202 G4int iz = G4int(z) - 1 ; 211 G4int iz = G4int(z) - 1 ;
203 if(iz < 0) iz = 0 ; 212 if(iz < 0) iz = 0 ;
204 else if(iz > 91) iz = 91 ; 213 else if(iz > 91) iz = 91 ;
205 vF = vFermi[iz] ; 214 vF = vFermi[iz] ;
>> 215 lF = lFactor[iz] ;
206 216
207 } else { 217 } else {
208 for (G4int iel=0; iel<NumberOfElements; ie 218 for (G4int iel=0; iel<NumberOfElements; iel++)
209 { 219 {
210 const G4Element* element = (*theElemen 220 const G4Element* element = (*theElementVector)[iel] ;
211 G4double z2 = element->GetZ() ; 221 G4double z2 = element->GetZ() ;
212 const G4double weight = theAtomicNumDe 222 const G4double weight = theAtomicNumDensityVector[iel] ;
213 norm += weight ; 223 norm += weight ;
214 z += z2 * weight ; 224 z += z2 * weight ;
215 G4int iz = G4int(z2) - 1 ; 225 G4int iz = G4int(z2) - 1 ;
216 if(iz < 0) iz = 0 ; 226 if(iz < 0) iz = 0 ;
217 else if(iz > 91) iz =91 ; 227 else if(iz > 91) iz =91 ;
218 vF += vFermi[iz] * weight ; 228 vF += vFermi[iz] * weight ;
>> 229 lF += lFactor[iz] * weight ;
219 } 230 }
220 if (norm > 0.0) { << 231 z /= norm ;
221 z /= norm ; << 232 vF /= norm ;
222 vF /= norm ; << 233 lF /= norm ;
223 } <<
224 } 234 }
225 235
226 // Helium ion case 236 // Helium ion case
227 if( ionCharge < 2.5 ) { 237 if( ionCharge < 2.5 ) {
228 238
229 G4double e = std::log(std::max(1.0, kineti << 239 G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ;
230 G4double x = c[0] ; 240 G4double x = c[0] ;
231 G4double y = 1.0 ; 241 G4double y = 1.0 ;
232 for (G4int i=1; i<6; i++) { 242 for (G4int i=1; i<6; i++) {
233 y *= e ; 243 y *= e ;
234 x += y * c[i] ; 244 x += y * c[i] ;
235 } 245 }
236 G4double q = 7.6 - e ; << 246 G4double q = 7.6 - e ;
237 q = 1.0 + ( 0.007 + 0.00005 * z ) * G4Exp( << 247 q = 1.0 + ( 0.007 + 0.00005 * z ) * std::exp( -q*q ) ;
238 return 4.0 * q * q * (1.0 - G4Exp(-x)) ; << 248 return 4.0 * q * q * (1.0 - std::exp(-x)) ;
239 249
240 // Heavy ion case 250 // Heavy ion case
241 } else { 251 } else {
242 252
243 // v1 is ion velocity in vF unit 253 // v1 is ion velocity in vF unit
244 G4double v1{0.0}, v2{0.0}; << 254 G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ;
245 if (vF > 0.0) { <<
246 v1 = std::sqrt( reducedEnergy / (25.0 * <<
247 v2 = 1.0/ (vF*vF); <<
248 } <<
249 G4double y ; 255 G4double y ;
250 G4double z13 = std::pow(ionCharge, 0.3333) 256 G4double z13 = std::pow(ionCharge, 0.3333) ;
251 257
252 // Faster than Fermi velocity 258 // Faster than Fermi velocity
253 if ( v1 > 1.0 ) { 259 if ( v1 > 1.0 ) {
254 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / 260 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
255 261
256 // Slower than Fermi velocity 262 // Slower than Fermi velocity
257 } else { 263 } else {
258 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + 264 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ;
259 } 265 }
260 266
261 G4double y3 = std::pow(y, 0.3) ; 267 G4double y3 = std::pow(y, 0.3) ;
262 G4double q = 1.0 - G4Exp( 0.803*y3 - 1.316 << 268 G4double q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 -
263 0.38157*y - 0.0089 << 269 0.38157*y - 0.008983*y*y ) ;
264 if( q < 0.0 ) q = 0.0 ; 270 if( q < 0.0 ) q = 0.0 ;
265 271
266 G4double sLocal = 7.6 - std::log(std::max << 272 G4double s = 7.6 - std::log(std::max(1.0, reducedEnergy/keV)) ;
267 sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * G4E << 273 s = 1.0 + ( 0.18 + 0.0015 * z ) * std::exp( -s*s )/ (ionCharge*ionCharge) ;
268 274
269 // Screen length according to 275 // Screen length according to
270 // J.F.Ziegler and J.M.Manoyan, The stoppi 276 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
271 // Nucl. Inst. & Meth. in Phys. Res. B35 ( 277 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
272 278
273 G4double lambda = 10.0 * vF * std::pow(1.0 279 G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ;
274 G4double qeff = ionCharge * sLocal * << 280 G4double qeff = ionCharge * s *
275 ( q + 0.5*(1.0-q) * std::log(1.0 + lambd << 281 ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) / (vF*vF) ) ;
276 if( 0.1 > qeff ) qeff = 0.1 ; << 282 if( 0.1 > qeff ) qeff = 0.1 ;
277 return qeff*qeff ; << 283 return qeff*qeff ;
278 } 284 }
279 } 285 }
>> 286
>> 287
280 288