#include <G4LivermorePolarizedComptonModel.hh>
Inheritance diagram for G4LivermorePolarizedComptonModel:
Public Member Functions | |
G4LivermorePolarizedComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedCompton") | |
virtual | ~G4LivermorePolarizedComptonModel () |
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 45 of file G4LivermorePolarizedComptonModel.hh.
G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "LivermorePolarizedCompton" | |||
) |
Definition at line 51 of file G4LivermorePolarizedComptonModel.cc.
References G4cout, G4endl, and G4VEmModel::SetHighEnergyLimit().
00053 :G4VEmModel(nam),fParticleChange(0),isInitialised(false), 00054 meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0) 00055 { 00056 lowEnergyLimit = 250 * eV; 00057 highEnergyLimit = 100 * GeV; 00058 //SetLowEnergyLimit(lowEnergyLimit); 00059 SetHighEnergyLimit(highEnergyLimit); 00060 00061 verboseLevel= 0; 00062 // Verbosity scale: 00063 // 0 = nothing 00064 // 1 = warning for energy non-conservation 00065 // 2 = details of energy budget 00066 // 3 = calculation of cross sections, file openings, sampling of atoms 00067 // 4 = entering in methods 00068 00069 if( verboseLevel>0 ) { 00070 G4cout << "Livermore Polarized Compton is constructed " << G4endl 00071 << "Energy range: " 00072 << lowEnergyLimit / eV << " eV - " 00073 << highEnergyLimit / GeV << " GeV" 00074 << G4endl; 00075 } 00076 }
G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel | ( | ) | [virtual] |
Definition at line 80 of file G4LivermorePolarizedComptonModel.cc.
00081 { 00082 if (meanFreePathTable) delete meanFreePathTable; 00083 if (crossSectionHandler) delete crossSectionHandler; 00084 if (scatterFunctionData) delete scatterFunctionData; 00085 }
G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0 , |
|||
G4double | cut = 0 , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 143 of file G4LivermorePolarizedComptonModel.cc.
References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.
00148 { 00149 if (verboseLevel > 3) 00150 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl; 00151 00152 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0; 00153 00154 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 00155 return cs; 00156 }
void G4LivermorePolarizedComptonModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 89 of file G4LivermorePolarizedComptonModel.cc.
References G4VCrossSectionHandler::BuildMeanFreePathForMaterials(), G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4ShellData::LoadData(), G4VEMDataSet::LoadData(), G4VCrossSectionHandler::LoadData(), G4VEmModel::LowEnergyLimit(), and G4ShellData::SetOccupancyData().
00091 { 00092 if (verboseLevel > 3) 00093 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl; 00094 00095 if (crossSectionHandler) 00096 { 00097 crossSectionHandler->Clear(); 00098 delete crossSectionHandler; 00099 } 00100 00101 // Reading of data files - all materials are read 00102 00103 crossSectionHandler = new G4CrossSectionHandler; 00104 crossSectionHandler->Clear(); 00105 G4String crossSectionFile = "comp/ce-cs-"; 00106 crossSectionHandler->LoadData(crossSectionFile); 00107 00108 meanFreePathTable = 0; 00109 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials(); 00110 00111 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation; 00112 G4String scatterFile = "comp/ce-sf-"; 00113 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.); 00114 scatterFunctionData->LoadData(scatterFile); 00115 00116 // For Doppler broadening 00117 shellData.SetOccupancyData(); 00118 G4String file = "/doppler/shell-doppler"; 00119 shellData.LoadData(file); 00120 00121 if (verboseLevel > 2) 00122 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl; 00123 00124 InitialiseElementSelectors(particle,cuts); 00125 00126 if( verboseLevel>0 ) { 00127 G4cout << "Livermore Polarized Compton model is initialized " << G4endl 00128 << "Energy range: " 00129 << LowEnergyLimit() / eV << " eV - " 00130 << HighEnergyLimit() / GeV << " GeV" 00131 << G4endl; 00132 } 00133 00134 // 00135 00136 if(isInitialised) return; 00137 fParticleChange = GetParticleChangeForGamma(); 00138 isInitialised = true; 00139 }
void G4LivermorePolarizedComptonModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 160 of file G4LivermorePolarizedComptonModel.cc.
References G4ShellData::BindingEnergy(), G4Electron::Electron(), G4VEMDataSet::FindValue(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetPolarization(), G4Element::GetZ(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4ParticleChangeForGamma::ProposePolarization(), G4VParticleChange::ProposeTrackStatus(), G4DopplerProfile::RandomSelectMomentum(), G4VEmModel::SelectRandomAtom(), G4ShellData::SelectRandomShell(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
00165 { 00166 // The scattered gamma energy is sampled according to Klein - Nishina formula. 00167 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15). 00168 // GEANT4 internal units 00169 // 00170 // Note : Effects due to binding of atomic electrons are negliged. 00171 00172 if (verboseLevel > 3) 00173 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl; 00174 00175 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy(); 00176 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization(); 00177 00178 // Protection: a polarisation parallel to the 00179 // direction causes problems; 00180 // in that case find a random polarization 00181 00182 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection(); 00183 00184 // Make sure that the polarization vector is perpendicular to the 00185 // gamma direction. If not 00186 00187 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0)) 00188 { // only for testing now 00189 gammaPolarization0 = GetRandomPolarization(gammaDirection0); 00190 } 00191 else 00192 { 00193 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0) 00194 { 00195 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0); 00196 } 00197 } 00198 00199 // End of Protection 00200 00201 // Within energy limit? 00202 00203 if(gammaEnergy0 <= lowEnergyLimit) 00204 { 00205 fParticleChange->ProposeTrackStatus(fStopAndKill); 00206 fParticleChange->SetProposedKineticEnergy(0.); 00207 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0); 00208 return; 00209 } 00210 00211 G4double E0_m = gammaEnergy0 / electron_mass_c2 ; 00212 00213 // Select randomly one element in the current material 00214 //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0); 00215 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 00216 const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0); 00217 G4int Z = (G4int)elm->GetZ(); 00218 00219 // Sample the energy and the polarization of the scattered photon 00220 00221 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ; 00222 00223 G4double epsilon0Local = 1./(1. + 2*E0_m); 00224 G4double epsilon0Sq = epsilon0Local*epsilon0Local; 00225 G4double alpha1 = - std::log(epsilon0Local); 00226 G4double alpha2 = 0.5*(1.- epsilon0Sq); 00227 00228 G4double wlGamma = h_Planck*c_light/gammaEnergy0; 00229 G4double gammaEnergy1; 00230 G4ThreeVector gammaDirection1; 00231 00232 do { 00233 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) 00234 { 00235 epsilon = std::exp(-alpha1*G4UniformRand()); 00236 epsilonSq = epsilon*epsilon; 00237 } 00238 else 00239 { 00240 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand(); 00241 epsilon = std::sqrt(epsilonSq); 00242 } 00243 00244 onecost = (1.- epsilon)/(epsilon*E0_m); 00245 sinThetaSqr = onecost*(2.-onecost); 00246 00247 // Protection 00248 if (sinThetaSqr > 1.) 00249 { 00250 G4cout 00251 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00252 << "sin(theta)**2 = " 00253 << sinThetaSqr 00254 << "; set to 1" 00255 << G4endl; 00256 sinThetaSqr = 1.; 00257 } 00258 if (sinThetaSqr < 0.) 00259 { 00260 G4cout 00261 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00262 << "sin(theta)**2 = " 00263 << sinThetaSqr 00264 << "; set to 0" 00265 << G4endl; 00266 sinThetaSqr = 0.; 00267 } 00268 // End protection 00269 00270 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);; 00271 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 00272 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction; 00273 00274 } while(greject < G4UniformRand()*Z); 00275 00276 00277 // **************************************************** 00278 // Phi determination 00279 // **************************************************** 00280 00281 G4double phi = SetPhi(epsilon,sinThetaSqr); 00282 00283 // 00284 // scattered gamma angles. ( Z - axis along the parent gamma) 00285 // 00286 00287 G4double cosTheta = 1. - onecost; 00288 00289 // Protection 00290 00291 if (cosTheta > 1.) 00292 { 00293 G4cout 00294 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00295 << "cosTheta = " 00296 << cosTheta 00297 << "; set to 1" 00298 << G4endl; 00299 cosTheta = 1.; 00300 } 00301 if (cosTheta < -1.) 00302 { 00303 G4cout 00304 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00305 << "cosTheta = " 00306 << cosTheta 00307 << "; set to -1" 00308 << G4endl; 00309 cosTheta = -1.; 00310 } 00311 // End protection 00312 00313 00314 G4double sinTheta = std::sqrt (sinThetaSqr); 00315 00316 // Protection 00317 if (sinTheta > 1.) 00318 { 00319 G4cout 00320 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00321 << "sinTheta = " 00322 << sinTheta 00323 << "; set to 1" 00324 << G4endl; 00325 sinTheta = 1.; 00326 } 00327 if (sinTheta < -1.) 00328 { 00329 G4cout 00330 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 00331 << "sinTheta = " 00332 << sinTheta 00333 << "; set to -1" 00334 << G4endl; 00335 sinTheta = -1.; 00336 } 00337 // End protection 00338 00339 00340 G4double dirx = sinTheta*std::cos(phi); 00341 G4double diry = sinTheta*std::sin(phi); 00342 G4double dirz = cosTheta ; 00343 00344 00345 // oneCosT , eom 00346 00347 // Doppler broadening - Method based on: 00348 // Y. Namito, S. Ban and H. Hirayama, 00349 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 00350 // NIM A 349, pp. 489-494, 1994 00351 00352 // Maximum number of sampling iterations 00353 00354 G4int maxDopplerIterations = 1000; 00355 G4double bindingE = 0.; 00356 G4double photonEoriginal = epsilon * gammaEnergy0; 00357 G4double photonE = -1.; 00358 G4int iteration = 0; 00359 G4double eMax = gammaEnergy0; 00360 00361 do 00362 { 00363 iteration++; 00364 // Select shell based on shell occupancy 00365 G4int shell = shellData.SelectRandomShell(Z); 00366 bindingE = shellData.BindingEnergy(Z,shell); 00367 00368 eMax = gammaEnergy0 - bindingE; 00369 00370 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) 00371 G4double pSample = profileData.RandomSelectMomentum(Z,shell); 00372 // Rescale from atomic units 00373 G4double pDoppler = pSample * fine_structure_const; 00374 G4double pDoppler2 = pDoppler * pDoppler; 00375 G4double var2 = 1. + onecost * E0_m; 00376 G4double var3 = var2*var2 - pDoppler2; 00377 G4double var4 = var2 - pDoppler2 * cosTheta; 00378 G4double var = var4*var4 - var3 + pDoppler2 * var3; 00379 if (var > 0.) 00380 { 00381 G4double varSqrt = std::sqrt(var); 00382 G4double scale = gammaEnergy0 / var3; 00383 // Random select either root 00384 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale; 00385 else photonE = (var4 + varSqrt) * scale; 00386 } 00387 else 00388 { 00389 photonE = -1.; 00390 } 00391 } while ( iteration <= maxDopplerIterations && 00392 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) ); 00393 00394 // End of recalculation of photon energy with Doppler broadening 00395 // Revert to original if maximum number of iterations threshold has been reached 00396 if (iteration >= maxDopplerIterations) 00397 { 00398 photonE = photonEoriginal; 00399 bindingE = 0.; 00400 } 00401 00402 gammaEnergy1 = photonE; 00403 00404 // 00405 // update G4VParticleChange for the scattered photon 00406 // 00407 00408 // gammaEnergy1 = epsilon*gammaEnergy0; 00409 00410 00411 // New polarization 00412 00413 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon, 00414 sinThetaSqr, 00415 phi, 00416 cosTheta); 00417 00418 // Set new direction 00419 G4ThreeVector tmpDirection1( dirx,diry,dirz ); 00420 gammaDirection1 = tmpDirection1; 00421 00422 // Change reference frame. 00423 00424 SystemOfRefChange(gammaDirection0,gammaDirection1, 00425 gammaPolarization0,gammaPolarization1); 00426 00427 if (gammaEnergy1 > 0.) 00428 { 00429 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ; 00430 fParticleChange->ProposeMomentumDirection( gammaDirection1 ); 00431 fParticleChange->ProposePolarization( gammaPolarization1 ); 00432 } 00433 else 00434 { 00435 gammaEnergy1 = 0.; 00436 fParticleChange->SetProposedKineticEnergy(0.) ; 00437 fParticleChange->ProposeTrackStatus(fStopAndKill); 00438 } 00439 00440 // 00441 // kinematic of the scattered electron 00442 // 00443 00444 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE; 00445 00446 // SI -protection against negative final energy: no e- is created 00447 // like in G4LivermoreComptonModel.cc 00448 if(ElecKineEnergy < 0.0) { 00449 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1); 00450 return; 00451 } 00452 00453 // SI - Removed range test 00454 00455 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2)); 00456 00457 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 - 00458 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum)); 00459 00460 fParticleChange->ProposeLocalEnergyDeposit(bindingE); 00461 00462 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ; 00463 fvect->push_back(dp); 00464 00465 }
Definition at line 73 of file G4LivermorePolarizedComptonModel.hh.
Referenced by Initialise(), and SampleSecondaries().