G4GammaConversionToMuons Class Reference

#include <G4GammaConversionToMuons.hh>

Inheritance diagram for G4GammaConversionToMuons:

G4VDiscreteProcess G4VProcess

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)
G4VParticleChangePostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
virtual G4double ComputeCrossSectionPerAtom (G4double GammaEnergy, G4double AtomicZ, G4double AtomicA)
G4double ComputeMeanFreePath (G4double GammaEnergy, G4Material *aMaterial)

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().

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

G4GammaConversionToMuons::~G4GammaConversionToMuons (  ) 

Definition at line 62 of file G4GammaConversionToMuons.cc.

00063 { }


Member Function Documentation

void G4GammaConversionToMuons::BuildPhysicsTable ( const G4ParticleDefinition  )  [virtual]

Reimplemented from G4VProcess.

Definition at line 75 of file G4GammaConversionToMuons.cc.

References PrintInfoDefinition().

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

G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom ( G4double  GammaEnergy,
G4double  AtomicZ,
G4double  AtomicA 
) [virtual]

Definition at line 142 of file G4GammaConversionToMuons.cc.

References G4ParticleDefinition::GetPDGMass(), and G4MuonPlus::MuonPlus().

Referenced by ComputeMeanFreePath(), and GetCrossSectionPerAtom().

00148 { static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
00149   static const G4double Mele=electron_mass_c2;
00150   static const G4double Rc=elm_coupling/Mmuon; // classical particle radius
00151   static const G4double sqrte=sqrt(exp(1.));
00152   static const G4double PowSat=-0.88;
00153 
00154   static G4double CrossSection = 0.0 ;
00155 
00156   if ( A < 1. ) return 0;
00157   if ( Egam < 4*Mmuon ) return 0 ; // below threshold return 0
00158 
00159   static G4double EgamLast=0,Zlast=0,PowThres,Ecor,B,Dn,Zthird,Winfty,WMedAppr,
00160       Wsatur,sigfac;
00161   
00162   if(Zlast==Z && Egam==EgamLast) return CrossSection; // already calculated
00163   EgamLast=Egam;
00164   
00165   if(Zlast!=Z) // new element
00166   { Zlast=Z;
00167     if(Z==1) // special case of Hydrogen
00168     { B=202.4;
00169       Dn=1.49;
00170     }
00171     else
00172     { B=183.;
00173       Dn=1.54*pow(A,0.27);
00174     }
00175     Zthird=pow(Z,-1./3.); // Z**(-1/3)
00176     Winfty=B*Zthird*Mmuon/(Dn*Mele);
00177     WMedAppr=1./(4.*Dn*sqrte*Mmuon);
00178     Wsatur=Winfty/WMedAppr;
00179     sigfac=4.*fine_structure_const*Z*Z*Rc*Rc;
00180     PowThres=1.479+0.00799*Dn;
00181     Ecor=-18.+4347./(B*Zthird);
00182   }
00183   G4double CorFuc=1.+.04*log(1.+Ecor/Egam);
00184   G4double Eg=pow(1.-4.*Mmuon/Egam,PowThres)*pow( pow(Wsatur,PowSat)+
00185               pow(Egam,PowSat),1./PowSat); // threshold and saturation
00186   CrossSection=7./9.*sigfac*log(1.+WMedAppr*CorFuc*Eg);
00187   CrossSection*=CrossSecFactor; // increase the CrossSection by  (by default 1)
00188   return CrossSection;
00189 }

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

Definition at line 104 of file G4GammaConversionToMuons.cc.

References ComputeCrossSectionPerAtom(), DBL_MAX, DBL_MIN, G4Material::GetElementVector(), G4Material::GetNumberOfElements(), and G4Material::GetVecNbOfAtomsPerVolume().

Referenced by GetMeanFreePath().

00108 {
00109   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00110   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
00111 
00112   G4double SIGMA = 0 ;
00113 
00114   for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ )
00115   {
00116     G4double AtomicZ = (*theElementVector)[i]->GetZ();
00117     G4double AtomicA = (*theElementVector)[i]->GetA()/(g/mole);
00118     SIGMA += NbOfAtomsPerVolume[i] *
00119       ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
00120   }
00121   return SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
00122 }

G4double G4GammaConversionToMuons::GetCrossSecFactor (  )  [inline]

Definition at line 86 of file G4GammaConversionToMuons.hh.

00086 { return CrossSecFactor;}

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

Definition at line 126 of file G4GammaConversionToMuons.cc.

References ComputeCrossSectionPerAtom(), G4Element::GetA(), G4DynamicParticle::GetKineticEnergy(), and G4Element::GetZ().

00131 {
00132    G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
00133    G4double AtomicZ = anElement->GetZ();
00134    G4double AtomicA = anElement->GetA()/(g/mole);
00135    G4double crossSection =
00136         ComputeCrossSectionPerAtom(GammaEnergy,AtomicZ,AtomicA);
00137    return crossSection;
00138 }

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().

00089 {
00090    const G4DynamicParticle* aDynamicGamma = aTrack.GetDynamicParticle();
00091    G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
00092    G4Material* aMaterial = aTrack.GetMaterial();
00093 
00094    if (GammaEnergy <  LowestEnergyLimit)
00095      MeanFreePath = DBL_MAX;
00096    else
00097      MeanFreePath = ComputeMeanFreePath(GammaEnergy,aMaterial);
00098 
00099    return MeanFreePath;
00100 }

G4bool G4GammaConversionToMuons::IsApplicable ( const G4ParticleDefinition  )  [virtual]

Reimplemented from G4VProcess.

Definition at line 67 of file G4GammaConversionToMuons.cc.

References G4Gamma::Gamma().

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

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, fStopAndKill, G4cout, G4endl, G4UniformRand, G4Element::GetA(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4Element::GetZ(), G4ParticleChange::Initialize(), G4MuonMinus::MuonMinus(), G4MuonPlus::MuonPlus(), G4INCL::Math::pi, G4VDiscreteProcess::PostStepDoIt(), G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), and G4VParticleChange::SetNumberOfSecondaries().

00208 {
00209   aParticleChange.Initialize(aTrack);
00210   G4Material* aMaterial = aTrack.GetMaterial();
00211 
00212   static const G4double Mmuon=G4MuonPlus::MuonPlus()->GetPDGMass();
00213   static const G4double Mele=electron_mass_c2;
00214   static const G4double sqrte=sqrt(exp(1.));
00215 
00216   // current Gamma energy and direction, return if energy too low
00217   const G4DynamicParticle *aDynamicGamma = aTrack.GetDynamicParticle();
00218   G4double Egam = aDynamicGamma->GetKineticEnergy();
00219   if (Egam < 4*Mmuon) return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
00220   G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
00221 
00222   // select randomly one element constituting the material
00223   const G4Element& anElement = *SelectRandomAtom(aDynamicGamma, aMaterial);
00224   G4double Z = anElement.GetZ();
00225   G4double A = anElement.GetA()/(g/mole);
00226 
00227   static G4double Zlast=0,B,Dn,Zthird,Winfty,A027,C1Num2,C2Term2;
00228   if(Zlast!=Z) // the element has changed
00229   { Zlast=Z;
00230     if(Z==1) // special case of Hydrogen
00231     { B=202.4;
00232       Dn=1.49;
00233     }
00234     else
00235     { B=183.;
00236       Dn=1.54*pow(A,0.27);
00237     }
00238     Zthird=pow(Z,-1./3.); // Z**(-1/3)
00239     Winfty=B*Zthird*Mmuon/(Dn*Mele);
00240     A027=pow(A,0.27);
00241     G4double C1Num=0.35*A027;
00242     C1Num2=C1Num*C1Num;
00243     C2Term2=Mele/(183.*Zthird*Mmuon);
00244   }
00245 
00246   G4double GammaMuonInv=Mmuon/Egam;
00247   G4double sqrtx=sqrt(.25-GammaMuonInv);
00248   G4double xmax=.5+sqrtx;
00249   G4double xmin=.5-sqrtx;
00250 
00251   // generate xPlus according to the differential cross section by rejection
00252   G4double Ds2=(Dn*sqrte-2.);
00253   G4double sBZ=sqrte*B*Zthird/Mele;
00254   G4double LogWmaxInv=1./log(Winfty*(1.+2.*Ds2*GammaMuonInv)
00255                              /(1.+2.*sBZ*Mmuon*GammaMuonInv));
00256   G4double xPlus,xMinus,xPM,result,W;
00257   do
00258   { xPlus=xmin+G4UniformRand()*(xmax-xmin);
00259     xMinus=1.-xPlus;
00260     xPM=xPlus*xMinus;
00261     G4double del=Mmuon*Mmuon/(2.*Egam*xPM);
00262     W=Winfty*(1.+Ds2*del/Mmuon)/(1.+sBZ*del);
00263     if(W<1.) W=1.; // to avoid negative cross section at xmin
00264     G4double xxp=1.-4./3.*xPM; // the main xPlus dependence
00265     result=xxp*log(W)*LogWmaxInv;
00266     if(result>1.) {
00267       G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
00268              << " in dSigxPlusGen, result=" << result << " > 1" << G4endl;
00269     }
00270   }
00271   while (G4UniformRand() > result);
00272 
00273   // now generate the angular variables via the auxilary variables t,psi,rho
00274   G4double t;
00275   G4double psi;
00276   G4double rho;
00277 
00278   G4double thetaPlus,thetaMinus,phiHalf; // final angular variables
00279 
00280   do      // t, psi, rho generation start  (while angle < pi)
00281   {
00282     //generate t by the rejection method
00283     G4double C1=C1Num2* GammaMuonInv/xPM;
00284     G4double f1_max=(1.-xPM) / (1.+C1);
00285     G4double f1; // the probability density
00286     do
00287     { t=G4UniformRand();
00288       f1=(1.-2.*xPM+4.*xPM*t*(1.-t)) / (1.+C1/(t*t));
00289       if(f1<0 || f1> f1_max) // should never happend
00290         {
00291           G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
00292                  << "outside allowed range f1=" << f1 << " is set to zero"
00293                  << G4endl;
00294           f1 = 0.0;
00295         }
00296     }
00297     while ( G4UniformRand()*f1_max > f1);
00298     // generate psi by the rejection method
00299     G4double f2_max=1.-2.*xPM*(1.-4.*t*(1.-t));
00300 
00301     // long version
00302     G4double f2;
00303     do
00304     { psi=2.*pi*G4UniformRand();
00305       f2=1.-2.*xPM+4.*xPM*t*(1.-t)*(1.+cos(2.*psi));
00306       if(f2<0 || f2> f2_max) // should never happend
00307         {
00308           G4cout << "G4GammaConversionToMuons::PostStepDoIt WARNING:"
00309                  << "outside allowed range f2=" << f2 << " is set to zero"
00310                  << G4endl;
00311           f2 = 0.0;
00312         }
00313     }
00314     while ( G4UniformRand()*f2_max > f2);
00315 
00316     // generate rho by direct transformation
00317     G4double C2Term1=GammaMuonInv/(2.*xPM*t);
00318     G4double C2=4./sqrt(xPM)*pow(C2Term1*C2Term1+C2Term2*C2Term2,2.);
00319     G4double rhomax=1.9/A027*(1./t-1.);
00320     G4double beta=log( (C2+pow(rhomax,4.))/C2 );
00321     rho=pow(C2 *( exp(beta*G4UniformRand())-1. ) ,0.25);
00322 
00323     //now get from t and psi the kinematical variables
00324     G4double u=sqrt(1./t-1.);
00325     G4double xiHalf=0.5*rho*cos(psi);
00326     phiHalf=0.5*rho/u*sin(psi);
00327 
00328     thetaPlus =GammaMuonInv*(u+xiHalf)/xPlus;
00329     thetaMinus=GammaMuonInv*(u-xiHalf)/xMinus;
00330 
00331   } while ( std::abs(thetaPlus)>pi || std::abs(thetaMinus) >pi);
00332 
00333   // now construct the vectors
00334   // azimuthal symmetry, take phi0 at random between 0 and 2 pi
00335   G4double phi0=2.*pi*G4UniformRand(); 
00336   G4double EPlus=xPlus*Egam;
00337   G4double EMinus=xMinus*Egam;
00338 
00339   // mu+ mu- directions for gamma in z-direction
00340   G4ThreeVector MuPlusDirection  ( sin(thetaPlus) *cos(phi0+phiHalf),
00341                    sin(thetaPlus)  *sin(phi0+phiHalf), cos(thetaPlus) );
00342   G4ThreeVector MuMinusDirection (-sin(thetaMinus)*cos(phi0-phiHalf),
00343                   -sin(thetaMinus) *sin(phi0-phiHalf), cos(thetaMinus) );
00344   // rotate to actual gamma direction
00345   MuPlusDirection.rotateUz(GammaDirection);
00346   MuMinusDirection.rotateUz(GammaDirection);
00347   aParticleChange.SetNumberOfSecondaries(2);
00348   // create G4DynamicParticle object for the particle1
00349   G4DynamicParticle* aParticle1= new G4DynamicParticle(
00350                            G4MuonPlus::MuonPlus(),MuPlusDirection,EPlus-Mmuon);
00351   aParticleChange.AddSecondary(aParticle1);
00352   // create G4DynamicParticle object for the particle2
00353   G4DynamicParticle* aParticle2= new G4DynamicParticle(
00354                        G4MuonMinus::MuonMinus(),MuMinusDirection,EMinus-Mmuon);
00355   aParticleChange.AddSecondary(aParticle2);
00356   //
00357   // Kill the incident photon
00358   //
00359   aParticleChange.ProposeMomentumDirection( 0., 0., 0. ) ;
00360   aParticleChange.ProposeEnergy( 0. ) ;
00361   aParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00362   //  Reset NbOfInteractionLengthLeft and return aParticleChange
00363   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
00364 }

void G4GammaConversionToMuons::PrintInfoDefinition (  ) 

Definition at line 396 of file G4GammaConversionToMuons.cc.

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

Referenced by BuildPhysicsTable().

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

void G4GammaConversionToMuons::SetCrossSecFactor ( G4double  fac  ) 

Definition at line 193 of file G4GammaConversionToMuons.cc.

References G4cout, and G4endl.

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


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:02 2013 for Geant4 by  doxygen 1.4.7