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