#include <G4WentzelOKandVIxSection.hh>

Inheritance diagram for G4WentzelOKandVIxSection:

Public Member Functions
	G4WentzelOKandVIxSection (G4bool comb=true)

virtual	~G4WentzelOKandVIxSection ()

void	Initialise (const G4ParticleDefinition *, G4double CosThetaLim)

void	SetupParticle (const G4ParticleDefinition *)

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 G4Material *	currentMaterial

G4NistManager *	fNistManager

G4Pow *	fG4pow

G4ScreeningMottCrossSection *	fMottXSection

G4ThreeVector	temp

G4double	numlimit

G4int	nwarnings

G4int	nwarnlimit

G4NuclearFormfactorType	fNucFormfactor

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	fMottFactor

G4double	gam0pcmp

G4double	pcmp2

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 73 of file G4WentzelOKandVIxSection.hh.

Constructor & Destructor Documentation

◆ G4WentzelOKandVIxSection() [1/2]

G4WentzelOKandVIxSection::G4WentzelOKandVIxSection ( G4bool comb = true )

explicit

Definition at line 70 of file G4WentzelOKandVIxSection.cc.

                                                              :
  temp(0.,0.,0.),
  numlimit(0.1),
  nwarnings(0),
  nwarnlimit(50),
  isCombined(comb),
  cosThetaMax(-1.0),
  alpha2(fine_structure_const*fine_structure_const)
{
  fNistManager = G4NistManager::Instance();
  fG4pow = G4Pow::GetInstance();
  fMottXSection = nullptr;
 
  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 = nullptr;
 
  fNucFormfactor = fExponentialNF;
 
  currentMaterial = nullptr;
  factB = factD = formfactA = screenZ = 0.0;
  cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = fMottFactor 
    = 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 = CLHEP::proton_mass_c2;
}

◆ ~G4WentzelOKandVIxSection()

G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection ( )

virtual

Definition at line 109 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 236 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));
}

Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeTransportXSectionPerVolume(), G4eCoulombScatteringModel::SampleSecondaries(), and G4hCoulombScatteringModel::SampleSecondaries().

◆ ComputeMaxElectronScattering()

void G4WentzelOKandVIxSection::ComputeMaxElectronScattering ( G4double cut )

protected

Definition at line 415 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 226 of file G4WentzelOKandVIxSection.hh.

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

Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeTransportXSectionPerVolume(), G4eCoulombScatteringModel::SampleSecondaries(), and G4hCoulombScatteringModel::SampleSecondaries().

◆ ComputeSecondTransportMoment()

G4double G4WentzelOKandVIxSection::ComputeSecondTransportMoment ( G4double CosThetaMax )

Definition at line 445 of file G4WentzelOKandVIxSection.cc.

{
  return 0.0;
}

◆ ComputeTransportCrossSectionPerAtom()

G4double G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom ( G4double CosThetaMax )

Definition at line 263 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; 
}

Referenced by G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), and G4WentzelVIModel::ComputeTransportXSectionPerVolume().

◆ FlatFormfactor()

G4double G4WentzelOKandVIxSection::FlatFormfactor ( G4double x )

inlineprotected

Definition at line 245 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 218 of file G4WentzelOKandVIxSection.hh.

{
  return cosTetMaxElec;
}

◆ GetCosThetaNuc()

G4double G4WentzelOKandVIxSection::GetCosThetaNuc ( ) const

inline

Definition at line 211 of file G4WentzelOKandVIxSection.hh.

{
  return cosTetMaxNuc;
}

◆ GetMomentumSquare()

G4double G4WentzelOKandVIxSection::GetMomentumSquare ( ) const

inline

Definition at line 204 of file G4WentzelOKandVIxSection.hh.

{
  return mom2;
}

Referenced by G4eCoulombScatteringModel::SampleSecondaries().

◆ Initialise()

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

Definition at line 116 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;
  */
}

Referenced by G4eCoulombScatteringModel::Initialise(), G4hCoulombScatteringModel::Initialise(), and G4WentzelVIModel::Initialise().

◆ InitialiseA()

void G4WentzelOKandVIxSection::InitialiseA ( )

protected

Definition at line 148 of file G4WentzelOKandVIxSection.cc.

{
  // Thomas-Fermi screening radii
  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
#ifdef G4MULTITHREADED
  G4MUTEXLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
  if(0.0 == ScreenRSquare[0]) {
#endif
    G4double a0 = electron_mass_c2/0.88534; 
    G4double constn = 6.937e-6/(MeV*MeV);
    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/(MeV*MeV);
 
    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;
    } 
#ifdef G4MULTITHREADED
  }
  G4MUTEXUNLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
#endif
  
  //G4cout << "G4WentzelOKandVIxSection::Initialise  mass= " << mass
  //         << "  " << p->GetParticleName() 
  //     << "  cosThetaMax= " << cosThetaMax << G4endl; 
  
}

Referenced by Initialise().

◆ operator=()

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

delete

◆ SampleSingleScattering()

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

Definition at line 352 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;
    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/MeV;
      G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
      fm = FlatFormfactor(x);
      fm *= FlatFormfactor(x*0.6
        *fG4pow->A13(fNistManager->GetAtomicMassAmu(targetZ)));
    }
    G4double grej;
    if(fMottXSection) {
      fMottXSection->SetupKinematic(tkin, targetZ);
      grej = fMottXSection->RatioMottRutherfordCosT(std::sqrt(z1))*fm*fm;
    } else {
      grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
      *fm*fm/(1.0 + z1*factD);
    }
    // G4cout << "SampleSingleScattering: E= " << tkin << " z1= " 
    //     << z1 << " grej= "<< grej << " mottFact= "<< fMottFactor<< G4endl;
    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;
}

Referenced by G4WentzelVIModel::SampleScattering(), G4eCoulombScatteringModel::SampleSecondaries(), and G4hCoulombScatteringModel::SampleSecondaries().

◆ SetTargetMass()

void G4WentzelOKandVIxSection::SetTargetMass ( G4double value )

inline

Definition at line 196 of file G4WentzelOKandVIxSection.hh.

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

Referenced by G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::DefineMaterial(), G4eCoulombScatteringModel::SampleSecondaries(), G4hCoulombScatteringModel::SampleSecondaries(), and SetupTarget().

◆ SetupKinematic()

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

Definition at line 203 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;
}

Referenced by G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeCrossSectionPerAtom(), G4WentzelVIRelModel::ComputeCrossSectionPerAtom(), G4WentzelVIModel::ComputeGeomPathLength(), G4LowEWentzelVIModel::ComputeTruePathLengthLimit(), G4WentzelVIModel::ComputeTruePathLengthLimit(), G4WentzelVIModel::ComputeTrueStepLength(), and G4eCoulombScatteringModel::SampleSecondaries().

◆ SetupParticle()

void G4WentzelOKandVIxSection::SetupParticle ( const G4ParticleDefinition * p )

Definition at line 186 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(), G4eCoulombScatteringModel::SetupParticle(), G4hCoulombScatteringModel::SetupParticle(), and G4WentzelVIModel::SetupParticle().

◆ SetupTarget()

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

Definition at line 221 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;
} 