00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
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
00065
00066 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name)
00067 : G4VLowEnergyModel(name),
00068 theHeMassAMU(4.0026)
00069 {;}
00070
00071
00072
00073 G4hIonEffChargeSquare::~G4hIonEffChargeSquare()
00074 {;}
00075
00076
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
00091
00092 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle,
00093 const G4Material* material,
00094 G4double kineticEnergy)
00095 {
00096
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
00106
00107 G4double G4hIonEffChargeSquare::HighEnergyLimit(
00108 const G4ParticleDefinition* ,
00109 const G4Material* ) const
00110 {
00111 return 1.0*TeV ;
00112 }
00113
00114
00115
00116 G4double G4hIonEffChargeSquare::LowEnergyLimit(
00117 const G4ParticleDefinition* ,
00118 const G4Material* ) const
00119 {
00120 return 0.0 ;
00121 }
00122
00123
00124
00125 G4double G4hIonEffChargeSquare::HighEnergyLimit(
00126 const G4ParticleDefinition* ) const
00127 {
00128 return 1.0*TeV ;
00129 }
00130
00131
00132
00133 G4double G4hIonEffChargeSquare::LowEnergyLimit(
00134 const G4ParticleDefinition* ) const
00135 {
00136 return 0.0 ;
00137 }
00138
00139
00140
00141 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* ,
00142 const G4Material* ) const
00143 {
00144 return true ;
00145 }
00146
00147
00148
00149 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* ,
00150 const G4Material* ) const
00151 {
00152 return true ;
00153 }
00154
00155
00156
00157 G4double G4hIonEffChargeSquare::IonEffChargeSquare(
00158 const G4Material* material,
00159 G4double kineticEnergy,
00160 G4double particleMass,
00161 G4double ionCharge) const
00162 {
00163
00164
00165
00166
00167
00168
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
00202 const G4ElementVector* theElementVector = material->GetElementVector() ;
00203 const G4double* theAtomicNumDensityVector =
00204 material->GetAtomicNumDensityVector() ;
00205 const G4int NumberOfElements = material->GetNumberOfElements() ;
00206
00207
00208
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
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
00253 } else {
00254
00255
00256 G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ;
00257 G4double y ;
00258 G4double z13 = std::pow(ionCharge, 0.3333) ;
00259
00260
00261 if ( v1 > 1.0 ) {
00262 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ;
00263
00264
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
00278
00279
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