Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ionEffectiveCharge.cc
Go to the documentation of this file.
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 // $Id$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4ionEffectiveCharge
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 07.05.2002
38 //
39 // Modifications:
40 // 12.09.2004 Set low energy limit to 1 keV (V.Ivanchenko)
41 // 25.01.2005 Add protection - min Charge 0.1 eplus (V.Ivanchenko)
42 // 28.04.2006 Set upper energy limit to 50 MeV (V.Ivanchenko)
43 // 23.05.2006 Set upper energy limit to Z*10 MeV (V.Ivanchenko)
44 // 15.08.2006 Add protection for not defined material (V.Ivanchenko)
45 // 27-09-2007 Use Fermi energy from material, optimazed formulas (V.Ivanchenko)
46 //
47 
48 // -------------------------------------------------------------------
49 //
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 
53 #include "G4ionEffectiveCharge.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4UnitsTable.hh"
57 #include "G4Material.hh"
58 #include "G4NistManager.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
63 {
64  chargeCorrection = 1.0;
65  energyHighLimit = 20.0*MeV;
66  energyLowLimit = 1.0*keV;
67  energyBohr = 25.*keV;
68  massFactor = amu_c2/(proton_mass_c2*keV);
69  minCharge = 1.0;
70  lastPart = 0;
71  lastMat = 0;
72  lastKinEnergy = 0.0;
73  effCharge = eplus;
74  nist = G4NistManager::Instance();
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
80 {}
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
85  const G4Material* material,
86  G4double kineticEnergy)
87 {
88  if(p == lastPart && material == lastMat && kineticEnergy == lastKinEnergy)
89  return effCharge;
90 
91  lastPart = p;
92  lastMat = material;
93  lastKinEnergy = kineticEnergy;
94 
95  G4double mass = p->GetPDGMass();
96  G4double charge = p->GetPDGCharge();
97  G4double Zi = charge/eplus;
98 
99  chargeCorrection = 1.0;
100  effCharge = charge;
101 
102  // The aproximation of ion effective charge from:
103  // J.F.Ziegler, J.P. Biersack, U. Littmark
104  // The Stopping and Range of Ions in Matter,
105  // Vol.1, Pergamon Press, 1985
106  // Fast ions or hadrons
107  G4double reducedEnergy = kineticEnergy * proton_mass_c2/mass ;
108 
109  //G4cout << "e= " << reducedEnergy << " Zi= " << Zi << " " << material->GetName() << G4endl;
110 
111  if( reducedEnergy > Zi*energyHighLimit || Zi < 1.5 || !material) return charge;
112 
113  G4double z = material->GetIonisation()->GetZeffective();
114  reducedEnergy = std::max(reducedEnergy,energyLowLimit);
115 
116  // Helium ion case
117  if( Zi < 2.5 ) {
118 
119  static G4double c[6] = {0.2865, 0.1266, -0.001429,
120  0.02402,-0.01135, 0.001475} ;
121 
122  G4double Q = std::max(0.0,std::log(reducedEnergy*massFactor));
123  G4double x = c[0];
124  G4double y = 1.0;
125  for (G4int i=1; i<6; i++) {
126  y *= Q;
127  x += y * c[i] ;
128  }
129  G4double ex;
130  if(x < 0.2) ex = x * (1 - 0.5*x);
131  else ex = 1. - std::exp(-x);
132 
133  G4double tq = 7.6 - Q;
134  G4double tq2= tq*tq;
135  G4double tt = ( 0.007 + 0.00005 * z );
136  if(tq2 < 0.2) tt *= (1.0 - tq2 + 0.5*tq2*tq2);
137  else tt *= std::exp(-tq2);
138 
139  effCharge = charge*(1.0 + tt) * std::sqrt(ex);
140 
141  // Heavy ion case
142  } else {
143 
144  G4double y;
145  // = nist->GetZ13(z);
146  //G4double z23 = y*y;
147  G4double zi13 = nist->GetZ13(Zi);
148  G4double zi23 = zi13*zi13;
149  // G4double e = std::max(reducedEnergy,energyBohr/z23);
150  //G4double e = reducedEnergy;
151 
152  // v1 is ion velocity in vF unit
153  G4double eF = material->GetIonisation()->GetFermiEnergy();
154  G4double v1sq = reducedEnergy/eF;
155  G4double vFsq = eF/energyBohr;
156  G4double vF = std::sqrt(eF/energyBohr);
157 
158  // Faster than Fermi velocity
159  if ( v1sq > 1.0 ) {
160  y = vF * std::sqrt(v1sq) * ( 1.0 + 0.2/v1sq ) / zi23 ;
161 
162  // Slower than Fermi velocity
163  } else {
164  y = 0.692308 * vF * (1.0 + 0.666666*v1sq + v1sq*v1sq/15.0) / zi23 ;
165  }
166 
167  G4double q;
168  G4double y3 = std::pow(y, 0.3) ;
169  // G4cout<<"y= "<<y<<" y3= "<<y3<<" v1= "<<v1<<" vF= "<<vF<<G4endl;
170  q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 - 0.38157*y - 0.008983*y*y ) ;
171 
172  //y *= 0.77;
173  //y *= (0.75 + 0.52/Zi);
174 
175  //if( y < 0.2 ) q = y*(1.0 - 0.5*y);
176  //else q = 1.0 - std::exp(-y);
177 
178  G4double qmin = minCharge/Zi;
179  if(q < qmin) q = qmin;
180 
181  effCharge = q*charge;
182 
183  /*
184  G4double x1 = 1.0*effCharge*(1.0 - 0.132*std::log(y))/(y*std::sqrt(z));
185  G4double x2 = 0.1*effCharge*effCharge*energyBohr/reducedEnergy;
186 
187  chargeCorrection = 1.0 + x1 - x2;
188 
189  G4cout << "x1= "<<x1<<" x2= "<< x2<<" corr= "<<chargeCorrection<<G4endl;
190  */
191 
192  G4double tq = 7.6 - std::log(reducedEnergy/keV);
193  G4double tq2= tq*tq;
194  G4double sq = ( 0.18 + 0.0015 * z ) / (Zi*Zi);
195  if(tq2 < 0.2) sq *= (1.0 - tq2 + 0.5*tq2*tq2);
196  else sq *= std::exp(-tq2);
197  sq += 1.0;
198  // G4cout << "sq= " << sq << G4endl;
199 
200  // Screen length according to
201  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
202  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
203 
204  G4double lambda = 10.0 * vF / (zi13 * (6.0 + q));
205  if(q < 0.2) lambda *= (1.0 - 0.66666667*q - q*q/9.0);
206  else lambda *= std::pow(1.0-q, 0.666666);
207 
208  G4double lambda2 = lambda*lambda;
209 
210  G4double xx = (0.5/q - 0.5)/vFsq;
211  if(lambda2 < 0.2) xx *= lambda2*(1.0 - 0.5*lambda2);
212  else xx *= std::log(1.0 + lambda2);
213 
214  chargeCorrection = sq * (1.0 + xx);
215 
216  }
217  // G4cout << "G4ionEffectiveCharge: charge= " << charge << " q= " << q
218  // << " chargeCor= " << chargeCorrection
219  // << " e(MeV)= " << kineticEnergy/MeV << G4endl;
220  return effCharge;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224 
225