| Geant4 Cross Reference |
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: G4hIonEffChargeSquare
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/06/2001 V.Ivanchenko Continuation for eff.charge (small change of y)
41 // 08/10/2002 V.Ivanchenko The charge of the nucleus is used not charge of
42 // DynamicParticle
43 //
44 // Class Description:
45 //
46 // Ion effective charge model
47 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
49 //
50 // Class Description: End
51 //
52 // -------------------------------------------------------------------
53 //
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56 #include "G4hIonEffChargeSquare.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4DynamicParticle.hh"
60 #include "G4ParticleDefinition.hh"
61 #include "G4Material.hh"
62 #include "G4Element.hh"
63 #include "G4Exp.hh"
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
67 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name)
68 : G4VLowEnergyModel(name),
69 theHeMassAMU(4.0026)
70 {;}
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74 G4hIonEffChargeSquare::~G4hIonEffChargeSquare()
75 {;}
76
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78
79 G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle,
80 const G4Material* material)
81 {
82 G4double energy = particle->GetKineticEnergy() ;
83 G4double particleMass = particle->GetMass() ;
84 G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ;
85
86 G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ;
87
88 return q ;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
93 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle,
94 const G4Material* material,
95 G4double kineticEnergy)
96 {
97 // SetRateMass(aParticle) ;
98 G4double particleMass = aParticle->GetPDGMass() ;
99 G4double charge = (aParticle->GetPDGCharge())/eplus ;
100
101 G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ;
102
103 return q ;
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107
108 G4double G4hIonEffChargeSquare::HighEnergyLimit(
109 const G4ParticleDefinition* ,
110 const G4Material* ) const
111 {
112 return 1.0*TeV ;
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
117 G4double G4hIonEffChargeSquare::LowEnergyLimit(
118 const G4ParticleDefinition* ,
119 const G4Material* ) const
120 {
121 return 0.0 ;
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126 G4double G4hIonEffChargeSquare::HighEnergyLimit(
127 const G4ParticleDefinition* ) const
128 {
129 return 1.0*TeV ;
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133
134 G4double G4hIonEffChargeSquare::LowEnergyLimit(
135 const G4ParticleDefinition* ) const
136 {
137 return 0.0 ;
138 }
139
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141
142 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* ,
143 const G4Material* ) const
144 {
145 return true ;
146 }
147
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149
150 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* ,
151 const G4Material* ) const
152 {
153 return true ;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158 G4double G4hIonEffChargeSquare::IonEffChargeSquare(
159 const G4Material* material,
160 G4double kineticEnergy,
161 G4double particleMass,
162 G4double ionCharge) const
163 {
164 // The aproximation of ion effective charge from:
165 // J.F.Ziegler, J.P. Biersack, U. Littmark
166 // The Stopping and Range of Ions in Matter,
167 // Vol.1, Pergamon Press, 1985
168
169 // Fast ions or hadrons
170 G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ;
171 if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV;
172 if( (reducedEnergy > ionCharge * 10.0 * MeV) ||
173 (ionCharge < 1.5) ) return ionCharge*ionCharge ;
174
175 static const G4double vFermi[92] = {
176 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.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
178 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.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
180 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.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
182 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, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
184 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} ;
186
187 static const G4double c[6] = {0.2865, 0.1266, -0.001429,
188 0.02402,-0.01135, 0.001475} ;
189
190 // get elements in the actual material,
191 const G4ElementVector* theElementVector = material->GetElementVector() ;
192 const G4double* theAtomicNumDensityVector =
193 material->GetAtomicNumDensityVector() ;
194 const G4int NumberOfElements = (G4int)material->GetNumberOfElements() ;
195
196 // loop for the elements in the material
197 // to find out average values Z, vF, lF
198 G4double z = 0.0, vF = 0.0, norm = 0.0 ;
199
200 if( 1 == NumberOfElements ) {
201 z = material->GetZ() ;
202 G4int iz = G4int(z) - 1 ;
203 if(iz < 0) iz = 0 ;
204 else if(iz > 91) iz = 91 ;
205 vF = vFermi[iz] ;
206
207 } else {
208 for (G4int iel=0; iel<NumberOfElements; iel++)
209 {
210 const G4Element* element = (*theElementVector)[iel] ;
211 G4double z2 = element->GetZ() ;
212 const G4double weight = theAtomicNumDensityVector[iel] ;
213 norm += weight ;
214 z += z2 * weight ;
215 G4int iz = G4int(z2) - 1 ;
216 if(iz < 0) iz = 0 ;
217 else if(iz > 91) iz =91 ;
218 vF += vFermi[iz] * weight ;
219 }
220 if (norm > 0.0) {
221 z /= norm ;
222 vF /= norm ;
223 }
224 }
225
226 // Helium ion case
227 if( ionCharge < 2.5 ) {
228
229 G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ;
230 G4double x = c[0] ;
231 G4double y = 1.0 ;
232 for (G4int i=1; i<6; i++) {
233 y *= e ;
234 x += y * c[i] ;
235 }
236 G4double q = 7.6 - e ;
237 q = 1.0 + ( 0.007 + 0.00005 * z ) * G4Exp( -q*q ) ;
238 return 4.0 * q * q * (1.0 - G4Exp(-x)) ;
239
240 // Heavy ion case
241 } else {
242
243 // v1 is ion velocity in vF unit
244 G4double v1{0.0}, v2{0.0};
245 if (vF > 0.0) {
246 v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF;
247 v2 = 1.0/ (vF*vF);
248 }
249 G4double y ;
250 G4double z13 = std::pow(ionCharge, 0.3333) ;
251
252 // Faster than Fermi velocity
253 if ( v1 > 1.0 ) {
254 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
255
256 // Slower than Fermi velocity
257 } else {
258 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ;
259 }
260
261 G4double y3 = std::pow(y, 0.3) ;
262 G4double q = 1.0 - G4Exp( 0.803*y3 - 1.3167*y3*y3 -
263 0.38157*y - 0.008983*y*y ) ;
264 if( q < 0.0 ) q = 0.0 ;
265
266 G4double sLocal = 7.6 - std::log(std::max(1.0, reducedEnergy/keV)) ;
267 sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * G4Exp( -sLocal*sLocal )/ (ionCharge*ionCharge) ;
268
269 // Screen length according to
270 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
271 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
272
273 G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ;
274 G4double qeff = ionCharge * sLocal *
275 ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) * v2) ;
276 if( 0.1 > qeff ) qeff = 0.1 ;
277 return qeff*qeff ;
278 }
279 }
280