G4LivermoreComptonModifiedModel Class Reference

#include <G4LivermoreComptonModifiedModel.hh>

Inheritance diagram for G4LivermoreComptonModifiedModel:

G4VEmModel

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

G4ParticleChangeForGammafParticleChange

Detailed Description

Definition at line 48 of file G4LivermoreComptonModifiedModel.hh.


Constructor & Destructor Documentation

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]

Definition at line 95 of file G4LivermoreComptonModifiedModel.cc.

00096 {  
00097   delete crossSectionHandler;
00098   delete scatterFunctionData;
00099 }


Member Function Documentation

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 }


Field Documentation

G4ParticleChangeForGamma* G4LivermoreComptonModifiedModel::fParticleChange [protected]

Definition at line 75 of file G4LivermoreComptonModifiedModel.hh.

Referenced by Initialise(), and SampleSecondaries().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:25 2013 for Geant4 by  doxygen 1.4.7