G4hIonEffChargeSquare.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // -------------------------------------------------------------------
00028 //
00029 // GEANT4 Class file
00030 //
00031 //
00032 // File name:     G4hIonEffChargeSquare
00033 //
00034 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
00035 // 
00036 // Creation date: 20 July 2000
00037 //
00038 // Modifications: 
00039 // 20/07/2000  V.Ivanchenko First implementation
00040 // 18/06/2001  V.Ivanchenko Continuation for eff.charge (small change of y)
00041 // 08/10/2002  V.Ivanchenko The charge of the nucleus is used not charge of 
00042 //                          DynamicParticle
00043 //
00044 // Class Description: 
00045 //
00046 // Ion effective charge model
00047 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
00048 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
00049 //
00050 // Class Description: End 
00051 //
00052 // -------------------------------------------------------------------
00053 //
00054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00055 
00056 #include "G4hIonEffChargeSquare.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4DynamicParticle.hh"
00060 #include "G4ParticleDefinition.hh"
00061 #include "G4Material.hh"
00062 #include "G4Element.hh"
00063 
00064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00065 
00066 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name)
00067   : G4VLowEnergyModel(name), 
00068     theHeMassAMU(4.0026)
00069 {;}
00070 
00071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00072 
00073 G4hIonEffChargeSquare::~G4hIonEffChargeSquare() 
00074 {;}
00075 
00076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00077 
00078 G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle,
00079                                          const G4Material* material) 
00080 {
00081   G4double energy = particle->GetKineticEnergy() ;
00082   G4double particleMass = particle->GetMass() ;
00083   G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ;
00084 
00085   G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ;
00086 
00087   return q ;
00088 }
00089 
00090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00091 
00092 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle,
00093                                          const G4Material* material,
00094                                          G4double kineticEnergy) 
00095 {
00096   //  SetRateMass(aParticle) ;
00097   G4double particleMass = aParticle->GetPDGMass() ;
00098   G4double charge = (aParticle->GetPDGCharge())/eplus ;
00099 
00100   G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ;
00101 
00102   return q ;
00103 }
00104 
00105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00106 
00107 G4double G4hIonEffChargeSquare::HighEnergyLimit(
00108                                                 const G4ParticleDefinition* ,
00109                                                 const G4Material* ) const
00110 {
00111   return 1.0*TeV ;
00112 }
00113 
00114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00115 
00116 G4double G4hIonEffChargeSquare::LowEnergyLimit(
00117                                                const G4ParticleDefinition* ,
00118                                                const G4Material* ) const
00119 {
00120   return 0.0 ;
00121 }
00122 
00123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00124 
00125 G4double G4hIonEffChargeSquare::HighEnergyLimit(
00126                                                 const G4ParticleDefinition* ) const
00127 {
00128   return 1.0*TeV ;
00129 }
00130 
00131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00132 
00133 G4double G4hIonEffChargeSquare::LowEnergyLimit(
00134                                                const G4ParticleDefinition* ) const
00135 {
00136   return 0.0 ;
00137 }
00138 
00139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00140  
00141 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* ,
00142                                          const G4Material* ) const
00143 {
00144   return true ;
00145 }
00146 
00147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00148  
00149 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* ,
00150                                          const G4Material* ) const
00151 {
00152   return true ;
00153 }
00154 
00155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00156 
00157 G4double G4hIonEffChargeSquare::IonEffChargeSquare(
00158                                                    const G4Material* material,
00159                                                    G4double kineticEnergy,
00160                                                    G4double particleMass,
00161                                                    G4double ionCharge) const
00162 {
00163   // The aproximation of ion effective charge from: 
00164   // J.F.Ziegler, J.P. Biersack, U. Littmark
00165   // The Stopping and Range of Ions in Matter,
00166   // Vol.1, Pergamon Press, 1985
00167 
00168   // Fast ions or hadrons
00169   G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ;
00170   if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV;
00171   if( (reducedEnergy > ionCharge * 10.0 * MeV) || 
00172       (ionCharge < 1.5) ) return ionCharge*ionCharge ;
00173 
00174   static G4double vFermi[92] = {
00175     1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
00176     0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
00177     0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
00178     0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
00179     1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
00180     0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
00181     0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
00182     0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
00183     0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
00184     1.9578,  1.0257} ;
00185 
00186   static G4double lFactor[92] = {
00187     1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
00188     0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
00189     0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
00190     0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
00191     0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
00192     0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9, 
00193     0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
00194     1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
00195     1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
00196     1.16, 1.16} ; 
00197 
00198   static G4double c[6] = {0.2865,  0.1266, -0.001429,
00199                           0.02402,-0.01135, 0.001475} ;
00200 
00201   // get elements in the actual material,
00202   const G4ElementVector* theElementVector = material->GetElementVector() ;
00203   const G4double* theAtomicNumDensityVector = 
00204                          material->GetAtomicNumDensityVector() ;
00205   const G4int NumberOfElements = material->GetNumberOfElements() ;
00206   
00207   //  loop for the elements in the material
00208   //  to find out average values Z, vF, lF
00209   G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ; 
00210 
00211   if( 1 == NumberOfElements ) {
00212     z = material->GetZ() ;
00213     G4int iz = G4int(z) - 1 ;
00214     if(iz < 0) iz = 0 ;
00215     else if(iz > 91) iz = 91 ;
00216     vF   = vFermi[iz] ;
00217     lF   = lFactor[iz] ;
00218 
00219   } else {
00220     for (G4int iel=0; iel<NumberOfElements; iel++)
00221       {
00222         const G4Element* element = (*theElementVector)[iel] ;
00223         G4double z2 = element->GetZ() ;
00224         const G4double weight = theAtomicNumDensityVector[iel] ;
00225         norm += weight ;
00226         z    += z2 * weight ;
00227         G4int iz = G4int(z2) - 1 ;
00228         if(iz < 0) iz = 0 ;
00229         else if(iz > 91) iz =91 ;
00230         vF   += vFermi[iz] * weight ;
00231         lF   += lFactor[iz] * weight ;
00232       }
00233     z  /= norm ;
00234     vF /= norm ;
00235     lF /= norm ;
00236   }
00237 
00238   // Helium ion case
00239   if( ionCharge < 2.5 ) {
00240 
00241     G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ; 
00242     G4double x = c[0] ;
00243     G4double y = 1.0 ;
00244     for (G4int i=1; i<6; i++) {
00245       y *= e ;
00246       x += y * c[i] ;
00247     }
00248     G4double q = 7.6 -  e ; 
00249     q = 1.0 + ( 0.007 + 0.00005 * z ) * std::exp( -q*q ) ;
00250     return  4.0 * q * q * (1.0 - std::exp(-x)) ;
00251 
00252     // Heavy ion case
00253   } else {
00254 
00255     // v1 is ion velocity in vF unit
00256     G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ;
00257     G4double y ;
00258     G4double z13 = std::pow(ionCharge, 0.3333) ;
00259 
00260     // Faster than Fermi velocity
00261     if ( v1 > 1.0 ) {
00262       y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
00263 
00264       // Slower than Fermi velocity
00265     } else {
00266       y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ;
00267     }
00268 
00269     G4double y3 = std::pow(y, 0.3) ;
00270     G4double q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 - 
00271                             0.38157*y - 0.008983*y*y ) ;     
00272     if( q < 0.0 ) q = 0.0 ;
00273 
00274     G4double sLocal = 7.6 -  std::log(std::max(1.0, reducedEnergy/keV)) ; 
00275     sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * std::exp( -sLocal*sLocal )/ (ionCharge*ionCharge) ;
00276 
00277     // Screen length according to
00278     // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
00279     // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
00280 
00281     G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ;
00282     G4double qeff   = ionCharge * sLocal *
00283       ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) / (vF*vF) ) ;
00284     if( 0.1 > qeff ) qeff = 0.1 ; 
00285     return qeff*qeff ;    
00286   }
00287 }
00288 
00289 

Generated on Mon May 27 17:48:32 2013 for Geant4 by  doxygen 1.4.7