G4GammaConversionToMuons.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4GammaConversionToMuons.cc 66996 2013-01-29 14:50:52Z gcosmo $
00028 //
00029 //         ------------ G4GammaConversionToMuons physics process ------
00030 //         by H.Burkhardt, S. Kelner and R. Kokoulin, April 2002
00031 //
00032 //
00033 // 07-08-02: missprint in OR condition in DoIt : f1<0 || f1>f1_max ..etc ...
00034 // 25-10-04: migrade to new interfaces of ParticleChange (vi)
00035 // ---------------------------------------------------------------------------
00036 
00037 #include "G4GammaConversionToMuons.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4UnitsTable.hh"
00041 #include "G4MuonPlus.hh"
00042 #include "G4MuonMinus.hh"
00043 
00044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00045 
00046 using namespace std;
00047 
00048 G4GammaConversionToMuons::G4GammaConversionToMuons(const G4String& processName,
00049     G4ProcessType type):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 }
00057 
00058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00059 
00060 // destructor
00061 
00062 G4GammaConversionToMuons::~G4GammaConversionToMuons() // (empty) destructor
00063 { }
00064 
00065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00066 
00067 G4bool G4GammaConversionToMuons::IsApplicable(
00068                                         const G4ParticleDefinition& particle)
00069 {
00070    return ( &particle == G4Gamma::Gamma() );
00071 }
00072 
00073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00074 
00075 void G4GammaConversionToMuons::BuildPhysicsTable(const G4ParticleDefinition&)
00076 // Build cross section and mean free path tables
00077 {  //here no tables, just calling PrintInfoDefinition
00078    PrintInfoDefinition();
00079 }
00080 
00081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00082 
00083 G4double G4GammaConversionToMuons::GetMeanFreePath(const G4Track& aTrack,
00084                                               G4double, G4ForceCondition*)
00085 
00086 // returns the photon mean free path in GEANT4 internal units
00087 // (MeanFreePath is a private member of the class)
00088 
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 }
00101 
00102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00103 
00104 G4double G4GammaConversionToMuons::ComputeMeanFreePath(G4double GammaEnergy,
00105                                                        G4Material* aMaterial)
00106 
00107 // computes and returns the photon mean free path in GEANT4 internal units
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 }
00123 
00124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00125 
00126 G4double G4GammaConversionToMuons::GetCrossSectionPerAtom(
00127                                    const G4DynamicParticle* aDynamicGamma,
00128                                          G4Element*         anElement)
00129 
00130 // gives the total cross section per atom in GEANT4 internal units
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 }
00139 
00140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00141 
00142 G4double G4GammaConversionToMuons::ComputeCrossSectionPerAtom(
00143                          G4double Egam, G4double Z, G4double A)
00144                          
00145 // Calculates the microscopic cross section in GEANT4 internal units.
00146 // Total cross section parametrisation from H.Burkhardt
00147 // It gives a good description at any energy (from 0 to 10**21 eV)
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 }
00190 
00191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00192 
00193 void G4GammaConversionToMuons::SetCrossSecFactor(G4double fac)
00194 // Set the factor to artificially increase the cross section
00195 { CrossSecFactor=fac;
00196   G4cout << "The cross section for GammaConversionToMuons is artificially "
00197          << "increased by the CrossSecFactor=" << CrossSecFactor << G4endl;
00198 }
00199 
00200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00201 
00202 G4VParticleChange* G4GammaConversionToMuons::PostStepDoIt(
00203                                                         const G4Track& aTrack,
00204                                                         const G4Step&  aStep)
00205 //
00206 // generation of gamma->mu+mu-
00207 //
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 }
00365 
00366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00367 
00368 G4Element* G4GammaConversionToMuons::SelectRandomAtom(
00369                                         const G4DynamicParticle* aDynamicGamma,
00370                                               G4Material* aMaterial)
00371 {
00372   // select randomly 1 element within the material, invoked by PostStepDoIt
00373 
00374   const G4int NumberOfElements            = aMaterial->GetNumberOfElements();
00375   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00376   if (NumberOfElements == 1) return (*theElementVector)[0];
00377 
00378   const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
00379 
00380   G4double PartialSumSigma = 0. ;
00381   G4double rval = G4UniformRand()/MeanFreePath;
00382 
00383 
00384   for ( G4int i=0 ; i < NumberOfElements ; i++ )
00385       { PartialSumSigma += NbOfAtomsPerVolume[i] *
00386                  GetCrossSectionPerAtom(aDynamicGamma, (*theElementVector)[i]);
00387         if (rval <= PartialSumSigma) return ((*theElementVector)[i]);
00388       }
00389   G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName()
00390        << "' has no elements, NULL pointer returned." << G4endl;
00391   return NULL;
00392 }
00393 
00394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
00395 
00396 void G4GammaConversionToMuons::PrintInfoDefinition()
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 }
00405 
00406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

Generated on Mon May 27 17:48:19 2013 for Geant4 by  doxygen 1.4.7