G4LivermorePolarizedPhotoElectricModel.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: G4LivermorePolarizedPhotoElectricModel.hh,v 1.2 2010-11-23 16:42:15 flongo Exp $
00027 // GEANT4 tag $Name: not supported by cvs2svn $
00028 //
00029 // Authors: G.Depaola & F.Longo
00030 //
00031 
00032 #include "G4LivermorePolarizedPhotoElectricModel.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4LossTableManager.hh"
00036 #include "G4VAtomDeexcitation.hh"
00037 #include "G4PhotoElectricAngularGeneratorPolarized.hh"
00038 #include "G4SauterGavrilaAngularDistribution.hh"
00039 #include "G4AtomicShell.hh"
00040 #include "G4ProductionCutsTable.hh"
00041 #include "G4Gamma.hh"
00042 
00043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00044 
00045 using namespace std;
00046 
00047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00048 
00049 G4LivermorePolarizedPhotoElectricModel::G4LivermorePolarizedPhotoElectricModel(
00050    const G4String& nam)
00051   :G4VEmModel(nam),fParticleChange(0),
00052    crossSectionHandler(0), shellCrossSectionHandler(0)
00053 {
00054   lowEnergyLimit  = 10 * eV; // SI - Could be 10 eV ?
00055   highEnergyLimit = 100 * GeV;
00056   
00057   verboseLevel= 0;
00058   // Verbosity scale:
00059   // 0 = nothing 
00060   // 1 = warning for energy non-conservation 
00061   // 2 = details of energy budget
00062   // 3 = calculation of cross sections, file openings, sampling of atoms
00063   // 4 = entering in methods
00064 
00065   theGamma    = G4Gamma::Gamma();
00066   theElectron = G4Electron::Electron();
00067 
00068   // use default generator
00069   SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
00070 
00071   // polarized generator needs fix
00072   //SetAngularDistribution(new G4PhotoElectricAngularGeneratorPolarized());
00073   SetDeexcitationFlag(true);
00074   fAtomDeexcitation = 0; 
00075   fDeexcitationActive = false; 
00076 
00077   if (verboseLevel > 1) {
00078     G4cout << "Livermore Polarized PhotoElectric is constructed " << G4endl
00079            << "Energy range: "
00080            << lowEnergyLimit / keV << " keV - "
00081            << highEnergyLimit / GeV << " GeV"
00082            << G4endl;
00083   }
00084 }
00085 
00086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00087 
00088 G4LivermorePolarizedPhotoElectricModel::~G4LivermorePolarizedPhotoElectricModel()
00089 {  
00090   delete crossSectionHandler;
00091   delete shellCrossSectionHandler;
00092 }
00093 
00094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00095 
00096 void 
00097 G4LivermorePolarizedPhotoElectricModel::Initialise(
00098                     const G4ParticleDefinition*,
00099                     const G4DataVector&)
00100 {
00101   if (verboseLevel > 3) {
00102     G4cout << "Calling G4LivermorePolarizedPhotoElectricModel::Initialise()" << G4endl;
00103   }
00104   if (crossSectionHandler)
00105   {
00106     crossSectionHandler->Clear();
00107     delete crossSectionHandler;
00108   }
00109 
00110   if (shellCrossSectionHandler)
00111     {
00112       shellCrossSectionHandler->Clear();
00113       delete shellCrossSectionHandler;
00114     }
00115 
00116   
00117   // Reading of data files - all materials are read  
00118   crossSectionHandler = new G4CrossSectionHandler;
00119   crossSectionHandler->Clear();
00120   G4String crossSectionFile = "phot/pe-cs-";
00121   crossSectionHandler->LoadData(crossSectionFile);
00122 
00123   shellCrossSectionHandler = new G4CrossSectionHandler();
00124   shellCrossSectionHandler->Clear();
00125   G4String shellCrossSectionFile = "phot/pe-ss-cs-";
00126   shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
00127 
00128   fParticleChange = GetParticleChangeForGamma();
00129 
00130   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00131   if(fAtomDeexcitation) { 
00132     fDeexcitationActive = fAtomDeexcitation->IsFluoActive(); 
00133   }
00134   if (verboseLevel > 1) {
00135     G4cout << "Livermore Polarized PhotoElectric model is initialized " 
00136            << G4endl
00137            << "Energy range: "
00138            << LowEnergyLimit() / keV << " keV - "
00139            << HighEnergyLimit() / GeV << " GeV"
00140            << G4endl;
00141   }
00142 }
00143 
00144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00145 
00146 G4double G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom(
00147                                        const G4ParticleDefinition*,
00148                                              G4double GammaEnergy,
00149                                              G4double Z, G4double,
00150                                              G4double, G4double)
00151 {
00152   G4double cs = crossSectionHandler->FindValue(G4lrint(Z), GammaEnergy);
00153 
00154   if (verboseLevel > 3) {
00155     G4cout << "G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom()" 
00156            << G4endl;
00157     G4cout << "  E(keV)= " << GammaEnergy/keV << "  Z= " << Z 
00158            << "  CrossSection(barn)= "
00159            << cs/barn << G4endl; 
00160   }
00161   return cs;
00162 }
00163 
00164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00165 
00166 void G4LivermorePolarizedPhotoElectricModel::SampleSecondaries(
00167                                 std::vector<G4DynamicParticle*>* fvect,
00168                                 const G4MaterialCutsCouple* couple,
00169                                 const G4DynamicParticle* aDynamicGamma,
00170                                 G4double,
00171                                 G4double)
00172 {
00173   if (verboseLevel > 3) {
00174     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricModel" 
00175            << G4endl;
00176   }
00177 
00178   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00179   
00180   // kill incident photon  
00181   fParticleChange->SetProposedKineticEnergy(0.);
00182   fParticleChange->ProposeTrackStatus(fStopAndKill);
00183   
00184   // low-energy gamma is absorpted by this process
00185 
00186   if (photonEnergy <= lowEnergyLimit)
00187     {
00188       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00189       return;
00190     }
00191 
00192   // Select randomly one element in the current material
00193   //G4cout << "Select random atom Egamma(keV)= " << photonEnergy/keV << G4endl;
00194   const G4Element* elm = 
00195     SelectRandomAtom(couple->GetMaterial(),theGamma,photonEnergy);
00196   G4int Z = G4lrint(elm->GetZ());
00197 
00198   // Select the ionised shell in the current atom according to shell cross sections
00199   //G4cout << "Select random shell Z= " << Z << G4endl;
00200   size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
00201   //G4cout << "Shell index= " << shellIndex 
00202   //     << " nShells= " << elm->GetNbOfAtomicShells() << G4endl;
00203   G4double bindingEnergy;
00204   const G4AtomicShell* shell = 0;
00205   if(fDeexcitationActive) {
00206     G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIndex);
00207     shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00208     bindingEnergy = shell->BindingEnergy();
00209   } else {
00210     G4int nshells = elm->GetNbOfAtomicShells() - 1;
00211     if(G4int(shellIndex) > nshells) { shellIndex = std::max(0, nshells); }
00212     bindingEnergy = elm->GetAtomicShell(shellIndex);
00213   }
00214 
00215   // There may be cases where the binding energy of the selected 
00216   // shell is > photon energy
00217   // In such cases do not generate secondaries
00218   if(photonEnergy < bindingEnergy) {
00219     fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00220     return;
00221   }
00222   //G4cout << "Z= " << Z <<  "  shellIndex= " << shellIndex 
00223   //     << " Ebind(keV)= " <<  bindingEnergy/keV << G4endl; 
00224 
00225 
00226   // Primary outcoming electron
00227   G4double eKineticEnergy = photonEnergy - bindingEnergy;
00228   G4double edep = bindingEnergy;
00229 
00230   // Calculate direction of the photoelectron
00231   G4ThreeVector electronDirection = 
00232     GetAngularDistribution()->SampleDirection(aDynamicGamma, 
00233                                               eKineticEnergy,
00234                                               shellIndex, 
00235                                               couple->GetMaterial());
00236 
00237   // The electron is created ...
00238   G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
00239                                                        electronDirection,
00240                                                        eKineticEnergy);
00241   fvect->push_back(electron);
00242 
00243   // Sample deexcitation
00244   if(fDeexcitationActive) {
00245     G4int index = couple->GetIndex();
00246     if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00247       size_t nbefore = fvect->size();
00248 
00249       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00250       size_t nafter = fvect->size();
00251       if(nafter > nbefore) {
00252         for (size_t j=nbefore; j<nafter; ++j) {
00253           edep -= ((*fvect)[j])->GetKineticEnergy();
00254         } 
00255       }
00256     }
00257   }
00258 
00259   // energy balance - excitation energy left
00260   if(edep > 0.0) {
00261     fParticleChange->ProposeLocalEnergyDeposit(edep);
00262   }
00263 }
00264 
00265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

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