#include <G4GammaConversionToMuons.hh>

Inheritance diagram for G4GammaConversionToMuons:

Collaboration diagram for G4GammaConversionToMuons:

Public Member Functions
	G4GammaConversionToMuons (const G4String &processName="GammaToMuPair", G4ProcessType type=fElectromagnetic)

	~G4GammaConversionToMuons ()

G4bool	IsApplicable (const G4ParticleDefinition &) override

void	BuildPhysicsTable (const G4ParticleDefinition &) override

void	PrintInfoDefinition ()

void	SetCrossSecFactor (G4double fac)

G4double	GetCrossSecFactor ()

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

G4double	GetCrossSectionPerAtom (const G4DynamicParticle aDynamicGamma, G4Element anElement)

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

virtual G4double	ComputeCrossSectionPerAtom (G4double GammaEnergy, G4double AtomicZ, G4double AtomicA)

G4double	ComputeMeanFreePath (G4double GammaEnergy, G4Material *aMaterial)

Public Member Functions inherited from G4VDiscreteProcess
	G4VDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)

	G4VDiscreteProcess (G4VDiscreteProcess &)

virtual	~G4VDiscreteProcess ()

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 &)

Public Member Functions inherited from G4VProcess
	G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)

	G4VProcess (const G4VProcess &right)

virtual	~G4VProcess ()

G4int	operator== (const G4VProcess &right) const

G4int	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 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

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 G4VProcess
void	SubtractNumberOfInteractionLengthLeft (G4double previousStepSize)

void	ClearNumberOfInteractionLengthLeft ()

Protected Attributes inherited from G4VProcess
const G4ProcessManager *	aProcessManager

G4VParticleChange *	pParticleChange

G4ParticleChange	aParticleChange

G4double	theNumberOfInteractionLengthLeft

G4double	currentInteractionLength

G4double	theInitialNumberOfInteractionLength

G4String	theProcessName

G4String	thePhysicsTableFileName

G4ProcessType	theProcessType

G4int	theProcessSubType

G4double	thePILfactor

G4bool	enableAtRestDoIt

G4bool	enableAlongStepDoIt

G4bool	enablePostStepDoIt

G4int	verboseLevel

Detailed Description

Definition at line 61 of file G4GammaConversionToMuons.hh.

Constructor & Destructor Documentation

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

explicit

Definition at line 55 of file G4GammaConversionToMuons.cc.

   : G4VDiscreteProcess (processName, type),
     Mmuon(G4MuonPlus::MuonPlus()->GetPDGMass()),
     Rc(elm_coupling/Mmuon),
     LowestEnergyLimit (4.*Mmuon), // 4*Mmuon
     HighestEnergyLimit(1e21*eV), // ok to 1e21eV=1e12GeV, then LPM suppression
     CrossSecFactor(1.)
 { 
   SetProcessSubType(fGammaConversionToMuMu);
   MeanFreePath = DBL_MAX;
 }

Here is the call graph for this function:

G4GammaConversionToMuons::~G4GammaConversionToMuons ( )

Definition at line 72 of file G4GammaConversionToMuons.cc.

73 {}

Member Function Documentation

void G4GammaConversionToMuons::BuildPhysicsTable ( const G4ParticleDefinition & )

overridevirtual

Reimplemented from G4VProcess.

Definition at line 85 of file G4GammaConversionToMuons.cc.

 {  //here no tables, just calling PrintInfoDefinition
    PrintInfoDefinition();
 }

Here is the call graph for this function:

G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom	(	G4double	GammaEnergy,
		G4double	AtomicZ,
		G4double	AtomicA
	)

virtual

Definition at line 152 of file G4GammaConversionToMuons.cc.

 { 
   if(Egam <= LowestEnergyLimit) return 0.0 ; // below threshold return 0
 
   G4int Z = G4lrint(ZZ);
   G4double CrossSection = 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=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
   //            pow(Egam,PowSat),1./PowSat); // threshold and saturation
   G4double Eg=G4Exp(G4Log(1.-4.*Mmuon/Egam)*PowThres)*
     G4Exp(G4Log( G4Exp(G4Log(Wsatur)*PowSat)+G4Exp(G4Log(Egam)*PowSat))/PowSat);
   CrossSection=7./9.*sigfac*G4Log(1.+WMedAppr*CorFuc*Eg);
   CrossSection*=CrossSecFactor; // increase the CrossSection by  (by default 1)
   return CrossSection;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 114 of file G4GammaConversionToMuons.cc.

 {
   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
 
   G4double SIGMA = 0 ;
 
   for ( size_t i=0 ; i < aMaterial->GetNumberOfElements(); ++i)
   {
     G4double AtomicZ = (*theElementVector)[i]->GetZ();
     G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole);
     SIGMA += NbOfAtomsPerVolume[i] *
       ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
   }
   return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4GammaConversionToMuons::GetCrossSecFactor ( )

inline

Definition at line 87 of file G4GammaConversionToMuons.hh.

87 { return CrossSecFactor;}

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

Definition at line 136 of file G4GammaConversionToMuons.cc.

 {
    G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
    G4double AtomicZ = anElement->GetZ();
    G4double AtomicA = anElement->GetN();
    G4double crossSection =
         ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
    return crossSection;
 }

Here is the call graph for this function:

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

overridevirtual

Implements G4VDiscreteProcess.

Definition at line 93 of file G4GammaConversionToMuons.cc.

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

Here is the call graph for this function:

G4bool G4GammaConversionToMuons::IsApplicable ( const G4ParticleDefinition & particle )

overridevirtual

Reimplemented from G4VProcess.

Definition at line 77 of file G4GammaConversionToMuons.cc.

 {
    return ( &particle == G4Gamma::Gamma() );
 }

Here is the call graph for this function:

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

overridevirtual

Reimplemented from G4VDiscreteProcess.

Definition at line 206 of file G4GammaConversionToMuons.cc.

 {
   aParticleChange.Initialize(aTrack);
   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);
   }
   G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
 
   // select randomly one element constituting the material
   const G4Element* anElement = SelectRandomAtom(aDynamicGamma, aMaterial);
   G4int Z = G4lrint(anElement->GetZ());
   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.35*A027;
   G4double C1Num2=C1Num*C1Num;
   G4double C2Term2=electron_mass_c2/(183.*Zthird*Mmuon);
 
   G4double GammaMuonInv=Mmuon/Egam;
   G4double sqrtx=sqrt(.25-GammaMuonInv);
   G4double xmax=.5+sqrtx;
   G4double xmin=.5-sqrtx;
 
   // generate xPlus according to the differential cross section by rejection
   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,xMinus,xPM,result,W;
   G4int nn = 0;
   const G4int nmax = 1000;
   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);
     if(W<=1. || nn > nmax) { break; } // to avoid negative cross section at xmin
     G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
     result=xxp*G4Log(W)*LogWmaxInv;
     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 thetaPlus,thetaMinus,phiHalf; // final angular variables
   nn = 0;
   do      // t, psi, rho generation start  (while angle < pi)
   {
     //generate t by the rejection method
     G4double C1=C1Num2* GammaMuonInv/xPM;
     G4double f1_max=(1.-xPM) / (1.+C1);
     G4double f1; // the probability density
     do
     { 
       ++nn;
       t=G4UniformRand();
       f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t));
       if(f1<0 || f1> f1_max) // should never happend
     {
       G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
          << "outside allowed range f1=" << f1 << " is set to zero"
          << 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=2.*pi*G4UniformRand();
       f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+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 C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.);
     G4double rhomax=1.9/A027*(1./t-1.);
     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=sqrt(1./t-1.);
     G4double xiHalf=0.5*rho*cos(psi);
     phiHalf=0.5*rho/u*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=2.*pi*G4UniformRand(); 
   G4double EPlus=xPlus*Egam;
   G4double EMinus=xMinus*Egam;
 
   // mu+ mu- directions for gamma in z-direction
   G4ThreeVector MuPlusDirection  ( sin(thetaPlus) *cos(phi0+phiHalf),
                    sin(thetaPlus)  *sin(phi0+phiHalf), cos(thetaPlus) );
   G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
                   -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
   // rotate to actual gamma direction
   MuPlusDirection.rotateUz(GammaDirection);
   MuMinusDirection.rotateUz(GammaDirection);
   aParticleChange.SetNumberOfSecondaries(2);
   // create G4DynamicParticle object for the particle1
   G4DynamicParticle* aParticle1= new G4DynamicParticle(
                            G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
   aParticleChange.AddSecondary(aParticle1);
   // create G4DynamicParticle object for the particle2
   G4DynamicParticle* aParticle2= new G4DynamicParticle(
                        G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
   aParticleChange.AddSecondary(aParticle2);
   //
   // Kill the incident photon
   //
   aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
   aParticleChange.ProposeEnergy( 0. ) ;
   aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
   //  Reset NbOfInteractionLengthLeft and return aParticleChange
   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
 }

Here is the call graph for this function:

void G4GammaConversionToMuons::PrintInfoDefinition ( )

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

Here is the call graph for this function:

Here is the caller graph for this function:

void G4GammaConversionToMuons::SetCrossSecFactor ( G4double fac )

Definition at line 196 of file G4GammaConversionToMuons.cc.

 { 
   CrossSecFactor=fac;
   G4cout << "The cross section for GammaConversionToMuons is artificially "
          << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
 }

Here is the caller graph for this function:

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

source/geant4.10.03.p03/source/processes/electromagnetic/highenergy/include/G4GammaConversionToMuons.hh
source/geant4.10.03.p03/source/processes/electromagnetic/highenergy/src/G4GammaConversionToMuons.cc

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation