#include <G4LowEPComptonModel.hh>
Inheritance diagram for G4LowEPComptonModel:
Public Member Functions | |
G4LowEPComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LowEPComptonModel") | |
virtual | ~G4LowEPComptonModel () |
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 78 of file G4LowEPComptonModel.hh.
G4LowEPComptonModel::G4LowEPComptonModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LowEPComptonModel" | |||
) |
Definition at line 85 of file G4LowEPComptonModel.cc.
References G4cout, G4endl, and G4VEmModel::SetDeexcitationFlag().
00087 :G4VEmModel(nam),fParticleChange(0),isInitialised(false), 00088 scatterFunctionData(0),crossSectionHandler(0),fAtomDeexcitation(0) 00089 { 00090 lowEnergyLimit = 250 * eV; 00091 highEnergyLimit = 100 * GeV; 00092 00093 verboseLevel=0 ; 00094 // Verbosity scale: 00095 // 0 = nothing 00096 // 1 = warning for energy non-conservation 00097 // 2 = details of energy budget 00098 // 3 = calculation of cross sections, file openings, sampling of atoms 00099 // 4 = entering in methods 00100 00101 if( verboseLevel>0 ) { 00102 G4cout << "Low energy photon Compton model is constructed " << G4endl 00103 << "Energy range: " 00104 << lowEnergyLimit / eV << " eV - " 00105 << highEnergyLimit / GeV << " GeV" 00106 << G4endl; 00107 } 00108 00109 //Mark this model as "applicable" for atomic deexcitation 00110 SetDeexcitationFlag(true); 00111 00112 }
G4LowEPComptonModel::~G4LowEPComptonModel | ( | ) | [virtual] |
G4double G4LowEPComptonModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0 , |
|||
G4double | cut = 0 , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 177 of file G4LowEPComptonModel.cc.
References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.
00182 { 00183 if (verboseLevel > 3) { 00184 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LowEPComptonModel" << G4endl; 00185 } 00186 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) { return 0.0; } 00187 00188 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 00189 return cs; 00190 }
void G4LowEPComptonModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 124 of file G4LowEPComptonModel.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().
00126 { 00127 if (verboseLevel > 2) { 00128 G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl; 00129 } 00130 00131 if (crossSectionHandler) 00132 { 00133 crossSectionHandler->Clear(); 00134 delete crossSectionHandler; 00135 } 00136 delete scatterFunctionData; 00137 00138 // Reading of data files - all materials are read 00139 crossSectionHandler = new G4CrossSectionHandler; 00140 G4String crossSectionFile = "comp/ce-cs-"; 00141 crossSectionHandler->LoadData(crossSectionFile); 00142 00143 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation; 00144 G4String scatterFile = "comp/ce-sf-"; 00145 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.); 00146 scatterFunctionData->LoadData(scatterFile); 00147 00148 // For Doppler broadening 00149 shellData.SetOccupancyData(); 00150 G4String file = "/doppler/shell-doppler"; 00151 shellData.LoadData(file); 00152 00153 InitialiseElementSelectors(particle,cuts); 00154 00155 if (verboseLevel > 2) { 00156 G4cout << "Loaded cross section files for low energy photon Compton model" << G4endl; 00157 } 00158 00159 if(isInitialised) { return; } 00160 isInitialised = true; 00161 00162 fParticleChange = GetParticleChangeForGamma(); 00163 00164 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 00165 00166 if( verboseLevel>0 ) { 00167 G4cout << "Low energy photon Compton model is initialized " << G4endl 00168 << "Energy range: " 00169 << LowEnergyLimit() / eV << " eV - " 00170 << HighEnergyLimit() / GeV << " GeV" 00171 << G4endl; 00172 } 00173 }
void G4LowEPComptonModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 199 of file G4LowEPComptonModel.cc.
References G4ShellData::BindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), G4Electron::Electron(), G4VEMDataSet::FindValue(), fParticleChange, fStopAndKill, G4cout, G4endl, G4lrint(), 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().
00203 { 00204 00205 // The scattered gamma energy is sampled according to Klein - Nishina formula. 00206 // then accepted or rejected depending on the Scattering Function multiplied 00207 // by factor from Klein - Nishina formula. 00208 // Expression of the angular distribution as Klein Nishina 00209 // angular and energy distribution and Scattering fuctions is taken from 00210 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth. 00211 // Phys. Res. B 101 (1995). Method of sampling with form factors is different 00212 // data are interpolated while in the article they are fitted. 00213 // Reference to the article is from J. Stepanek New Photon, Positron 00214 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10 00215 // TeV (draft). 00216 // The random number techniques of Butcher & Messel are used 00217 // (Nucl Phys 20(1960),15). 00218 00219 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV; 00220 00221 if (verboseLevel > 3) { 00222 G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= " 00223 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName() 00224 << G4endl; 00225 } 00226 00227 // low-energy gamma is absorpted by this process 00228 if (photonEnergy0 <= lowEnergyLimit) 00229 { 00230 fParticleChange->ProposeTrackStatus(fStopAndKill); 00231 fParticleChange->SetProposedKineticEnergy(0.); 00232 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); 00233 return ; 00234 } 00235 00236 G4double e0m = photonEnergy0 / electron_mass_c2 ; 00237 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 00238 00239 // Select randomly one element in the current material 00240 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 00241 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 00242 G4int Z = (G4int)elm->GetZ(); 00243 00244 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m); 00245 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0; 00246 G4double alpha1 = -std::log(LowEPCepsilon0); 00247 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq); 00248 00249 G4double wlPhoton = h_Planck*c_light/photonEnergy0; 00250 00251 // Sample the energy of the scattered photon 00252 G4double LowEPCepsilon; 00253 G4double LowEPCepsilonSq; 00254 G4double oneCosT; 00255 G4double sinT2; 00256 G4double gReject; 00257 00258 do 00259 { 00260 if ( alpha1/(alpha1+alpha2) > G4UniformRand()) 00261 { 00262 LowEPCepsilon = std::exp(-alpha1 * G4UniformRand()); 00263 LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon; 00264 } 00265 else 00266 { 00267 LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand(); 00268 LowEPCepsilon = std::sqrt(LowEPCepsilonSq); 00269 } 00270 00271 oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m); 00272 sinT2 = oneCosT * (2. - oneCosT); 00273 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm); 00274 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 00275 gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction; 00276 00277 } while(gReject < G4UniformRand()*Z); 00278 00279 G4double cosTheta = 1. - oneCosT; 00280 G4double sinTheta = std::sqrt(sinT2); 00281 G4double phi = twopi * G4UniformRand(); 00282 G4double dirx = sinTheta * std::cos(phi); 00283 G4double diry = sinTheta * std::sin(phi); 00284 G4double dirz = cosTheta ; 00285 00286 00287 // Scatter photon energy and Compton electron direction - Method based on: 00288 // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin' 00289 // "A low energy bound atomic electron Compton scattering model for Geant4" 00290 // TNS ISSUE, PG, 2012 00291 00292 // Set constants and initialize scattering parameters 00293 00294 G4double vel_c = 299792458; 00295 G4double momentum_au_to_nat = (pi/2.0)*1.992851740*std::pow(10.,-24.); 00296 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.); 00297 00298 G4int maxDopplerIterations = 1000; 00299 G4double bindingE = 0.; 00300 G4double pEIncident = photonEnergy0 ; 00301 G4double pERecoil = -1.; 00302 G4double eERecoil = -1.; 00303 G4double e_alpha =0.; 00304 G4double e_beta = 0.; 00305 00306 G4double CE_emission_flag = 0.; 00307 G4double ePAU = -1; 00308 //G4double Alpha=0; 00309 G4int shellIdx = 0; 00310 G4double u_temp = 0; 00311 G4double cosPhiE =0; 00312 G4double sinThetaE =0; 00313 G4double cosThetaE =0; 00314 G4int iteration = 0; 00315 do{ 00316 00317 00318 // ****************************************** 00319 // | Determine scatter photon energy | 00320 // ****************************************** 00321 00322 do 00323 { 00324 iteration++; 00325 00326 00327 // ******************************************** 00328 // | Sample bound electron information | 00329 // ******************************************** 00330 00331 // Select shell based on shell occupancy 00332 00333 shellIdx = shellData.SelectRandomShell(Z); 00334 bindingE = shellData.BindingEnergy(Z,shellIdx)/MeV; 00335 00336 //G4cout << "New sample" << G4endl; 00337 00338 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) 00339 ePAU = profileData.RandomSelectMomentum(Z,shellIdx); 00340 00341 // Convert to SI units 00342 00343 G4double ePSI = ePAU * momentum_au_to_nat; 00344 00345 //Calculate bound electron velocity and normalise to natural units 00346 00347 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c; 00348 00349 // Sample incident electron direction, amorphous material, to scattering photon scattering plane 00350 00351 e_alpha = pi*G4UniformRand(); 00352 e_beta = twopi*G4UniformRand(); 00353 00354 // Total energy of system 00355 00356 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp)); 00357 G4double systemE = eEIncident + pEIncident; 00358 00359 00360 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp)); 00361 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha)); 00362 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha); 00363 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta); 00364 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident); 00365 pERecoil = (numerator/denominator); 00366 eERecoil = systemE - pERecoil; 00367 CE_emission_flag = pEIncident - pERecoil; 00368 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE)); 00369 00370 00371 00372 // End of recalculation of photon energy with Doppler broadening 00373 00374 00375 00376 // ******************************************************* 00377 // | Determine ejected Compton electron direction | 00378 // ******************************************************* 00379 00380 // Calculate velocity of ejected Compton electron 00381 00382 G4double a_temp = eERecoil / electron_mass_c2; 00383 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp))); 00384 00385 // Coefficients and terms from simulatenous equations 00386 00387 G4double sinAlpha = std::sin(e_alpha); 00388 G4double cosAlpha = std::cos(e_alpha); 00389 G4double sinBeta = std::sin(e_beta); 00390 G4double cosBeta = std::cos(e_beta); 00391 00392 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp)); 00393 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp)); 00394 00395 G4double var_A = pERecoil*u_p_temp*sinTheta; 00396 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident); 00397 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta); 00398 00399 G4double var_D1 = gamma*electron_mass_c2*pERecoil; 00400 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha)); 00401 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil); 00402 G4double var_D = var_D1*var_D2 + var_D3; 00403 00404 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha); 00405 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta; 00406 G4double var_E = var_E1 - var_E2; 00407 00408 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha); 00409 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta); 00410 G4double var_F = var_F1 - var_F2; 00411 00412 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha; 00413 00414 // Two equations form a quadratic form of Wx^2 + Yx + Z = 0 00415 // Coefficents and solution to quadratic 00416 00417 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A); 00418 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B); 00419 G4double var_W = var_W1 + var_W2; 00420 00421 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C)); 00422 00423 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C); 00424 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A); 00425 G4double var_Z = var_Z1 + var_Z2; 00426 00427 G4double diff = ((var_Y*var_Y)-4*var_W*var_Z); 00428 00429 00430 // Check if diff has FPE, if so set var_Y*var_Y and 4*var_W*var_Z to six significant figures to fix 00431 if (diff < 0.0) 00432 { 00433 00434 G4double funorder = 0.0; 00435 G4double sf = 6.0; 00436 00437 G4double diff1 = var_Y*var_Y; 00438 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(diff1)))+std::numeric_limits<double>::epsilon()))+sf; 00439 diff1 = G4lrint(diff1*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder)); 00440 00441 G4double diff2 = 4*var_W*var_Z; 00442 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(diff2)))+std::numeric_limits<double>::epsilon()))+sf; 00443 diff2 = G4lrint(diff2*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder)); 00444 00445 diff = diff1 -diff2; 00446 00447 } 00448 00449 00450 // Plus and minus of quadratic 00451 00452 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W); 00453 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W); 00454 00455 00456 // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron 00457 G4double ThetaE = 0.; 00458 G4double sol_select = G4UniformRand(); 00459 00460 if (sol_select < 0.5) 00461 { 00462 ThetaE = std::acos(X_p); 00463 } 00464 if (sol_select > 0.5) 00465 { 00466 ThetaE = std::acos(X_m); 00467 } 00468 00469 cosThetaE = std::cos(ThetaE); 00470 sinThetaE = std::sin(ThetaE); 00471 G4double Theta = std::acos(cosTheta); 00472 00473 //Calculate electron Phi 00474 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE)); 00475 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta)); 00476 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp); 00477 // Trigs 00478 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE); 00479 00480 00481 // Check if cosPhiE has FPE, if so set to four significant figures to fix 00482 if (abs(cosPhiE) > 1.0) 00483 { 00484 G4double funorder = 0.0; 00485 G4double sf = 4.0; 00486 funorder = abs(G4lrint(std::log10(static_cast<double>(abs(cosPhiE)))+std::numeric_limits<double>::epsilon()))+sf; 00487 cosPhiE = G4lrint(cosPhiE*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder)); 00488 00489 } 00490 // End of calculation of ejection Compton electron direction 00491 00492 //Fix for floating point errors 00493 00494 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1)); 00495 00496 // Revert to original if maximum number of iterations threshold has been reached 00497 00498 00499 if (iteration >= maxDopplerIterations) 00500 { 00501 pERecoil = photonEnergy0 ; 00502 bindingE = 0.; 00503 dirx=0.0; 00504 diry=0.0; 00505 dirz=1.0; 00506 } 00507 00508 // Set "scattered" photon direction and energy 00509 00510 G4ThreeVector photonDirection1(dirx,diry,dirz); 00511 photonDirection1.rotateUz(photonDirection0); 00512 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 00513 00514 if (pERecoil > 0.) 00515 { 00516 fParticleChange->SetProposedKineticEnergy(pERecoil) ; 00517 00518 // Set ejected Compton electron direction and energy 00519 G4double PhiE = std::acos(cosPhiE); 00520 G4double eDirX = sinThetaE * std::cos(phi+PhiE); 00521 G4double eDirY = sinThetaE * std::sin(phi+PhiE); 00522 G4double eDirZ = cosThetaE; 00523 00524 G4double eKineticEnergy = pEIncident - pERecoil - bindingE; 00525 00526 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 00527 eDirection.rotateUz(photonDirection0); 00528 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(), 00529 eDirection,eKineticEnergy) ; 00530 fvect->push_back(dp); 00531 00532 } 00533 else 00534 { 00535 fParticleChange->SetProposedKineticEnergy(0.); 00536 fParticleChange->ProposeTrackStatus(fStopAndKill); 00537 } 00538 00539 00540 00541 00542 // sample deexcitation 00543 00544 if(fAtomDeexcitation && iteration < maxDopplerIterations) { 00545 G4int index = couple->GetIndex(); 00546 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 00547 size_t nbefore = fvect->size(); 00548 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx); 00549 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 00550 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 00551 size_t nafter = fvect->size(); 00552 if(nafter > nbefore) { 00553 for (size_t i=nbefore; i<nafter; ++i) { 00554 bindingE -= ((*fvect)[i])->GetKineticEnergy(); 00555 } 00556 } 00557 } 00558 } 00559 if(bindingE < 0.0) { bindingE = 0.0; } 00560 fParticleChange->ProposeLocalEnergyDeposit(bindingE); 00561 00562 }
Definition at line 105 of file G4LowEPComptonModel.hh.
Referenced by Initialise(), and SampleSecondaries().