G4KleinNishinaModel Class Reference

#include <G4KleinNishinaModel.hh>

Inheritance diagram for G4KleinNishinaModel:

G4VEmModel

Public Member Functions

 G4KleinNishinaModel (const G4String &nam="KleinNishina")
virtual ~G4KleinNishinaModel ()
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double cut, G4double emax)
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)

Protected Attributes

G4ParticleDefinition * theGamma
G4ParticleDefinition * theElectron
G4ParticleChangeForGamma * fParticleChange
G4double lowestGammaEnergy

Detailed Description

Definition at line 59 of file G4KleinNishinaModel.hh.

Constructor & Destructor Documentation

G4KleinNishinaModel::G4KleinNishinaModel ( const G4String &  nam = "KleinNishina"  ) 

Definition at line 65 of file G4KleinNishinaModel.cc.

References G4Electron::Electron(), fParticleChange, G4Gamma::Gamma(), lowestGammaEnergy, G4VEmModel::SetDeexcitationFlag(), theElectron, and theGamma.

00066   : G4VEmModel(nam)
00067 {
00068   theGamma = G4Gamma::Gamma();
00069   theElectron = G4Electron::Electron();
00070   lowestGammaEnergy = 1.0*eV;
00071   limitFactor       = 4;
00072   fProbabilities.resize(9,0.0);
00073   SetDeexcitationFlag(true);
00074   fParticleChange = 0;
00075   fAtomDeexcitation = 0;
00076 }

G4KleinNishinaModel::~G4KleinNishinaModel (  )  [virtual]

Definition at line 80 of file G4KleinNishinaModel.cc.

00081 {}

Member Function Documentation

G4double G4KleinNishinaModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  ,
G4double  kinEnergy,
G4double  Z,
G4double  A,
G4double  cut,
G4double  emax 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 96 of file G4KleinNishinaModel.cc.

00100 {
00101   G4double xSection = 0.0 ;
00102   if ( Z < 0.9999 || GammaEnergy < 0.1*keV) { return xSection; }
00103 
00104   static const G4double a = 20.0 , b = 230.0 , c = 440.0;
00105   
00106   static const G4double
00107     d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527   *barn, d4=-1.9798e+1*barn,
00108     e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
00109     f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
00110      
00111   G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
00112            p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
00113 
00114   G4double T0  = 15.0*keV; 
00115   if (Z < 1.5) { T0 = 40.0*keV; } 
00116 
00117   G4double X   = max(GammaEnergy, T0) / electron_mass_c2;
00118   xSection = p1Z*std::log(1.+2.*X)/X
00119                + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00120                 
00121   //  modification for low energy. (special case for Hydrogen)
00122   if (GammaEnergy < T0) {
00123     G4double dT0 = keV;
00124     X = (T0+dT0) / electron_mass_c2 ;
00125     G4double sigma = p1Z*log(1.+2*X)/X
00126                     + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00127     G4double   c1 = -T0*(sigma-xSection)/(xSection*dT0);             
00128     G4double   c2 = 0.150; 
00129     if (Z > 1.5) { c2 = 0.375-0.0556*log(Z); }
00130     G4double    y = log(GammaEnergy/T0);
00131     xSection *= exp(-y*(c1+c2*y));          
00132   }
00133 
00134   if(xSection < 0.0) { xSection = 0.0; }
00135   //  G4cout << "e= " << GammaEnergy << " Z= " << Z 
00136   //  << " cross= " << xSection << G4endl;
00137   return xSection;
00138 }

void G4KleinNishinaModel::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
) [virtual]

Implements G4VEmModel.

Definition at line 85 of file G4KleinNishinaModel.cc.

References G4LossTableManager::AtomDeexcitation(), fParticleChange, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::InitialiseElementSelectors(), and G4LossTableManager::Instance().

00087 {
00088   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00089   InitialiseElementSelectors(p, cuts);
00090   if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
00091 }

void G4KleinNishinaModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  ,
const G4MaterialCutsCouple *  ,
const G4DynamicParticle *  ,
G4double  tmin,
G4double  maxEnergy 
) [virtual]

Implements G4VEmModel.

Definition at line 142 of file G4KleinNishinaModel.cc.

References G4InuclSpecialFunctions::bindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), fParticleChange, fStopAndKill, G4lrint(), G4UniformRand, G4VAtomDeexcitation::GenerateParticles(), G4VAtomDeexcitation::GetAtomicShell(), G4Element::GetAtomicShell(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4Element::GetNbOfAtomicShells(), G4Element::GetNbOfShellElectrons(), G4Element::GetZ(), lowestGammaEnergy, G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), G4ParticleChangeForGamma::SetProposedKineticEnergy(), theElectron, and theGamma.

00148 {
00149   // primary gamma
00150   G4double energy = aDynamicGamma->GetKineticEnergy();
00151   G4ThreeVector direction = aDynamicGamma->GetMomentumDirection();
00152 
00153   // select atom
00154   const G4Element* elm = SelectRandomAtom(couple, theGamma, energy);
00155 
00156   // select shell first
00157   G4int nShells = elm->GetNbOfAtomicShells();
00158   if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
00159   G4double totprob = 0.0;
00160   G4int i;
00161   for(i=0; i<nShells; ++i) {
00162     //G4double bindingEnergy = elm->GetAtomicShell(i);
00163     totprob += elm->GetNbOfShellElectrons(i);
00164     //totprob += elm->GetNbOfShellElectrons(i)/(bindingEnergy*bindingEnergy);
00165     fProbabilities[i] = totprob; 
00166   }
00167   //if(totprob == 0.0) { return; }
00168 
00169   // Loop on sampling
00170   //  const G4int nlooplim = 100;
00171   //G4int nloop = 0;
00172 
00173   G4double bindingEnergy, ePotEnergy, eKinEnergy;
00174   G4double gamEnergy0, gamEnergy1;
00175 
00176   //static const G4double eminus2 =  1.0 - exp(-2.0);
00177 
00178   do {
00179     //++nloop;
00180     G4double xprob = totprob*G4UniformRand();
00181 
00182     // select shell
00183     for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
00184    
00185     bindingEnergy = elm->GetAtomicShell(i);
00186     //    ePotEnergy    = bindingEnergy;
00187     //    gamEnergy0 = energy;
00188     lv1.set(0.0,0.0,energy,energy);
00189 
00190     //G4cout << "nShells= " << nShells << " i= " << i 
00191     //   << " Egamma= " << energy << " Ebind= " << bindingEnergy
00192     //   << " Elim= " << limitEnergy 
00193     //   << G4endl;
00194 
00195     // for rest frame of the electron
00196     G4double x = -log(G4UniformRand());
00197     eKinEnergy = bindingEnergy*x;
00198     ePotEnergy = bindingEnergy*(1.0 + x);
00199 
00200     // for rest frame of the electron
00201     G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
00202     G4double phi = G4UniformRand()*twopi;
00203     G4double costet = 2*G4UniformRand() - 1;
00204     G4double sintet = sqrt((1 - costet)*(1 + costet));
00205     lv2.set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
00206             eTotMomentum*costet,eKinEnergy + electron_mass_c2);
00207     bst = lv2.boostVector();
00208     lv1.boost(-bst);
00209 
00210     gamEnergy0 = lv1.e();
00211    
00212     // In the rest frame of the electron
00213     // The scattered gamma energy is sampled according to Klein-Nishina formula
00214     // The random number techniques of Butcher & Messel are used 
00215     // (Nuc Phys 20(1960),15). 
00216     G4double E0_m = gamEnergy0/electron_mass_c2;
00217 
00218     //
00219     // sample the energy rate of the scattered gamma 
00220     //
00221 
00222     G4double epsilon, epsilonsq, onecost, sint2, greject ;
00223 
00224     G4double eps0       = 1./(1 + 2*E0_m);
00225     G4double epsilon0sq = eps0*eps0;
00226     G4double alpha1     = - log(eps0);
00227     G4double alpha2     = 0.5*(1 - epsilon0sq);
00228 
00229     do {
00230       if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
00231         epsilon   = exp(-alpha1*G4UniformRand());   // epsilon0**r
00232         epsilonsq = epsilon*epsilon; 
00233 
00234       } else {
00235         epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
00236         epsilon   = sqrt(epsilonsq);
00237       }
00238 
00239       onecost = (1.- epsilon)/(epsilon*E0_m);
00240       sint2   = onecost*(2.-onecost);
00241       greject = 1. - epsilon*sint2/(1.+ epsilonsq);
00242 
00243     } while (greject < G4UniformRand());
00244     gamEnergy1 = epsilon*gamEnergy0;
00245  
00246     // before scattering total 4-momentum in e- system
00247     lv2.set(0.0,0.0,0.0,electron_mass_c2);
00248     lv2 += lv1;
00249  
00250     //
00251     // scattered gamma angles. ( Z - axis along the parent gamma)
00252     //
00253     if(sint2 < 0.0) { sint2 = 0.0; }
00254     costet = 1. - onecost; 
00255     sintet = sqrt(sint2);
00256     phi  = twopi * G4UniformRand();
00257 
00258     // e- recoil
00259     //
00260     // in  rest frame of the electron
00261     G4ThreeVector gamDir = lv1.vect().unit();
00262     G4ThreeVector v = G4ThreeVector(sintet*cos(phi),sintet*sin(phi),costet);
00263     v.rotateUz(gamDir);
00264     lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),gamEnergy1*v.z(),gamEnergy1);
00265     lv2 -= lv1;
00266     //G4cout<<"Egam= "<<lv1.e()<<" Ee= "<< lv2.e()-electron_mass_c2 << G4endl;
00267     lv2.boost(bst);
00268     eKinEnergy = lv2.e() - electron_mass_c2 - ePotEnergy;   
00269     //G4cout << "eKinEnergy= " << eKinEnergy << G4endl;
00270 
00271   } while ( eKinEnergy < 0.0 );
00272 
00273   //
00274   // update G4VParticleChange for the scattered gamma
00275   //
00276    
00277   lv1.boost(bst);
00278   gamEnergy1 = lv1.e();
00279   if(gamEnergy1 > lowestGammaEnergy) {
00280     G4ThreeVector gamDirection1 = lv1.vect().unit();
00281     gamDirection1.rotateUz(direction);
00282     fParticleChange->ProposeMomentumDirection(gamDirection1);
00283   } else { 
00284     fParticleChange->ProposeTrackStatus(fStopAndKill);
00285     gamEnergy1 = 0.0;
00286   }
00287   fParticleChange->SetProposedKineticEnergy(gamEnergy1);
00288 
00289   //
00290   // kinematic of the scattered electron
00291   //
00292 
00293   if(eKinEnergy > lowestGammaEnergy) {
00294     G4ThreeVector eDirection = lv2.vect().unit();
00295     eDirection.rotateUz(direction);
00296     G4DynamicParticle* dp = 
00297       new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
00298     fvect->push_back(dp);
00299   } else { eKinEnergy = 0.0; }
00300 
00301   G4double edep = energy - gamEnergy1 - eKinEnergy;
00302   
00303   // sample deexcitation
00304   //
00305   if(fAtomDeexcitation) {
00306     G4int index = couple->GetIndex();
00307     if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00308       G4int Z = G4lrint(elm->GetZ());
00309       G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
00310       const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00311       size_t nbefore = fvect->size();
00312       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00313       size_t nafter = fvect->size();
00314       if(nafter > nbefore) {
00315         for (size_t j=nbefore; j<nafter; ++j) {
00316           edep -= ((*fvect)[j])->GetKineticEnergy();
00317         } 
00318       }
00319     }
00320   }
00321   // energy balance
00322   if(edep < 0.0) { edep = 0.0; }
00323   fParticleChange->ProposeLocalEnergyDeposit(edep);
00324 }

Field Documentation

G4ParticleChangeForGamma* G4KleinNishinaModel::fParticleChange [protected]

Definition at line 88 of file G4KleinNishinaModel.hh.

Referenced by G4KleinNishinaModel(), Initialise(), and SampleSecondaries().

G4double G4KleinNishinaModel::lowestGammaEnergy [protected]

Definition at line 89 of file G4KleinNishinaModel.hh.

Referenced by G4KleinNishinaModel(), and SampleSecondaries().

G4ParticleDefinition* G4KleinNishinaModel::theElectron [protected]

Definition at line 87 of file G4KleinNishinaModel.hh.

Referenced by G4KleinNishinaModel(), and SampleSecondaries().

G4ParticleDefinition* G4KleinNishinaModel::theGamma [protected]

Definition at line 86 of file G4KleinNishinaModel.hh.

Referenced by G4KleinNishinaModel(), and SampleSecondaries().

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