Geant4_10
G4hIonEffChargeSquare.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 //
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 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
67  : G4VLowEnergyModel(name),
68  theHeMassAMU(4.0026)
69 {;}
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
74 {;}
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79  const G4Material* material)
80 {
81  G4double energy = particle->GetKineticEnergy() ;
82  G4double particleMass = particle->GetMass() ;
83  G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ;
84 
85  G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ;
86 
87  return q ;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
93  const G4Material* material,
94  G4double kineticEnergy)
95 {
96  // SetRateMass(aParticle) ;
97  G4double particleMass = aParticle->GetPDGMass() ;
98  G4double charge = (aParticle->GetPDGCharge())/eplus ;
99 
100  G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ;
101 
102  return q ;
103 }
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
108  const G4ParticleDefinition* ,
109  const G4Material* ) const
110 {
111  return 1.0*TeV ;
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
117  const G4ParticleDefinition* ,
118  const G4Material* ) const
119 {
120  return 0.0 ;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
126  const G4ParticleDefinition* ) const
127 {
128  return 1.0*TeV ;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 
134  const G4ParticleDefinition* ) const
135 {
136  return 0.0 ;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 
142  const G4Material* ) const
143 {
144  return true ;
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
150  const G4Material* ) const
151 {
152  return true ;
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
157 G4double G4hIonEffChargeSquare::IonEffChargeSquare(
158  const G4Material* material,
159  G4double kineticEnergy,
160  G4double particleMass,
161  G4double ionCharge) const
162 {
163  // The aproximation of ion effective charge from:
164  // J.F.Ziegler, J.P. Biersack, U. Littmark
165  // The Stopping and Range of Ions in Matter,
166  // Vol.1, Pergamon Press, 1985
167 
168  // Fast ions or hadrons
169  G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ;
170  if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV;
171  if( (reducedEnergy > ionCharge * 10.0 * MeV) ||
172  (ionCharge < 1.5) ) return ionCharge*ionCharge ;
173 
174  static const G4double vFermi[92] = {
175  1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424,
176  0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
177  0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411,
178  0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
179  1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054,
180  0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
181  0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
182  0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
183  0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
184  1.9578, 1.0257} ;
185 
186  static const G4double lFactor[92] = {
187  1.0, 1.0, 1.1, 1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
188  0.82, 0.81, 0.83, 0.88, 1.0, 0.95, 0.97, 0.99, 0.98, 0.97,
189  0.98, 0.97, 0.96, 0.93, 0.91, 0.9, 0.88, 0.9, 0.9, 0.9,
190  0.9, 0.85, 0.9, 0.9, 0.91, 0.92, 0.9, 0.9, 0.9, 0.9,
191  0.9, 0.88, 0.9, 0.88, 0.88, 0.9, 0.9, 0.88, 0.9, 0.9,
192  0.9, 0.9, 0.96, 1.2, 0.9, 0.88, 0.88, 0.85, 0.9, 0.9,
193  0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1, 1.08, 1.08,
194  1.08, 1.08, 1.09, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15,
195  1.17, 1.2, 1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
196  1.16, 1.16} ;
197 
198  static const G4double c[6] = {0.2865, 0.1266, -0.001429,
199  0.02402,-0.01135, 0.001475} ;
200 
201  // get elements in the actual material,
202  const G4ElementVector* theElementVector = material->GetElementVector() ;
203  const G4double* theAtomicNumDensityVector =
204  material->GetAtomicNumDensityVector() ;
205  const G4int NumberOfElements = material->GetNumberOfElements() ;
206 
207  // loop for the elements in the material
208  // to find out average values Z, vF, lF
209  G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
210 
211  if( 1 == NumberOfElements ) {
212  z = material->GetZ() ;
213  G4int iz = G4int(z) - 1 ;
214  if(iz < 0) iz = 0 ;
215  else if(iz > 91) iz = 91 ;
216  vF = vFermi[iz] ;
217  lF = lFactor[iz] ;
218 
219  } else {
220  for (G4int iel=0; iel<NumberOfElements; iel++)
221  {
222  const G4Element* element = (*theElementVector)[iel] ;
223  G4double z2 = element->GetZ() ;
224  const G4double weight = theAtomicNumDensityVector[iel] ;
225  norm += weight ;
226  z += z2 * weight ;
227  G4int iz = G4int(z2) - 1 ;
228  if(iz < 0) iz = 0 ;
229  else if(iz > 91) iz =91 ;
230  vF += vFermi[iz] * weight ;
231  lF += lFactor[iz] * weight ;
232  }
233  z /= norm ;
234  vF /= norm ;
235  lF /= norm ;
236  }
237 
238  // Helium ion case
239  if( ionCharge < 2.5 ) {
240 
241  G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ;
242  G4double x = c[0] ;
243  G4double y = 1.0 ;
244  for (G4int i=1; i<6; i++) {
245  y *= e ;
246  x += y * c[i] ;
247  }
248  G4double q = 7.6 - e ;
249  q = 1.0 + ( 0.007 + 0.00005 * z ) * std::exp( -q*q ) ;
250  return 4.0 * q * q * (1.0 - std::exp(-x)) ;
251 
252  // Heavy ion case
253  } else {
254 
255  // v1 is ion velocity in vF unit
256  G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ;
257  G4double y ;
258  G4double z13 = std::pow(ionCharge, 0.3333) ;
259 
260  // Faster than Fermi velocity
261  if ( v1 > 1.0 ) {
262  y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
263 
264  // Slower than Fermi velocity
265  } else {
266  y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ;
267  }
268 
269  G4double y3 = std::pow(y, 0.3) ;
270  G4double q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 -
271  0.38157*y - 0.008983*y*y ) ;
272  if( q < 0.0 ) q = 0.0 ;
273 
274  G4double sLocal = 7.6 - std::log(std::max(1.0, reducedEnergy/keV)) ;
275  sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * std::exp( -sLocal*sLocal )/ (ionCharge*ionCharge) ;
276 
277  // Screen length according to
278  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
279  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
280 
281  G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ;
282  G4double qeff = ionCharge * sLocal *
283  ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) / (vF*vF) ) ;
284  if( 0.1 > qeff ) qeff = 0.1 ;
285  return qeff*qeff ;
286  }
287 }
288 
289 
G4double GetZ() const
Definition: G4Material.cc:606
std::vector< G4Element * > G4ElementVector
G4hIonEffChargeSquare(const G4String &name)
G4double GetKineticEnergy() const
Float_t norm
G4double GetZ() const
Definition: G4Element.hh:131
const XML_Char * name
Definition: expat.h:151
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
double weight
Definition: plottest35.C:25
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
Double_t y
Definition: plot.C:279
string material
Definition: eplot.py:19
G4double TheValue(const G4DynamicParticle *particle, const G4Material *material)
double energy
Definition: plottest35.C:25
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
G4double iz
Definition: TRTMaterials.hh:39
float proton_mass_c2
Definition: hepunit.py:275
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:214
G4bool IsInCharge(const G4DynamicParticle *particle, const G4Material *material) const
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
tuple z
Definition: test.py:28
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double GetPDGCharge() const
G4double HighEnergyLimit(const G4ParticleDefinition *aParticle, const G4Material *material) const
G4double LowEnergyLimit(const G4ParticleDefinition *aParticle, const G4Material *material) const