G4PenelopeGammaConversionModel.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 // Author: Luciano Pandola
00029 //
00030 // History:
00031 // --------
00032 // 13 Jan 2010   L Pandola    First implementation (updated to Penelope08)
00033 // 24 May 2011   L Pandola    Renamed (make v2008 as default Penelope)
00034 //
00035 
00036 #include "G4PenelopeGammaConversionModel.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4MaterialCutsCouple.hh"
00041 #include "G4ProductionCutsTable.hh"
00042 #include "G4DynamicParticle.hh"
00043 #include "G4Element.hh"
00044 #include "G4Gamma.hh"
00045 #include "G4Electron.hh"
00046 #include "G4Positron.hh"
00047 #include "G4PhysicsFreeVector.hh"
00048 
00049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00050 
00051 
00052 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition*,
00053                                                                const G4String& nam)
00054   :G4VEmModel(nam),fParticleChange(0),logAtomicCrossSection(0),
00055    fEffectiveCharge(0),fMaterialInvScreeningRadius(0),
00056    fScreeningFunction(0),isInitialised(false)
00057 {
00058   fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2;
00059   fIntrinsicHighEnergyLimit = 100.0*GeV;
00060   fSmallEnergy = 1.1*MeV;
00061   InitializeScreeningRadii();
00062 
00063   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
00064   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00065   //
00066   verboseLevel= 0;
00067   // Verbosity scale:
00068   // 0 = nothing 
00069   // 1 = warning for energy non-conservation 
00070   // 2 = details of energy budget
00071   // 3 = calculation of cross sections, file openings, sampling of atoms
00072   // 4 = entering in methods
00073 }
00074 
00075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00076 
00077 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel()
00078 {
00079   std::map <G4int,G4PhysicsFreeVector*>::iterator i;
00080   if (logAtomicCrossSection)
00081     {
00082       for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++)
00083         if (i->second) delete i->second;
00084       delete logAtomicCrossSection;
00085     }
00086   if (fEffectiveCharge)
00087     delete fEffectiveCharge;
00088   if (fMaterialInvScreeningRadius)
00089     delete fMaterialInvScreeningRadius;
00090   if (fScreeningFunction)
00091     delete fScreeningFunction;
00092 }
00093 
00094 
00095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00096 
00097 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition*,
00098                                                 const G4DataVector&)
00099 {
00100   if (verboseLevel > 3)
00101     G4cout << "Calling  G4PenelopeGammaConversionModel::Initialise()" << G4endl;
00102 
00103   // logAtomicCrossSection is created only once, since it is  never cleared
00104   if (!logAtomicCrossSection)
00105     logAtomicCrossSection =  new std::map<G4int,G4PhysicsFreeVector*>;
00106 
00107   //delete old material data...
00108   if (fEffectiveCharge)
00109     {
00110       delete fEffectiveCharge;
00111       fEffectiveCharge = 0;
00112     }
00113   if (fMaterialInvScreeningRadius)
00114     {
00115       delete fMaterialInvScreeningRadius;
00116       fMaterialInvScreeningRadius = 0;
00117     }
00118   if (fScreeningFunction)
00119     {
00120       delete fScreeningFunction;
00121       fScreeningFunction = 0;
00122     }      
00123   //and create new ones
00124   fEffectiveCharge = new std::map<const G4Material*,G4double>;
00125   fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
00126   fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
00127 
00128   if (verboseLevel > 0) { 
00129     G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl
00130            << "Energy range: "
00131            << LowEnergyLimit() / MeV << " MeV - "
00132            << HighEnergyLimit() / GeV << " GeV"
00133            << G4endl;
00134   }
00135 
00136   if(isInitialised) return;
00137   fParticleChange = GetParticleChangeForGamma();
00138   isInitialised = true;
00139 }
00140 
00141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00142 
00143 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom(
00144                                                                     const G4ParticleDefinition*,
00145                                                                     G4double energy,
00146                                                                     G4double Z, G4double,
00147                                                                     G4double, G4double)
00148 {
00149   //
00150   // Penelope model v2008.
00151   // Cross section (including triplet production) read from database and managed 
00152   // through the G4CrossSectionHandler utility. Cross section data are from
00153   // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598
00154   //
00155   
00156   if (energy < fIntrinsicLowEnergyLimit)
00157     return 0;
00158 
00159   G4int iZ = (G4int) Z;
00160 
00161  //read data files
00162   if (!logAtomicCrossSection->count(iZ))
00163     ReadDataFile(iZ);
00164   //now it should be ok
00165   if (!logAtomicCrossSection->count(iZ))
00166      {
00167        G4ExceptionDescription ed;
00168        ed << "Unable to retrieve cross section table for Z=" << iZ << G4endl;
00169        G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()",
00170                    "em2018",FatalException,ed);
00171      }
00172 
00173   G4double cs = 0;
00174   G4double logene = std::log(energy);
00175   G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second;
00176 
00177   G4double logXS = theVec->Value(logene);
00178   cs = std::exp(logXS);
00179 
00180   if (verboseLevel > 2)
00181     G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z << 
00182       " = " << cs/barn << " barn" << G4endl;
00183   return cs;
00184 }
00185 
00186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00187 
00188 void 
00189 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00190                                                   const G4MaterialCutsCouple* couple,
00191                                                   const G4DynamicParticle* aDynamicGamma,
00192                                                   G4double,
00193                                                   G4double)
00194 {
00195   //
00196   // Penelope model v2008.
00197   // Final state is sampled according to the Bethe-Heitler model with Coulomb
00198   // corrections, according to the semi-empirical model of
00199   //  J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531.
00200   //
00201   // The model uses the high energy Coulomb correction from 
00202   //  H. Davies et al., Phys. Rev. 93 (1954) 788
00203   // and atomic screening radii tabulated from 
00204   //  J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023
00205   // for Z= 1 to 92. 
00206   //
00207   if (verboseLevel > 3)
00208     G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl;
00209 
00210   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00211 
00212   // Always kill primary
00213   fParticleChange->ProposeTrackStatus(fStopAndKill);
00214   fParticleChange->SetProposedKineticEnergy(0.);
00215 
00216   if (photonEnergy <= fIntrinsicLowEnergyLimit)
00217     {
00218       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00219       return ;
00220     }
00221 
00222   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00223   const G4Material* mat = couple->GetMaterial();
00224 
00225   //check if material data are available
00226   if (!fEffectiveCharge->count(mat))
00227     InitializeScreeningFunctions(mat); 
00228   if (!fEffectiveCharge->count(mat))
00229     {
00230       G4ExceptionDescription ed;      
00231       ed << "Unable to allocate the EffectiveCharge data for " << 
00232         mat->GetName() << G4endl;
00233       G4Exception("G4PenelopeGammaConversion::SampleSecondaries()",
00234                   "em2019",FatalException,ed);            
00235     }
00236 
00237   // eps is the fraction of the photon energy assigned to e- (including rest mass)
00238   G4double eps = 0;
00239   G4double eki = electron_mass_c2/photonEnergy;
00240 
00241   //Do it fast for photon energy < 1.1 MeV (close to threshold)
00242   if (photonEnergy < fSmallEnergy)
00243     eps = eki + (1.0-2.0*eki)*G4UniformRand();
00244   else
00245     {
00246       //Complete calculation
00247       G4double effC = fEffectiveCharge->find(mat)->second;
00248       G4double alz = effC*fine_structure_const;
00249       G4double T = std::sqrt(2.0*eki);
00250       G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
00251          +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
00252          -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
00253         +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
00254      
00255       G4double F0b = fScreeningFunction->find(mat)->second.second;
00256       G4double g0 = F0b + F00;
00257       G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
00258       G4double bmin = 4.0*eki/invRad;
00259       std::pair<G4double,G4double> scree =  GetScreeningFunctions(bmin);
00260       G4double g1 = scree.first;
00261       G4double g2 = scree.second;
00262       G4double g1min = g1+g0;
00263       G4double g2min = g2+g0;
00264       G4double xr = 0.5-eki;
00265       G4double a1 = 2.*g1min*xr*xr/3.;      
00266       G4double p1 = a1/(a1+g2min);
00267 
00268       G4bool loopAgain = false;
00269       //Random sampling of eps
00270       do{
00271         loopAgain = false;
00272         if (G4UniformRand() <= p1)
00273           {
00274             G4double  ru2m1 = 2.0*G4UniformRand()-1.0;
00275             if (ru2m1 < 0)
00276               eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
00277             else
00278               eps = 0.5+xr*std::pow(ru2m1,1./3.);
00279             G4double B = eki/(invRad*eps*(1.0-eps));
00280             scree =  GetScreeningFunctions(B);
00281             g1 = scree.first;
00282             g1 = std::max(g1+g0,0.);
00283             if (G4UniformRand()*g1min > g1) 
00284               loopAgain = true;
00285           }
00286         else
00287           {
00288             eps = eki+2.0*xr*G4UniformRand();
00289             G4double B = eki/(invRad*eps*(1.0-eps));
00290             scree =  GetScreeningFunctions(B);
00291             g2 = scree.second;
00292             g2 = std::max(g2+g0,0.);
00293             if (G4UniformRand()*g2min > g2)
00294               loopAgain = true;
00295           }      
00296       }while(loopAgain);
00297       
00298     }
00299   if (verboseLevel > 4)
00300     G4cout << "Sampled eps = " << eps << G4endl;
00301 
00302   G4double electronTotEnergy = eps*photonEnergy;
00303   G4double positronTotEnergy = (1.0-eps)*photonEnergy;
00304   
00305   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
00306 
00307   //electron kinematics
00308   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 
00309   G4double costheta_el = G4UniformRand()*2.0-1.0;
00310   G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
00311   costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
00312   G4double phi_el  = twopi * G4UniformRand() ;
00313   G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
00314   G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
00315   G4double dirZ_el = costheta_el;
00316 
00317   //positron kinematics
00318   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00319   G4double costheta_po = G4UniformRand()*2.0-1.0;
00320   kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
00321   costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
00322   G4double phi_po  = twopi * G4UniformRand() ;
00323   G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
00324   G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
00325   G4double dirZ_po = costheta_po;
00326 
00327   // Kinematics of the created pair:
00328   // the electron and positron are assumed to have a symetric angular
00329   // distribution with respect to the Z axis along the parent photon
00330   G4double localEnergyDeposit = 0. ;
00331  
00332   if (electronKineEnergy > 0.0)
00333     {
00334       G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
00335       electronDirection.rotateUz(photonDirection);
00336       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
00337                                                            electronDirection,
00338                                                            electronKineEnergy);
00339       fvect->push_back(electron);
00340     }
00341   else
00342     {
00343       localEnergyDeposit += electronKineEnergy;
00344       electronKineEnergy = 0;
00345     }
00346 
00347   //Generate the positron. Real particle in any case, because it will annihilate. If below
00348   //threshold, produce it at rest
00349   if (positronKineEnergy < 0.0)
00350     {
00351       localEnergyDeposit += positronKineEnergy;
00352       positronKineEnergy = 0; //produce it at rest
00353     }
00354   G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po);
00355   positronDirection.rotateUz(photonDirection);
00356   G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(),
00357                                                       positronDirection, positronKineEnergy);
00358   fvect->push_back(positron);
00359 
00360   //Add rest of energy to the local energy deposit
00361   fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
00362   
00363   if (verboseLevel > 1)
00364     {
00365       G4cout << "-----------------------------------------------------------" << G4endl;
00366       G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl;
00367       G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
00368       G4cout << "-----------------------------------------------------------" << G4endl;
00369       if (electronKineEnergy)
00370         G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" 
00371                << G4endl;
00372       if (positronKineEnergy)
00373         G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl;
00374       G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl;
00375       if (localEnergyDeposit)
00376         G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
00377       G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+
00378                                           localEnergyDeposit+2.0*electron_mass_c2)/keV <<
00379         " keV" << G4endl;
00380       G4cout << "-----------------------------------------------------------" << G4endl;
00381     }
00382  if (verboseLevel > 0)
00383     {
00384       G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
00385                                       localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
00386       if (energyDiff > 0.05*keV)
00387         G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " 
00388                << (electronKineEnergy+positronKineEnergy+
00389                    localEnergyDeposit+2.0*electron_mass_c2)/keV 
00390                << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl;
00391     } 
00392 }
00393 
00394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00395 
00396 void G4PenelopeGammaConversionModel::ReadDataFile(const G4int Z)
00397 {
00398  if (verboseLevel > 2)
00399     {
00400       G4cout << "G4PenelopeGammaConversionModel::ReadDataFile()" << G4endl;
00401       G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl;
00402     }
00403  
00404   char* path = getenv("G4LEDATA");
00405   if (!path)
00406     {
00407       G4String excep = 
00408         "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!";
00409       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
00410                   "em0006",FatalException,excep);
00411       return;
00412     }
00413  
00414   /*
00415     Read the cross section file
00416   */
00417   std::ostringstream ost;
00418   if (Z>9)
00419     ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08";
00420   else
00421     ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08";
00422   std::ifstream file(ost.str().c_str());
00423   if (!file.is_open())
00424     {
00425       G4String excep = "G4PenelopeGammaConversionModel - data file " + 
00426         G4String(ost.str()) + " not found!";
00427       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
00428                   "em0003",FatalException,excep);
00429     }
00430 
00431   //I have to know in advance how many points are in the data list
00432   //to initialize the G4PhysicsFreeVector()
00433   size_t ndata=0;
00434   G4String line;
00435   while( getline(file, line) )
00436     ndata++;
00437   ndata -= 1; //remove one header line
00438   //G4cout << "Found: " << ndata << " lines" << G4endl;
00439 
00440   file.clear();
00441   file.close();
00442   file.open(ost.str().c_str());
00443   G4int readZ =0;
00444   file >> readZ; 
00445 
00446   if (verboseLevel > 3)
00447     G4cout << "Element Z=" << Z << G4endl;
00448 
00449   //check the right file is opened.
00450   if (readZ != Z)
00451     {
00452       G4ExceptionDescription ed;
00453       ed << "Corrupted data file for Z=" << Z << G4endl;      
00454       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
00455                   "em0005",FatalException,ed);
00456     }
00457 
00458   G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata);
00459   G4double ene=0,xs=0;
00460   for (size_t i=0;i<ndata;i++)
00461     {
00462       file >> ene >> xs;
00463       //dimensional quantities
00464       ene *= eV;
00465       xs *= barn;
00466       if (xs < 1e-40*cm2) //protection against log(0)
00467         xs = 1e-40*cm2;
00468       theVec->PutValue(i,std::log(ene),std::log(xs));      
00469     }
00470   file.close();
00471 
00472   if (!logAtomicCrossSection)
00473     {
00474       G4ExceptionDescription ed;
00475       ed << "Problem with allocation of logAtomicCrossSection data table " << G4endl;
00476       G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()",
00477                   "em2020",FatalException,ed);
00478       delete theVec;
00479       return;
00480     }
00481   logAtomicCrossSection->insert(std::make_pair(Z,theVec));
00482 
00483   return;
00484 
00485 }
00486 
00487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00488 
00489 void G4PenelopeGammaConversionModel::InitializeScreeningRadii()
00490 {
00491   G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01,
00492                        6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01,
00493                        4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01,
00494                        4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01,
00495                        3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01,
00496                        3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01,
00497                        3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01,
00498                        3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01,
00499                        2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01,
00500                        2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01,
00501                        2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01,
00502                        2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01,
00503                        2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01,
00504                        2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01,
00505                        2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01,
00506                        2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01,
00507                        2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01};
00508 
00509   //copy temporary vector in class data member
00510   for (G4int i=0;i<99;i++)
00511     fAtomicScreeningRadius[i] = temp[i];
00512 }
00513 
00514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00515 
00516 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(const G4Material* material)
00517 {
00518   // This is subroutine GPPa0 of Penelope
00519   //
00520   // 1) calculate the effective Z for the purpose
00521   //
00522   G4double zeff = 0;
00523   G4int intZ = 0;
00524   G4int nElements = material->GetNumberOfElements();
00525   const G4ElementVector* elementVector = material->GetElementVector();
00526 
00527   //avoid calculations if only one building element!
00528   if (nElements == 1)
00529     {
00530       zeff = (*elementVector)[0]->GetZ();
00531       intZ = (G4int) zeff;
00532     }
00533   else // many elements...let's do the calculation
00534     {
00535       const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume();
00536       
00537       G4double atot = 0;
00538       for (G4int i=0;i<nElements;i++)
00539         {
00540           G4double Zelement = (*elementVector)[i]->GetZ();
00541           G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu();
00542           atot += Aelement*fractionVector[i];
00543           zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei
00544         }
00545       atot /= material->GetTotNbOfAtomsPerVolume();
00546       zeff /= (material->GetTotNbOfAtomsPerVolume()*atot);
00547       
00548       intZ = (G4int) (zeff+0.25);
00549       if (intZ <= 0)
00550         intZ = 1;
00551       if (intZ > 99)
00552         intZ = 99;
00553     }
00554 
00555   if (fEffectiveCharge)
00556     fEffectiveCharge->insert(std::make_pair(material,zeff));
00557 
00558   //
00559   // 2) Calculate Coulomb Correction
00560   //
00561   G4double alz = fine_structure_const*zeff;
00562   G4double alzSquared = alz*alz;
00563   G4double fc =  alzSquared*(0.202059-alzSquared*
00564                              (0.03693-alzSquared*
00565                               (0.00835-alzSquared*(0.00201-alzSquared*
00566                                                    (0.00049-alzSquared*
00567                                                     (0.00012-alzSquared*0.00003)))))
00568                              +1.0/(alzSquared+1.0));
00569   //
00570   // 3) Screening functions and low-energy corrections
00571   //
00572   G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1];
00573   if (fMaterialInvScreeningRadius)
00574     fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius));
00575 
00576   std::pair<G4double,G4double> myPair(0,0);
00577   G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]);
00578   G4double f0b = f0a - 4.0*fc;
00579   myPair.first = f0a;
00580   myPair.second = f0b;
00581 
00582   if (fScreeningFunction)
00583     fScreeningFunction->insert(std::make_pair(material,myPair));
00584 
00585   if (verboseLevel > 2)
00586     {
00587       G4cout << "Average Z for material " << material->GetName() << " = " << 
00588         zeff << G4endl;
00589       G4cout << "Effective radius for material " << material->GetName() << " = " << 
00590         fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " << 
00591         matRadius << G4endl;
00592       G4cout << "Screening parameters F0 for material " << material->GetName() << " = " << 
00593         f0a << "," << f0b << G4endl;
00594     }
00595   return;
00596 }
00597 
00598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00599 
00600 std::pair<G4double,G4double> 
00601 G4PenelopeGammaConversionModel::GetScreeningFunctions(G4double B)
00602 {
00603   // This is subroutine SCHIFF of Penelope
00604   //
00605   // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross 
00606   // section for pair production
00607   //
00608   std::pair<G4double,G4double> result(0.,0.);
00609   G4double BSquared = B*B;
00610   G4double f1 = 2.0-2.0*std::log(1.0+BSquared);
00611   G4double f2 = f1 - 6.66666666e-1; // (-2/3)
00612   if (B < 1.0e-10)
00613     f1 = f1-twopi*B;
00614   else
00615     {
00616       G4double a0 = 4.0*B*std::atan(1./B);
00617       f1 = f1 - a0;
00618       f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared));
00619     }
00620   G4double g1 = 0.5*(3.0*f1-f2);
00621   G4double g2 = 0.25*(3.0*f1+f2);
00622 
00623   result.first = g1;
00624   result.second = g2;
00625 
00626   return result;
00627 }

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