G4GammaConversionToMuons Class Reference

#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 G4VParticleChange *	PostStepDoIt (const G4Track &, const G4Step &)
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
virtual G4VParticleChange *	PostStepDoIt (const G4Track &track, const G4Step &stepData)=0
virtual G4VParticleChange *	AlongStepDoIt (const G4Track &track, const G4Step &stepData)=0
virtual G4VParticleChange *	AtRestDoIt (const G4Track &track, const G4Step &stepData)=0
virtual G4double	AlongStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)=0
virtual G4double	AtRestGetPhysicalInteractionLength (const G4Track &track, G4ForceCondition *condition)=0
virtual G4double	PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)=0
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 G4bool	IsApplicable (const G4ParticleDefinition &)
virtual void	BuildPhysicsTable (const G4ParticleDefinition &)
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

Additional Inherited Members
Static Public Member Functions inherited from G4VProcess
static const G4String &	GetProcessTypeName (G4ProcessType)
virtual G4double	GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)=0
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::G4GammaConversionToMuons	(	const G4String &	processName = "GammaToMuPair",
		G4ProcessType	type = fElectromagnetic
	)

Definition at line 48 of file G4GammaConversionToMuons.cc.

                       :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::~G4GammaConversionToMuons

(

)

Definition at line 62 of file G4GammaConversionToMuons.cc.

{ }

Member Function Documentation

BuildPhysicsTable()

void G4GammaConversionToMuons::BuildPhysicsTable ( const G4ParticleDefinition &  )

virtual

Reimplemented from G4VProcess.

Definition at line 75 of file G4GammaConversionToMuons.cc.

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

ComputeCrossSectionPerAtom()

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

virtual

Definition at line 142 of file G4GammaConversionToMuons.cc.

{ 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 G4double CrossSection = 0.0 ;
 
  if ( A < 1. ) return 0;
  if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0
 
  static 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;
}

Referenced by ComputeMeanFreePath(), and GetCrossSectionPerAtom().

ComputeMeanFreePath()

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

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

Referenced by GetMeanFreePath().

GetCrossSecFactor()

G4double G4GammaConversionToMuons::GetCrossSecFactor ( )

inline

Definition at line 86 of file G4GammaConversionToMuons.hh.

{ return CrossSecFactor;}

GetCrossSectionPerAtom()

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

Definition at line 126 of file G4GammaConversionToMuons.cc.

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

GetMeanFreePath()

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

virtual

Implements G4VDiscreteProcess.

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

IsApplicable()

G4bool G4GammaConversionToMuons::IsApplicable ( const G4ParticleDefinition &  particle )

virtual

Reimplemented from G4VProcess.

Definition at line 67 of file G4GammaConversionToMuons.cc.

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

PostStepDoIt()

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

virtual

Reimplemented from G4VDiscreteProcess.

Definition at line 202 of file G4GammaConversionToMuons.cc.

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

PrintInfoDefinition()

void G4GammaConversionToMuons::PrintInfoDefinition

(

)

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

Referenced by BuildPhysicsTable().

SetCrossSecFactor()

void G4GammaConversionToMuons::SetCrossSecFactor

(

G4double

fac

)

Definition at line 193 of file G4GammaConversionToMuons.cc.

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

• G4GammaConversionToMuons.hh
• G4GammaConversionToMuons.cc