G4LivermoreComptonModel Class Reference

#include <G4LivermoreComptonModel.hh>

Inheritance diagram for G4LivermoreComptonModel:

G4VEmModel

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

G4ParticleChangeForGammafParticleChange

Detailed Description

Definition at line 48 of file G4LivermoreComptonModel.hh.


Constructor & Destructor Documentation

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]

Definition at line 95 of file G4LivermoreComptonModel.cc.

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


Member Function Documentation

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 }


Field Documentation

G4ParticleChangeForGamma* G4LivermoreComptonModel::fParticleChange [protected]

Definition at line 75 of file G4LivermoreComptonModel.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