#include <G4WentzelOKandVIxSection.hh>

Inheritance diagram for G4WentzelOKandVIxSection:

Public Member Functions
	G4WentzelOKandVIxSection (G4bool combined=true)

virtual	~G4WentzelOKandVIxSection ()

void	Initialise (const G4ParticleDefinition *, G4double CosThetaLim)

void	SetupParticle (const G4ParticleDefinition *)

virtual G4double	SetupKinematic (G4double kinEnergy, const G4Material *mat)

G4double	SetupTarget (G4int Z, G4double cut)

G4double	ComputeTransportCrossSectionPerAtom (G4double CosThetaMax)

G4ThreeVector &	SampleSingleScattering (G4double CosThetaMin, G4double CosThetaMax, G4double elecRatio)

G4double	ComputeSecondTransportMoment (G4double CosThetaMax)

G4double	ComputeNuclearCrossSection (G4double CosThetaMin, G4double CosThetaMax)

G4double	ComputeElectronCrossSection (G4double CosThetaMin, G4double CosThetaMax)

void	SetTargetMass (G4double value)

G4double	GetMomentumSquare () const

G4double	GetCosThetaNuc () const

G4double	GetCosThetaElec () const

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

	G4WentzelOKandVIxSection (const G4WentzelOKandVIxSection &)=delete

Protected Member Functions
void	ComputeMaxElectronScattering (G4double cut)

void	InitialiseA ()

G4double	FlatFormfactor (G4double x)

Protected Attributes
const G4ParticleDefinition *	theProton

const G4ParticleDefinition *	theElectron

const G4ParticleDefinition *	thePositron

const G4ParticleDefinition *	particle = nullptr

const G4Material *	currentMaterial = nullptr

G4NistManager *	fNistManager

G4Pow *	fG4pow

G4ScreeningMottCrossSection *	fMottXSection = nullptr

G4ThreeVector	temp

G4double	coeff

G4double	cosTetMaxElec = 1.0

G4double	cosTetMaxNuc = 1.0

G4double	cosThetaMax = 1.0

G4double	chargeSquare = 0.0

G4double	charge3 = 0.0

G4double	spin = 0.0

G4double	mass = 0.0

G4double	tkin = 0.0

G4double	mom2 = 0.0

G4double	momCM2 = 0.0

G4double	invbeta2 = 1.0

G4double	kinFactor = 1.0

G4double	etag = DBL_MAX

G4double	ecut = DBL_MAX

G4double	targetMass

G4double	screenZ = 0.0

G4double	formfactA = 0.0

G4double	factorA2 = 0.0

G4double	factB = 0.0

G4double	factD = 0.0

G4double	fMottFactor = 1.0

G4double	gam0pcmp = 1.0

G4double	pcmp2 = 1.0

G4int	targetZ = 0

G4int	nwarnings = 0

G4NuclearFormfactorType	fNucFormfactor = fExponentialNF

G4bool	isCombined

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

static G4double	ScreenRSquare [100] = {0.0}

static G4double	FormFactor [100] = {0.0}

Detailed Description

Definition at line 72 of file G4WentzelOKandVIxSection.hh.

Constructor & Destructor Documentation

◆ G4WentzelOKandVIxSection() [1/2]

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( G4bool combined = true )

explicit

Definition at line 77 of file G4WentzelOKandVIxSection.cc.

                                                              :
  temp(0.,0.,0.),
  isCombined(comb)
{
  fNistManager = G4NistManager::Instance();
  fG4pow = G4Pow::GetInstance();
 
  theElectron = G4Electron::Electron();
  thePositron = G4Positron::Positron();
  theProton   = G4Proton::Proton();
 
  G4double p0 = CLHEP::electron_mass_c2*CLHEP::classic_electr_radius;
  coeff = CLHEP::twopi*p0*p0;
  targetMass = CLHEP::proton_mass_c2;
}

Referenced by G4WentzelOKandVIxSection(), and operator=().

◆ ~G4WentzelOKandVIxSection()

G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection ( )

virtual

Definition at line 95 of file G4WentzelOKandVIxSection.cc.

{
  delete fMottXSection;
}

◆ G4WentzelOKandVIxSection() [2/2]

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( const G4WentzelOKandVIxSection & )

delete

Member Function Documentation

◆ ComputeElectronCrossSection()

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

inline

Definition at line 225 of file G4WentzelOKandVIxSection.hh.

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

◆ ComputeMaxElectronScattering()

void G4WentzelOKandVIxSection::ComputeMaxElectronScattering ( G4double cut )

protected

Definition at line 392 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);
    cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
  } else {
 
    G4double tmax = (particle == theElectron) ? 0.5*tkin : tkin;
    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; 
      }
    }
  }
}

Referenced by SetupTarget().

◆ ComputeNuclearCrossSection()

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

inline

Definition at line 215 of file G4WentzelOKandVIxSection.hh.

{
  return targetZ*kinFactor*fMottFactor*(cosTMin - cosTMax)/
    ((1.0 - cosTMin + screenZ)*(1.0 - cosTMax + screenZ));
}

◆ ComputeSecondTransportMoment()

G4double G4WentzelOKandVIxSection::ComputeSecondTransportMoment ( G4double CosThetaMax )

Definition at line 422 of file G4WentzelOKandVIxSection.cc.

{
  return 0.0;
}

◆ ComputeTransportCrossSectionPerAtom()

G4double G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom ( G4double CosThetaMax )

Definition at line 242 of file G4WentzelOKandVIxSection.cc.

{
  G4double xSection = 0.0;
  if(cosTMax >= 1.0) { return xSection; }
 
  G4double costm = std::max(cosTMax,cosTetMaxElec); 
  G4double fb = screenZ*factB;
 
  // scattering off electrons
  if(costm < 1.0) {
    G4double x = (1.0 - costm)/screenZ;
    if(x < numlimit) { 
      G4double x2 = 0.5*x*x;
      xSection = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x)); 
    } else { 
      G4double x1= x/(1 + x);
      G4double xlog = G4Log(1.0 + x);  
      xSection = xlog - x1 - fb*(x + x1 - 2*xlog);
    }
 
    if(xSection < 0.0) {
      ++nwarnings;
      if(nwarnings < nwarnlimit) {
        G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
               << " scattering on e- <0"
               << G4endl;
        G4cout << "cross= " << xSection
               << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
               << " Z= " << targetZ << "  " 
               << particle->GetParticleName() << G4endl;
        G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
               << " x= " << x << G4endl;
      }
      xSection = 0.0;
    }
  }
  /*  
      G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom: \n"
      << " Z= " << targetZ 
      << " e(MeV)= " << tkin/MeV << " XSel= " << xSection  
      << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
      << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
      << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
  */
  // scattering off nucleus
  if(cosTMax < 1.0) {
    G4double x = (1.0 - cosTMax)/screenZ;
    G4double y;
    if(x < numlimit) { 
      G4double x2 = 0.5*x*x;
      y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x)); 
    } else { 
      G4double x1= x/(1 + x);
      G4double xlog = G4Log(1.0 + x);  
      y = xlog - x1 - 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= " << x <<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; 
}

◆ FlatFormfactor()

G4double G4WentzelOKandVIxSection::FlatFormfactor ( G4double x )

inlineprotected

Definition at line 234 of file G4WentzelOKandVIxSection.hh.

{
  return 3.0*(std::sin(x) - x*std::cos(x))/(x*x*x);
}

Referenced by SampleSingleScattering().

◆ GetCosThetaElec()

G4double G4WentzelOKandVIxSection::GetCosThetaElec ( ) const

inline

Definition at line 207 of file G4WentzelOKandVIxSection.hh.

{
  return cosTetMaxElec;
}

◆ GetCosThetaNuc()

G4double G4WentzelOKandVIxSection::GetCosThetaNuc ( ) const

inline

Definition at line 200 of file G4WentzelOKandVIxSection.hh.

{
  return cosTetMaxNuc;
}

◆ GetMomentumSquare()

G4double G4WentzelOKandVIxSection::GetMomentumSquare ( ) const

inline

Definition at line 193 of file G4WentzelOKandVIxSection.hh.

{
  return mom2;
}

◆ Initialise()

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

Definition at line 102 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; } 
  G4EmParameters* param = G4EmParameters::Instance();
  G4double a = param->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
  factorA2 = 0.5*a*a;
  currentMaterial = nullptr;
 
  fNucFormfactor = param->NuclearFormfactorType();
  if(0.0 == ScreenRSquare[0]) { InitialiseA(); }
 
  // Mott corrections always added
  if((p == theElectron || p == thePositron) && !fMottXSection) {
    fMottXSection = new G4ScreeningMottCrossSection();
    fMottXSection->Initialise(p, 1.0);
  }
  /*
  G4cout << "G4WentzelOKandVIxSection::Initialise for " 
     << p->GetParticleName() << " cosThetaMax= " << cosThetaMax 
     << "  " << ScreenRSquare[0] << " coeff= " << coeff << G4endl;
  */
}

◆ InitialiseA()

void G4WentzelOKandVIxSection::InitialiseA ( )

protected

Definition at line 134 of file G4WentzelOKandVIxSection.cc.

{
  // Thomas-Fermi screening radii
  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
  if(0.0 != ScreenRSquare[0]) { return; }
  G4AutoLock l(&theWOKVIMutex);
  if(0.0 == ScreenRSquare[0]) {
    const G4double invmev2 = 1./(CLHEP::MeV*CLHEP::MeV);
    G4double a0 = CLHEP::electron_mass_c2/0.88534; 
    G4double constn = 6.937e-6*invmev2;
    G4double fct = G4EmParameters::Instance()->ScreeningFactor();
 
    G4double afact = 0.5*fct*alpha2*a0*a0;
    ScreenRSquare[0] = afact;
    ScreenRSquare[1] = afact;
    ScreenRSquareElec[1] = afact; 
    FormFactor[1] = 3.097e-6*invmev2;
 
    for(G4int j=2; j<100; ++j) {
      G4double x = fG4pow->Z13(j);
      ScreenRSquare[j] = afact*(1 + G4Exp(-j*j*0.001))*x*x;
      ScreenRSquareElec[j] = afact*x*x;
      x = fNistManager->GetA27(j);
      FormFactor[j] = constn*x*x;
    } 
  }  
  l.unlock();
}

Referenced by Initialise().

◆ operator=()

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

delete

◆ SampleSingleScattering()

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

Definition at line 331 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;
    G4double w2 = 1. - cost2;
    G4double w3 = rndmEngineMod->flat()*(w2 - w1);
    G4double z1 = ((w2 - w3)*screenZ + w1*w2)/(screenZ + w1 + w3);
    G4double fm = 1.0;
 
    if(fNucFormfactor == fExponentialNF) {
      fm += formf*z1;
      fm = 1.0/(fm*fm);
    } else if(fNucFormfactor == fGaussianNF) {
      fm = G4Exp(-2*formf*z1);
    } else if(fNucFormfactor == fFlatNF) {
      static const G4double ccoef = 0.00508/CLHEP::MeV;
      G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
      fm = FlatFormfactor(x);
      fm *= FlatFormfactor(x*0.6*fG4pow->A13(fNistManager->GetAtomicMassAmu(targetZ)));
    }
    // G4cout << " fm=" << fm << "  " << fMottXSection << G4endl;
    G4double grej;
    if(nullptr != fMottXSection) {
      fMottXSection->SetupKinematic(tkin, targetZ);
      grej = fMottXSection->RatioMottRutherfordCosT(std::sqrt(z1))*fm;
    } else {
      grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
      *fm/(1.0 + z1*factD);
    }
    if(fMottFactor*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;
}

◆ SetTargetMass()

void G4WentzelOKandVIxSection::SetTargetMass ( G4double value )

inline

Definition at line 185 of file G4WentzelOKandVIxSection.hh.

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

Referenced by G4WentzelVIRelModel::DefineMaterial(), and SetupTarget().

◆ SetupKinematic()

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

virtual

Reimplemented in G4WentzelVIRelXSection.

Definition at line 182 of file G4WentzelOKandVIxSection.cc.

{
  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 = isCombined ? 
      std::max(cosThetaMax, 1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2)
      : cosThetaMax;
  } 
  return cosTetMaxNuc;
}

◆ SetupParticle()

void G4WentzelOKandVIxSection::SetupParticle ( const G4ParticleDefinition * p )

Definition at line 165 of file G4WentzelOKandVIxSection.cc.

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

Referenced by Initialise().

◆ SetupTarget()

G4double G4WentzelOKandVIxSection::SetupTarget	(	G4int	Z,
		G4double	cut )

Definition at line 200 of file G4WentzelOKandVIxSection.cc.

{
  G4double cosTetMaxNuc2 = cosTetMaxNuc;
  if(Z != targetZ || tkin != etag) {
    etag    = tkin; 
    targetZ = std::min(Z, 99);
    G4double massT = (1 == Z) ? CLHEP::proton_mass_c2 :
      fNistManager->GetAtomicMassAmu(Z)*CLHEP::amu_c2;
    SetTargetMass(massT);
 
    kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
    if(particle == theElectron && fMottXSection) {
      fMottFactor = (1.0 + 2.0e-4*Z*Z);
    }
 
    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 && particle == theProton && cosTetMaxNuc2 < 0.0) {
      cosTetMaxNuc2 = 0.0;
    }
    formfactA = FormFactor[targetZ]*mom2;
 
    cosTetMaxElec = 1.0;
    ComputeMaxElectronScattering(cut); 
  }
  //G4cout << "SetupTarget:  Z= " << targetZ << " kinFactor= " << kinFactor
  //     << " fMottFactor= " << fMottFactor << " screenZ= " << screenZ <<G4endl;
  return cosTetMaxNuc2;
} 