G4LivermoreGammaConversionModelRC Class Reference

#include <G4LivermoreGammaConversionModelRC.hh>

Inheritance diagram for G4LivermoreGammaConversionModelRC:

G4VEmModel

Public Member Functions

 G4LivermoreGammaConversionModelRC (const G4ParticleDefinition *p=0, const G4String &nam="LivermoreConversion")
virtual ~G4LivermoreGammaConversionModelRC ()
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 Member Functions

G4double GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)

Protected Attributes

G4ParticleChangeForGamma * fParticleChange

Detailed Description

Definition at line 44 of file G4LivermoreGammaConversionModelRC.hh.

Constructor & Destructor Documentation

G4LivermoreGammaConversionModelRC::G4LivermoreGammaConversionModelRC ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "LivermoreConversion" 
)

Definition at line 50 of file G4LivermoreGammaConversionModelRC.cc.

References G4cout, G4endl, and G4VEmModel::SetHighEnergyLimit().

00052   :G4VEmModel(nam),fParticleChange(0),smallEnergy(2.*MeV),isInitialised(false),
00053    crossSectionHandler(0),meanFreePathTable(0)
00054 {
00055   lowEnergyLimit = 2.0*electron_mass_c2;
00056   highEnergyLimit = 100 * GeV;
00057   SetHighEnergyLimit(highEnergyLimit);
00058          
00059   verboseLevel= 0;
00060   // Verbosity scale:
00061   // 0 = nothing 
00062   // 1 = warning for energy non-conservation 
00063   // 2 = details of energy budget
00064   // 3 = calculation of cross sections, file openings, sampling of atoms
00065   // 4 = entering in methods
00066 
00067   if(verboseLevel > 0) {
00068     G4cout << "Livermore Gamma conversion is constructed " << G4endl
00069            << "Energy range: "
00070            << lowEnergyLimit / MeV << " MeV - "
00071            << highEnergyLimit / GeV << " GeV"
00072            << G4endl;
00073   }
00074 }

G4LivermoreGammaConversionModelRC::~G4LivermoreGammaConversionModelRC (  )  [virtual]

Definition at line 78 of file G4LivermoreGammaConversionModelRC.cc.

00079 {  
00080   if (crossSectionHandler) delete crossSectionHandler;
00081 }

Member Function Documentation

G4double G4LivermoreGammaConversionModelRC::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  ,
G4double  kinEnergy,
G4double  Z,
G4double  A = 0,
G4double  cut = 0,
G4double  emax = DBL_MAX 
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 126 of file G4LivermoreGammaConversionModelRC.cc.

References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.

00130 {
00131   if (verboseLevel > 3) {
00132     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC" 
00133            << G4endl;
00134   }
00135   if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
00136 
00137   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00138   return cs;
00139 }

G4double G4LivermoreGammaConversionModelRC::GetMeanFreePath ( const G4Track &  aTrack,
G4double  previousStepSize,
G4ForceCondition *  condition 
) [protected]

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

Implements G4VEmModel.

Definition at line 86 of file G4LivermoreGammaConversionModelRC.cc.

References G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VCrossSectionHandler::Initialise(), G4VCrossSectionHandler::LoadData(), and G4VEmModel::LowEnergyLimit().

00088 {
00089   if (verboseLevel > 3)
00090     G4cout << "Calling G4LivermoreGammaConversionModelRC::Initialise()" << G4endl;
00091 
00092   if (crossSectionHandler)
00093   {
00094     crossSectionHandler->Clear();
00095     delete crossSectionHandler;
00096   }
00097 
00098   // Read data tables for all materials
00099   
00100   crossSectionHandler = new G4CrossSectionHandler();
00101   crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
00102   G4String crossSectionFile = "pair/pp-cs-";
00103   crossSectionHandler->LoadData(crossSectionFile);
00104 
00105   //
00106   
00107   if (verboseLevel > 2) 
00108     G4cout << "Loaded cross section files for Livermore Gamma Conversion model RC" << G4endl;
00109 
00110   if (verboseLevel > 0) { 
00111     G4cout << "Livermore Gamma Conversion model is initialized " << G4endl
00112            << "Energy range: "
00113            << LowEnergyLimit() / MeV << " MeV - "
00114            << HighEnergyLimit() / GeV << " GeV"
00115            << G4endl;
00116   }
00117 
00118   if(isInitialised) return;
00119   fParticleChange = GetParticleChangeForGamma();
00120   isInitialised = true;
00121 }

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

Implements G4VEmModel.

Definition at line 143 of file G4LivermoreGammaConversionModelRC.cc.

References G4Electron::Electron(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetDefinition(), G4Element::GetfCoulomb(), G4Element::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamElm::GetlogZ3(), G4DynamicParticle::GetMomentumDirection(), G4IonisParamElm::GetZ3(), G4Positron::Positron(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().

00148 {
00149 
00150 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler
00151 // cross sections with Coulomb correction. A modified version of the random
00152 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15).
00153 
00154 // Note 1 : Effects due to the breakdown of the Born approximation at low
00155 // energy are ignored.
00156 // Note 2 : The differential cross section implicitly takes account of
00157 // pair creation in both nuclear and atomic electron fields. However triplet
00158 // prodution is not generated.
00159 
00160   if (verboseLevel > 3)
00161     G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC" << G4endl;
00162 
00163   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00164   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00165 
00166   G4double epsilon ;
00167   G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
00168   G4double electronTotEnergy;
00169   G4double positronTotEnergy;
00170 
00171 
00172   // Do it fast if photon energy < 2. MeV
00173   if (photonEnergy < smallEnergy )
00174     {
00175       epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
00176  
00177      if (G4int(2*G4UniformRand()))
00178       {
00179         electronTotEnergy = (1. - epsilon) * photonEnergy;
00180         positronTotEnergy = epsilon * photonEnergy;
00181       }
00182      else
00183       {
00184          positronTotEnergy = (1. - epsilon) * photonEnergy;
00185          electronTotEnergy = epsilon * photonEnergy;
00186        }
00187     }
00188   else
00189     {
00190       // Select randomly one element in the current material
00191       //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
00192       const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
00193       const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
00194       G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries" << G4endl;
00195 
00196       if (element == 0)
00197         {
00198           G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0" 
00199                  << G4endl;
00200           return;
00201         }
00202       G4IonisParamElm* ionisation = element->GetIonisation();
00203       if (ionisation == 0)
00204         {
00205           G4cout << "G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0" 
00206                  << G4endl;
00207           return;
00208         }
00209 
00210       // Extract Coulomb factor for this Element
00211       G4double fZ = 8. * (ionisation->GetlogZ3());
00212       if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
00213 
00214       // Limits of the screening variable
00215       G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
00216       G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
00217       G4double screenMin = std::min(4.*screenFactor,screenMax) ;
00218 
00219       // Limits of the energy sampling
00220       G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
00221       G4double epsilonMin = std::max(epsilon0Local,epsilon1);
00222       G4double epsilonRange = 0.5 - epsilonMin ;
00223 
00224       // Sample the energy rate of the created electron (or positron)
00225       G4double screen;
00226       G4double gReject ;
00227 
00228       G4double f10 = ScreenFunction1(screenMin) - fZ;
00229       G4double f20 = ScreenFunction2(screenMin) - fZ;
00230       G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
00231       G4double normF2 = std::max(1.5 * f20,0.);
00232       G4double a=393.3750918, b=115.3070201, c=810.6428451, d=19.96497475, e=1016.874592, f=1.936685510,
00233                gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
00234       G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287, ff=-117602.840,
00235                gg=-721455.467, hh=693957.8635, ii=156266.1085, jj=533209.9347;                            
00236       G4double Rechazo = 0.;
00237       G4double logepsMin = log(epsilonMin);
00238       G4double NormaRC = a + b*logepsMin + c/logepsMin + d*pow(logepsMin,2.) + e/pow(logepsMin,2.) + f*pow(logepsMin,3.) +
00239                             gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
00240                             k/pow(logepsMin,5.);
00241  
00242       do {
00243          do {
00244            if (normF1 / (normF1 + normF2) > G4UniformRand() )
00245              {
00246                epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.3333) ;
00247                screen = screenFactor / (epsilon * (1. - epsilon));
00248                gReject = (ScreenFunction1(screen) - fZ) / f10 ;
00249              }
00250            else
00251              {
00252                epsilon = epsilonMin + epsilonRange * G4UniformRand();
00253                screen = screenFactor / (epsilon * (1 - epsilon));
00254                gReject = (ScreenFunction2(screen) - fZ) / f20 ;
00255              }
00256          } while ( gReject < G4UniformRand() );
00257                   
00258          if (G4int(2*G4UniformRand())) epsilon = (1. - epsilon); // Extención de Epsilon hasta 1.
00259        
00260          G4double logepsilon = log(epsilon);
00261          G4double deltaP_R1 = 1. + (a + b*logepsilon + c/logepsilon + d*pow(logepsilon,2.) + e/pow(logepsilon,2.) + 
00262                               f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) + 
00263                               j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
00264          G4double deltaP_R2 = 1.+((aa + cc*logepsilon +  ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
00265                              / (1. + bb*logepsilon + dd*pow(logepsilon,2.) + ff*pow(logepsilon,3.) + hh*pow(logepsilon,4.) 
00266                              + jj*pow(logepsilon,5.) ))/100.;
00267        
00268          if (epsilon <= 0.5)
00269             {
00270             Rechazo = deltaP_R1/NormaRC;
00271             }
00272           else
00273             {
00274             Rechazo = deltaP_R2/NormaRC;
00275             }       
00276             G4cout << Rechazo << " " << NormaRC << " " << epsilon << G4endl;
00277       } while (Rechazo < G4UniformRand() );
00278       
00279       electronTotEnergy = (1. - epsilon) * photonEnergy;
00280       positronTotEnergy = epsilon * photonEnergy;
00281 
00282     }   //  End of epsilon sampling
00283     
00284   // Fix charges randomly
00285 
00286   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
00287   // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
00288   // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
00289 
00290   G4double u;
00291   const G4double a1 = 0.625;
00292   G4double a2 = 3. * a1;
00293   //  G4double d = 27. ;
00294 
00295   //  if (9. / (9. + d) > G4UniformRand())
00296   if (0.25 > G4UniformRand())
00297     {
00298       u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ;
00299     }
00300   else
00301     {
00302       u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ;
00303     }
00304 
00305   G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
00306   G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
00307   G4double phi  = twopi * G4UniformRand();
00308 
00309   G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00310   G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00311   
00312   
00313   // Kinematics of the created pair:
00314   // the electron and positron are assumed to have a symetric angular 
00315   // distribution with respect to the Z axis along the parent photon
00316   
00317   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00318   
00319   // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
00320 
00321   G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00322   electronDirection.rotateUz(photonDirection);
00323       
00324   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00325                                                             electronDirection,
00326                                                             electronKineEnergy);
00327 
00328   // The e+ is always created (even with kinetic energy = 0) for further annihilation
00329   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00330 
00331   // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
00332 
00333   G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00334   positronDirection.rotateUz(photonDirection);   
00335   
00336   // Create G4DynamicParticle object for the particle2 
00337   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00338                                                        positronDirection, positronKineEnergy);
00339   // Fill output vector
00340 //  G4cout  << "Cree el e+ " <<  epsilon << G4endl;
00341   fvect->push_back(particle1);
00342   fvect->push_back(particle2);
00343 
00344   // kill incident photon
00345   fParticleChange->SetProposedKineticEnergy(0.);
00346   fParticleChange->ProposeTrackStatus(fStopAndKill);   
00347 
00348 }

Field Documentation

G4ParticleChangeForGamma* G4LivermoreGammaConversionModelRC::fParticleChange [protected]

Definition at line 72 of file G4LivermoreGammaConversionModelRC.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