#include <G4WentzelOKandVIxSection.hh>

Collaboration diagram for G4WentzelOKandVIxSection:

Public Member Functions
	G4WentzelOKandVIxSection (G4bool combined=true)

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	ComputeSecondTransportMoment (G4double CosThetaMax)

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

Private Member Functions
void	ComputeMaxElectronScattering (G4double cut)

G4WentzelOKandVIxSection &	operator= (const G4WentzelOKandVIxSection &right)

	G4WentzelOKandVIxSection (const G4WentzelOKandVIxSection &)

Private Attributes
const G4ParticleDefinition *	theProton

const G4ParticleDefinition *	theElectron

const G4ParticleDefinition *	thePositron

const G4Material *	currentMaterial

G4NistManager *	fNistManager

G4Pow *	fG4pow

G4ThreeVector	temp

G4double	numlimit

G4int	nwarnings

G4int	nwarnlimit

G4bool	isCombined

G4double	coeff

G4double	cosTetMaxElec

G4double	cosTetMaxNuc

G4double	cosThetaMax

G4double	alpha2

const G4ParticleDefinition *	particle

G4double	chargeSquare

G4double	charge3

G4double	spin

G4double	mass

G4double	tkin

G4double	mom2

G4double	momCM2

G4double	invbeta2

G4double	kinFactor

G4double	etag

G4double	ecut

G4double	lowEnergyLimit

G4int	targetZ

G4double	targetMass

G4double	screenZ

G4double	formfactA

G4double	factorA2

G4double	factB

G4double	factB1

G4double	factD

G4double	gam0pcmp

G4double	pcmp2

Static Private Attributes
static G4double	ScreenRSquareElec [100] = {0.0}

static G4double	ScreenRSquare [100] = {0.0}

static G4double	FormFactor [100] = {0.0}

Detailed Description

Definition at line 71 of file G4WentzelOKandVIxSection.hh.

Constructor & Destructor Documentation

◆ G4WentzelOKandVIxSection() [1/2]

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( G4bool combined = true )

Definition at line 66 of file G4WentzelOKandVIxSection.cc.

                                                                   :
   temp(0.,0.,0.),
   numlimit(0.1),
   nwarnings(0),
   nwarnlimit(50),
   isCombined(combined),
   cosThetaMax(-1.0),
   alpha2(fine_structure_const*fine_structure_const)
 {
   fNistManager = G4NistManager::Instance();
   fG4pow = G4Pow::GetInstance();
   theElectron = G4Electron::Electron();
   thePositron = G4Positron::Positron();
   theProton   = G4Proton::Proton();
   lowEnergyLimit = 1.0*eV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff = twopi*p0*p0;
   particle = 0;
 
   // Thomas-Fermi screening radii
   // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
 
   if(0.0 == ScreenRSquare[0]) {
     G4double a0 = electron_mass_c2/0.88534; 
     G4double constn = 6.937e-6/(MeV*MeV);
 
     ScreenRSquare[0] = alpha2*a0*a0;
     ScreenRSquareElec[0] = ScreenRSquare[0];
     for(G4int j=1; j<100; ++j) {
       G4double x = a0*fG4pow->Z13(j);
       if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
       else {
         ScreenRSquare[j] = 0.5*(1 + G4Exp(-j*j*0.001))*alpha2*x*x;
         ScreenRSquareElec[j] = 0.5*alpha2*x*x;
       }
       x = fNistManager->GetA27(j);
       FormFactor[j] = constn*x*x;
     } 
   }
   currentMaterial = 0;
   factB = factD = formfactA = screenZ = 0.0;
   cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
 
   factB1= 0.5*CLHEP::pi*fine_structure_const;
 
   tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
   ecut = etag = DBL_MAX;
   targetZ = 0;
   targetMass = proton_mass_c2;
   particle = 0;
 }

Here is the call graph for this function:

◆ ~G4WentzelOKandVIxSection()

G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection ( )

virtual

Definition at line 120 of file G4WentzelOKandVIxSection.cc.

121 {}

◆ G4WentzelOKandVIxSection() [2/2]

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( const G4WentzelOKandVIxSection & )

private

Member Function Documentation

◆ ComputeElectronCrossSection()

G4double G4WentzelOKandVIxSection::ComputeElectronCrossSection	(	G4double	CosThetaMin,
		G4double	CosThetaMax
	)

inline

Definition at line 244 of file G4WentzelOKandVIxSection.hh.

 {
   G4double xsec = 0.0;
   G4double cost1 = std::max(cosTMin,cosTetMaxElec);
   G4double cost2 = std::max(cosTMax,cosTetMaxElec);
   if(cost1 > cost2) {
     xsec = kinFactor*(cost1 - cost2)/
       ((1.0 - cost1 + screenZ)*(1.0 - cost2 + screenZ));
   }
   return xsec;
 }

Here is the caller graph for this function:

◆ ComputeMaxElectronScattering()

void G4WentzelOKandVIxSection::ComputeMaxElectronScattering ( G4double cut )

private

Definition at line 348 of file G4WentzelOKandVIxSection.cc.

 {
   if(mass > MeV) {
     G4double ratio = electron_mass_c2/mass;
     G4double tau = tkin/mass;
     G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
       (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
     //tmax = std::min(tmax, targetZ*targetZ*10*eV); 
     cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
   } else {
 
     G4double tmax = tkin;
     if(particle == theElectron) { tmax *= 0.5; }
     //tmax = std::min(tmax, targetZ*targetZ*10*eV); 
     G4double t = std::min(cutEnergy, tmax);
     G4double mom21 = t*(t + 2.0*electron_mass_c2);
     G4double t1 = tkin - t;
     //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= " 
     //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
     if(t1 > 0.0) {
       G4double mom22 = t1*(t1 + 2.0*mass);
       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
       if(ctm <  1.0) { cosTetMaxElec = ctm; }
       if(particle == theElectron && cosTetMaxElec < 0.0) { 
         cosTetMaxElec = 0.0; 
       }
     }
   }
 }

Here is the caller graph for this function:

◆ ComputeNuclearCrossSection()

G4double G4WentzelOKandVIxSection::ComputeNuclearCrossSection	(	G4double	CosThetaMin,
		G4double	CosThetaMax
	)

inline

Definition at line 230 of file G4WentzelOKandVIxSection.hh.

 {
   G4double xsec = 0.0;
   if(cosTMax < cosTMin) {
     xsec = targetZ*kinFactor*(cosTMin - cosTMax)/
       ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
   }
   return xsec;
 }

Here is the caller graph for this function:

◆ ComputeSecondTransportMoment()

G4double G4WentzelOKandVIxSection::ComputeSecondTransportMoment ( G4double CosThetaMax )

Definition at line 381 of file G4WentzelOKandVIxSection.cc.

 {
   return 0.0;
 }

Here is the caller graph for this function:

◆ ComputeTransportCrossSectionPerAtom()

G4double G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom ( G4double CosThetaMax )

Definition at line 206 of file G4WentzelOKandVIxSection.cc.

 {
   G4double xSection = 0.0;
   if(cosTMax >= 1.0) { return xSection; }
  
   G4double x = 0; 
   G4double y = 0;
   G4double x1= 0;
   G4double x2= 0;
   G4double xlog = 0.0;
 
   G4double costm = std::max(cosTMax,cosTetMaxElec); 
   G4double fb = screenZ*factB;
 
   // scattering off electrons
   if(costm < 1.0) {
     x = (1.0 - costm)/screenZ;
     if(x < numlimit) { 
       x2 = 0.5*x*x;
       y  = x2*(1.0 - 1.3333333*x + 3*x2);
       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
     } else { 
       x1= x/(1 + x);
       xlog = G4Log(1.0 + x);  
       y = xlog - x1; 
       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
     }
 
     if(y < 0.0) {
       ++nwarnings;
       if(nwarnings < nwarnlimit) {
         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
                << " scattering on e- <0"
                << G4endl;
         G4cout << "y= " << y 
                << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
                << " Z= " << targetZ << "  " 
                << particle->GetParticleName() << G4endl;
         G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
                << " x= " << x << G4endl;
       }
       y = 0.0;
     }
     xSection = y;
   }
   /* 
        G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
        << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
          << " cut(MeV)= " << ecut/MeV  
            << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
          << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
          << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
   */
   // scattering off nucleus
   if(cosTMax < 1.0) {
     x = (1.0 - cosTMax)/screenZ;
     if(x < numlimit) { 
       x2 = 0.5*x*x;
       y  = x2*(1.0 - 1.3333333*x + 3*x2); 
       if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
     } else { 
       x1= x/(1 + x);
       xlog = G4Log(1.0 + x);  
       y = xlog - x1; 
       if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
     }
 
     if(y < 0.0) {
       ++nwarnings;
       if(nwarnings < nwarnlimit) {
         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
                << " scattering on nucleus <0"
                << G4endl;
         G4cout << "y= " << y 
                << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
                << particle->GetParticleName() << G4endl;
         G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 
                << " x= " << " x1= " << x1 <<G4endl;
       }
       y = 0.0;
     }
     xSection += y*targetZ; 
   }
   xSection *= kinFactor;
  
   /*
   G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
          << " screenZ= " << screenZ << " formF= " << formfactA 
          << " for " << particle->GetParticleName() 
            << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
          << " x= " << x 
          << G4endl;
   */
   return xSection; 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ GetCosThetaElec()

G4double G4WentzelOKandVIxSection::GetCosThetaElec ( ) const

inline

Definition at line 222 of file G4WentzelOKandVIxSection.hh.

 {
   return cosTetMaxElec;
 }

◆ GetCosThetaNuc()

G4double G4WentzelOKandVIxSection::GetCosThetaNuc ( ) const

inline

Definition at line 215 of file G4WentzelOKandVIxSection.hh.

 {
   return cosTetMaxNuc;
 }

◆ GetMomentumSquare()

G4double G4WentzelOKandVIxSection::GetMomentumSquare ( ) const

inline

Definition at line 208 of file G4WentzelOKandVIxSection.hh.

 {
   return mom2;
 }

Here is the caller graph for this function:

◆ Initialise()

void G4WentzelOKandVIxSection::Initialise	(	const G4ParticleDefinition *	p,
		G4double	CosThetaLim
	)

Definition at line 125 of file G4WentzelOKandVIxSection.cc.

 {
   SetupParticle(p);
   tkin = mom2 = momCM2 = 0.0;
   ecut = etag = DBL_MAX;
   targetZ = 0;
 
   // cosThetaMax is below 1.0 only when MSC is combined with SS
   if(isCombined) { cosThetaMax = cosThetaLim; }
   
   G4double a = G4EmParameters::Instance()->FactorForAngleLimit()
     *CLHEP::hbarc/CLHEP::fermi;
   factorA2 = 0.5*a*a;
   currentMaterial = 0;
   
   //G4cout << "G4WentzelOKandVIxSection::Initialise  mass= " << mass
   //         << "  " << p->GetParticleName() 
   //     << "  cosThetaMax= " << cosThetaMax << G4endl; 
   
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ operator=()

G4WentzelOKandVIxSection& G4WentzelOKandVIxSection::operator= ( const G4WentzelOKandVIxSection & right )

private

◆ SampleSingleScattering()

G4ThreeVector & G4WentzelOKandVIxSection::SampleSingleScattering	(	G4double	CosThetaMin,
		G4double	CosThetaMax,
		G4double	elecRatio = `0.0`
	)

Definition at line 305 of file G4WentzelOKandVIxSection.cc.

 {
   temp.set(0.0,0.0,1.0);
   CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
  
   G4double formf = formfactA;
   G4double cost1 = cosTMin;
   G4double cost2 = cosTMax;
   if(elecRatio > 0.0) {
     if(rndmEngineMod->flat() <= elecRatio) {
       formf = 0.0;
       cost1 = std::max(cost1,cosTetMaxElec);
       cost2 = std::max(cost2,cosTetMaxElec);
     }
   }
   if(cost1 > cost2) {
     G4double w1 = 1. - cost1 + screenZ;
     G4double w2 = 1. - cost2 + screenZ;
     G4double z1 = w1*w2/(w1 + rndmEngineMod->flat()*(w2 - w1)) - screenZ;
 
     G4double fm = 1.0 + formf*z1;
     G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))
       /((1.0 + z1*factD)*fm*fm);
 
     if(rndmEngineMod->flat() <= grej) {
       // exclude "false" scattering due to formfactor and spin effect
       G4double cost = 1.0 - z1;
       if(cost > 1.0)       { cost = 1.0; }
       else if(cost < -1.0) { cost =-1.0; }
       G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
       //G4cout << "sint= " << sint << G4endl;
       G4double phi  = twopi*rndmEngineMod->flat();
       temp.set(sint*cos(phi),sint*sin(phi),cost);
     }
   }
   return temp;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ SetTargetMass()

void G4WentzelOKandVIxSection::SetTargetMass ( G4double value )

inline

Definition at line 200 of file G4WentzelOKandVIxSection.hh.

 {
   targetMass = value;
   factD = std::sqrt(mom2)/value;
 }

Here is the caller graph for this function:

◆ SetupKinematic()

G4double G4WentzelOKandVIxSection::SetupKinematic	(	G4double	kinEnergy,
		const G4Material *	mat
	)

inline

Definition at line 181 of file G4WentzelOKandVIxSection.hh.

 {
   if(ekin != tkin || mat != currentMaterial) { 
     currentMaterial = mat;
     tkin  = ekin;
     mom2  = tkin*(tkin + 2.0*mass);
     invbeta2 = 1.0 +  mass*mass/mom2;
     factB = spin/invbeta2; 
     cosTetMaxNuc = cosThetaMax;
     if(isCombined) {
       cosTetMaxNuc = std::max(cosTetMaxNuc, 
     1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2);
     }
   } 
   return cosTetMaxNuc;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ SetupParticle()

void G4WentzelOKandVIxSection::SetupParticle ( const G4ParticleDefinition * p )

Definition at line 149 of file G4WentzelOKandVIxSection.cc.

 {
   particle = p;
   mass = particle->GetPDGMass();
   spin = particle->GetPDGSpin();
   if(0.0 != spin) { spin = 0.5; }
   G4double q = std::fabs(particle->GetPDGCharge()/eplus);
   chargeSquare = q*q;
   charge3 = chargeSquare*q;
   tkin = 0.0;
   currentMaterial = 0;
   targetZ = 0;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ SetupTarget()

G4double G4WentzelOKandVIxSection::SetupTarget	(	G4int	Z,
		G4double	cut = `DBL_MAX`
	)

Definition at line 166 of file G4WentzelOKandVIxSection.cc.

 {
   G4double cosTetMaxNuc2 = cosTetMaxNuc;
   if(Z != targetZ || tkin != etag) {
     etag    = tkin; 
     targetZ = Z;
     if(targetZ > 99) { targetZ = 99; }
     G4double massT = proton_mass_c2;
     if(targetZ > 1) { 
       massT = fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2;
     }
     SetTargetMass(massT);
 
     kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
 
     if(1 == Z) {
       screenZ = ScreenRSquare[targetZ]/mom2;
     } else if(mass > MeV) {
       screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
         ScreenRSquare[targetZ]/mom2;
     } else {
       G4double tau = tkin/mass;
       screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
           *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
         ScreenRSquareElec[targetZ]/mom2;
     }
     if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
       cosTetMaxNuc2 = 0.0;
     }
     formfactA = FormFactor[targetZ]*mom2;
 
     cosTetMaxElec = 1.0;
     ComputeMaxElectronScattering(cut); 
   }
   return cosTetMaxNuc2;
 } 