#include <G4PenelopeGammaConversionModel.hh>
Inheritance diagram for G4PenelopeGammaConversionModel:
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 | |
G4ParticleChangeForGamma * | fParticleChange |
Definition at line 56 of file G4PenelopeGammaConversionModel.hh.
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 }
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] |
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 83 of file G4PenelopeGammaConversionModel.hh.
Referenced by Initialise(), and SampleSecondaries().