#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 G4MaterialCutsCouple , const G4DynamicParticle , 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)

void	SetMasterThread (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)

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

static LPMFuncs	gLPMFuncs

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 120 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());
}

◆ ~G4PairProductionRelModel()

G4PairProductionRelModel::~G4PairProductionRelModel ( )

override

Definition at line 136 of file G4PairProductionRelModel.cc.

{
  if (isFirstInstance) {
    // clear ElementData container
    for (auto const & ptr : gElementData) { delete ptr; }
    gElementData.clear();
    // clear LPMFunctions (if any)
    if (fIsUseLPMCorrection) {
      gLPMFuncs.fLPMFuncG.clear();
      gLPMFuncs.fLPMFuncPhi.clear();
      gLPMFuncs.fIsInitialized = false;
    }
  }
}

◆ 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 320 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 237 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 657 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 201 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 286 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 182 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 151 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();
    if (fIsUseLPMCorrection) {
      InitLPMFunctions();
    }
    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 176 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 357 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 215 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 229 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 222 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 350 of file G4PairProductionRelModel.cc.

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

Member Data Documentation

◆ fCoulombCorrectionThreshold

G4double G4PairProductionRelModel::fCoulombCorrectionThreshold

protected

Definition at line 176 of file G4PairProductionRelModel.hh.

Referenced by G4PairProductionRelModel(), and SampleSecondaries().

◆ fG4Calc

G4Pow* G4PairProductionRelModel::fG4Calc

protected

Definition at line 178 of file G4PairProductionRelModel.hh.

Referenced by SampleSecondaries().

◆ fIsUseCompleteScreening

G4bool G4PairProductionRelModel::fIsUseCompleteScreening

protected

Definition at line 171 of file G4PairProductionRelModel.hh.

Referenced by ComputeDXSectionPerAtom(), and ComputeRelDXSectionPerAtom().

◆ fIsUseLPMCorrection

G4bool G4PairProductionRelModel::fIsUseLPMCorrection

protected

Definition at line 170 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom(), Initialise(), LPMflag(), SampleSecondaries(), SetLPMflag(), and ~G4PairProductionRelModel().

◆ fLPMEnergy

G4double G4PairProductionRelModel::fLPMEnergy

protected

Definition at line 173 of file G4PairProductionRelModel.hh.

Referenced by SampleSecondaries(), and SetupForMaterial().

◆ fParametrizedXSectionThreshold

G4double G4PairProductionRelModel::fParametrizedXSectionThreshold

protected

Definition at line 175 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom(), and G4PairProductionRelModel().

◆ fParticleChange

G4ParticleChangeForGamma* G4PairProductionRelModel::fParticleChange

protected

Definition at line 182 of file G4PairProductionRelModel.hh.

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

◆ fTheElectron

G4ParticleDefinition* G4PairProductionRelModel::fTheElectron

protected

Definition at line 180 of file G4PairProductionRelModel.hh.

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

◆ fTheGamma

G4ParticleDefinition* G4PairProductionRelModel::fTheGamma

protected

Definition at line 179 of file G4PairProductionRelModel.hh.

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

◆ fThePositron

G4ParticleDefinition* G4PairProductionRelModel::fThePositron

protected

Definition at line 181 of file G4PairProductionRelModel.hh.

Referenced by SampleSecondaries().

◆ gEgLPMActivation

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

staticprotected

Definition at line 164 of file G4PairProductionRelModel.hh.

Referenced by ComputeXSectionPerAtom(), and SampleSecondaries().

◆ gElementData

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

staticprotected

Definition at line 166 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 92 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 95 of file G4PairProductionRelModel.hh.

96 { fIsUseLPMCorrection = val; }

97 inline G4bool LPMflag() const { return fIsUseLPMCorrection; }

G4PairProductionRelModel::LPMflag

G4bool LPMflag() const

Definition G4PairProductionRelModel.hh:97

◆ 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 156 of file G4PairProductionRelModel.hh.

Referenced by SampleSecondaries(), and SetupForMaterial().

◆ gLPMFuncs

G4PairProductionRelModel::LPMFuncs G4PairProductionRelModel::gLPMFuncs

staticprotected

Definition at line 167 of file G4PairProductionRelModel.hh.

Referenced by ~G4PairProductionRelModel().

◆ gMaxZet

const G4int G4PairProductionRelModel::gMaxZet = 120

staticprotected

Definition at line 154 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 85 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 81 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 163 of file G4PairProductionRelModel.hh.

Referenced by ComputeCrossSectionPerAtom().

◆ isFirstInstance

G4bool G4PairProductionRelModel::isFirstInstance {false}

protected

Definition at line 169 of file G4PairProductionRelModel.hh.

169{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

◆ gLPMFuncs

◆ gMaxZet

◆ gWGL

◆ gXGL

◆ gXSecFactor

◆ isFirstInstance