G4GammaConversionToMuons Class Reference

#include <G4GammaConversionToMuons.hh>

Inheritance diagram for G4GammaConversionToMuons:

Public Member Functions
	G4GammaConversionToMuons (const G4String &processName="GammaToMuPair", G4ProcessType type=fElectromagnetic)
	~G4GammaConversionToMuons () override
G4bool	IsApplicable (const G4ParticleDefinition &) override
void	BuildPhysicsTable (const G4ParticleDefinition &) override
void	PrintInfoDefinition ()
void	SetCrossSecFactor (G4double fac)
G4double	GetCrossSecFactor () const
G4double	GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
G4double	GetCrossSectionPerAtom (const G4DynamicParticle *aDynamicGamma, const G4Element *anElement)
G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep) override
G4double	ComputeCrossSectionPerAtom (G4double GammaEnergy, G4int Z)
G4double	ComputeMeanFreePath (G4double GammaEnergy, const G4Material *aMaterial)
G4GammaConversionToMuons &	operator= (const G4GammaConversionToMuons &right)=delete
	G4GammaConversionToMuons (const G4GammaConversionToMuons &)=delete
Public Member Functions inherited from G4VDiscreteProcess
	G4VDiscreteProcess (const G4String &aName, G4ProcessType aType=fNotDefined)
	G4VDiscreteProcess (G4VDiscreteProcess &)
virtual	~G4VDiscreteProcess ()
G4VDiscreteProcess &	operator= (const G4VDiscreteProcess &)=delete
virtual G4double	PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual G4double	AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
virtual G4double	AtRestGetPhysicalInteractionLength (const G4Track &, G4ForceCondition *)
virtual G4VParticleChange *	AtRestDoIt (const G4Track &, const G4Step &)
virtual G4VParticleChange *	AlongStepDoIt (const G4Track &, const G4Step &)
virtual G4double	GetCrossSection (const G4double, const G4MaterialCutsCouple *)
virtual G4double	MinPrimaryEnergy (const G4ParticleDefinition *, const G4Material *)
Public Member Functions inherited from G4VProcess
	G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)
	G4VProcess (const G4VProcess &right)
virtual	~G4VProcess ()
G4VProcess &	operator= (const G4VProcess &)=delete
G4bool	operator== (const G4VProcess &right) const
G4bool	operator!= (const G4VProcess &right) const
G4double	GetCurrentInteractionLength () const
void	SetPILfactor (G4double value)
G4double	GetPILfactor () const
G4double	AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)
G4double	AtRestGPIL (const G4Track &track, G4ForceCondition *condition)
G4double	PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
virtual void	PreparePhysicsTable (const G4ParticleDefinition &)
virtual G4bool	StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
virtual G4bool	RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)
const G4String &	GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
const G4String &	GetProcessName () const
G4ProcessType	GetProcessType () const
void	SetProcessType (G4ProcessType)
G4int	GetProcessSubType () const
void	SetProcessSubType (G4int)
virtual const G4VProcess *	GetCreatorProcess () const
virtual void	StartTracking (G4Track *)
virtual void	EndTracking ()
virtual void	SetProcessManager (const G4ProcessManager *)
virtual const G4ProcessManager *	GetProcessManager ()
virtual void	ResetNumberOfInteractionLengthLeft ()
G4double	GetNumberOfInteractionLengthLeft () const
G4double	GetTotalNumberOfInteractionLengthTraversed () const
G4bool	isAtRestDoItIsEnabled () const
G4bool	isAlongStepDoItIsEnabled () const
G4bool	isPostStepDoItIsEnabled () const
virtual void	DumpInfo () const
virtual void	ProcessDescription (std::ostream &outfile) const
void	SetVerboseLevel (G4int value)
G4int	GetVerboseLevel () const
virtual void	SetMasterProcess (G4VProcess *masterP)
const G4VProcess *	GetMasterProcess () const
virtual void	BuildWorkerPhysicsTable (const G4ParticleDefinition &part)
virtual void	PrepareWorkerPhysicsTable (const G4ParticleDefinition &)

Additional Inherited Members
Static Public Member Functions inherited from G4VProcess
static const G4String &	GetProcessTypeName (G4ProcessType)
Protected Member Functions inherited from G4VDiscreteProcess
Protected Member Functions inherited from G4VProcess
void	SubtractNumberOfInteractionLengthLeft (G4double prevStepSize)
void	ClearNumberOfInteractionLengthLeft ()
Protected Attributes inherited from G4VProcess
const G4ProcessManager *	aProcessManager = nullptr
G4VParticleChange *	pParticleChange = nullptr
G4ParticleChange	aParticleChange
G4double	theNumberOfInteractionLengthLeft = -1.0
G4double	currentInteractionLength = -1.0
G4double	theInitialNumberOfInteractionLength = -1.0
G4String	theProcessName
G4String	thePhysicsTableFileName
G4ProcessType	theProcessType = fNotDefined
G4int	theProcessSubType = -1
G4double	thePILfactor = 1.0
G4int	verboseLevel = 0
G4bool	enableAtRestDoIt = true
G4bool	enableAlongStepDoIt = true
G4bool	enablePostStepDoIt = true

Detailed Description

Definition at line 61 of file G4GammaConversionToMuons.hh.

Constructor & Destructor Documentation

G4GammaConversionToMuons() [1/2]

G4GammaConversionToMuons::G4GammaConversionToMuons ( const G4String & processName = "GammaToMuPair", G4ProcessType type = fElectromagnetic )

explicit

Definition at line 59 of file G4GammaConversionToMuons.cc.

  : G4VDiscreteProcess (processName, type),
    Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
    Rc(CLHEP::elm_coupling / Mmuon),
    LimitEnergy(5. * Mmuon),
    LowestEnergyLimit(2. * Mmuon),
    HighestEnergyLimit(1e12 * CLHEP::GeV),  // ok to 1e12GeV, then LPM suppression
    theGamma(G4Gamma::Gamma()),
    theMuonPlus(G4MuonPlus::MuonPlus()),
    theMuonMinus(G4MuonMinus::MuonMinus())
{
  SetProcessSubType(fGammaConversionToMuMu);
  fManager = G4LossTableManager::Instance();
  fManager->Register(this);
}

~G4GammaConversionToMuons()

G4GammaConversionToMuons::~G4GammaConversionToMuons ( )

override

Definition at line 78 of file G4GammaConversionToMuons.cc.

{
  fManager->DeRegister(this);
}

G4GammaConversionToMuons() [2/2]

G4GammaConversionToMuons::G4GammaConversionToMuons ( const G4GammaConversionToMuons & )

delete

Member Function Documentation

BuildPhysicsTable()

void G4GammaConversionToMuons::BuildPhysicsTable ( const G4ParticleDefinition & p )

overridevirtual

Reimplemented from G4VProcess.

Definition at line 92 of file G4GammaConversionToMuons.cc.

{
  Energy5DLimit = G4EmParameters::Instance()->MaxEnergyFor5DMuPair();
 
  auto table = G4Material::GetMaterialTable();
  std::size_t nelm = 0;
  for (auto const& mat : *table) {
    std::size_t n = mat->GetNumberOfElements();
    nelm = std::max(nelm, n);
  }
  temp.resize(nelm, 0);
 
  if (Energy5DLimit > 0.0 && nullptr != f5Dmodel) {
    f5Dmodel = new G4BetheHeitler5DModel();
    f5Dmodel->SetLeptonPair(theMuonPlus, theMuonMinus);
    const std::size_t numElems = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
    const G4DataVector cuts(numElems);
    f5Dmodel->Initialise(&p, cuts);
  }
  PrintInfoDefinition();
}

Referenced by G4GammaGeneralProcess::BuildPhysicsTable().

ComputeCrossSectionPerAtom()

G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom	(	G4double	GammaEnergy,
		G4int	Z )

Definition at line 170 of file G4GammaConversionToMuons.cc.

{ 
  if(Egam <= LowestEnergyLimit) { return 0.0; }
 
  G4NistManager* nist = G4NistManager::Instance();
 
  G4double PowThres, Ecor, B, Dn, Zthird, Winfty, WMedAppr, Wsatur, sigfac;
 
  if (Z == 1) {  // special case of Hydrogen
    B = 202.4;
    Dn = 1.49;
  }
  else {
    B = 183.;
    Dn = 1.54 * nist->GetA27(Z);
  }
  Zthird = 1. / nist->GetZ13(Z);  // Z**(-1/3)
  Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
  WMedAppr = 1. / (4. * Dn * sqrte * Mmuon);
  Wsatur = Winfty / WMedAppr;
  sigfac = 4. * fine_structure_const * Z * Z * Rc * Rc;
  PowThres = 1.479 + 0.00799 * Dn;
  Ecor = -18. + 4347. / (B * Zthird);
 
  G4double CorFuc = 1. + .04 * G4Log(1. + Ecor / Egam);
  G4double Eg =
    G4Exp(G4Log(1. - 4. * Mmuon / Egam) * PowThres)
    * G4Exp(G4Log(G4Exp(G4Log(Wsatur) * PowSat) + G4Exp(G4Log(Egam) * PowSat)) / PowSat);
 
  G4double CrossSection = 7. / 9. * sigfac * G4Log(1. + WMedAppr * CorFuc * Eg);
  CrossSection *= CrossSecFactor;  // increase the CrossSection by  (by default 1)
  return CrossSection;
}

Referenced by ComputeMeanFreePath(), and GetCrossSectionPerAtom().

ComputeMeanFreePath()

G4double G4GammaConversionToMuons::ComputeMeanFreePath	(	G4double	GammaEnergy,
		const G4Material *	aMaterial )

Definition at line 129 of file G4GammaConversionToMuons.cc.

{
  if(GammaEnergy <= LowestEnergyLimit) { return DBL_MAX; }
  const G4ElementVector* theElementVector = aMaterial->GetElementVector();
  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
 
  G4double SIGMA = 0.0;
  G4double fact  = 1.0;
  G4double e = GammaEnergy;
  // low energy approximation as in Bethe-Heitler model
  if(e < LimitEnergy) {
    G4double y = (e - LowestEnergyLimit)/(LimitEnergy - LowestEnergyLimit);
    fact = y*y;
    e = LimitEnergy;
  } 
 
  for ( std::size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
  {
    SIGMA += NbOfAtomsPerVolume[i] * fact *
      ComputeCrossSectionPerAtom(e, (*theElementVector)[i]->GetZasInt());
  }
  return (SIGMA > 0.0) ? 1./SIGMA : DBL_MAX;
}

Referenced by G4GammaGeneralProcess::BuildPhysicsTable(), and GetMeanFreePath().

GetCrossSecFactor()

G4double G4GammaConversionToMuons::GetCrossSecFactor ( ) const

inline

Definition at line 86 of file G4GammaConversionToMuons.hh.

{ return CrossSecFactor;}

GetCrossSectionPerAtom()

G4double G4GammaConversionToMuons::GetCrossSectionPerAtom	(	const G4DynamicParticle *	aDynamicGamma,
		const G4Element *	anElement )

Definition at line 158 of file G4GammaConversionToMuons.cc.

{
   return ComputeCrossSectionPerAtom(aDynamicGamma->GetKineticEnergy(),
                                     anElement->GetZasInt());
}

GetMeanFreePath()

G4double G4GammaConversionToMuons::GetMeanFreePath ( const G4Track & aTrack, G4double previousStepSize, G4ForceCondition * condition )

overridevirtual

Implements G4VDiscreteProcess.

Definition at line 116 of file G4GammaConversionToMuons.cc.

{
  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
  const G4Material* aMaterial = aTrack.GetMaterial();
  return ComputeMeanFreePath(GammaEnergy, aMaterial);
}

IsApplicable()

G4bool G4GammaConversionToMuons::IsApplicable ( const G4ParticleDefinition & part )

overridevirtual

Reimplemented from G4VProcess.

Definition at line 85 of file G4GammaConversionToMuons.cc.

{
  return (&part == theGamma);
}

operator=()

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

delete

PostStepDoIt()

G4VParticleChange * G4GammaConversionToMuons::PostStepDoIt ( const G4Track & aTrack, const G4Step & aStep )

overridevirtual

Reimplemented from G4VDiscreteProcess.

Definition at line 224 of file G4GammaConversionToMuons.cc.

{
  aParticleChange.Initialize(aTrack);
  const G4Material* aMaterial = aTrack.GetMaterial();
 
  // current Gamma energy and direction, return if energy too low
  const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
  G4double Egam = aDynamicGamma->GetKineticEnergy();
  if (Egam <= LowestEnergyLimit) {
    return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
  }
  //
  // Kill the incident photon
  //
  aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
  aParticleChange.ProposeEnergy( 0. ) ;
  aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
 
  if (Egam <= Energy5DLimit) {
    std::vector<G4DynamicParticle*> fvect;
    f5Dmodel->SampleSecondaries(&fvect, aTrack.GetMaterialCutsCouple(), 
                aTrack.GetDynamicParticle(), 0.0, DBL_MAX);
    for(auto dp : fvect) { aParticleChange.AddSecondary(dp); }
    return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
  }  
 
  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
 
  // select randomly one element constituting the material
  const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
  G4int Z = anElement->GetZasInt();
  G4NistManager* nist = G4NistManager::Instance();
 
  G4double B, Dn;
  G4double A027 = nist->GetA27(Z);
 
  if (Z == 1) {  // special case of Hydrogen
    B = 202.4;
    Dn = 1.49;
  }
  else {
    B = 183.;
    Dn = 1.54 * A027;
  }
  G4double Zthird = 1. / nist->GetZ13(Z);  // Z**(-1/3)
  G4double Winfty = B * Zthird * Mmuon / (Dn * electron_mass_c2);
 
  G4double C1Num = 0.138 * A027;
  G4double C1Num2 = C1Num * C1Num;
  G4double C2Term2 = electron_mass_c2 / (183. * Zthird * Mmuon);
 
  G4double GammaMuonInv = Mmuon / Egam;
 
  // generate xPlus according to the differential cross section by rejection
  G4double xmin = (Egam <= LimitEnergy) ? 0.5 : 0.5 - std::sqrt(0.25 - GammaMuonInv);
  G4double xmax = 1. - xmin;
 
  G4double Ds2 = (Dn * sqrte - 2.);
  G4double sBZ = sqrte * B * Zthird / electron_mass_c2;
  G4double LogWmaxInv =
    1. / G4Log(Winfty * (1. + 2. * Ds2 * GammaMuonInv) / (1. + 2. * sBZ * Mmuon * GammaMuonInv));
  G4double xPlus = 0.5;
  G4double xMinus = 0.5;
  G4double xPM = 0.25;
 
  G4int nn = 0;
  const G4int nmax = 1000;
 
  // sampling for Egam > LimitEnergy
  if (xmin < 0.5) {
    G4double result, W;
    do {
      xPlus = xmin + G4UniformRand() * (xmax - xmin);
      xMinus = 1. - xPlus;
      xPM = xPlus * xMinus;
      G4double del = Mmuon * Mmuon / (2. * Egam * xPM);
      W = Winfty * (1. + Ds2 * del / Mmuon) / (1. + sBZ * del);
      G4double xxp = 1. - 4. / 3. * xPM;  // the main xPlus dependence
      result = (xxp > 0.) ? xxp * G4Log(W) * LogWmaxInv : 0.0;
      if (result > 1.) {
        G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
               << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
      }
      ++nn;
      if(nn >= nmax) { break; }
    }
    // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    while (G4UniformRand() > result);
  }
 
  // now generate the angular variables via the auxilary variables t,psi,rho
  G4double t;
  G4double psi;
  G4double rho;
 
  G4double a3 = (GammaMuonInv / (2. * xPM));
  G4double a33 = a3 * a3;
  G4double f1;
  G4double b1  = 1./(4.*C1Num2);
  G4double b3  = b1*b1*b1;
  G4double a21 = a33 + b1;
  
  G4double f1_max=-(1.-xPM)*(2.*b1+(a21+a33)*G4Log(a33/a21))/(2*b3);  
 
  G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
  nn = 0;
  // t, psi, rho generation start  (while angle < pi)
  do {
    //generate t by the rejection method
    do { 
      ++nn;
      t=G4UniformRand();
      G4double a34=a33/(t*t);
      G4double a22 = a34 + b1;
      if(std::abs(b1)<0.0001*a34) {
        // special case of a34=a22 because of logarithm accuracy
        f1=(1.-2.*xPM+4.*xPM*t*(1.-t))/(12.*a34*a34*a34*a34);
      }
      else {
        f1=-(1.-2.*xPM+4.*xPM*t*(1.-t))*(2.*b1+(a22+a34)*G4Log(a34/a22))/(2*b3);
      }
      if (f1 < 0.0 || f1 > f1_max) {  // should never happend
        G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
        << "outside allowed range f1=" << f1
        << " is set to zero, a34 = "<< a34 << " a22 = "<<a22<<"."
        << G4endl;
        f1 = 0.0;
      }
      if(nn > nmax) { break; }
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko  
    } while ( G4UniformRand()*f1_max > f1);
    // generate psi by the rejection method
    G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
    // long version
    G4double f2;
    do { 
      ++nn;
      psi=twopi*G4UniformRand();
      f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+std::cos(2.*psi));
      if(f2<0 || f2> f2_max) { // should never happend
        G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
               << "outside allowed range f2=" << f2 << " is set to zero" << G4endl;
        f2 = 0.0;
      }
      if(nn >= nmax) { break; }
      // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
    } while ( G4UniformRand()*f2_max > f2);
 
    // generate rho by direct transformation
    G4double C2Term1=GammaMuonInv/(2.*xPM*t);
    G4double C22 = C2Term1*C2Term1+C2Term2*C2Term2;
    G4double C2=4.*C22*C22/std::sqrt(xPM);
    G4double rhomax=(1./t-1.)*1.9/A027;
    G4double beta=G4Log( (C2+rhomax*rhomax*rhomax*rhomax)/C2 );
    rho=G4Exp(G4Log(C2 *( G4Exp(beta*G4UniformRand())-1. ))*0.25);
 
    //now get from t and psi the kinematical variables
    G4double u=std::sqrt(1./t-1.);
    G4double xiHalf=0.5*rho*std::cos(psi);
    phiHalf=0.5*rho/u*std::sin(psi);
 
    thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
    thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
 
    // protection against infinite loop
    if(nn > nmax) {
      if(std::abs(thetaPlus)>pi) { thetaPlus = 0.0; }
      if(std::abs(thetaMinus)>pi) { thetaMinus = 0.0; }
    }
 
    // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
  } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
 
  // now construct the vectors
  // azimuthal symmetry, take phi0 at random between 0 and 2 pi
  G4double phi0=twopi*G4UniformRand(); 
  G4double EPlus=xPlus*Egam;
  G4double EMinus=xMinus*Egam;
 
  // mu+ mu- directions for gamma in z-direction
  G4ThreeVector MuPlusDirection  ( std::sin(thetaPlus) *std::cos(phi0+phiHalf),
                   std::sin(thetaPlus)  *std::sin(phi0+phiHalf), std::cos(thetaPlus) );
  G4ThreeVector MuMinusDirection (-std::sin(thetaMinus)*std::cos(phi0-phiHalf),
                  -std::sin(thetaMinus) *std::sin(phi0-phiHalf), std::cos(thetaMinus) );
  // rotate to actual gamma direction
  MuPlusDirection.rotateUz(GammaDirection);
  MuMinusDirection.rotateUz(GammaDirection);
 
  // create G4DynamicParticle object for the particle1
  auto aParticle1 = new G4DynamicParticle(theMuonPlus, MuPlusDirection, EPlus - Mmuon);
  aParticleChange.AddSecondary(aParticle1);
  // create G4DynamicParticle object for the particle2
  auto aParticle2 = new G4DynamicParticle(theMuonMinus, MuMinusDirection, EMinus - Mmuon);
  aParticleChange.AddSecondary(aParticle2);
  //  Reset NbOfInteractionLengthLeft and return aParticleChange
  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
}

PrintInfoDefinition()

void G4GammaConversionToMuons::PrintInfoDefinition

(

)

Definition at line 462 of file G4GammaConversionToMuons.cc.

{
  G4String comments = "gamma->mu+mu- Bethe Heitler process, SubType= ";
  G4cout << G4endl << GetProcessName() << ":  " << comments << GetProcessSubType() << G4endl;
  G4cout << "        good cross section parametrization from "
         << G4BestUnit(LowestEnergyLimit, "Energy") << " to " << HighestEnergyLimit / GeV
         << " GeV for all Z." << G4endl;
  G4cout << "        cross section factor: " << CrossSecFactor << G4endl;
}

Referenced by BuildPhysicsTable().

SetCrossSecFactor()

void G4GammaConversionToMuons::SetCrossSecFactor

(

G4double

fac

)

Definition at line 211 of file G4GammaConversionToMuons.cc.

{
  if (fac < 0.0) return;
  CrossSecFactor = fac;
  if (verboseLevel > 1) {
    G4cout << "The cross section for GammaConversionToMuons is artificially "
           << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
  }
}

Referenced by G4EmExtraPhysics::ConstructProcess().

---

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

• G4GammaConversionToMuons.hh
• G4GammaConversionToMuons.cc