source: trunk/source/processes/electromagnetic/lowenergy/src/G4hIonEffChargeSquare.cc@ 1192

Last change on this file since 1192 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 10.2 KB
RevLine 
[819]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 "G4DynamicParticle.hh"
58#include "G4ParticleDefinition.hh"
59#include "G4Material.hh"
60#include "G4Element.hh"
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name)
65 : G4VLowEnergyModel(name),
66 theHeMassAMU(4.0026)
67{;}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71G4hIonEffChargeSquare::~G4hIonEffChargeSquare()
72{;}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle,
77 const G4Material* material)
78{
79 G4double energy = particle->GetKineticEnergy() ;
80 G4double particleMass = particle->GetMass() ;
81 G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ;
82
83 G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ;
84
85 return q ;
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle,
91 const G4Material* material,
92 G4double kineticEnergy)
93{
94 // SetRateMass(aParticle) ;
95 G4double particleMass = aParticle->GetPDGMass() ;
96 G4double charge = (aParticle->GetPDGCharge())/eplus ;
97
98 G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ;
99
100 return q ;
101}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105G4double G4hIonEffChargeSquare::HighEnergyLimit(
106 const G4ParticleDefinition* ,
107 const G4Material* ) const
108{
109 return 1.0*TeV ;
110}
111
112//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113
114G4double G4hIonEffChargeSquare::LowEnergyLimit(
115 const G4ParticleDefinition* ,
116 const G4Material* ) const
117{
118 return 0.0 ;
119}
120
121//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122
123G4double G4hIonEffChargeSquare::HighEnergyLimit(
124 const G4ParticleDefinition* ) const
125{
126 return 1.0*TeV ;
127}
128
129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130
131G4double G4hIonEffChargeSquare::LowEnergyLimit(
132 const G4ParticleDefinition* ) const
133{
134 return 0.0 ;
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* ,
140 const G4Material* ) const
141{
142 return true ;
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146
147G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* ,
148 const G4Material* ) const
149{
150 return true ;
151}
152
153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154
155G4double G4hIonEffChargeSquare::IonEffChargeSquare(
156 const G4Material* material,
157 G4double kineticEnergy,
158 G4double particleMass,
159 G4double ionCharge) const
160{
161 // The aproximation of ion effective charge from:
162 // J.F.Ziegler, J.P. Biersack, U. Littmark
163 // The Stopping and Range of Ions in Matter,
164 // Vol.1, Pergamon Press, 1985
165
166 // Fast ions or hadrons
167 G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ;
168 if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV;
169 if( (reducedEnergy > ionCharge * 10.0 * MeV) ||
170 (ionCharge < 1.5) ) return ionCharge*ionCharge ;
171
172 static G4double vFermi[92] = {
173 1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424,
174 0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
175 0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411,
176 0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
177 1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054,
178 0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
179 0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
180 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
181 0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
182 1.9578, 1.0257} ;
183
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,
197 0.02402,-0.01135, 0.001475} ;
198
199 // get elements in the actual material,
200 const G4ElementVector* theElementVector = material->GetElementVector() ;
201 const G4double* theAtomicNumDensityVector =
202 material->GetAtomicNumDensityVector() ;
203 const G4int NumberOfElements = material->GetNumberOfElements() ;
204
205 // loop for the elements in the material
206 // to find out average values Z, vF, lF
207 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
208
209 if( 1 == NumberOfElements ) {
210 z = material->GetZ() ;
211 G4int iz = G4int(z) - 1 ;
212 if(iz < 0) iz = 0 ;
213 else if(iz > 91) iz = 91 ;
214 vF = vFermi[iz] ;
215 lF = lFactor[iz] ;
216
217 } else {
218 for (G4int iel=0; iel<NumberOfElements; iel++)
219 {
220 const G4Element* element = (*theElementVector)[iel] ;
221 G4double z2 = element->GetZ() ;
222 const G4double weight = theAtomicNumDensityVector[iel] ;
223 norm += weight ;
224 z += z2 * weight ;
225 G4int iz = G4int(z2) - 1 ;
226 if(iz < 0) iz = 0 ;
227 else if(iz > 91) iz =91 ;
228 vF += vFermi[iz] * weight ;
229 lF += lFactor[iz] * weight ;
230 }
231 z /= norm ;
232 vF /= norm ;
233 lF /= norm ;
234 }
235
236 // Helium ion case
237 if( ionCharge < 2.5 ) {
238
239 G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ;
240 G4double x = c[0] ;
241 G4double y = 1.0 ;
242 for (G4int i=1; i<6; i++) {
243 y *= e ;
244 x += y * c[i] ;
245 }
246 G4double q = 7.6 - e ;
247 q = 1.0 + ( 0.007 + 0.00005 * z ) * std::exp( -q*q ) ;
248 return 4.0 * q * q * (1.0 - std::exp(-x)) ;
249
250 // Heavy ion case
251 } else {
252
253 // v1 is ion velocity in vF unit
254 G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ;
255 G4double y ;
256 G4double z13 = std::pow(ionCharge, 0.3333) ;
257
258 // Faster than Fermi velocity
259 if ( v1 > 1.0 ) {
260 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
261
262 // Slower than Fermi velocity
263 } else {
264 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ;
265 }
266
267 G4double y3 = std::pow(y, 0.3) ;
268 G4double q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 -
269 0.38157*y - 0.008983*y*y ) ;
270 if( q < 0.0 ) q = 0.0 ;
271
272 G4double s = 7.6 - std::log(std::max(1.0, reducedEnergy/keV)) ;
273 s = 1.0 + ( 0.18 + 0.0015 * z ) * std::exp( -s*s )/ (ionCharge*ionCharge) ;
274
275 // Screen length according to
276 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
277 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
278
279 G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ;
280 G4double qeff = ionCharge * s *
281 ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) / (vF*vF) ) ;
282 if( 0.1 > qeff ) qeff = 0.1 ;
283 return qeff*qeff ;
284 }
285}
286
287
Note: See TracBrowser for help on using the repository browser.