G4PenelopeGammaConversionModel Class Reference

#include <G4PenelopeGammaConversionModel.hh>

Inheritance diagram for G4PenelopeGammaConversionModel:

G4VEmModel

Public Member Functions

 G4PenelopeGammaConversionModel (const G4ParticleDefinition *p=0, const G4String &processName="PenConversion")
virtual ~G4PenelopeGammaConversionModel ()
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void SetVerbosityLevel (G4int lev)
G4int GetVerbosityLevel ()

Protected Attributes

G4ParticleChangeForGammafParticleChange

Detailed Description

Definition at line 56 of file G4PenelopeGammaConversionModel.hh.


Constructor & Destructor Documentation

G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel ( const G4ParticleDefinition p = 0,
const G4String processName = "PenConversion" 
)

Definition at line 52 of file G4PenelopeGammaConversionModel.cc.

References G4VEmModel::SetHighEnergyLimit().

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 }

G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel (  )  [virtual]

Definition at line 77 of file G4PenelopeGammaConversionModel.cc.

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 }


Member Function Documentation

G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0,
G4double  cut = 0,
G4double  emax = DBL_MAX 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 143 of file G4PenelopeGammaConversionModel.cc.

References FatalException, G4cout, G4endl, G4Exception(), and G4PhysicsVector::Value().

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 }

G4int G4PenelopeGammaConversionModel::GetVerbosityLevel (  )  [inline]

Definition at line 83 of file G4PenelopeGammaConversionModel.hh.

00083 {return verboseLevel;};

void G4PenelopeGammaConversionModel::Initialise ( const G4ParticleDefinition ,
const G4DataVector  
) [virtual]

Implements G4VEmModel.

Definition at line 97 of file G4PenelopeGammaConversionModel.cc.

References fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), and G4VEmModel::LowEnergyLimit().

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 }

void G4PenelopeGammaConversionModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  ,
const G4MaterialCutsCouple ,
const G4DynamicParticle ,
G4double  tmin,
G4double  maxEnergy 
) [virtual]

Implements G4VEmModel.

Definition at line 189 of file G4PenelopeGammaConversionModel.cc.

References G4Electron::Electron(), F00, FatalException, fParticleChange, fStopAndKill, G4cout, G4endl, G4Exception(), G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Positron::Positron(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4VParticleChange::ProposeTrackStatus(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().

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 }

void G4PenelopeGammaConversionModel::SetVerbosityLevel ( G4int  lev  )  [inline]

Definition at line 82 of file G4PenelopeGammaConversionModel.hh.

00082 {verboseLevel = lev;};


Field Documentation

G4ParticleChangeForGamma* G4PenelopeGammaConversionModel::fParticleChange [protected]

Definition at line 83 of file G4PenelopeGammaConversionModel.hh.

Referenced by Initialise(), and SampleSecondaries().


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