#include <G4LivermoreComptonModel.hh>
Inheritance diagram for G4LivermoreComptonModel:
Public Member Functions | |
G4LivermoreComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreCompton") | |
virtual | ~G4LivermoreComptonModel () |
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 G4LivermoreComptonModel.hh.
G4LivermoreComptonModel::G4LivermoreComptonModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LivermoreCompton" | |||
) |
Definition at line 63 of file G4LivermoreComptonModel.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 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 }
G4LivermoreComptonModel::~G4LivermoreComptonModel | ( | ) | [virtual] |
G4double G4LivermoreComptonModel::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 G4LivermoreComptonModel.cc.
References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.
00161 { 00162 if (verboseLevel > 3) { 00163 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreComptonModel" << 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 G4LivermoreComptonModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 103 of file G4LivermoreComptonModel.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 G4LivermoreComptonModel::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 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 G4LivermoreComptonModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 173 of file G4LivermoreComptonModel.cc.
References 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(), 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 << "G4LivermoreComptonModel::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 eMax = photonEnergy0; 00273 00274 G4int shellIdx = 0; 00275 00276 do 00277 { 00278 ++iteration; 00279 // Select shell based on shell occupancy 00280 shellIdx = shellData.SelectRandomShell(Z); 00281 bindingE = shellData.BindingEnergy(Z,shellIdx); 00282 00283 eMax = photonEnergy0 - bindingE; 00284 00285 // Randomly sample bound electron momentum 00286 // (memento: the data set is in Atomic Units) 00287 G4double pSample = profileData.RandomSelectMomentum(Z,shellIdx); 00288 // Rescale from atomic units 00289 G4double pDoppler = pSample * fine_structure_const; 00290 G4double pDoppler2 = pDoppler * pDoppler; 00291 G4double var2 = 1. + oneCosT * e0m; 00292 G4double var3 = var2*var2 - pDoppler2; 00293 G4double var4 = var2 - pDoppler2 * cosTheta; 00294 G4double var = var4*var4 - var3 + pDoppler2 * var3; 00295 if (var > 0.) 00296 { 00297 G4double varSqrt = std::sqrt(var); 00298 G4double scale = photonEnergy0 / var3; 00299 // Random select either root 00300 if (G4UniformRand() < 0.5) { photonE = (var4 - varSqrt) * scale; } 00301 else { photonE = (var4 + varSqrt) * scale; } 00302 } 00303 else 00304 { 00305 photonE = -1.; 00306 } 00307 } while ( iteration <= maxDopplerIterations && 00308 00309 // JB : corrected the following condition 00310 // (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) ); 00311 00312 (photonE > eMax ) ); 00313 00314 // End of recalculation of photon energy with Doppler broadening 00315 // Revert to original if maximum number of iterations threshold has been reached 00316 00317 if (iteration >= maxDopplerIterations) 00318 { 00319 photonE = photonEoriginal; 00320 bindingE = 0.; 00321 } 00322 00323 // Update G4VParticleChange for the scattered photon 00324 00325 G4ThreeVector photonDirection1(dirx,diry,dirz); 00326 photonDirection1.rotateUz(photonDirection0); 00327 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 00328 00329 G4double photonEnergy1 = photonE; 00330 00331 if (photonEnergy1 > 0.) 00332 { 00333 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ; 00334 } 00335 else 00336 { 00337 photonEnergy1 = 0.; 00338 fParticleChange->SetProposedKineticEnergy(0.) ; 00339 fParticleChange->ProposeTrackStatus(fStopAndKill); 00340 } 00341 00342 // Kinematics of the scattered electron 00343 G4double eKineticEnergy = photonEnergy0 - photonEnergy1 - bindingE; 00344 00345 // protection against negative final energy: no e- is created 00346 if(eKineticEnergy < 0.0) { 00347 bindingE = photonEnergy0 - photonEnergy1; 00348 00349 } else { 00350 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2; 00351 00352 G4double electronE = photonEnergy0 * (1. - epsilon) + electron_mass_c2; 00353 G4double electronP2 = electronE*electronE - electron_mass_c2*electron_mass_c2; 00354 G4double sinThetaE = -1.; 00355 G4double cosThetaE = 0.; 00356 if (electronP2 > 0.) 00357 { 00358 cosThetaE = (eTotalEnergy + photonEnergy1 )* (1. - epsilon) / std::sqrt(electronP2); 00359 sinThetaE = -sqrt((1. - cosThetaE) * (1. + cosThetaE)); 00360 } 00361 00362 G4double eDirX = sinThetaE * std::cos(phi); 00363 G4double eDirY = sinThetaE * std::sin(phi); 00364 G4double eDirZ = cosThetaE; 00365 00366 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 00367 eDirection.rotateUz(photonDirection0); 00368 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(), 00369 eDirection,eKineticEnergy) ; 00370 fvect->push_back(dp); 00371 } 00372 00373 // sample deexcitation 00374 // 00375 if(fAtomDeexcitation && iteration < maxDopplerIterations) { 00376 G4int index = couple->GetIndex(); 00377 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 00378 size_t nbefore = fvect->size(); 00379 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx); 00380 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 00381 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 00382 size_t nafter = fvect->size(); 00383 if(nafter > nbefore) { 00384 for (size_t i=nbefore; i<nafter; ++i) { 00385 bindingE -= ((*fvect)[i])->GetKineticEnergy(); 00386 } 00387 } 00388 } 00389 } 00390 if(bindingE < 0.0) { bindingE = 0.0; } 00391 fParticleChange->ProposeLocalEnergyDeposit(bindingE); 00392 }
Definition at line 75 of file G4LivermoreComptonModel.hh.
Referenced by Initialise(), and SampleSecondaries().