#include <G4GammaConversionToMuons.hh>

Inheritance diagram for G4GammaConversionToMuons:

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

	~G4GammaConversionToMuons ()

G4bool	IsApplicable (const G4ParticleDefinition &)

void	BuildPhysicsTable (const G4ParticleDefinition &)

void	PrintInfoDefinition ()

void	SetCrossSecFactor (G4double fac)

G4double	GetCrossSecFactor ()

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

G4double	GetCrossSectionPerAtom (const G4DynamicParticle aDynamicGamma, G4Element anElement)

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

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

Definition at line 48 of file G4GammaConversionToMuons.cc.

References DBL_MAX, and G4VProcess::SetProcessSubType().

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

G4GammaConversionToMuons::~G4GammaConversionToMuons ( )

Definition at line 62 of file G4GammaConversionToMuons.cc.

63 { }

Member Function Documentation

void G4GammaConversionToMuons::BuildPhysicsTable ( const G4ParticleDefinition & )

virtual

Reimplemented from G4VProcess.

Definition at line 75 of file G4GammaConversionToMuons.cc.

References PrintInfoDefinition().

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

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

virtual

Definition at line 142 of file G4GammaConversionToMuons.cc.

References python.hepunit::electron_mass_c2, python.hepunit::elm_coupling, python.hepunit::fine_structure_const, G4ThreadLocal, G4ParticleDefinition::GetPDGMass(), and G4MuonPlus::MuonPlus().

Referenced by ComputeMeanFreePath(), and GetCrossSectionPerAtom().

 { static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
   static const G4double Mele=electron_mass_c2;
   static const G4double Rc=elm_coupling/Mmuon; // classical particle radius
   static const G4double sqrte=sqrt(exp(1.));
   static const G4double PowSat=-0.88;
 
   static G4ThreadLocal G4double CrossSection = 0.0 ;
 
   if ( A < 1. ) return 0;
   if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0
 
   static G4ThreadLocal G4double EgamLast=0,Zlast=0,PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
       Wsatur,sigfac;
   
   if(Zlast==Z && Egam==EgamLast) return CrossSection; // already calculated
   EgamLast=Egam;
   
   if(Zlast!=Z) // new element
   { Zlast=Z;
     if(Z==1) // special case of Hydrogen
     { B=202.4;
       Dn=1.49;
     }
     else
     { B=183.;
       Dn=1.54*pow(A,0.27);
     }
     Zthird=pow(Z,-1./3.); // Z**(-1/3)
     Winfty=B*Zthird*Mmuon/(Dn*Mele);
     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*log(1.+Ecor/Egam);
   G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
               pow(Egam,PowSat),1./PowSat); // threshold and saturation
   CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg);
   CrossSection*=CrossSecFactor; // increase the CrossSection by  (by default 1)
   return CrossSection;
 }

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

Definition at line 104 of file G4GammaConversionToMuons.cc.

References ComputeCrossSectionPerAtom(), DBL_MAX, DBL_MIN, g(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4Material::GetVecNbOfAtomsPerVolume(), and python.hepunit::mole.

Referenced by GetMeanFreePath().

 {
   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;
 }

G4double G4GammaConversionToMuons::GetCrossSecFactor ( )

inline

Definition at line 86 of file G4GammaConversionToMuons.hh.

86 { return CrossSecFactor;}

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

Definition at line 126 of file G4GammaConversionToMuons.cc.

References ComputeCrossSectionPerAtom(), g(), G4Element::GetA(), G4DynamicParticle::GetKineticEnergy(), G4Element::GetZ(), and python.hepunit::mole.

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

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

virtual

Implements G4VDiscreteProcess.

Definition at line 83 of file G4GammaConversionToMuons.cc.

References ComputeMeanFreePath(), DBL_MAX, G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), and G4Track::GetMaterial().

 {
    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;
 }

G4bool G4GammaConversionToMuons::IsApplicable ( const G4ParticleDefinition & particle )

virtual

Reimplemented from G4VProcess.

Definition at line 67 of file G4GammaConversionToMuons.cc.

References G4Gamma::Gamma().

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

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

virtual

Reimplemented from G4VDiscreteProcess.

Definition at line 202 of file G4GammaConversionToMuons.cc.

References G4ParticleChange::AddSecondary(), G4VProcess::aParticleChange, C1, python.hepunit::electron_mass_c2, fStopAndKill, g(), G4cout, G4endl, G4ThreadLocal, G4UniformRand, G4Element::GetA(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4Element::GetZ(), G4ParticleChange::Initialize(), python.hepunit::mole, G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), python.hepunit::pi, G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), CLHEP::Hep3Vector::rotateUz(), and G4VParticleChange::SetNumberOfSecondaries().

 {
   aParticleChange.Initialize(aTrack);
   G4Material* aMaterial = aTrack.GetMaterial();
 
   static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
   static const G4double Mele=electron_mass_c2;
   static const G4double sqrte=sqrt(exp(1.));
 
   // current Gamma energy and direction, return if energy too low
   const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle();
   G4double Egam = aDynamicGamma->GetKineticEnergy();
   if (Egam < 4*Mmuon) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
   G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
 
   // select randomly one element constituting the material
   const G4Element& anElement = *SelectRandomAtom(aDynamicGamma, aMaterial);
   G4double Z = anElement.GetZ();
   G4double A = anElement.GetA()/(g/mole);
 
   static G4ThreadLocal G4double Zlast=0,B,Dn,Zthird,Winfty,A027,C1Num2,C2Term2;
   if(Zlast!=Z) // the element has changed
   { Zlast=Z;
     if(Z==1) // special case of Hydrogen
     { B=202.4;
       Dn=1.49;
     }
     else
     { B=183.;
       Dn=1.54*pow(A,0.27);
     }
     Zthird=pow(Z,-1./3.); // Z**(-1/3)
     Winfty=B*Zthird*Mmuon/(Dn*Mele);
     A027=pow(A,0.27);
     G4double C1Num=0.35*A027;
     C1Num2=C1Num*C1Num;
     C2Term2=Mele/(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/Mele;
   G4double LogWmaxInv=1./log(Winfty*(1.+2.*Ds2*GammaMuonInv)
                              /(1.+2.*sBZ*Mmuon*GammaMuonInv));
   G4double xPlus,xMinus,xPM,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);
     if(W<1.) W=1.; // to avoid negative cross section at xmin
     G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
     result=xxp*log(W)*LogWmaxInv;
     if(result>1.) {
       G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
          << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
     }
   }
   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
 
   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
     { 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;
     }
     }
     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
     { 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;
     }
     }
     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=log( (C2+pow(rhomax,4.))/C2 );
     rho=pow(C2 *( exp(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;
 
   } 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 );
 }

void G4GammaConversionToMuons::PrintInfoDefinition ( )

Definition at line 396 of file G4GammaConversionToMuons.cc.

References G4BestUnit, G4cout, G4endl, G4VProcess::GetProcessName(), G4VProcess::GetProcessSubType(), and python.hepunit::GeV.

Referenced by BuildPhysicsTable().

 {
   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;
 }

void G4GammaConversionToMuons::SetCrossSecFactor ( G4double fac )

Definition at line 193 of file G4GammaConversionToMuons.cc.

References G4cout, and G4endl.

Referenced by PhysicsList::SetGammaToMuPairFac().

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

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

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation