#include <G4LivermoreComptonModifiedModel.hh>
Inheritance diagram for G4LivermoreComptonModifiedModel:
Public Member Functions | |
G4LivermoreComptonModifiedModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreModifiedCompton") | |
virtual | ~G4LivermoreComptonModifiedModel () |
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 Attributes | |
G4ParticleChangeForGamma * | fParticleChange |
Definition at line 48 of file G4LivermoreComptonModifiedModel.hh.
G4LivermoreComptonModifiedModel::G4LivermoreComptonModifiedModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LivermoreModifiedCompton" | |||
) |
Definition at line 63 of file G4LivermoreComptonModifiedModel.cc.
References G4cout, G4endl, and G4VEmModel::SetDeexcitationFlag().
00065 :G4VEmModel(nam),fParticleChange(0),isInitialised(false), 00066 scatterFunctionData(0), 00067 crossSectionHandler(0),fAtomDeexcitation(0) 00068 { 00069 lowEnergyLimit = 250 * eV; 00070 highEnergyLimit = 100 * GeV; 00071 00072 verboseLevel=0 ; 00073 // Verbosity scale: 00074 // 0 = nothing 00075 // 1 = warning for energy non-conservation 00076 // 2 = details of energy budget 00077 // 3 = calculation of cross sections, file openings, sampling of atoms 00078 // 4 = entering in methods 00079 00080 if( verboseLevel>0 ) { 00081 G4cout << "Livermore Modified Compton model is constructed " << G4endl 00082 << "Energy range: " 00083 << lowEnergyLimit / eV << " eV - " 00084 << highEnergyLimit / GeV << " GeV" 00085 << G4endl; 00086 } 00087 00088 //Mark this model as "applicable" for atomic deexcitation 00089 SetDeexcitationFlag(true); 00090 00091 }
G4LivermoreComptonModifiedModel::~G4LivermoreComptonModifiedModel | ( | ) | [virtual] |
G4double G4LivermoreComptonModifiedModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0 , |
|||
G4double | cut = 0 , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 156 of file G4LivermoreComptonModifiedModel.cc.
References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.
00161 { 00162 if (verboseLevel > 3) { 00163 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModifiedModel" << G4endl; 00164 } 00165 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; } 00166 00167 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 00168 return cs; 00169 }
void G4LivermoreComptonModifiedModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 103 of file G4LivermoreComptonModifiedModel.cc.
References G4LossTableManager::AtomDeexcitation(), G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4LossTableManager::Instance(), G4ShellData::LoadData(), G4VEMDataSet::LoadData(), G4VCrossSectionHandler::LoadData(), G4VEmModel::LowEnergyLimit(), and G4ShellData::SetOccupancyData().
00105 { 00106 if (verboseLevel > 2) { 00107 G4cout << "Calling G4LivermoreComptonModifiedModel::Initialise()" << G4endl; 00108 } 00109 00110 if (crossSectionHandler) 00111 { 00112 crossSectionHandler->Clear(); 00113 delete crossSectionHandler; 00114 } 00115 delete scatterFunctionData; 00116 00117 // Reading of data files - all materials are read 00118 crossSectionHandler = new G4CrossSectionHandler; 00119 G4String crossSectionFile = "comp/ce-cs-"; 00120 crossSectionHandler->LoadData(crossSectionFile); 00121 00122 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation; 00123 G4String scatterFile = "comp/ce-sf-"; 00124 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.); 00125 scatterFunctionData->LoadData(scatterFile); 00126 00127 // For Doppler broadening 00128 shellData.SetOccupancyData(); 00129 G4String file = "/doppler/shell-doppler"; 00130 shellData.LoadData(file); 00131 00132 InitialiseElementSelectors(particle,cuts); 00133 00134 if (verboseLevel > 2) { 00135 G4cout << "Loaded cross section files for Livermore Modified Compton model" << G4endl; 00136 } 00137 00138 if(isInitialised) { return; } 00139 isInitialised = true; 00140 00141 fParticleChange = GetParticleChangeForGamma(); 00142 00143 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 00144 00145 if( verboseLevel>0 ) { 00146 G4cout << "Livermore Compton model is initialized " << G4endl 00147 << "Energy range: " 00148 << LowEnergyLimit() / eV << " eV - " 00149 << HighEnergyLimit() / GeV << " GeV" 00150 << G4endl; 00151 } 00152 }
void G4LivermoreComptonModifiedModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 173 of file G4LivermoreComptonModifiedModel.cc.
References Alpha, G4ShellData::BindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), G4Electron::Electron(), G4VEMDataSet::FindValue(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4VAtomDeexcitation::GenerateParticles(), G4VAtomDeexcitation::GetAtomicShell(), G4DynamicParticle::GetDefinition(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Material::GetName(), G4Element::GetZ(), G4INCL::Math::pi, G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4DopplerProfile::RandomSelectMomentum(), G4VEmModel::SelectRandomAtom(), G4ShellData::SelectRandomShell(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
00177 { 00178 00179 // The scattered gamma energy is sampled according to Klein - Nishina formula. 00180 // then accepted or rejected depending on the Scattering Function multiplied 00181 // by factor from Klein - Nishina formula. 00182 // Expression of the angular distribution as Klein Nishina 00183 // angular and energy distribution and Scattering fuctions is taken from 00184 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth. 00185 // Phys. Res. B 101 (1995). Method of sampling with form factors is different 00186 // data are interpolated while in the article they are fitted. 00187 // Reference to the article is from J. Stepanek New Photon, Positron 00188 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10 00189 // TeV (draft). 00190 // The random number techniques of Butcher & Messel are used 00191 // (Nucl Phys 20(1960),15). 00192 00193 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 00194 00195 if (verboseLevel > 3) { 00196 G4cout << "G4LivermoreComptonModifiedModel::SampleSecondaries() E(MeV)= " 00197 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName() 00198 << G4endl; 00199 } 00200 00201 // low-energy gamma is absorpted by this process 00202 if (photonEnergy0 <= lowEnergyLimit) 00203 { 00204 fParticleChange->ProposeTrackStatus(fStopAndKill); 00205 fParticleChange->SetProposedKineticEnergy(0.); 00206 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); 00207 return ; 00208 } 00209 00210 G4double e0m = photonEnergy0 / electron_mass_c2 ; 00211 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 00212 00213 // Select randomly one element in the current material 00214 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 00215 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 00216 G4int Z = (G4int)elm->GetZ(); 00217 00218 G4double epsilon0Local = 1. / (1. + 2. * e0m); 00219 G4double epsilon0Sq = epsilon0Local * epsilon0Local; 00220 G4double alpha1 = -std::log(epsilon0Local); 00221 G4double alpha2 = 0.5 * (1. - epsilon0Sq); 00222 00223 G4double wlPhoton = h_Planck*c_light/photonEnergy0; 00224 00225 // Sample the energy of the scattered photon 00226 G4double epsilon; 00227 G4double epsilonSq; 00228 G4double oneCosT; 00229 G4double sinT2; 00230 G4double gReject; 00231 00232 do 00233 { 00234 if ( alpha1/(alpha1+alpha2) > G4UniformRand()) 00235 { 00236 // std::pow(epsilon0Local,G4UniformRand()) 00237 epsilon = std::exp(-alpha1 * G4UniformRand()); 00238 epsilonSq = epsilon * epsilon; 00239 } 00240 else 00241 { 00242 epsilonSq = epsilon0Sq + (1. - epsilon0Sq) * G4UniformRand(); 00243 epsilon = std::sqrt(epsilonSq); 00244 } 00245 00246 oneCosT = (1. - epsilon) / ( epsilon * e0m); 00247 sinT2 = oneCosT * (2. - oneCosT); 00248 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm); 00249 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 00250 gReject = (1. - epsilon * sinT2 / (1. + epsilonSq)) * scatteringFunction; 00251 00252 } while(gReject < G4UniformRand()*Z); 00253 00254 G4double cosTheta = 1. - oneCosT; 00255 G4double sinTheta = std::sqrt (sinT2); 00256 G4double phi = twopi * G4UniformRand() ; 00257 G4double dirx = sinTheta * std::cos(phi); 00258 G4double diry = sinTheta * std::sin(phi); 00259 G4double dirz = cosTheta ; 00260 00261 // Doppler broadening - Method based on: 00262 // Y. Namito, S. Ban and H. Hirayama, 00263 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon 00264 // into the EGS4 Code", NIM A 349, pp. 489-494, 1994 00265 00266 // Maximum number of sampling iterations 00267 G4int maxDopplerIterations = 1000; 00268 G4double bindingE = 0.; 00269 G4double photonEoriginal = epsilon * photonEnergy0; 00270 G4double photonE = -1.; 00271 G4int iteration = 0; 00272 G4double systemE = 0; 00273 G4double ePAU = -1; 00274 G4int shellIdx = 0; 00275 G4double vel_c = 299792458; 00276 G4double momentum_au_to_nat = 1.992851740*std::pow(10.,-24.); 00277 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.); 00278 G4double eMax = -1; 00279 G4double Alpha=0; 00280 do 00281 { 00282 ++iteration; 00283 // Select shell based on shell occupancy 00284 shellIdx = shellData.SelectRandomShell(Z); 00285 bindingE = shellData.BindingEnergy(Z,shellIdx); 00286 00287 00288 00289 // Randomly sample bound electron momentum 00290 // (memento: the data set is in Atomic Units) 00291 G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx); 00292 // Rescale from atomic units 00293 00294 00295 //Kinetic energy of target electron 00296 00297 00298 // Reverse vector projection onto scattering vector 00299 00300 do { 00301 Alpha = G4UniformRand()*pi/2.0; 00302 } while(Alpha >= (pi/2.0)); 00303 00304 ePAU = pSample / std::cos(Alpha); 00305 00306 // Convert to SI and the calculate electron energy in natural units 00307 00308 G4double ePSI = ePAU * momentum_au_to_nat; 00309 G4double u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c; 00310 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp)); 00311 00312 //Total energy of the system 00313 systemE = eEIncident+photonEnergy0; 00314 00315 eMax = systemE - bindingE - electron_mass_c2; 00316 G4double pDoppler = pSample * fine_structure_const; 00317 G4double pDoppler2 = pDoppler * pDoppler; 00318 G4double var2 = 1. + oneCosT * e0m; 00319 G4double var3 = var2*var2 - pDoppler2; 00320 G4double var4 = var2 - pDoppler2 * cosTheta; 00321 G4double var = var4*var4 - var3 + pDoppler2 * var3; 00322 if (var > 0.) 00323 { 00324 G4double varSqrt = std::sqrt(var); 00325 G4double scale = photonEnergy0 / var3; 00326 // Random select either root 00327 if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; } 00328 else { photonE = (var4 + varSqrt) * scale; } 00329 } 00330 else 00331 { 00332 photonE = -1.; 00333 } 00334 } while ( iteration <= maxDopplerIterations && 00335 (photonE < 0. || photonE > eMax ) ); 00336 00337 // End of recalculation of photon energy with Doppler broadening 00338 // Kinematics of the scattered electron 00339 G4double eKineticEnergy = systemE - photonE - bindingE - electron_mass_c2; 00340 00341 // protection against negative final energy: no e- is created 00342 G4double eDirX = 0.0; 00343 G4double eDirY = 0.0; 00344 G4double eDirZ = 1.0; 00345 00346 if(eKineticEnergy < 0.0) { 00347 G4cout << "Error, kinetic energy of electron less than zero" << G4endl; 00348 } 00349 00350 else{ 00351 // Estimation of Compton electron polar angle taken from: 00352 // The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport 00353 // Eqn 2.2.25 Pg 42, NRCC Report PIRS-701 00354 G4double E_num = photonEnergy0 - photonE*cosTheta; 00355 G4double E_dom = sqrt(photonEnergy0*photonEnergy0 + photonE*photonE -2*photonEnergy0*photonE*cosTheta); 00356 G4double cosThetaE = E_num / E_dom; 00357 G4double sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE)); 00358 00359 eDirX = sinThetaE * std::cos(phi); 00360 eDirY = sinThetaE * std::sin(phi); 00361 eDirZ = cosThetaE; 00362 00363 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 00364 eDirection.rotateUz(photonDirection0); 00365 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(), 00366 eDirection,eKineticEnergy) ; 00367 fvect->push_back(dp); 00368 } 00369 00370 00371 // Revert to original if maximum number of iterations threshold has been reached 00372 00373 if (iteration >= maxDopplerIterations) 00374 { 00375 photonE = photonEoriginal; 00376 bindingE = 0.; 00377 } 00378 00379 // Update G4VParticleChange for the scattered photon 00380 00381 G4ThreeVector photonDirection1(dirx,diry,dirz); 00382 photonDirection1.rotateUz(photonDirection0); 00383 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 00384 00385 G4double photonEnergy1 = photonE; 00386 00387 if (photonEnergy1 > 0.) 00388 { 00389 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ; 00390 00391 if (iteration < maxDopplerIterations) 00392 { 00393 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 00394 eDirection.rotateUz(photonDirection0); 00395 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(), 00396 eDirection,eKineticEnergy) ; 00397 fvect->push_back(dp); 00398 } 00399 } 00400 else 00401 { 00402 photonEnergy1 = 0.; 00403 fParticleChange->SetProposedKineticEnergy(0.) ; 00404 fParticleChange->ProposeTrackStatus(fStopAndKill); 00405 } 00406 00407 // sample deexcitation 00408 // 00409 if(fAtomDeexcitation && iteration < maxDopplerIterations) { 00410 G4int index = couple->GetIndex(); 00411 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 00412 size_t nbefore = fvect->size(); 00413 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx); 00414 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 00415 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 00416 size_t nafter = fvect->size(); 00417 if(nafter > nbefore) { 00418 for (size_t i=nbefore; i<nafter; ++i) { 00419 bindingE -= ((*fvect)[i])->GetKineticEnergy(); 00420 } 00421 } 00422 } 00423 } 00424 if(bindingE < 0.0) { bindingE = 0.0; } 00425 fParticleChange->ProposeLocalEnergyDeposit(bindingE); 00426 }
Definition at line 75 of file G4LivermoreComptonModifiedModel.hh.
Referenced by Initialise(), and SampleSecondaries().