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

Last change on this file since 1007 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 10.2 KB
Line 
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.