#include <G4LivermorePolarizedGammaConversionModel.hh>
Inheritance diagram for G4LivermorePolarizedGammaConversionModel:
Public Member Functions | |
G4LivermorePolarizedGammaConversionModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedGammaConversion") | |
virtual | ~G4LivermorePolarizedGammaConversionModel () |
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) |
Protected Member Functions | |
G4double | GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) |
Protected Attributes | |
G4ParticleChangeForGamma * | fParticleChange |
Definition at line 45 of file G4LivermorePolarizedGammaConversionModel.hh.
G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LivermorePolarizedGammaConversion" | |||
) |
Definition at line 43 of file G4LivermorePolarizedGammaConversionModel.cc.
References G4cout, G4endl, G4VCrossSectionHandler::Initialise(), G4VEmModel::SetHighEnergyLimit(), and G4VEmModel::SetLowEnergyLimit().
00045 :G4VEmModel(nam),fParticleChange(0), 00046 isInitialised(false),meanFreePathTable(0),crossSectionHandler(0) 00047 { 00048 lowEnergyLimit = 2*electron_mass_c2; 00049 highEnergyLimit = 100 * GeV; 00050 SetLowEnergyLimit(lowEnergyLimit); 00051 SetHighEnergyLimit(highEnergyLimit); 00052 smallEnergy = 2.*MeV; 00053 00054 Phi=0.; 00055 Psi=0.; 00056 00057 verboseLevel= 0; 00058 // Verbosity scale: 00059 // 0 = nothing 00060 // 1 = warning for energy non-conservation 00061 // 2 = details of energy budget 00062 // 3 = calculation of cross sections, file openings, samping of atoms 00063 // 4 = entering in methods 00064 00065 if(verboseLevel > 0) { 00066 G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl 00067 << "Energy range: " 00068 << lowEnergyLimit / keV << " keV - " 00069 << highEnergyLimit / GeV << " GeV" 00070 << G4endl; 00071 } 00072 00073 crossSectionHandler = new G4CrossSectionHandler(); 00074 crossSectionHandler->Initialise(0,lowEnergyLimit,highEnergyLimit,400); 00075 }
G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel | ( | ) | [virtual] |
G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0 , |
|||
G4double | cut = 0 , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 149 of file G4LivermorePolarizedGammaConversionModel.cc.
References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.
00154 { 00155 if (verboseLevel > 3) { 00156 G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()" 00157 << G4endl; 00158 } 00159 if(Z < 0.9 || GammaEnergy <= lowEnergyLimit) { return 0.0; } 00160 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 00161 return cs; 00162 }
G4double G4LivermorePolarizedGammaConversionModel::GetMeanFreePath | ( | const G4Track & | aTrack, | |
G4double | previousStepSize, | |||
G4ForceCondition * | condition | |||
) | [protected] |
void G4LivermorePolarizedGammaConversionModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 88 of file G4LivermorePolarizedGammaConversionModel.cc.
References G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4VCrossSectionHandler::LoadData(), and G4VEmModel::LowEnergyLimit().
00090 { 00091 if (verboseLevel > 3) 00092 G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl; 00093 00094 if (crossSectionHandler) 00095 { 00096 crossSectionHandler->Clear(); 00097 delete crossSectionHandler; 00098 } 00099 00100 // Energy limits 00101 /* 00102 // V.Ivanchenko: this was meanless check 00103 if (LowEnergyLimit() < lowEnergyLimit) 00104 { 00105 G4cout << "G4LivermorePolarizedGammaConversionModel: low energy limit increased from " << 00106 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl; 00107 // SetLowEnergyLimit(lowEnergyLimit); 00108 } 00109 */ 00110 if (HighEnergyLimit() > highEnergyLimit) 00111 { 00112 G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " << 00113 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl; 00114 // V.Ivanchenko: this is forbidden 00115 // SetHighEnergyLimit(highEnergyLimit); 00116 } 00117 00118 // Reading of data files - all materials are read 00119 00120 crossSectionHandler = new G4CrossSectionHandler; 00121 crossSectionHandler->Clear(); 00122 G4String crossSectionFile = "pair/pp-cs-"; 00123 crossSectionHandler->LoadData(crossSectionFile); 00124 00125 // 00126 if (verboseLevel > 2) { 00127 G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model" 00128 << G4endl; 00129 } 00130 InitialiseElementSelectors(particle,cuts); 00131 00132 if(verboseLevel > 0) { 00133 G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl 00134 << "Energy range: " 00135 << LowEnergyLimit() / keV << " keV - " 00136 << HighEnergyLimit() / GeV << " GeV" 00137 << G4endl; 00138 } 00139 00140 // 00141 if(!isInitialised) { 00142 isInitialised = true; 00143 fParticleChange = GetParticleChangeForGamma(); 00144 } 00145 }
void G4LivermorePolarizedGammaConversionModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 167 of file G4LivermorePolarizedGammaConversionModel.cc.
References G4Electron::Electron(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetDefinition(), G4Element::GetfCoulomb(), G4Element::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamElm::GetlogZ3(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetPolarization(), G4IonisParamElm::GetZ3(), G4Positron::Positron(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
00172 { 00173 00174 // Fluorescence generated according to: 00175 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic 00176 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides", 00177 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997) 00178 00179 if (verboseLevel > 3) 00180 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl; 00181 00182 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 00183 // Within energy limit? 00184 00185 if(photonEnergy <= lowEnergyLimit) 00186 { 00187 fParticleChange->ProposeTrackStatus(fStopAndKill); 00188 fParticleChange->SetProposedKineticEnergy(0.); 00189 return; 00190 } 00191 00192 00193 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization(); 00194 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection(); 00195 00196 // Make sure that the polarization vector is perpendicular to the 00197 // gamma direction. If not 00198 00199 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0)) 00200 { // only for testing now 00201 gammaPolarization0 = GetRandomPolarization(gammaDirection0); 00202 } 00203 else 00204 { 00205 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0) 00206 { 00207 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0); 00208 } 00209 } 00210 00211 // End of Protection 00212 00213 00214 G4double epsilon ; 00215 G4double epsilon0Local = electron_mass_c2 / photonEnergy ; 00216 00217 // Do it fast if photon energy < 2. MeV 00218 00219 if (photonEnergy < smallEnergy ) 00220 { 00221 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand(); 00222 } 00223 else 00224 { 00225 00226 // Select randomly one element in the current material 00227 00228 // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy); 00229 //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy); 00230 00231 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 00232 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy); 00233 00234 /* 00235 if (element == 0) 00236 { 00237 G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - element = 0" << G4endl; 00238 } 00239 */ 00240 00241 G4IonisParamElm* ionisation = element->GetIonisation(); 00242 00243 /* 00244 if (ionisation == 0) 00245 { 00246 G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - ionisation = 0" << G4endl; 00247 } 00248 */ 00249 00250 00251 // Extract Coulomb factor for this Element 00252 00253 G4double fZ = 8. * (ionisation->GetlogZ3()); 00254 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb()); 00255 00256 // Limits of the screening variable 00257 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ; 00258 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ; 00259 G4double screenMin = std::min(4.*screenFactor,screenMax) ; 00260 00261 // Limits of the energy sampling 00262 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ; 00263 G4double epsilonMin = std::max(epsilon0Local,epsilon1); 00264 G4double epsilonRange = 0.5 - epsilonMin ; 00265 00266 // Sample the energy rate of the created electron (or positron) 00267 G4double screen; 00268 G4double gReject ; 00269 00270 G4double f10 = ScreenFunction1(screenMin) - fZ; 00271 G4double f20 = ScreenFunction2(screenMin) - fZ; 00272 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.); 00273 G4double normF2 = std::max(1.5 * f20,0.); 00274 00275 do { 00276 if (normF1 / (normF1 + normF2) > G4UniformRand() ) 00277 { 00278 epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ; 00279 screen = screenFactor / (epsilon * (1. - epsilon)); 00280 gReject = (ScreenFunction1(screen) - fZ) / f10 ; 00281 } 00282 else 00283 { 00284 epsilon = epsilonMin + epsilonRange * G4UniformRand(); 00285 screen = screenFactor / (epsilon * (1 - epsilon)); 00286 gReject = (ScreenFunction2(screen) - fZ) / f20 ; 00287 00288 00289 } 00290 } while ( gReject < G4UniformRand() ); 00291 00292 } // End of epsilon sampling 00293 00294 // Fix charges randomly 00295 00296 G4double electronTotEnergy; 00297 G4double positronTotEnergy; 00298 00299 00300 if (G4int(2*G4UniformRand())) 00301 { 00302 electronTotEnergy = (1. - epsilon) * photonEnergy; 00303 positronTotEnergy = epsilon * photonEnergy; 00304 } 00305 else 00306 { 00307 positronTotEnergy = (1. - epsilon) * photonEnergy; 00308 electronTotEnergy = epsilon * photonEnergy; 00309 } 00310 00311 // Scattered electron (positron) angles. ( Z - axis along the parent photon) 00312 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211), 00313 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977) 00314 00315 /* 00316 G4double u; 00317 const G4double a1 = 0.625; 00318 G4double a2 = 3. * a1; 00319 00320 if (0.25 > G4UniformRand()) 00321 { 00322 u = - log(G4UniformRand() * G4UniformRand()) / a1 ; 00323 } 00324 else 00325 { 00326 u = - log(G4UniformRand() * G4UniformRand()) / a2 ; 00327 } 00328 */ 00329 00330 G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy 00331 00332 G4double cosTheta = 0.; 00333 G4double sinTheta = 0.; 00334 00335 SetTheta(&cosTheta,&sinTheta,Ene); 00336 00337 // G4double theta = u * electron_mass_c2 / photonEnergy ; 00338 // G4double phi = twopi * G4UniformRand() ; 00339 00340 G4double phi,psi=0.; 00341 00342 //corrected e+ e- angular angular distribution //preliminary! 00343 00344 // if(photonEnergy>50*MeV) 00345 // { 00346 phi = SetPhi(photonEnergy); 00347 psi = SetPsi(photonEnergy,phi); 00348 // } 00349 //else 00350 // { 00351 //psi = G4UniformRand()*2.*pi; 00352 //phi = pi; // coplanar 00353 // } 00354 00355 Psi = psi; 00356 Phi = phi; 00357 //G4cout << "PHI " << phi << G4endl; 00358 //G4cout << "PSI " << psi << G4endl; 00359 00360 G4double phie = psi; //azimuthal angle for the electron 00361 00362 G4double dirX = sinTheta*cos(phie); 00363 G4double dirY = sinTheta*sin(phie); 00364 G4double dirZ = cosTheta; 00365 G4ThreeVector electronDirection(dirX,dirY,dirZ); 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 localEnergyDeposit = 0. ; 00371 00372 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 00373 00374 SystemOfRefChange(gammaDirection0,electronDirection, 00375 gammaPolarization0); 00376 00377 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(), 00378 electronDirection, 00379 electronKineEnergy); 00380 00381 // The e+ is always created (even with kinetic energy = 0) for further annihilation 00382 00383 Ene = positronTotEnergy/electron_mass_c2; // Normalized energy 00384 00385 cosTheta = 0.; 00386 sinTheta = 0.; 00387 00388 SetTheta(&cosTheta,&sinTheta,Ene); 00389 G4double phip = phie+phi; //azimuthal angle for the positron 00390 00391 dirX = sinTheta*cos(phip); 00392 dirY = sinTheta*sin(phip); 00393 dirZ = cosTheta; 00394 G4ThreeVector positronDirection(dirX,dirY,dirZ); 00395 00396 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; 00397 SystemOfRefChange(gammaDirection0,positronDirection, 00398 gammaPolarization0); 00399 00400 // Create G4DynamicParticle object for the particle2 00401 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(), 00402 positronDirection, positronKineEnergy); 00403 00404 00405 fvect->push_back(particle1); 00406 fvect->push_back(particle2); 00407 00408 00409 00410 // Kill the incident photon 00411 00412 00413 00414 // Create lists of pointers to DynamicParticles (photons and electrons) 00415 // (Is the electron vector necessary? To be checked) 00416 // std::vector<G4DynamicParticle*>* photonVector = 0; 00417 //std::vector<G4DynamicParticle*> electronVector; 00418 00419 fParticleChange->ProposeMomentumDirection( 0., 0., 0. ); 00420 fParticleChange->SetProposedKineticEnergy(0.); 00421 fParticleChange->ProposeTrackStatus(fStopAndKill); 00422 00423 }
Definition at line 73 of file G4LivermorePolarizedGammaConversionModel.hh.
Referenced by Initialise(), and SampleSecondaries().