G4WentzelVIRelXSection Class Reference

#include <G4WentzelVIRelXSection.hh>


Public Member Functions

 G4WentzelVIRelXSection ()
virtual ~G4WentzelVIRelXSection ()
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)
void SetRelativisticMass (G4double value)
G4double GetMomentumSquare () const
G4double GetCosThetaNuc () const
G4double GetCosThetaElec () const

Detailed Description

Definition at line 71 of file G4WentzelVIRelXSection.hh.

Constructor & Destructor Documentation

G4WentzelVIRelXSection::G4WentzelVIRelXSection (  ) 

Definition at line 65 of file G4WentzelVIRelXSection.cc.

References G4Electron::Electron(), G4NistManager::GetA27(), G4Pow::GetInstance(), Initialise(), G4NistManager::Instance(), G4INCL::Math::pi, G4Positron::Positron(), G4Proton::Proton(), and G4Pow::Z13().

00065                                                :
00066   numlimit(0.1),
00067   nwarnings(0),
00068   nwarnlimit(50),
00069   alpha2(fine_structure_const*fine_structure_const)
00070 {
00071   fNistManager = G4NistManager::Instance();
00072   fG4pow = G4Pow::GetInstance();
00073   theElectron = G4Electron::Electron();
00074   thePositron = G4Positron::Positron();
00075   theProton   = G4Proton::Proton();
00076   lowEnergyLimit = 1.0*eV;
00077   G4double p0 = electron_mass_c2*classic_electr_radius;
00078   coeff = twopi*p0*p0;
00079   particle = 0;
00080 
00081   // Thomas-Fermi screening radii
00082   // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
00083 
00084   if(0.0 == ScreenRSquare[0]) {
00085     G4double a0 = electron_mass_c2/0.88534; 
00086     G4double constn = 6.937e-6/(MeV*MeV);
00087 
00088     ScreenRSquare[0] = alpha2*a0*a0;
00089     for(G4int j=1; j<100; ++j) {
00090       G4double x = a0*fG4pow->Z13(j);
00091       //ScreenRSquare[j] = 0.5*(1 + exp(-j*j*0.001))*alpha2*x*x;
00092       ScreenRSquare[j] = 0.5*alpha2*x*x;
00093       x = fNistManager->GetA27(j);
00094       FormFactor[j] = constn*x*x;
00095     } 
00096   }
00097   currentMaterial = 0;
00098   elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
00099   cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
00100 
00101   factB1= 0.5*CLHEP::pi*fine_structure_const;
00102 
00103   Initialise(theElectron, 1.0);
00104   targetMass = proton_mass_c2;
00105 }

G4WentzelVIRelXSection::~G4WentzelVIRelXSection (  )  [virtual]

Definition at line 109 of file G4WentzelVIRelXSection.cc.

00110 {}

Member Function Documentation

G4double G4WentzelVIRelXSection::ComputeElectronCrossSection ( G4double  CosThetaMin,
G4double  CosThetaMax 
) [inline]

Definition at line 246 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom().

00248 {
00249   G4double xsec = 0.0;
00250   G4double cost1 = std::max(cosTMin,cosTetMaxElec);
00251   G4double cost2 = std::max(cosTMax,cosTetMaxElec);
00252   if(cost1 > cost2) {
00253     xsec = kinFactor*(cost1 - cost2)/((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
00254   }
00255   return xsec;
00256 }

G4double G4WentzelVIRelXSection::ComputeNuclearCrossSection ( G4double  CosThetaMin,
G4double  CosThetaMax 
) [inline]

Definition at line 232 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom().

00234 {
00235   G4double xsec = 0.0;
00236   if(cosTMax < cosTMin) {
00237     xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
00238       ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
00239   }
00240   return xsec;
00241 }

G4double G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom ( G4double  CosThetaMax  ) 

Definition at line 192 of file G4WentzelVIRelXSection.cc.

References G4cout, G4endl, and G4ParticleDefinition::GetParticleName().

Referenced by G4WentzelVIRelModel::ComputeCrossSectionPerAtom().

00193 {
00194   G4double xsec = 0.0;
00195   if(cosTMax >= 1.0) { return xsec; }
00196  
00197   G4double xSection = 0.0;
00198   G4double x = 0; 
00199   G4double y = 0;
00200   G4double x1= 0;
00201   G4double x2= 0;
00202   G4double xlog = 0.0;
00203 
00204   G4double costm = std::max(cosTMax,cosTetMaxElec); 
00205   G4double fb = screenZ*factB;
00206 
00207   // scattering off electrons
00208   if(costm < 1.0) {
00209     x = (1.0 - costm)/screenZ;
00210     if(x < numlimit) { 
00211       x2 = 0.5*x*x;
00212       y  = x2*(1.0 - 1.3333333*x + 3*x2);
00213       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00214     } else { 
00215       x1= x/(1 + x);
00216       xlog = log(1.0 + x);  
00217       y = xlog - x1; 
00218       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00219     }
00220 
00221     if(y < 0.0) {
00222       ++nwarnings;
00223       if(nwarnings < nwarnlimit) {
00224         G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00225                << G4endl;
00226         G4cout << "y= " << y 
00227                << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
00228                << " Z= " << targetZ << "  " 
00229                << particle->GetParticleName() << G4endl;
00230         G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
00231                << " x= " << x << G4endl;
00232       }
00233       y = 0.0;
00234     }
00235     xSection = y;
00236   }
00237   /* 
00238   G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
00239          << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
00240          << " cut(MeV)= " << ecut/MeV  
00241          << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
00242          << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
00243          << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
00244   */
00245   // scattering off nucleus
00246   if(cosTMax < 1.0) {
00247     x = (1.0 - cosTMax)/screenZ;
00248     if(x < numlimit) { 
00249       x2 = 0.5*x*x;
00250       y  = x2*(1.0 - 1.3333333*x + 3*x2); 
00251       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
00252     } else { 
00253       x1= x/(1 + x);
00254       xlog = log(1.0 + x);  
00255       y = xlog - x1; 
00256       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
00257     }
00258 
00259     if(y < 0.0) {
00260       ++nwarnings;
00261       if(nwarnings < nwarnlimit) {
00262         G4cout << "G4WentzelVIRelXSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
00263                << G4endl;
00264         G4cout << "y= " << y 
00265                << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
00266                << particle->GetParticleName() << G4endl;
00267         G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 
00268                << " x= " << " x1= " << x1 <<G4endl;
00269       }
00270       y = 0.0;
00271     }
00272     xSection += y*targetZ; 
00273   }
00274   xSection *= kinFactor;
00275   /*
00276   G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
00277          << " screenZ= " << screenZ << " formF= " << formfactA 
00278          << " for " << particle->GetParticleName() 
00279          << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
00280          << " x= " << x 
00281          << G4endl;
00282   */
00283   return xSection; 
00284 }

G4double G4WentzelVIRelXSection::GetCosThetaElec (  )  const [inline]

Definition at line 224 of file G4WentzelVIRelXSection.hh.

00225 {
00226   return cosTetMaxElec;
00227 }

G4double G4WentzelVIRelXSection::GetCosThetaNuc (  )  const [inline]

Definition at line 217 of file G4WentzelVIRelXSection.hh.

00218 {
00219   return cosTetMaxNuc;
00220 }

G4double G4WentzelVIRelXSection::GetMomentumSquare (  )  const [inline]

Definition at line 210 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::SampleSecondaries().

00211 {
00212   return mom2;
00213 }

void G4WentzelVIRelXSection::Initialise ( const G4ParticleDefinition *  ,
G4double  CosThetaLim 
)

Definition at line 114 of file G4WentzelVIRelXSection.cc.

References DBL_MAX, G4LossTableManager::FactorForAngleLimit(), G4LossTableManager::Instance(), and SetupParticle().

Referenced by G4WentzelVIRelXSection(), G4WentzelVIRelModel::Initialise(), and G4hCoulombScatteringModel::Initialise().

00116 {
00117   SetupParticle(p);
00118   tkin = mom2 = momCM2 = 0.0;
00119   ecut = etag = DBL_MAX;
00120   targetZ = 0;
00121   cosThetaMax = CosThetaLim;
00122   G4double a = 
00123     G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
00124   factorA2 = 0.5*a*a;
00125   currentMaterial = 0;
00126 }

G4ThreeVector G4WentzelVIRelXSection::SampleSingleScattering ( G4double  CosThetaMin,
G4double  CosThetaMax,
G4double  elecRatio = 0.0 
)

Definition at line 289 of file G4WentzelVIRelXSection.cc.

References G4UniformRand.

Referenced by G4WentzelVIRelModel::SampleScattering(), and G4hCoulombScatteringModel::SampleSecondaries().

00292 {
00293   G4ThreeVector v(0.0,0.0,1.0);
00294  
00295   G4double formf = formfactA;
00296   G4double cost1 = cosTMin;
00297   G4double cost2 = cosTMax;
00298   if(elecRatio > 0.0) {
00299     if(G4UniformRand() <= elecRatio) {
00300       formf = 0.0;
00301       cost1 = std::max(cost1,cosTetMaxElec);
00302       cost2 = std::max(cost2,cosTetMaxElec);
00303     }
00304   }
00305   if(cost1 < cost2) { return v; }
00306 
00307   G4double w1 = 1. - cost1 + screenZ;
00308   G4double w2 = 1. - cost2 + screenZ;
00309   G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
00310 
00311   //G4double fm =  1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
00312   G4double fm =  1.0 + formf*z1;
00313   //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
00314   G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
00315   // "false" scattering
00316   if( G4UniformRand() > grej ) { return v; }
00317     // } 
00318   G4double cost = 1.0 - z1;
00319 
00320   if(cost > 1.0)       { cost = 1.0; }
00321   else if(cost < -1.0) { cost =-1.0; }
00322   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00323   //G4cout << "sint= " << sint << G4endl;
00324   G4double phi  = twopi*G4UniformRand();
00325   G4double vx1 = sint*cos(phi);
00326   G4double vy1 = sint*sin(phi);
00327 
00328   // only direction is changed
00329   v.set(vx1,vy1,cost);
00330   return v;
00331 }

void G4WentzelVIRelXSection::SetRelativisticMass ( G4double  value  )  [inline]

Definition at line 203 of file G4WentzelVIRelXSection.hh.

00204 {
00205   mass = value;
00206 }

void G4WentzelVIRelXSection::SetTargetMass ( G4double  value  )  [inline]

Definition at line 195 of file G4WentzelVIRelXSection.hh.

Referenced by G4hCoulombScatteringModel::SampleSecondaries(), and SetupTarget().

00196 {
00197   targetMass = value;
00198   factD = std::sqrt(mom2)/value;
00199 }

G4double G4WentzelVIRelXSection::SetupKinematic ( G4double  kinEnergy,
const G4Material *  mat 
) [inline]

Definition at line 179 of file G4WentzelVIRelXSection.hh.

References G4IonisParamMat::GetInvA23(), and G4Material::GetIonisation().

Referenced by G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeGeomPathLength(), G4WentzelVIRelModel::ComputeTruePathLengthLimit(), and G4WentzelVIRelModel::ComputeTrueStepLength().

00180 {
00181   if(ekin != tkin || mat != currentMaterial) { 
00182     currentMaterial = mat;
00183     tkin  = ekin;
00184     mom2  = tkin*(tkin + 2.0*mass);
00185     invbeta2 = 1.0 +  mass*mass/mom2;
00186     factB = spin/invbeta2;
00187     cosTetMaxNuc = 
00188         std::max(cosThetaMax,1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2);
00189   } 
00190   return cosTetMaxNuc;
00191 }

void G4WentzelVIRelXSection::SetupParticle ( const G4ParticleDefinition *   ) 

Definition at line 130 of file G4WentzelVIRelXSection.cc.

References G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and G4ParticleDefinition::GetPDGSpin().

Referenced by Initialise(), and G4hCoulombScatteringModel::SetupParticle().

00131 {
00132   particle = p;
00133   mass = particle->GetPDGMass();
00134   spin = particle->GetPDGSpin();
00135   if(0.0 != spin) { spin = 0.5; }
00136   G4double q = std::fabs(particle->GetPDGCharge()/eplus);
00137   chargeSquare = q*q;
00138   charge3 = chargeSquare*q;
00139   tkin = 0.0;
00140   currentMaterial = 0;
00141   targetZ = 0;
00142 }

G4double G4WentzelVIRelXSection::SetupTarget ( G4int  Z,
G4double  cut = DBL_MAX 
)

Definition at line 147 of file G4WentzelVIRelXSection.cc.

References G4NistManager::GetAtomicMassAmu(), and SetTargetMass().

Referenced by G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), and G4WentzelVIRelModel::SampleScattering().

00148 {
00149   G4double cosTetMaxNuc2 = cosTetMaxNuc;
00150   if(Z != targetZ || tkin != etag) {
00151     etag    = tkin; 
00152     targetZ = Z;
00153     if(targetZ > 99) { targetZ = 99; }
00154     SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
00155     //G4double tmass2 = targetMass*targetMass;
00156     //G4double etot = tkin + mass;
00157     //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
00158     //momCM2 = mom2*tmass2/invmass2;
00159     //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
00160     //pcmp2    = tmass2/invmass2;
00161 
00162     kinFactor = coeff*targetZ*chargeSquare*invbeta2/mom2;
00163 
00164     screenZ = ScreenRSquare[targetZ]/mom2;
00165     if(Z > 1) {
00166       screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
00167       /*
00168       if(mass > MeV) {
00169         screenZ *= std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare);
00170       } else {
00171         G4double tau = tkin/mass;
00172         //      screenZ *= std::min(Z*invbeta2,
00173         screenZ *= std::min(Z*1.13,
00174         (1.13 +3.76*Z*Z*invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(Z)))));
00175       }
00176       */
00177     }
00178     if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
00179       cosTetMaxNuc2 = 0.0;
00180     }
00181     formfactA = FormFactor[targetZ]*mom2;
00182 
00183     cosTetMaxElec = 1.0;
00184     ComputeMaxElectronScattering(cut); 
00185   }
00186   return cosTetMaxNuc2;
00187 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:59 2013 for Geant4 by  doxygen 1.4.7