G4PairProductionRelModel.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 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4PairProductionRelModel
00034 //
00035 // Author:        Andreas Schaelicke
00036 //
00037 // Creation date: 02.04.2009
00038 //
00039 // Modifications:
00040 //
00041 // Class Description:
00042 //
00043 // Main References:
00044 //  J.W.Motz et.al., Rev. Mod. Phys. 41 (1969) 581.
00045 //  S.Klein,  Rev. Mod. Phys. 71 (1999) 1501.
00046 //  T.Stanev et.al., Phys. Rev. D25 (1982) 1291.
00047 //  M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media, 
00048 //                     Wiley, 1972.
00049 //
00050 // -------------------------------------------------------------------
00051 //
00052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00054 
00055 #include "G4PairProductionRelModel.hh"
00056 #include "G4PhysicalConstants.hh"
00057 #include "G4SystemOfUnits.hh"
00058 #include "G4Gamma.hh"
00059 #include "G4Electron.hh"
00060 #include "G4Positron.hh"
00061 
00062 #include "G4ParticleChangeForGamma.hh"
00063 #include "G4LossTableManager.hh"
00064 
00065 
00066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00067 
00068 using namespace std;
00069 
00070 
00071 const G4double G4PairProductionRelModel::facFel = log(184.15);
00072 const G4double G4PairProductionRelModel::facFinel = log(1194.); // 1440.
00073 
00074 const G4double G4PairProductionRelModel::preS1 = 1./(184.15*184.15);
00075 const G4double G4PairProductionRelModel::logTwo = log(2.);
00076 
00077 const G4double G4PairProductionRelModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
00078                                                  0.5917, 0.7628, 0.8983, 0.9801 };
00079 const G4double G4PairProductionRelModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
00080                                                  0.1813, 0.1569, 0.1112, 0.0506 };
00081 const G4double G4PairProductionRelModel::Fel_light[]  = {0., 5.31  , 4.79  , 4.74 ,  4.71};
00082 const G4double G4PairProductionRelModel::Finel_light[] = {0., 6.144 , 5.621 , 5.805 , 5.924};
00083 
00084 
00085 
00086 G4PairProductionRelModel::G4PairProductionRelModel(const G4ParticleDefinition*,
00087                                          const G4String& nam)
00088   : G4VEmModel(nam),
00089     fLPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)*0.5),
00090     fLPMflag(true),
00091     lpmEnergy(0.),
00092     use_completescreening(false)
00093 {
00094   fParticleChange = 0;
00095   theGamma    = G4Gamma::Gamma();
00096   thePositron = G4Positron::Positron();
00097   theElectron = G4Electron::Electron();
00098 
00099   nist = G4NistManager::Instance();  
00100 
00101   currentZ = z13 = z23 = lnZ = Fel = Finel = fCoulomb = phiLPM = gLPM = xiLPM = 0;
00102 
00103 }
00104 
00105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00106 
00107 G4PairProductionRelModel::~G4PairProductionRelModel()
00108 {}
00109 
00110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00111 
00112 void G4PairProductionRelModel::Initialise(const G4ParticleDefinition* p,
00113                                           const G4DataVector& cuts)
00114 {
00115   if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
00116   InitialiseElementSelectors(p, cuts);
00117 }
00118 
00119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00120 
00121 G4double G4PairProductionRelModel::ComputeXSectionPerAtom(G4double totalEnergy, G4double Z)
00122 {
00123   G4double cross = 0.0;
00124 
00125   // number of intervals and integration step 
00126   G4double vcut = electron_mass_c2/totalEnergy ;
00127 
00128   // limits by the screening variable
00129   G4double dmax = DeltaMax();
00130   G4double dmin = min(DeltaMin(totalEnergy),dmax);
00131   G4double vcut1 = 0.5 - 0.5*sqrt(1. - dmin/dmax) ;
00132   vcut = max(vcut, vcut1);
00133 
00134 
00135   G4double vmax = 0.5;
00136   G4int n = 1;  // needs optimisation
00137 
00138   G4double delta = (vmax - vcut)*totalEnergy/G4double(n);
00139 
00140   G4double e0 = vcut*totalEnergy;
00141   G4double xs; 
00142 
00143   // simple integration
00144   for(G4int l=0; l<n; l++,e0 += delta) {
00145     for(G4int i=0; i<8; i++) {
00146 
00147       G4double eg = (e0 + xgi[i]*delta);
00148       if (fLPMflag && totalEnergy>100.*GeV) 
00149         xs = ComputeRelDXSectionPerAtom(eg,totalEnergy,Z);
00150       else
00151         xs = ComputeDXSectionPerAtom(eg,totalEnergy,Z);
00152       cross += wgi[i]*xs;
00153 
00154     }
00155   }
00156 
00157   cross *= delta*2.;
00158 
00159   return cross;
00160 }
00161 
00162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00163 
00164 G4double 
00165 G4PairProductionRelModel::ComputeDXSectionPerAtom(G4double eplusEnergy, 
00166                                                   G4double totalEnergy, 
00167                                                   G4double /*Z*/)
00168 {
00169   // most simple case - complete screening:
00170 
00171   //  dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
00172   //     * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3)  + 1/9 * y*(1-y) ]
00173   // y = E+/k
00174   G4double yp=eplusEnergy/totalEnergy;
00175   G4double ym=1.-yp;
00176 
00177   G4double cross = 0.;
00178   if (use_completescreening) 
00179     cross = (yp*yp + ym*ym + 2./3.*ym*yp)*(Fel - fCoulomb) + 1./9.*yp*ym;
00180   else {
00181     G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
00182     cross = (yp*yp + ym*ym)*(0.25*Phi1(delta) - lnZ/3. - fCoulomb)
00183       + 2./3.*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
00184   }
00185   return cross/totalEnergy;
00186 
00187 }
00188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00189 
00190 G4double 
00191 G4PairProductionRelModel::ComputeRelDXSectionPerAtom(G4double eplusEnergy, 
00192                                                      G4double totalEnergy, 
00193                                                      G4double /*Z*/)
00194 {
00195   // most simple case - complete screening:
00196 
00197   //  dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
00198   //     * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3)  + 1/9 * y*(1-y) ]
00199   // y = E+/k
00200   G4double yp=eplusEnergy/totalEnergy;
00201   G4double ym=1.-yp;
00202 
00203   CalcLPMFunctions(totalEnergy,eplusEnergy); // gamma
00204 
00205   G4double cross = 0.;
00206   if (use_completescreening) 
00207     cross = xiLPM*(2./3.*phiLPM*(yp*yp + ym*ym) + gLPM)*(Fel - fCoulomb);
00208   else {
00209     G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
00210     cross = (1./3.*gLPM + 2./3.*phiLPM)*(yp*yp + ym*ym)
00211                              *(0.25*Phi1(delta) - lnZ/3. - fCoulomb) 
00212            + 2./3.*gLPM*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
00213     cross *= xiLPM;
00214   }
00215   return cross/totalEnergy;
00216 
00217 }
00218 
00219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00220 
00221 void  
00222 G4PairProductionRelModel::CalcLPMFunctions(G4double k, G4double eplusEnergy)
00223 {
00224   // *** calculate lpm variable s & sprime ***
00225   // Klein eqs. (78) & (79)
00226   G4double sprime = sqrt(0.125*k*lpmEnergy/(eplusEnergy*(k-eplusEnergy)));
00227 
00228   G4double s1 = preS1*z23;
00229   G4double logS1 = 2./3.*lnZ-2.*facFel;
00230   G4double logTS1 = logTwo+logS1;
00231 
00232   xiLPM = 2.;
00233 
00234   if (sprime>1) 
00235     xiLPM = 1.;
00236   else if (sprime>sqrt(2.)*s1) {
00237     G4double h  = log(sprime)/logTS1;
00238     xiLPM = 1+h-0.08*(1-h)*(1-sqr(1-h))/logTS1;
00239   }
00240 
00241   G4double s0 = sprime/sqrt(xiLPM); 
00242   //   G4cout<<"k="<<k<<" y="<<eplusEnergy/k<<G4endl;
00243   //   G4cout<<"s0="<<s0<<G4endl;
00244   
00245   // *** calculate supression functions phi and G ***
00246   // Klein eqs. (77)
00247   G4double s2=s0*s0;
00248   G4double s3=s0*s2;
00249   G4double s4=s2*s2;
00250 
00251   if (s0<0.1) {
00252     // high suppression limit
00253     phiLPM = 6.*s0 - 18.84955592153876*s2 + 39.47841760435743*s3 
00254       - 57.69873135166053*s4;
00255     gLPM = 37.69911184307752*s2 - 236.8705056261446*s3 + 807.7822389*s4;
00256   }
00257   else if (s0<1.9516) {
00258     // intermediate suppression
00259     // using eq.77 approxim. valid s0<2.      
00260     phiLPM = 1.-exp(-6.*s0*(1.+(3.-pi)*s0)
00261                 +s3/(0.623+0.795*s0+0.658*s2));
00262     if (s0<0.415827397755) {
00263       // using eq.77 approxim. valid 0.07<s<2
00264       G4double psiLPM = 1-exp(-4*s0-8*s2/(1+3.936*s0+4.97*s2-0.05*s3+7.50*s4));
00265       gLPM = 3*psiLPM-2*phiLPM;
00266     }
00267     else {
00268       // using alternative parametrisiation
00269       G4double pre = -0.16072300849123999 + s0*3.7550300067531581 + s2*-1.7981383069010097 
00270         + s3*0.67282686077812381 + s4*-0.1207722909879257;
00271       gLPM = tanh(pre);
00272     }
00273   }
00274   else {
00275     // low suppression limit valid s>2.
00276     phiLPM = 1. - 0.0119048/s4;
00277     gLPM = 1. - 0.0230655/s4;
00278   }
00279 
00280   // *** make sure suppression is smaller than 1 ***
00281   // *** caused by Migdal approximation in xi    ***
00282   if (xiLPM*phiLPM>1. || s0>0.57)  xiLPM=1./phiLPM;
00283 }
00284 
00285 
00286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00287 
00288 G4double 
00289 G4PairProductionRelModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00290                                                      G4double gammaEnergy, G4double Z,
00291                                                      G4double, G4double, G4double)
00292 {
00293   //  static const G4double gammaEnergyLimit = 1.5*MeV;
00294   G4double crossSection = 0.0 ;
00295   if ( Z < 0.9 ) return crossSection;
00296   if ( gammaEnergy <= 2.0*electron_mass_c2 ) return crossSection;
00297 
00298   SetCurrentElement(Z);
00299   // choose calculator according to parameters and switches
00300   // in the moment only one calculator:
00301   crossSection=ComputeXSectionPerAtom(gammaEnergy,Z);
00302 
00303   G4double xi = Finel/(Fel - fCoulomb); // inelastic contribution
00304   crossSection*=4.*Z*(Z+xi)*fine_structure_const*classic_electr_radius*classic_electr_radius;
00305 
00306   return crossSection;
00307 }
00308 
00309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00310 
00311 void 
00312 G4PairProductionRelModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00313                                             const G4MaterialCutsCouple* couple,
00314                                             const G4DynamicParticle* aDynamicGamma,
00315                                             G4double,
00316                                             G4double)
00317 // The secondaries e+e- energies are sampled using the Bethe - Heitler
00318 // cross sections with Coulomb correction.
00319 // A modified version of the random number techniques of Butcher & Messel
00320 // is used (Nuc Phys 20(1960),15).
00321 //
00322 // GEANT4 internal units.
00323 //
00324 // Note 1 : Effects due to the breakdown of the Born approximation at
00325 //          low energy are ignored.
00326 // Note 2 : The differential cross section implicitly takes account of 
00327 //          pair creation in both nuclear and atomic electron fields.
00328 //          However triplet prodution is not generated.
00329 {
00330   const G4Material* aMaterial = couple->GetMaterial();
00331 
00332   G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
00333   G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
00334 
00335   G4double epsil ;
00336   G4double epsil0 = electron_mass_c2/GammaEnergy ;
00337   if(epsil0 > 1.0) { return; }
00338 
00339   // do it fast if GammaEnergy < 2. MeV
00340   static const G4double Egsmall=2.*MeV;
00341 
00342   // select randomly one element constituing the material
00343   const G4Element* anElement = SelectRandomAtom(aMaterial, theGamma, GammaEnergy);
00344 
00345   if (GammaEnergy < Egsmall) {
00346 
00347     epsil = epsil0 + (0.5-epsil0)*G4UniformRand();
00348 
00349   } else {
00350     // now comes the case with GammaEnergy >= 2. MeV
00351 
00352     // Extract Coulomb factor for this Element
00353     G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
00354     if (GammaEnergy > 50.*MeV) FZ += 8.*(anElement->GetfCoulomb());
00355 
00356     // limits of the screening variable
00357     G4double screenfac = 136.*epsil0/(anElement->GetIonisation()->GetZ3());
00358     G4double screenmax = exp ((42.24 - FZ)/8.368) - 0.952 ;
00359     G4double screenmin = min(4.*screenfac,screenmax);
00360 
00361     // limits of the energy sampling
00362     G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ;
00363     G4double epsilmin = max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin;
00364 
00365     //
00366     // sample the energy rate of the created electron (or positron)
00367     //
00368     //G4double epsil, screenvar, greject ;
00369     G4double  screenvar, greject ;
00370 
00371     G4double F10 = ScreenFunction1(screenmin) - FZ;
00372     G4double F20 = ScreenFunction2(screenmin) - FZ;
00373     G4double NormF1 = max(F10*epsilrange*epsilrange,0.); 
00374     G4double NormF2 = max(1.5*F20,0.);
00375 
00376     do {
00377       if ( NormF1/(NormF1+NormF2) > G4UniformRand() ) {
00378         epsil = 0.5 - epsilrange*pow(G4UniformRand(), 0.333333);
00379         screenvar = screenfac/(epsil*(1-epsil));
00380         if (fLPMflag && GammaEnergy>100.*GeV) {
00381           CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
00382           greject = xiLPM*((gLPM+2.*phiLPM)*Phi1(screenvar) - gLPM*Phi2(screenvar) - phiLPM*FZ)/F10;
00383         }
00384         else {
00385           greject = (ScreenFunction1(screenvar) - FZ)/F10;
00386         }
00387               
00388       } else { 
00389         epsil = epsilmin + epsilrange*G4UniformRand();
00390         screenvar = screenfac/(epsil*(1-epsil));
00391         if (fLPMflag && GammaEnergy>100.*GeV) {
00392           CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
00393           greject = xiLPM*((0.5*gLPM+phiLPM)*Phi1(screenvar) + 0.5*gLPM*Phi2(screenvar) - 0.5*(gLPM+phiLPM)*FZ)/F20;
00394         }
00395         else {
00396           greject = (ScreenFunction2(screenvar) - FZ)/F20;
00397         }
00398       }
00399 
00400     } while( greject < G4UniformRand() );
00401 
00402   }   //  end of epsil sampling
00403    
00404   //
00405   // fixe charges randomly
00406   //
00407 
00408   G4double ElectTotEnergy, PositTotEnergy;
00409   if (G4UniformRand() > 0.5) {
00410 
00411     ElectTotEnergy = (1.-epsil)*GammaEnergy;
00412     PositTotEnergy = epsil*GammaEnergy;
00413      
00414   } else {
00415     
00416     PositTotEnergy = (1.-epsil)*GammaEnergy;
00417     ElectTotEnergy = epsil*GammaEnergy;
00418   }
00419 
00420   //
00421   // scattered electron (positron) angles. ( Z - axis along the parent photon)
00422   //
00423   //  universal distribution suggested by L. Urban 
00424   // (Geant3 manual (1993) Phys211),
00425   //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
00426 
00427   G4double u;
00428   const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00429 
00430   if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
00431   else                            u= - log(G4UniformRand()*G4UniformRand())/a2;
00432 
00433   G4double TetEl = u*electron_mass_c2/ElectTotEnergy;
00434   G4double TetPo = u*electron_mass_c2/PositTotEnergy;
00435   G4double Phi  = twopi * G4UniformRand();
00436   G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
00437   G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
00438    
00439   //
00440   // kinematic of the created pair
00441   //
00442   // the electron and positron are assumed to have a symetric
00443   // angular distribution with respect to the Z axis along the parent photon.
00444 
00445   G4double ElectKineEnergy = max(0.,ElectTotEnergy - electron_mass_c2);
00446 
00447   G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
00448   ElectDirection.rotateUz(GammaDirection);   
00449 
00450   // create G4DynamicParticle object for the particle1  
00451   G4DynamicParticle* aParticle1= new G4DynamicParticle(
00452                      theElectron,ElectDirection,ElectKineEnergy);
00453   
00454   // the e+ is always created (even with Ekine=0) for further annihilation.
00455 
00456   G4double PositKineEnergy = max(0.,PositTotEnergy - electron_mass_c2);
00457 
00458   G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
00459   PositDirection.rotateUz(GammaDirection);   
00460 
00461   // create G4DynamicParticle object for the particle2 
00462   G4DynamicParticle* aParticle2= new G4DynamicParticle(
00463                       thePositron,PositDirection,PositKineEnergy);
00464 
00465   // Fill output vector
00466   fvect->push_back(aParticle1);
00467   fvect->push_back(aParticle2);
00468 
00469   // kill incident photon
00470   fParticleChange->SetProposedKineticEnergy(0.);
00471   fParticleChange->ProposeTrackStatus(fStopAndKill);   
00472 }
00473 
00474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00475 
00476 
00477 void G4PairProductionRelModel::SetupForMaterial(const G4ParticleDefinition*,
00478                                                 const G4Material* mat, G4double)
00479 {
00480   lpmEnergy = mat->GetRadlen()*fLPMconstant;
00481   //  G4cout<<" lpmEnergy="<<lpmEnergy<<G4endl;
00482 }
00483 
00484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

Generated on Mon May 27 17:49:14 2013 for Geant4 by  doxygen 1.4.7