G4LivermoreGammaConversionModelRC.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Author: Francesco Longo & Gerardo Depaola
00029 //         on base of G4LivermoreGammaConversionModel
00030 //
00031 // History:
00032 // --------
00033 // 12 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
00034 //                  - apply internal high-energy limit only in constructor 
00035 //                  - do not apply low-energy limit (default is 0)
00036 //                  - use CLHEP electron mass for low-enegry limit
00037 //                  - remove MeanFreePath method and table
00038 
00039 
00040 #include "G4LivermoreGammaConversionModelRC.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 
00044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00045 
00046 using namespace std;
00047 
00048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00049 
00050 G4LivermoreGammaConversionModelRC::G4LivermoreGammaConversionModelRC(const G4ParticleDefinition*,
00051                                                                  const G4String& nam)
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 }
00075 
00076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00077 
00078 G4LivermoreGammaConversionModelRC::~G4LivermoreGammaConversionModelRC()
00079 {  
00080   if (crossSectionHandler) delete crossSectionHandler;
00081 }
00082 
00083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00084 
00085 void 
00086 G4LivermoreGammaConversionModelRC::Initialise(const G4ParticleDefinition*,
00087                                             const G4DataVector&)
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 }
00122 
00123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00124 
00125 G4double 
00126 G4LivermoreGammaConversionModelRC::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00127                                                             G4double GammaEnergy,
00128                                                             G4double Z, G4double,
00129                                                             G4double, G4double)
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 }
00140 
00141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00142 
00143 void G4LivermoreGammaConversionModelRC::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00144                                               const G4MaterialCutsCouple* couple,
00145                                               const G4DynamicParticle* aDynamicGamma,
00146                                               G4double,
00147                                               G4double)
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 }
00349 
00350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00351 
00352 G4double G4LivermoreGammaConversionModelRC::ScreenFunction1(G4double screenVariable)
00353 {
00354   // Compute the value of the screening function 3*phi1 - phi2
00355 
00356   G4double value;
00357   
00358   if (screenVariable > 1.)
00359     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00360   else
00361     value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
00362   
00363   return value;
00364 } 
00365 
00366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00367 
00368 G4double G4LivermoreGammaConversionModelRC::ScreenFunction2(G4double screenVariable)
00369 {
00370   // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
00371   
00372   G4double value;
00373   
00374   if (screenVariable > 1.)
00375     value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
00376   else
00377     value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
00378   
00379   return value;
00380 } 
00381 

Generated on Mon May 27 17:48:48 2013 for Geant4 by  doxygen 1.4.7