#include <G4WentzelOKandVIxSection.hh>
Public Member Functions | |
G4WentzelOKandVIxSection () | |
virtual | ~G4WentzelOKandVIxSection () |
void | Initialise (const G4ParticleDefinition *, G4double CosThetaLim) |
void | SetupParticle (const G4ParticleDefinition *) |
G4double | SetupTarget (G4int Z, G4double cut=DBL_MAX) |
G4double | ComputeTransportCrossSectionPerAtom (G4double CosThetaMax) |
G4ThreeVector | SampleSingleScattering (G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio=0.0) |
G4double | ComputeNuclearCrossSection (G4double CosThetaMin, G4double CosThetaMax) |
G4double | ComputeElectronCrossSection (G4double CosThetaMin, G4double CosThetaMax) |
G4double | SetupKinematic (G4double kinEnergy, const G4Material *mat) |
void | SetTargetMass (G4double value) |
G4double | GetMomentumSquare () const |
G4double | GetCosThetaNuc () const |
G4double | GetCosThetaElec () const |
Definition at line 71 of file G4WentzelOKandVIxSection.hh.
G4WentzelOKandVIxSection::G4WentzelOKandVIxSection | ( | ) |
Definition at line 66 of file G4WentzelOKandVIxSection.cc.
References DBL_MAX, G4Electron::Electron(), G4NistManager::GetA27(), G4Pow::GetInstance(), G4NistManager::Instance(), G4INCL::Math::pi, G4Positron::Positron(), G4Proton::Proton(), and G4Pow::Z13().
00066 : 00067 numlimit(0.1), 00068 nwarnings(0), 00069 nwarnlimit(50), 00070 alpha2(fine_structure_const*fine_structure_const) 00071 { 00072 fNistManager = G4NistManager::Instance(); 00073 fG4pow = G4Pow::GetInstance(); 00074 theElectron = G4Electron::Electron(); 00075 thePositron = G4Positron::Positron(); 00076 theProton = G4Proton::Proton(); 00077 lowEnergyLimit = 1.0*eV; 00078 G4double p0 = electron_mass_c2*classic_electr_radius; 00079 coeff = twopi*p0*p0; 00080 particle = 0; 00081 00082 // Thomas-Fermi screening radii 00083 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282 00084 00085 if(0.0 == ScreenRSquare[0]) { 00086 G4double a0 = electron_mass_c2/0.88534; 00087 G4double constn = 6.937e-6/(MeV*MeV); 00088 00089 ScreenRSquare[0] = alpha2*a0*a0; 00090 ScreenRSquareElec[0] = ScreenRSquare[0]; 00091 for(G4int j=1; j<100; ++j) { 00092 G4double x = a0*fG4pow->Z13(j); 00093 if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; } 00094 else { 00095 ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x; 00096 ScreenRSquareElec[j] = 0.5*alpha2*x*x; 00097 } 00098 x = fNistManager->GetA27(j); 00099 FormFactor[j] = constn*x*x; 00100 } 00101 } 00102 currentMaterial = 0; 00103 elecXSRatio = factB = factD = formfactA = screenZ = 0.0; 00104 cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0; 00105 00106 factB1= 0.5*CLHEP::pi*fine_structure_const; 00107 00108 tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0; 00109 ecut = etag = DBL_MAX; 00110 targetZ = 0; 00111 cosThetaMax = 1.0; 00112 targetMass = proton_mass_c2; 00113 particle = 0; 00114 }
G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection | ( | ) | [virtual] |
G4double G4WentzelOKandVIxSection::ComputeElectronCrossSection | ( | G4double | CosThetaMin, | |
G4double | CosThetaMax | |||
) | [inline] |
Definition at line 237 of file G4WentzelOKandVIxSection.hh.
Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom().
00239 { 00240 G4double xsec = 0.0; 00241 G4double cost1 = std::max(cosTMin,cosTetMaxElec); 00242 G4double cost2 = std::max(cosTMax,cosTetMaxElec); 00243 if(cost1 > cost2) { 00244 xsec = kinFactor*(cost1 - cost2)/((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ)); 00245 } 00246 return xsec; 00247 }
G4double G4WentzelOKandVIxSection::ComputeNuclearCrossSection | ( | G4double | CosThetaMin, | |
G4double | CosThetaMax | |||
) | [inline] |
Definition at line 223 of file G4WentzelOKandVIxSection.hh.
Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom().
00225 { 00226 G4double xsec = 0.0; 00227 if(cosTMax < cosTMin) { 00228 xsec = targetZ*kinFactor*(cosTMin - cosTMax)/ 00229 ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ)); 00230 } 00231 return xsec; 00232 }
Definition at line 202 of file G4WentzelOKandVIxSection.cc.
References G4cout, G4endl, and G4ParticleDefinition::GetParticleName().
Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom().
00203 { 00204 G4double xsec = 0.0; 00205 if(cosTMax >= 1.0) { return xsec; } 00206 00207 G4double xSection = 0.0; 00208 G4double x = 0; 00209 G4double y = 0; 00210 G4double x1= 0; 00211 G4double x2= 0; 00212 G4double xlog = 0.0; 00213 00214 G4double costm = std::max(cosTMax,cosTetMaxElec); 00215 G4double fb = screenZ*factB; 00216 00217 // scattering off electrons 00218 if(costm < 1.0) { 00219 x = (1.0 - costm)/screenZ; 00220 if(x < numlimit) { 00221 x2 = 0.5*x*x; 00222 y = x2*(1.0 - 1.3333333*x + 3*x2); 00223 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } 00224 } else { 00225 x1= x/(1 + x); 00226 xlog = log(1.0 + x); 00227 y = xlog - x1; 00228 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } 00229 } 00230 00231 if(y < 0.0) { 00232 ++nwarnings; 00233 if(nwarnings < nwarnlimit) { 00234 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" 00235 << G4endl; 00236 G4cout << "y= " << y 00237 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 00238 << " Z= " << targetZ << " " 00239 << particle->GetParticleName() << G4endl; 00240 G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 00241 << " x= " << x << G4endl; 00242 } 00243 y = 0.0; 00244 } 00245 xSection = y; 00246 } 00247 /* 00248 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 00249 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection 00250 << " cut(MeV)= " << ecut/MeV 00251 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 00252 << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 00253 << " 1-costm= " << 1.0 - cosThetaMax << G4endl; 00254 */ 00255 // scattering off nucleus 00256 if(cosTMax < 1.0) { 00257 x = (1.0 - cosTMax)/screenZ; 00258 if(x < numlimit) { 00259 x2 = 0.5*x*x; 00260 y = x2*(1.0 - 1.3333333*x + 3*x2); 00261 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); } 00262 } else { 00263 x1= x/(1 + x); 00264 xlog = log(1.0 + x); 00265 y = xlog - x1; 00266 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); } 00267 } 00268 00269 if(y < 0.0) { 00270 ++nwarnings; 00271 if(nwarnings < nwarnlimit) { 00272 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0" 00273 << G4endl; 00274 G4cout << "y= " << y 00275 << " e(MeV)= " << tkin << " Z= " << targetZ << " " 00276 << particle->GetParticleName() << G4endl; 00277 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 00278 << " x= " << " x1= " << x1 <<G4endl; 00279 } 00280 y = 0.0; 00281 } 00282 xSection += y*targetZ; 00283 } 00284 xSection *= kinFactor; 00285 /* 00286 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 00287 << " screenZ= " << screenZ << " formF= " << formfactA 00288 << " for " << particle->GetParticleName() 00289 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2) 00290 << " x= " << x 00291 << G4endl; 00292 */ 00293 return xSection; 00294 }
G4double G4WentzelOKandVIxSection::GetCosThetaElec | ( | ) | const [inline] |
G4double G4WentzelOKandVIxSection::GetCosThetaNuc | ( | ) | const [inline] |
G4double G4WentzelOKandVIxSection::GetMomentumSquare | ( | ) | const [inline] |
Definition at line 201 of file G4WentzelOKandVIxSection.hh.
Referenced by G4eCoulombScatteringModel::SampleSecondaries().
void G4WentzelOKandVIxSection::Initialise | ( | const G4ParticleDefinition * | , | |
G4double | CosThetaLim | |||
) |
Definition at line 123 of file G4WentzelOKandVIxSection.cc.
References DBL_MAX, G4LossTableManager::FactorForAngleLimit(), G4LossTableManager::Instance(), and SetupParticle().
Referenced by G4WentzelVIModel::Initialise(), and G4eCoulombScatteringModel::Initialise().
00125 { 00126 SetupParticle(p); 00127 tkin = mom2 = momCM2 = 0.0; 00128 ecut = etag = DBL_MAX; 00129 targetZ = 0; 00130 cosThetaMax = CosThetaLim; 00131 G4double a = G4LossTableManager::Instance()->FactorForAngleLimit() 00132 *CLHEP::hbarc/CLHEP::fermi; 00133 factorA2 = 0.5*a*a; 00134 currentMaterial = 0; 00135 00136 //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass 00137 // << " " << p->GetParticleName() << G4endl; 00138 00139 }
G4ThreeVector G4WentzelOKandVIxSection::SampleSingleScattering | ( | G4double | CosThetaMin, | |
G4double | CosThetaMax, | |||
G4double | elecRatio = 0.0 | |||
) |
Definition at line 299 of file G4WentzelOKandVIxSection.cc.
References G4UniformRand.
Referenced by G4WentzelVIModel::SampleScattering(), and G4eCoulombScatteringModel::SampleSecondaries().
00302 { 00303 G4ThreeVector v(0.0,0.0,1.0); 00304 00305 G4double formf = formfactA; 00306 G4double cost1 = cosTMin; 00307 G4double cost2 = cosTMax; 00308 if(elecRatio > 0.0) { 00309 if(G4UniformRand() <= elecRatio) { 00310 formf = 0.0; 00311 cost1 = std::max(cost1,cosTetMaxElec); 00312 cost2 = std::max(cost2,cosTetMaxElec); 00313 } 00314 } 00315 if(cost1 < cost2) { return v; } 00316 00317 G4double w1 = 1. - cost1 + screenZ; 00318 G4double w2 = 1. - cost2 + screenZ; 00319 G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ; 00320 00321 //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass); 00322 G4double fm = 1.0 + formf*z1; 00323 //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm ); 00324 G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm ); 00325 // "false" scattering 00326 if( G4UniformRand() > grej ) { return v; } 00327 // } 00328 G4double cost = 1.0 - z1; 00329 00330 if(cost > 1.0) { cost = 1.0; } 00331 else if(cost < -1.0) { cost =-1.0; } 00332 G4double sint = sqrt((1.0 - cost)*(1.0 + cost)); 00333 //G4cout << "sint= " << sint << G4endl; 00334 G4double phi = twopi*G4UniformRand(); 00335 G4double vx1 = sint*cos(phi); 00336 G4double vy1 = sint*sin(phi); 00337 00338 // only direction is changed 00339 v.set(vx1,vy1,cost); 00340 return v; 00341 }
void G4WentzelOKandVIxSection::SetTargetMass | ( | G4double | value | ) | [inline] |
Definition at line 193 of file G4WentzelOKandVIxSection.hh.
Referenced by G4eCoulombScatteringModel::SampleSecondaries(), and SetupTarget().
G4double G4WentzelOKandVIxSection::SetupKinematic | ( | G4double | kinEnergy, | |
const G4Material * | mat | |||
) | [inline] |
Definition at line 177 of file G4WentzelOKandVIxSection.hh.
References G4IonisParamMat::GetInvA23(), and G4Material::GetIonisation().
Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeGeomPathLength(), G4WentzelVIModel::ComputeTruePathLengthLimit(), and G4WentzelVIModel::ComputeTrueStepLength().
00178 { 00179 if(ekin != tkin || mat != currentMaterial) { 00180 currentMaterial = mat; 00181 tkin = ekin; 00182 mom2 = tkin*(tkin + 2.0*mass); 00183 invbeta2 = 1.0 + mass*mass/mom2; 00184 factB = spin/invbeta2; 00185 cosTetMaxNuc = 00186 std::max(cosThetaMax,1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2); 00187 } 00188 return cosTetMaxNuc; 00189 }
void G4WentzelOKandVIxSection::SetupParticle | ( | const G4ParticleDefinition * | ) |
Definition at line 143 of file G4WentzelOKandVIxSection.cc.
References G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().
Referenced by Initialise(), and G4eCoulombScatteringModel::SetupParticle().
00144 { 00145 particle = p; 00146 mass = particle->GetPDGMass(); 00147 spin = particle->GetPDGSpin(); 00148 if(0.0 != spin) { spin = 0.5; } 00149 G4double q = std::fabs(particle->GetPDGCharge()/eplus); 00150 chargeSquare = q*q; 00151 charge3 = chargeSquare*q; 00152 tkin = 0.0; 00153 currentMaterial = 0; 00154 targetZ = 0; 00155 }
Definition at line 160 of file G4WentzelOKandVIxSection.cc.
References G4NistManager::GetAtomicMassAmu(), SetTargetMass(), and G4Pow::Z23().
Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), and G4WentzelVIModel::SampleScattering().
00161 { 00162 G4double cosTetMaxNuc2 = cosTetMaxNuc; 00163 if(Z != targetZ || tkin != etag) { 00164 etag = tkin; 00165 targetZ = Z; 00166 if(targetZ > 99) { targetZ = 99; } 00167 SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2); 00168 //G4double tmass2 = targetMass*targetMass; 00169 //G4double etot = tkin + mass; 00170 //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass; 00171 //momCM2 = mom2*tmass2/invmass2; 00172 //gam0pcmp = (etot + targetMass)*targetMass/invmass2; 00173 //pcmp2 = tmass2/invmass2; 00174 00175 kinFactor = coeff*Z*chargeSquare*invbeta2/mom2; 00176 00177 if(1 == Z) { 00178 screenZ = ScreenRSquare[targetZ]/mom2; 00179 } else if(mass > MeV) { 00180 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)* 00181 ScreenRSquare[targetZ]/mom2; 00182 } else { 00183 G4double tau = tkin/mass; 00184 screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z 00185 *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))* 00186 ScreenRSquareElec[targetZ]/mom2; 00187 } 00188 if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) { 00189 cosTetMaxNuc2 = 0.0; 00190 } 00191 formfactA = FormFactor[targetZ]*mom2; 00192 00193 cosTetMaxElec = 1.0; 00194 ComputeMaxElectronScattering(cut); 00195 } 00196 return cosTetMaxNuc2; 00197 }