#include <G4LivermoreGammaConversionModel.hh>
Inheritance diagram for G4LivermoreGammaConversionModel:
Public Member Functions | |
G4LivermoreGammaConversionModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreConversion") | |
virtual | ~G4LivermoreGammaConversionModel () |
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) |
Definition at line 40 of file G4LivermoreGammaConversionModel.hh.
G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LivermoreConversion" | |||
) |
Definition at line 40 of file G4LivermoreGammaConversionModel.cc.
References G4cout, and G4endl.
00042 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),maxZ(99) 00043 { 00044 fParticleChange = 0; 00045 00046 lowEnergyLimit = 2.0*electron_mass_c2; 00047 data.resize(maxZ+1,0); 00048 00049 verboseLevel= 0; 00050 // Verbosity scale for debugging purposes: 00051 // 0 = nothing 00052 // 1 = calculation of cross sections, file openings... 00053 // 2 = entering in methods 00054 00055 if(verboseLevel > 0) 00056 { 00057 G4cout << "G4LivermoreGammaConversionModel is constructed " << G4endl; 00058 } 00059 }
G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel | ( | ) | [virtual] |
Definition at line 63 of file G4LivermoreGammaConversionModel.cc.
00064 { 00065 for(G4int i=0; i<=maxZ; ++i) { delete data[i]; } 00066 }
G4double G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0 , |
|||
G4double | cut = 0 , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 182 of file G4LivermoreGammaConversionModel.cc.
References G4cout, G4endl, G4PhysicsVector::GetVectorLength(), CLHEP::detail::n, and G4PhysicsVector::Value().
00186 { 00187 if (verboseLevel > 1) 00188 { 00189 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel" 00190 << G4endl; 00191 } 00192 00193 if (GammaEnergy < lowEnergyLimit) { return 0.0; } 00194 00195 G4double xs = 0.0; 00196 00197 G4int intZ=G4int(Z); 00198 00199 if(intZ < 1 || intZ > maxZ) { return xs; } 00200 00201 G4LPhysicsFreeVector* pv = data[intZ]; 00202 00203 // element was not initialised 00204 if(!pv) 00205 { 00206 char* path = getenv("G4LEDATA"); 00207 ReadData(intZ, path); 00208 pv = data[intZ]; 00209 if(!pv) { return xs; } 00210 } 00211 // x-section is taken from the table 00212 xs = pv->Value(GammaEnergy); 00213 00214 if(verboseLevel > 0) 00215 { 00216 G4int n = pv->GetVectorLength() - 1; 00217 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << GammaEnergy/MeV << G4endl; 00218 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; 00219 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl; 00220 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl; 00221 G4cout << "*********************************************************" << G4endl; 00222 } 00223 00224 return xs; 00225 00226 }
void G4LivermoreGammaConversionModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 71 of file G4LivermoreGammaConversionModel.cc.
References G4cout, G4endl, G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4VEmModel::GetParticleChangeForGamma(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), and G4VEmModel::LowEnergyLimit().
00073 { 00074 if (verboseLevel > 1) 00075 { 00076 G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel." << G4endl 00077 << "Energy range: " 00078 << LowEnergyLimit() / MeV << " MeV - " 00079 << HighEnergyLimit() / GeV << " GeV" 00080 << G4endl; 00081 } 00082 00083 // Initialise element selector 00084 00085 InitialiseElementSelectors(particle, cuts); 00086 00087 // Access to elements 00088 00089 char* path = getenv("G4LEDATA"); 00090 00091 G4ProductionCutsTable* theCoupleTable = 00092 G4ProductionCutsTable::GetProductionCutsTable(); 00093 G4int numOfCouples = theCoupleTable->GetTableSize(); 00094 00095 for(G4int i=0; i<numOfCouples; ++i) 00096 { 00097 const G4Material* material = 00098 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 00099 const G4ElementVector* theElementVector = material->GetElementVector(); 00100 G4int nelm = material->GetNumberOfElements(); 00101 00102 for (G4int j=0; j<nelm; ++j) 00103 { 00104 00105 G4int Z = (G4int)(*theElementVector)[j]->GetZ(); 00106 if(Z < 1) { Z = 1; } 00107 else if(Z > maxZ) { Z = maxZ; } 00108 if(!data[Z]) { ReadData(Z, path); } 00109 } 00110 } 00111 // 00112 00113 if(isInitialised) { return; } 00114 fParticleChange = GetParticleChangeForGamma(); 00115 isInitialised = true; 00116 }
void G4LivermoreGammaConversionModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 230 of file G4LivermoreGammaConversionModel.cc.
References G4Electron::Electron(), fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetDefinition(), G4Element::GetfCoulomb(), G4Element::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamElm::GetlogZ3(), G4DynamicParticle::GetMomentumDirection(), G4IonisParamElm::GetZ3(), G4Positron::Positron(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
00235 { 00236 00237 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler 00238 // cross sections with Coulomb correction. A modified version of the random 00239 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15). 00240 00241 // Note 1 : Effects due to the breakdown of the Born approximation at low 00242 // energy are ignored. 00243 // Note 2 : The differential cross section implicitly takes account of 00244 // pair creation in both nuclear and atomic electron fields. However triplet 00245 // prodution is not generated. 00246 00247 if (verboseLevel > 1) 00248 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl; 00249 00250 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 00251 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); 00252 00253 G4double epsilon ; 00254 G4double epsilon0Local = electron_mass_c2 / photonEnergy ; 00255 00256 // Do it fast if photon energy < 2. MeV 00257 if (photonEnergy < smallEnergy ) 00258 { 00259 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand(); 00260 } 00261 else 00262 { 00263 // Select randomly one element in the current material 00264 00265 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 00266 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy); 00267 00268 if (element == 0) 00269 { 00270 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0" 00271 << G4endl; 00272 return; 00273 } 00274 G4IonisParamElm* ionisation = element->GetIonisation(); 00275 if (ionisation == 0) 00276 { 00277 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0" 00278 << G4endl; 00279 return; 00280 } 00281 00282 // Extract Coulomb factor for this Element 00283 G4double fZ = 8. * (ionisation->GetlogZ3()); 00284 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb()); 00285 00286 // Limits of the screening variable 00287 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ; 00288 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ; 00289 G4double screenMin = std::min(4.*screenFactor,screenMax) ; 00290 00291 // Limits of the energy sampling 00292 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ; 00293 G4double epsilonMin = std::max(epsilon0Local,epsilon1); 00294 G4double epsilonRange = 0.5 - epsilonMin ; 00295 00296 // Sample the energy rate of the created electron (or positron) 00297 G4double screen; 00298 G4double gReject ; 00299 00300 G4double f10 = ScreenFunction1(screenMin) - fZ; 00301 G4double f20 = ScreenFunction2(screenMin) - fZ; 00302 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.); 00303 G4double normF2 = std::max(1.5 * f20,0.); 00304 00305 do 00306 { 00307 if (normF1 / (normF1 + normF2) > G4UniformRand() ) 00308 { 00309 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ; 00310 screen = screenFactor / (epsilon * (1. - epsilon)); 00311 gReject = (ScreenFunction1(screen) - fZ) / f10 ; 00312 } 00313 else 00314 { 00315 epsilon = epsilonMin + epsilonRange * G4UniformRand(); 00316 screen = screenFactor / (epsilon * (1 - epsilon)); 00317 gReject = (ScreenFunction2(screen) - fZ) / f20 ; 00318 } 00319 } while ( gReject < G4UniformRand() ); 00320 00321 } // End of epsilon sampling 00322 00323 // Fix charges randomly 00324 00325 G4double electronTotEnergy; 00326 G4double positronTotEnergy; 00327 00328 if (G4UniformRand() > 0.5) 00329 { 00330 electronTotEnergy = (1. - epsilon) * photonEnergy; 00331 positronTotEnergy = epsilon * photonEnergy; 00332 } 00333 else 00334 { 00335 positronTotEnergy = (1. - epsilon) * photonEnergy; 00336 electronTotEnergy = epsilon * photonEnergy; 00337 } 00338 00339 // Scattered electron (positron) angles. ( Z - axis along the parent photon) 00340 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211), 00341 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977) 00342 00343 G4double u; 00344 const G4double a1 = 0.625; 00345 G4double a2 = 3. * a1; 00346 // G4double d = 27. ; 00347 00348 // if (9. / (9. + d) > G4UniformRand()) 00349 if (0.25 > G4UniformRand()) 00350 { 00351 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ; 00352 } 00353 else 00354 { 00355 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ; 00356 } 00357 00358 G4double thetaEle = u*electron_mass_c2/electronTotEnergy; 00359 G4double thetaPos = u*electron_mass_c2/positronTotEnergy; 00360 G4double phi = twopi * G4UniformRand(); 00361 00362 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle); 00363 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos); 00364 00365 00366 // Kinematics of the created pair: 00367 // the electron and positron are assumed to have a symetric angular 00368 // distribution with respect to the Z axis along the parent photon 00369 00370 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 00371 00372 G4ThreeVector electronDirection (dxEle, dyEle, dzEle); 00373 electronDirection.rotateUz(photonDirection); 00374 00375 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(), 00376 electronDirection, 00377 electronKineEnergy); 00378 00379 // The e+ is always created 00380 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; 00381 00382 G4ThreeVector positronDirection (dxPos, dyPos, dzPos); 00383 positronDirection.rotateUz(photonDirection); 00384 00385 // Create G4DynamicParticle object for the particle2 00386 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(), 00387 positronDirection, 00388 positronKineEnergy); 00389 // Fill output vector 00390 fvect->push_back(particle1); 00391 fvect->push_back(particle2); 00392 00393 // kill incident photon 00394 fParticleChange->SetProposedKineticEnergy(0.); 00395 fParticleChange->ProposeTrackStatus(fStopAndKill); 00396 00397 }