G4LivermoreComptonModifiedModel.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 //
00029 // Author: Sebastien Incerti
00030 //         30 October 2008
00031 //         on base of G4LowEnergyCompton developed by A.Forti and M.G.Pia
00032 //
00033 // History:
00034 // --------
00035 // 18 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
00036 //                  - apply internal high-energy limit only in constructor 
00037 //                  - do not apply low-energy limit (default is 0)
00038 //                  - remove GetMeanFreePath method and table
00039 //                  - added protection against numerical problem in energy sampling 
00040 //                  - use G4ElementSelector
00041 // 26 Dec 2010   V Ivanchenko Load data tables only once to avoid memory leak
00042 // 30 May 2011   V Ivanchenko Migration to model design for deexcitation
00043 
00044 #include "G4LivermoreComptonModifiedModel.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4Electron.hh"
00048 #include "G4ParticleChangeForGamma.hh"
00049 #include "G4LossTableManager.hh"
00050 #include "G4VAtomDeexcitation.hh"
00051 #include "G4AtomicShell.hh"
00052 #include "G4CrossSectionHandler.hh"
00053 #include "G4CompositeEMDataSet.hh"
00054 #include "G4LogLogInterpolation.hh"
00055 #include "G4Gamma.hh"
00056 
00057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00058 
00059 using namespace std;
00060 
00061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00062 
00063 G4LivermoreComptonModifiedModel::G4LivermoreComptonModifiedModel(const G4ParticleDefinition*,
00064                                                  const G4String& nam)
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 }
00092 
00093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00094 
00095 G4LivermoreComptonModifiedModel::~G4LivermoreComptonModifiedModel()
00096 {  
00097   delete crossSectionHandler;
00098   delete scatterFunctionData;
00099 }
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00102 
00103 void G4LivermoreComptonModifiedModel::Initialise(const G4ParticleDefinition* particle,
00104                                          const G4DataVector& cuts)
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 }
00153 
00154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00155 
00156 G4double G4LivermoreComptonModifiedModel::ComputeCrossSectionPerAtom(
00157                                        const G4ParticleDefinition*,
00158                                              G4double GammaEnergy,
00159                                              G4double Z, G4double,
00160                                              G4double, G4double)
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 }
00170 
00171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00172 
00173 void G4LivermoreComptonModifiedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00174                                                 const G4MaterialCutsCouple* couple,
00175                                                 const G4DynamicParticle* aDynamicGamma,
00176                                                 G4double, G4double)
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 }
00427 

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