#include <G4PairProductionRelModel.hh>

Inheritance diagram for G4PairProductionRelModel:

Public Member Functions
	G4PairProductionRelModel (const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheHeitlerLPM")
	~G4PairProductionRelModel () override
void	Initialise (const G4ParticleDefinition *, const G4DataVector &) override
void	InitialiseLocal (const G4ParticleDefinition , G4VEmModel masterModel) override
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) override
void	SampleSecondaries (std::vector< G4DynamicParticle * > , const G4MaterialCutsCouple , const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
void	SetupForMaterial (const G4ParticleDefinition , const G4Material , G4double) override
void	SetLPMflag (G4bool val)
G4bool	LPMflag () const
G4PairProductionRelModel &	operator= (const G4PairProductionRelModel &right)=delete
	G4PairProductionRelModel (const G4PairProductionRelModel &)=delete
Public Member Functions inherited from G4VEmModel
	G4VEmModel (const G4String &nam)
virtual	~G4VEmModel ()
virtual void	InitialiseForMaterial (const G4ParticleDefinition , const G4Material )
virtual void	InitialiseForElement (const G4ParticleDefinition *, G4int Z)
virtual G4double	ComputeDEDXPerVolume (const G4Material , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
virtual G4double	CrossSectionPerVolume (const G4Material , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double	GetPartialCrossSection (const G4Material , G4int level, const G4ParticleDefinition , G4double kineticEnergy)
virtual G4double	ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double	ChargeSquareRatio (const G4Track &)
virtual G4double	GetChargeSquareRatio (const G4ParticleDefinition , const G4Material , G4double kineticEnergy)
virtual G4double	GetParticleCharge (const G4ParticleDefinition , const G4Material , G4double kineticEnergy)
virtual void	StartTracking (G4Track *)
virtual void	CorrectionsAlongStep (const G4Material , const G4ParticleDefinition , const G4double kinEnergy, const G4double cutEnergy, const G4double &length, G4double &eloss)
virtual G4double	Value (const G4MaterialCutsCouple , const G4ParticleDefinition , G4double kineticEnergy)
virtual G4double	MinPrimaryEnergy (const G4Material , const G4ParticleDefinition , G4double cut=0.0)
virtual G4double	MinEnergyCut (const G4ParticleDefinition , const G4MaterialCutsCouple )
virtual void	DefineForRegion (const G4Region *)
virtual void	FillNumberOfSecondaries (G4int &numberOfTriplets, G4int &numberOfRecoil)
virtual void	ModelDescription (std::ostream &outFile) const
void	InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
std::vector< G4EmElementSelector * > *	GetElementSelectors ()
void	SetElementSelectors (std::vector< G4EmElementSelector * > *)
G4double	ComputeDEDX (const G4MaterialCutsCouple , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double	CrossSection (const G4MaterialCutsCouple , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeMeanFreePath (const G4ParticleDefinition , G4double kineticEnergy, const G4Material , G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition , const G4Element , G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectRandomAtom (const G4MaterialCutsCouple , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectTargetAtom (const G4MaterialCutsCouple , const G4ParticleDefinition , G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectRandomAtom (const G4Material , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	GetCurrentElement (const G4Material *mat=nullptr) const
G4int	SelectRandomAtomNumber (const G4Material *) const
const G4Isotope *	GetCurrentIsotope (const G4Element *elm=nullptr) const
G4int	SelectIsotopeNumber (const G4Element *) const
void	SetParticleChange (G4VParticleChange , G4VEmFluctuationModel f=nullptr)
void	SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
G4ElementData *	GetElementData ()
G4PhysicsTable *	GetCrossSectionTable ()
G4VEmFluctuationModel *	GetModelOfFluctuations ()
G4VEmAngularDistribution *	GetAngularDistribution ()
G4VEmModel *	GetTripletModel ()
void	SetTripletModel (G4VEmModel *)
void	SetAngularDistribution (G4VEmAngularDistribution *)
G4double	HighEnergyLimit () const
G4double	LowEnergyLimit () const
G4double	HighEnergyActivationLimit () const
G4double	LowEnergyActivationLimit () const
G4double	PolarAngleLimit () const
G4double	SecondaryThreshold () const
G4bool	DeexcitationFlag () const
G4bool	ForceBuildTableFlag () const
G4bool	UseAngularGeneratorFlag () const
void	SetAngularGeneratorFlag (G4bool)
void	SetHighEnergyLimit (G4double)
void	SetLowEnergyLimit (G4double)
void	SetActivationHighEnergyLimit (G4double)
void	SetActivationLowEnergyLimit (G4double)
G4bool	IsActive (G4double kinEnergy) const
void	SetPolarAngleLimit (G4double)
void	SetSecondaryThreshold (G4double)
void	SetDeexcitationFlag (G4bool val)
void	SetForceBuildTable (G4bool val)
void	SetFluctuationFlag (G4bool val)
G4bool	IsMaster () const
void	SetUseBaseMaterials (G4bool val)
G4bool	UseBaseMaterials () const
G4double	MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4String &	GetName () const
void	SetCurrentCouple (const G4MaterialCutsCouple *)
G4bool	IsLocked () const
void	SetLocked (G4bool)
void	SetLPMFlag (G4bool)
void	SetMasterThread (G4bool)
G4VEmModel &	operator= (const G4VEmModel &right)=delete
	G4VEmModel (const G4VEmModel &)=delete

Protected Member Functions
void	ComputePhi12 (const G4double delta, G4double &phi1, G4double &phi2)
G4double	ScreenFunction1 (const G4double delta)
G4double	ScreenFunction2 (const G4double delta)
void	ScreenFunction12 (const G4double delta, G4double &f1, G4double &f2)
G4double	ComputeParametrizedXSectionPerAtom (G4double gammaEnergy, G4double Z)
G4double	ComputeXSectionPerAtom (G4double gammaEnergy, G4double Z)
G4double	ComputeDXSectionPerAtom (G4double eplusEnergy, G4double gammaEnergy, G4double Z)
G4double	ComputeRelDXSectionPerAtom (G4double eplusEnergy, G4double gammaEnergy, G4double Z)
Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss *	GetParticleChangeForLoss ()
G4ParticleChangeForGamma *	GetParticleChangeForGamma ()
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCouple *	CurrentCouple () const
void	SetCurrentElement (const G4Element *)

Protected Attributes
G4bool	isFirstInstance {false}
G4bool	fIsUseLPMCorrection
G4bool	fIsUseCompleteScreening
G4double	fLPMEnergy
G4double	fParametrizedXSectionThreshold
G4double	fCoulombCorrectionThreshold
G4Pow *	fG4Calc
G4ParticleDefinition *	fTheGamma
G4ParticleDefinition *	fTheElectron
G4ParticleDefinition *	fThePositron
G4ParticleChangeForGamma *	fParticleChange
Protected Attributes inherited from G4VEmModel
G4ElementData *	fElementData = nullptr
G4VParticleChange *	pParticleChange = nullptr
G4PhysicsTable *	xSectionTable = nullptr
const G4Material *	pBaseMaterial = nullptr
const std::vector< G4double > *	theDensityFactor = nullptr
const std::vector< G4int > *	theDensityIdx = nullptr
G4double	inveplus
G4double	pFactor = 1.0
std::size_t	currentCoupleIndex = 0
std::size_t	basedCoupleIndex = 0
G4bool	lossFlucFlag = true

Static Protected Attributes
static const G4int	gMaxZet = 120
static const G4double	gLPMconstant
static const G4double	gXGL [8]
static const G4double	gWGL [8]
static const G4double	gFelLowZet [8]
static const G4double	gFinelLowZet [8]
static const G4double	gXSecFactor
static const G4double	gEgLPMActivation = 100.*CLHEP::GeV
static std::vector< ElementData * >	gElementData

Detailed Description

Definition at line 65 of file G4PairProductionRelModel.hh.

Constructor & Destructor Documentation

◆ G4PairProductionRelModel() [1/2]

G4PairProductionRelModel::G4PairProductionRelModel	(	const G4ParticleDefinition *	p = nullptr,
		const G4String &	nam = "BetheHeitlerLPM" )

explicit

Definition at line 118 of file G4PairProductionRelModel.cc.

  : G4VEmModel(nam), fIsUseLPMCorrection(true), fIsUseCompleteScreening(false),
  fLPMEnergy(0.), fG4Calc(G4Pow::GetInstance()), fTheGamma(G4Gamma::Gamma()),
  fTheElectron(G4Electron::Electron()), fThePositron(G4Positron::Positron()),
  fParticleChange(nullptr)
{
  // gamma energy below which the parametrized atomic x-section is used (30 GeV)
  fParametrizedXSectionThreshold = 30.0*CLHEP::GeV;
  // gamma energy below the Coulomb correction is turned off (50 MeV)
  fCoulombCorrectionThreshold    = 50.0*CLHEP::MeV;
  // set angular generator used in the final state kinematics computation
  SetAngularDistribution(new G4ModifiedTsai());
}

Referenced by G4BetheHeitler5DModel::G4BetheHeitler5DModel(), G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(), G4PairProductionRelModel(), and operator=().

◆ ~G4PairProductionRelModel()

G4PairProductionRelModel::~G4PairProductionRelModel ( )

override

Definition at line 134 of file G4PairProductionRelModel.cc.

{
  if (isFirstInstance) {
    // clear ElementData container
    for (auto const & ptr : gElementData) { delete ptr; }
    gElementData.clear();
  }
}

◆ G4PairProductionRelModel() [2/2]

G4PairProductionRelModel::G4PairProductionRelModel ( const G4PairProductionRelModel & )

delete

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

G4double G4PairProductionRelModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	kinEnergy,
		G4double	Z,
		G4double	A = 0.,
		G4double	cut = 0.,
		G4double	emax = DBL_MAX )

overridevirtual

Reimplemented from G4VEmModel.

Definition at line 309 of file G4PairProductionRelModel.cc.

{
  G4double crossSection = 0.0 ;
  // check kinematical limit
  if ( gammaEnergy <= 2.0*electron_mass_c2 ) { return crossSection; }
  // compute the atomic cross section either by using x-section parametrization
  // or by numerically integrationg the DCS (with or without LPM)
  if ( gammaEnergy < fParametrizedXSectionThreshold) {
    // using the parametrized cross sections (max up to 80 GeV)
    crossSection = ComputeParametrizedXSectionPerAtom(gammaEnergy, Z);
  } else {
    // by numerical integration of the DCS:
    // Computes the cross section with or without LPM suppression depending on
    // settings (by default with if the gamma energy is above a given threshold)
    // and using or not using complete sreening approximation (by default not).
    // Only the dependent part is computed in the numerical integration of the DCS
    // i.e. the result must be multiplied here with 4 \alpha r_0^2 Z(Z+\eta(Z))
    crossSection = ComputeXSectionPerAtom(gammaEnergy, Z);
    // apply the constant factors:
    // - eta(Z) is a correction to account interaction in the field of e-
    // - gXSecFactor = 4 \alpha r_0^2
    const G4int iz     = std::min(gMaxZet, G4lrint(Z));
    const G4double eta = gElementData[iz]->fEtaValue;
    crossSection      *= gXSecFactor*Z*(Z+eta);
  }
  // final protection
  return std::max(crossSection, 0.);
}

◆ ComputeDXSectionPerAtom()

G4double G4PairProductionRelModel::ComputeDXSectionPerAtom	(	G4double	eplusEnergy,
		G4double	gammaEnergy,
		G4double	Z )

protected

Definition at line 226 of file G4PairProductionRelModel.cc.

{
  G4double xSection = 0.;
  const G4int    iz   = std::min(gMaxZet, G4lrint(Z));
  const G4double eps  = pEnergy/gammaEnergy;
  const G4double epsm = 1.-eps;
  const G4double dum  = eps*epsm;
  if (fIsUseCompleteScreening) {
    // complete screening:
    const G4double Lel   = gElementData[iz]->fLradEl;
    const G4double fc    = gElementData[iz]->fCoulomb;
    xSection = (eps*eps + epsm*epsm + 2.*dum/3.)*(Lel-fc) - dum/9.;
  } else {
    // normal case:
    const G4double eps0  = CLHEP::electron_mass_c2/gammaEnergy;
    const G4double fc    = gElementData[iz]->fCoulomb;
    const G4double lnZ13 = gElementData[iz]->fLogZ13;
    const G4double delta = gElementData[iz]->fDeltaFactor*eps0/dum;
    G4double phi1, phi2;
    ComputePhi12(delta, phi1, phi2);
    xSection =  (eps*eps + epsm*epsm)*(0.25*phi1-lnZ13-fc)
               + 2.*dum*(0.25*phi2-lnZ13-fc)/3.;
  }
  // non-const. part of the DCS differential in total energy transfer not in eps
  // ds/dEt=ds/deps deps/dEt with deps/dEt=1/Eg
  return std::max(xSection, 0.0)/gammaEnergy;
}

Referenced by ComputeXSectionPerAtom().

◆ ComputeParametrizedXSectionPerAtom()

G4double G4PairProductionRelModel::ComputeParametrizedXSectionPerAtom	(	G4double	gammaEnergy,
		G4double	Z )

protected

Definition at line 572 of file G4PairProductionRelModel.cc.

{
  G4double xSection = 0.0 ;
  // short versions
  static const G4double kMC2  = CLHEP::electron_mass_c2;
  // zero cross section below the kinematical limit: Eg<2mc^2
  if (Z < 0.9 || gammaE <= 2.0*kMC2) { return xSection; }
  //
  static const G4double gammaEnergyLimit = 1.5*CLHEP::MeV;
  // set coefficients a, b c
  static const G4double a0 =  8.7842e+2*CLHEP::microbarn;
  static const G4double a1 = -1.9625e+3*CLHEP::microbarn;
  static const G4double a2 =  1.2949e+3*CLHEP::microbarn;
  static const G4double a3 = -2.0028e+2*CLHEP::microbarn;
  static const G4double a4 =  1.2575e+1*CLHEP::microbarn;
  static const G4double a5 = -2.8333e-1*CLHEP::microbarn;
 
  static const G4double b0 = -1.0342e+1*CLHEP::microbarn;
  static const G4double b1 =  1.7692e+1*CLHEP::microbarn;
  static const G4double b2 = -8.2381   *CLHEP::microbarn;
  static const G4double b3 =  1.3063   *CLHEP::microbarn;
  static const G4double b4 = -9.0815e-2*CLHEP::microbarn;
  static const G4double b5 =  2.3586e-3*CLHEP::microbarn;
 
  static const G4double c0 = -4.5263e+2*CLHEP::microbarn;
  static const G4double c1 =  1.1161e+3*CLHEP::microbarn;
  static const G4double c2 = -8.6749e+2*CLHEP::microbarn;
  static const G4double c3 =  2.1773e+2*CLHEP::microbarn;
  static const G4double c4 = -2.0467e+1*CLHEP::microbarn;
  static const G4double c5 =  6.5372e-1*CLHEP::microbarn;
  // check low energy limit of the approximation (1.5 MeV)
  G4double gammaEnergyOrg = gammaE;
  if (gammaE < gammaEnergyLimit) { gammaE = gammaEnergyLimit; }
  // compute gamma energy variables
  const G4double x  = G4Log(gammaE/kMC2);
  const G4double x2 = x *x;
  const G4double x3 = x2*x;
  const G4double x4 = x3*x;
  const G4double x5 = x4*x;
  //
  const G4double F1 = a0 + a1*x + a2*x2 + a3*x3 + a4*x4 + a5*x5;
  const G4double F2 = b0 + b1*x + b2*x2 + b3*x3 + b4*x4 + b5*x5;
  const G4double F3 = c0 + c1*x + c2*x2 + c3*x3 + c4*x4 + c5*x5;
  // compute the approximated cross section
  xSection = (Z + 1.)*(F1*Z + F2*Z*Z + F3);
  // check if we are below the limit of the approximation and apply correction
  if (gammaEnergyOrg < gammaEnergyLimit) {
    const G4double dum = (gammaEnergyOrg-2.*kMC2)/(gammaEnergyLimit-2.*kMC2);
    xSection *= dum*dum;
  }
  return xSection;
}

Referenced by ComputeCrossSectionPerAtom().

◆ ComputePhi12()

void G4PairProductionRelModel::ComputePhi12	(	const G4double	delta,
		G4double &	phi1,
		G4double &	phi2 )

inlineprotected

Definition at line 186 of file G4PairProductionRelModel.hh.

{
    if (delta > 1.4) {
      phi1 = 21.0190 - 4.145*G4Log(delta + 0.958);
      phi2 = phi1;
    } else {
      phi1 = 20.806 - delta*(3.190 - 0.5710*delta);
      phi2 = 20.234 - delta*(2.126 - 0.0903*delta);
    }
}

Referenced by ComputeDXSectionPerAtom(), ComputeRelDXSectionPerAtom(), and SampleSecondaries().

◆ ComputeRelDXSectionPerAtom()

G4double G4PairProductionRelModel::ComputeRelDXSectionPerAtom	(	G4double	eplusEnergy,
		G4double	gammaEnergy,
		G4double	Z )

protected

Definition at line 275 of file G4PairProductionRelModel.cc.

{
  G4double xSection = 0.;
  const G4int    iz   = std::min(gMaxZet, G4lrint(Z));
  const G4double eps  = pEnergy/gammaEnergy;
  const G4double epsm = 1.-eps;
  const G4double dum  = eps*epsm;
  // evaluate LPM suppression functions
  G4double fXiS, fGS, fPhiS;
  ComputeLPMfunctions(fXiS, fGS, fPhiS, eps, gammaEnergy, iz);
  if (fIsUseCompleteScreening) {
    // complete screening:
    const G4double Lel   = gElementData[iz]->fLradEl;
    const G4double fc    = gElementData[iz]->fCoulomb;
    xSection = (Lel-fc)*((eps*eps+epsm*epsm)*2.*fPhiS + fGS)/3. - dum*fGS/9.;
  } else {
    // normal case:
    const G4double eps0  = CLHEP::electron_mass_c2/gammaEnergy;
    const G4double fc    = gElementData[iz]->fCoulomb;
    const G4double lnZ13 = gElementData[iz]->fLogZ13;
    const G4double delta = gElementData[iz]->fDeltaFactor*eps0/dum;
    G4double phi1, phi2;
    ComputePhi12(delta, phi1, phi2);
    xSection =  (eps*eps + epsm*epsm)*(2.*fPhiS+fGS)*(0.25*phi1-lnZ13-fc)/3.
               + 2.*dum*fGS*(0.25*phi2-lnZ13-fc)/3.;
  }
  // non-const. part of the DCS differential in total energy transfer not in eps
  // ds/dEt=ds/deps deps/dEt with deps/dEt=1/Eg
  return std::max(fXiS*xSection, 0.0)/gammaEnergy;
}

Referenced by ComputeXSectionPerAtom().

◆ ComputeXSectionPerAtom()

G4double G4PairProductionRelModel::ComputeXSectionPerAtom	(	G4double	gammaEnergy,
		G4double	Z )

protected

Definition at line 171 of file G4PairProductionRelModel.cc.

{
  G4double xSection = 0.0;
  // check if LPM suppression needs to be used
  const G4bool   isLPM  = (fIsUseLPMCorrection && gammaEnergy>gEgLPMActivation);
  // determine the kinematical limits (taken into account the correction due to
  // the way in which the Coulomb correction is applied i.e. avoid negative DCS)
  const G4int    iz     = std::min(gMaxZet, G4lrint(Z));
  const G4double eps0   = CLHEP::electron_mass_c2/gammaEnergy;
  // Coulomb correction is always included in the DCS even below 50 MeV (note:
  // that this DCS is only used to get the integrated x-section)
  const G4double dmax   = gElementData[iz]->fDeltaMaxHigh;
  const G4double dmin   = 4.*eps0*gElementData[iz]->fDeltaFactor;
  const G4double eps1   = 0.5 - 0.5*std::sqrt(1.-dmin/dmax);
  const G4double epsMin = std::max(eps0, eps1);
  const G4double epsMax = 0.5; // DCS is symmetric around eps=0.5
  // let Et be the total energy transferred to the e- or to the e+
  // the [Et-min, Et-max] interval will be divided into i=1,2,..,n subintervals
  // with width of dInterv = (Et-max - Et-min)/n and numerical integration will
  // be done in each sub-inteval using the xi = (Et - Et_i-min)/dInterv variable
  // that is in [0,1]. The 8-point GL q. is used for the integration on [0,1].
  const G4int    numSub = 2;
  const G4double dInterv= (epsMax - epsMin)*gammaEnergy/G4double(numSub);
  G4double minEti       = epsMin*gammaEnergy;  // Et-min i.e. Et_0-min
  for (G4int i = 0; i < numSub; ++i) {
    for (G4int ngl = 0; ngl < 8; ++ngl) {
      const G4double Et = (minEti + gXGL[ngl]*dInterv);
      const G4double xs = isLPM ? ComputeRelDXSectionPerAtom(Et, gammaEnergy, Z)
                                : ComputeDXSectionPerAtom(Et, gammaEnergy, Z);
      xSection += gWGL[ngl]*xs;
    }
    // update minimum Et of the sub-inteval
    minEti += dInterv;
  }
  // apply corrections of variable transformation and half interval integration
  xSection = std::max(2.*xSection*dInterv, 0.);
  return xSection;
}

Referenced by ComputeCrossSectionPerAtom().

◆ Initialise()

void G4PairProductionRelModel::Initialise	(	const G4ParticleDefinition *	p,
		const G4DataVector &	cuts )

overridevirtual

Implements G4VEmModel.

Definition at line 143 of file G4PairProductionRelModel.cc.

{
  if(nullptr == fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
 
  if (isFirstInstance || gElementData.empty()) {
    // init element data and LPM funcs
    G4AutoLock l(&thePairProdRelMutex);
    if (gElementData.empty()) {
      isFirstInstance = true;
      gElementData.resize(gMaxZet+1, nullptr);
    }
    // static data should be initialised only in the one instance
    InitialiseElementData();
    l.unlock();
  }
  // element selectors should be initialised in the master thread
  if (IsMaster()) {
    InitialiseElementSelectors(p, cuts);
  }
}

Referenced by G4BetheHeitler5DModel::Initialise(), and G4LivermoreGammaConversionModel::Initialise().

◆ InitialiseLocal()

void G4PairProductionRelModel::InitialiseLocal	(	const G4ParticleDefinition *	,
		G4VEmModel *	masterModel )

overridevirtual

Reimplemented from G4VEmModel.

Definition at line 165 of file G4PairProductionRelModel.cc.

{
  SetElementSelectors(masterModel->GetElementSelectors());
}

◆ LPMflag()

G4bool G4PairProductionRelModel::LPMflag ( ) const

inline

Definition at line 97 of file G4PairProductionRelModel.hh.

97{ return fIsUseLPMCorrection; }

◆ operator=()

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

delete

◆ SampleSecondaries()

void G4PairProductionRelModel::SampleSecondaries	(	std::vector< G4DynamicParticle * > *	fvect,
		const G4MaterialCutsCouple *	couple,
		const G4DynamicParticle *	aDynamicGamma,
		G4double	tmin,
		G4double	maxEnergy )

overridevirtual

Implements G4VEmModel.

Definition at line 346 of file G4PairProductionRelModel.cc.

{
  const G4Material* mat         = couple->GetMaterial();
  const G4double    gammaEnergy = aDynamicGamma->GetKineticEnergy();
  const G4double    eps0        = CLHEP::electron_mass_c2/gammaEnergy ;
  //
  // check kinematical limit: gamma energy(Eg) must be at least 2 e- rest mass
  // (but the model should be used at higher energies above 100 MeV)
  if (eps0 > 0.5) { return; }
  //
  // select target atom of the material
  const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, gammaEnergy,
                                         aDynamicGamma->GetLogKineticEnergy());
  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
  //
  // 'eps' is the total energy transferred to one of the e-/e+ pair in initial
  // gamma energy units Eg. Since the corresponding DCS is symmetric on eps=0.5,
  // the kinematical limits for eps0=mc^2/Eg <= eps <= 0.5
  // 1. 'eps' is sampled uniformly on the [eps0, 0.5] inteval if Eg<Egsmall
  // 2. otherwise, on the [eps_min, 0.5] interval according to the DCS (case 2.)
  G4double eps;
  // case 1.
  static const G4double Egsmall = 2.*CLHEP::MeV;
  if (gammaEnergy < Egsmall) {
    eps = eps0 + (0.5-eps0)*rndmEngine->flat();
  } else {
  // case 2.
    // get the Coulomb factor for the target element (Z) and gamma energy (Eg)
    // F(Z) = 8*ln(Z)/3           if Eg <= 50 [MeV] => no Coulomb correction
    // F(Z) = 8*ln(Z)/3 + 8*fc(Z) if Eg  > 50 [MeV] => fc(Z) is the Coulomb cor.
    //
    // The screening variable 'delta(eps)' = 136*Z^{-1/3}*eps0/[eps(1-eps)]
    // Due to the Coulomb correction, the DCS can go below zero even at
    // kinematicaly allowed eps > eps0 values. In order to exclude this eps
    // range with negative DCS, the minimum eps value will be set to eps_min =
    // max[eps0, epsp] with epsp is the solution of SF(delta(epsp)) - F(Z)/2 = 0
    // with SF being the screening function (SF1=SF2 at high value of delta).
    // The solution is epsp = 0.5 - 0.5*sqrt[ 1 - 4*136*Z^{-1/3}eps0/deltap]
    // with deltap = Exp[(42.038-F(Z))/8.29]-0.958. So the limits are:
    // - when eps=eps_max = 0.5            => delta_min = 136*Z^{-1/3}*eps0/4
    // - epsp = 0.5 - 0.5*sqrt[ 1 - delta_min/deltap]
    // - and eps_min = max[eps0, epsp]
    const G4int    iZet        = std::min(gMaxZet, anElement->GetZasInt());
    const G4double deltaFactor = gElementData[iZet]->fDeltaFactor*eps0;
    const G4double deltaMin    = 4.*deltaFactor;
    G4double       deltaMax    = gElementData[iZet]->fDeltaMaxLow;
    G4double       FZ          = 8.*gElementData[iZet]->fLogZ13;
    if ( gammaEnergy > fCoulombCorrectionThreshold ) {   // Eg > 50 MeV ?
      FZ      += 8.*gElementData[iZet]->fCoulomb;
      deltaMax = gElementData[iZet]->fDeltaMaxHigh;
    }
    // compute the limits of eps
    const G4double epsp        = 0.5 - 0.5*std::sqrt(1. - deltaMin/deltaMax) ;
    const G4double epsMin      = std::max(eps0,epsp);
    const G4double epsRange    = 0.5 - epsMin;
    //
    // sample the energy rate (eps) of the created electron (or positron)
    G4double F10, F20;
    ScreenFunction12(deltaMin, F10, F20);
    F10 -= FZ;
    F20 -= FZ;
    const G4double NormF1   = std::max(F10 * epsRange * epsRange, 0.);
    const G4double NormF2   = std::max(1.5 * F20                , 0.);
    const G4double NormCond = NormF1/(NormF1 + NormF2);
    // check if LPM correction is active
    const G4bool isLPM = (fIsUseLPMCorrection && gammaEnergy>gEgLPMActivation);
    fLPMEnergy = mat->GetRadlen()*gLPMconstant;
    // we will need 3 uniform random number for each trial of sampling
    G4double rndmv[3];
    G4double greject = 0.;
    do {
      rndmEngine->flatArray(3, rndmv);
      if (NormCond > rndmv[0]) {
        eps = 0.5 - epsRange * fG4Calc->A13(rndmv[1]);
        const G4double delta = deltaFactor/(eps*(1.-eps));
        if (isLPM) {
          G4double lpmXiS, lpmGS, lpmPhiS, phi1, phi2;
          ComputePhi12(delta, phi1, phi2);
          ComputeLPMfunctions(lpmXiS, lpmGS, lpmPhiS, eps, gammaEnergy, iZet);
          greject = lpmXiS*((2.*lpmPhiS+lpmGS)*phi1-lpmGS*phi2-lpmPhiS*FZ)/F10;
        } else {
          greject = (ScreenFunction1(delta)-FZ)/F10;
        }
      } else {
        eps = epsMin + epsRange*rndmv[1];
        const G4double delta = deltaFactor/(eps*(1.-eps));
        if (isLPM) {
          G4double lpmXiS, lpmGS, lpmPhiS, phi1, phi2;
          ComputePhi12(delta, phi1, phi2);
          ComputeLPMfunctions(lpmXiS, lpmGS, lpmPhiS, eps, gammaEnergy, iZet);
          greject = lpmXiS*( (lpmPhiS+0.5*lpmGS)*phi1 + 0.5*lpmGS*phi2
                             -0.5*(lpmGS+lpmPhiS)*FZ )/F20;
        } else {
          greject = (ScreenFunction2(delta)-FZ)/F20;
        }
      }
      // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
    } while (greject < rndmv[2]);
    //  end of eps sampling
  }
  //
  // select charges randomly
  G4double eTotEnergy, pTotEnergy;
  if (rndmEngine->flat() > 0.5) {
    eTotEnergy = (1.-eps)*gammaEnergy;
    pTotEnergy = eps*gammaEnergy;
  } else {
    pTotEnergy = (1.-eps)*gammaEnergy;
    eTotEnergy = eps*gammaEnergy;
  }
  //
  // sample pair kinematics
  //
  const G4double eKinEnergy = std::max(0.,eTotEnergy - CLHEP::electron_mass_c2);
  const G4double pKinEnergy = std::max(0.,pTotEnergy - CLHEP::electron_mass_c2);
  //
  G4ThreeVector eDirection, pDirection;
  //
  GetAngularDistribution()->SamplePairDirections(aDynamicGamma,
                         eKinEnergy, pKinEnergy, eDirection, pDirection);
  // create G4DynamicParticle object for the particle1
  auto aParticle1 = new G4DynamicParticle(fTheElectron,eDirection,eKinEnergy);
 
  // create G4DynamicParticle object for the particle2
  auto aParticle2 = new G4DynamicParticle(fThePositron,pDirection,pKinEnergy);
  // Fill output vector
  fvect->push_back(aParticle1);
  fvect->push_back(aParticle2);
  // kill incident photon
  fParticleChange->SetProposedKineticEnergy(0.);
  fParticleChange->ProposeTrackStatus(fStopAndKill);
}

◆ ScreenFunction1()

G4double G4PairProductionRelModel::ScreenFunction1 ( const G4double delta )

inlineprotected

Definition at line 200 of file G4PairProductionRelModel.hh.

{
  return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
                       : 42.184 - delta*(7.444 - 1.623*delta);
}

Referenced by SampleSecondaries().

◆ ScreenFunction12()

void G4PairProductionRelModel::ScreenFunction12	(	const G4double	delta,
		G4double &	f1,
		G4double &	f2 )

inlineprotected

Definition at line 214 of file G4PairProductionRelModel.hh.

{
  if (delta > 1.4) {
    f1 = 42.038 - 8.29*G4Log(delta + 0.958);
    f2 = f1;
  } else {
    f1 = 42.184 - delta*(7.444 - 1.623*delta);
    f2 = 41.326 - delta*(5.848 - 0.902*delta);
  }
}

Referenced by SampleSecondaries().

◆ ScreenFunction2()

G4double G4PairProductionRelModel::ScreenFunction2 ( const G4double delta )

inlineprotected

Definition at line 207 of file G4PairProductionRelModel.hh.

{
  return (delta > 1.4) ? 42.038 - 8.29*G4Log(delta + 0.958)
                       : 41.326 - delta*(5.848 - 0.902*delta);
}

Referenced by SampleSecondaries().

◆ SetLPMflag()

void G4PairProductionRelModel::SetLPMflag ( G4bool val )

inline

Definition at line 96 of file G4PairProductionRelModel.hh.

96{ fIsUseLPMCorrection = val; }

◆ SetupForMaterial()

void G4PairProductionRelModel::SetupForMaterial	(	const G4ParticleDefinition *	,
		const G4Material *	mat,
		G4double	)

overridevirtual

Reimplemented from G4VEmModel.

Definition at line 339 of file G4PairProductionRelModel.cc.

{
  fLPMEnergy = mat->GetRadlen()*gLPMconstant;
}

Member Data Documentation

◆ fCoulombCorrectionThreshold

G4double G4PairProductionRelModel::fCoulombCorrectionThreshold

protected

Definition at line 161 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), and SampleSecondaries().

◆ fG4Calc

G4Pow* G4PairProductionRelModel::fG4Calc

protected

Definition at line 163 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), and SampleSecondaries().

◆ fIsUseCompleteScreening

G4bool G4PairProductionRelModel::fIsUseCompleteScreening

protected

Definition at line 156 of file G4PairProductionRelModel.hh.

Referenced by ComputeDXSectionPerAtom(), ComputeRelDXSectionPerAtom(), and G4PairProductionRelModel().

◆ fIsUseLPMCorrection

G4bool G4PairProductionRelModel::fIsUseLPMCorrection

protected

Definition at line 155 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom(), G4PairProductionRelModel(), LPMflag(), and SampleSecondaries().

◆ fLPMEnergy

G4double G4PairProductionRelModel::fLPMEnergy

protected

Definition at line 158 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), SampleSecondaries(), and SetupForMaterial().

◆ fParametrizedXSectionThreshold

G4double G4PairProductionRelModel::fParametrizedXSectionThreshold

protected

Definition at line 160 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom(), and G4PairProductionRelModel().

◆ fParticleChange

G4ParticleChangeForGamma* G4PairProductionRelModel::fParticleChange

protected

Definition at line 167 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), Initialise(), G4BetheHeitler5DModel::SampleSecondaries(), and SampleSecondaries().

◆ fTheElectron

G4ParticleDefinition* G4PairProductionRelModel::fTheElectron

protected

Definition at line 165 of file G4PairProductionRelModel.hh.

Referenced by G4BetheHeitler5DModel::G4BetheHeitler5DModel(), G4PairProductionRelModel(), G4BetheHeitler5DModel::SampleSecondaries(), and SampleSecondaries().

◆ fTheGamma

G4ParticleDefinition* G4PairProductionRelModel::fTheGamma

protected

Definition at line 164 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), G4BetheHeitler5DModel::SampleSecondaries(), and SampleSecondaries().

◆ fThePositron

G4ParticleDefinition* G4PairProductionRelModel::fThePositron

protected

Definition at line 166 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), and SampleSecondaries().

◆ gEgLPMActivation

const G4double G4PairProductionRelModel::gEgLPMActivation = 100.*CLHEP::GeV

staticprotected

Definition at line 150 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom(), and SampleSecondaries().

◆ gElementData

std::vector< G4PairProductionRelModel::ElementData * > G4PairProductionRelModel::gElementData

staticprotected

Definition at line 152 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom(), ComputeDXSectionPerAtom(), ComputeRelDXSectionPerAtom(), ComputeXSectionPerAtom(), Initialise(), SampleSecondaries(), and ~G4PairProductionRelModel().

◆ gFelLowZet

const G4double G4PairProductionRelModel::gFelLowZet

staticprotected

Initial value:

                                                      = {
  0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 4.6134, 4.5520
}

Definition at line 93 of file G4PairProductionRelModel.hh.

◆ gFinelLowZet

const G4double G4PairProductionRelModel::gFinelLowZet

staticprotected

Initial value:

                                                      = {
  0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 5.3688, 5.3236
}

Definition at line 96 of file G4PairProductionRelModel.hh.

                                       { fIsUseLPMCorrection = val;  }
  inline G4bool LPMflag() const        { return fIsUseLPMCorrection; }
 

◆ gLPMconstant

const G4double G4PairProductionRelModel::gLPMconstant

staticprotected

Initial value:

                                                    =
CLHEP::fine_structure_const*CLHEP::electron_mass_c2*CLHEP::electron_mass_c2
/(4.*CLHEP::pi*CLHEP::hbarc)

Definition at line 142 of file G4PairProductionRelModel.hh.

Referenced by SampleSecondaries(), and SetupForMaterial().

◆ gMaxZet

const G4int G4PairProductionRelModel::gMaxZet = 120

staticprotected

Definition at line 140 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom(), ComputeDXSectionPerAtom(), ComputeRelDXSectionPerAtom(), ComputeXSectionPerAtom(), Initialise(), and SampleSecondaries().

◆ gWGL

const G4double G4PairProductionRelModel::gWGL

staticprotected

Initial value:

                                              = {
  5.06142681e-02, 1.11190517e-01, 1.56853323e-01, 1.81341892e-01,
  1.81341892e-01, 1.56853323e-01, 1.11190517e-01, 5.06142681e-02
}

Definition at line 86 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom().

◆ gXGL

const G4double G4PairProductionRelModel::gXGL

staticprotected

Initial value:

                                              = {
  1.98550718e-02, 1.01666761e-01, 2.37233795e-01, 4.08282679e-01,
  5.91717321e-01, 7.62766205e-01, 8.98333239e-01, 9.80144928e-01
}

Definition at line 82 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom().

◆ gXSecFactor

const G4double G4PairProductionRelModel::gXSecFactor

staticprotected

Initial value:

                                                   =
4.*CLHEP::fine_structure_const*CLHEP::classic_electr_radius
*CLHEP::classic_electr_radius

Definition at line 149 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom().

◆ isFirstInstance

G4bool G4PairProductionRelModel::isFirstInstance {false}

protected

Definition at line 154 of file G4PairProductionRelModel.hh.

154{false};

Referenced by Initialise(), and ~G4PairProductionRelModel().

The documentation for this class was generated from the following files:

Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4PairProductionRelModel() [1/2]

◆ ~G4PairProductionRelModel()

◆ G4PairProductionRelModel() [2/2]

Member Function Documentation

◆ ComputeCrossSectionPerAtom()

◆ ComputeDXSectionPerAtom()

◆ ComputeParametrizedXSectionPerAtom()

◆ ComputePhi12()

◆ ComputeRelDXSectionPerAtom()

◆ ComputeXSectionPerAtom()

◆ Initialise()

◆ InitialiseLocal()

◆ LPMflag()

◆ operator=()

◆ SampleSecondaries()

◆ ScreenFunction1()

◆ ScreenFunction12()

◆ ScreenFunction2()

◆ SetLPMflag()

◆ SetupForMaterial()

Member Data Documentation

◆ fCoulombCorrectionThreshold

◆ fG4Calc

◆ fIsUseCompleteScreening

◆ fIsUseLPMCorrection

◆ fLPMEnergy

◆ fParametrizedXSectionThreshold

◆ fParticleChange

◆ fTheElectron

◆ fTheGamma

◆ fThePositron

◆ gEgLPMActivation

◆ gElementData

◆ gFelLowZet

◆ gFinelLowZet

◆ gLPMconstant

◆ gMaxZet

◆ gWGL

◆ gXGL

◆ gXSecFactor

◆ isFirstInstance