G4LivermoreComptonModel.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 "G4LivermoreComptonModel.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 G4LivermoreComptonModel::G4LivermoreComptonModel(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 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 G4LivermoreComptonModel::~G4LivermoreComptonModel()
00096 {  
00097   delete crossSectionHandler;
00098   delete scatterFunctionData;
00099 }
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00102 
00103 void G4LivermoreComptonModel::Initialise(const G4ParticleDefinition* particle,
00104                                          const G4DataVector& cuts)
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 }
00153 
00154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00155 
00156 G4double G4LivermoreComptonModel::ComputeCrossSectionPerAtom(
00157                                        const G4ParticleDefinition*,
00158                                              G4double GammaEnergy,
00159                                              G4double Z, G4double,
00160                                              G4double, G4double)
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 }
00170 
00171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00172 
00173 void G4LivermoreComptonModel::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 << "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 }
00393 

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