physics/src/G4LowEnergyPhotoElectric.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author: A. Forti
0029 //         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0030 //
0031 // History:
0032 // --------
0033 // October 1998 - low energy modifications by Alessandra Forti
0034 // Added Livermore data table construction methods A. Forti
0035 // Modified BuildMeanFreePath to read new data tables A. Forti
0036 // Added EnergySampling method A. Forti
0037 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
0038 // Added SelectRandomAtom A. Forti
0039 // Added map of the elements A. Forti
0040 //   10.04.2000 VL
0041 // - Correcting Fluorescence transition probabilities in order to take into account
0042 //   non-radiative transitions. No Auger electron simulated yet: energy is locally deposited.
0043 // 17.02.2000 Veronique Lefebure
0044 // - bugs corrected in fluorescence simulation:
0045 //   . when final use of binding energy: no photon was ever created
0046 //   . no Fluorescence was simulated when the photo-electron energy
0047 //     was below production threshold.
0048 //
0049 // 07-09-99,  if no e- emitted: edep=photon energy, mma
0050 // 24.04.01   V.Ivanchenko     remove RogueWave
0051 // 12.08.2001 MGP              Revised according to a design iteration
0052 // 16.09.2001 E. Guardincerri  Added fluorescence generation
0053 // 06.10.2001 MGP              Added protection to avoid negative electron energies
0054 //                             when binding energy of selected shell > photon energy
0055 // 18.04.2001 V.Ivanchenko     Fix problem with low energy gammas from fluorescence
0056 //                             MeanFreePath is calculated by crosSectionHandler directly
0057 // 31.05.2002 V.Ivanchenko     Add path of Fluo + Auger cuts to AtomicDeexcitation
0058 // 14.06.2002 V.Ivanchenko     By default do not cheak range of e-
0059 // 21.01.2003 V.Ivanchenko     Cut per region
0060 // 10.05.2004 P.Rodrigues      Changes to accommodate new angular generators
0061 // 20.01.2006 A.Trindade       Changes to accommodate polarized angular generator
0062 //
0063 // --------------------------------------------------------------
0064
0065 #include "G4LowEnergyPhotoElectric.hh"
0066
0067 #include "G4RDVPhotoElectricAngularDistribution.hh"
0068 #include "G4RDPhotoElectricAngularGeneratorSimple.hh"
0069 #include "G4RDPhotoElectricAngularGeneratorSauterGavrila.hh"
0070 #include "G4RDPhotoElectricAngularGeneratorPolarized.hh"
0071
0072 #include "G4PhysicalConstants.hh"
0073 #include "G4SystemOfUnits.hh"
0074 #include "G4ParticleDefinition.hh"
0075 #include "G4Track.hh"
0076 #include "G4Step.hh"
0077 #include "G4ForceCondition.hh"
0078 #include "G4Gamma.hh"
0079 #include "G4Electron.hh"
0080 #include "G4DynamicParticle.hh"
0081 #include "G4VParticleChange.hh"
0082 #include "G4ThreeVector.hh"
0083 #include "G4RDVCrossSectionHandler.hh"
0084 #include "G4RDCrossSectionHandler.hh"
0085 #include "G4RDVEMDataSet.hh"
0086 #include "G4RDCompositeEMDataSet.hh"
0087 #include "G4RDVDataSetAlgorithm.hh"
0088 #include "G4RDLogLogInterpolation.hh"
0089 #include "G4RDVRangeTest.hh"
0090 #include "G4RDRangeNoTest.hh"
0091 #include "G4RDAtomicTransitionManager.hh"
0092 #include "G4RDAtomicShell.hh"
0093 #include "G4ProductionCutsTable.hh"
0094
0095 G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric(const G4String& processName)
0096   : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
0097     intrinsicLowEnergyLimit(10*eV),
0098     intrinsicHighEnergyLimit(100*GeV),
0099     cutForLowEnergySecondaryPhotons(250.*eV),
0100     cutForLowEnergySecondaryElectrons(250.*eV)
0101 {
0102   if (lowEnergyLimit < intrinsicLowEnergyLimit ||
0103       highEnergyLimit > intrinsicHighEnergyLimit)
0104     {
0105       G4Exception("G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric()",
0106                   "OutOfRange", FatalException,
0107                   "Energy limit outside intrinsic process validity range!");
0108     }
0109
0110   crossSectionHandler = new G4RDCrossSectionHandler();
0111   shellCrossSectionHandler = new G4RDCrossSectionHandler();
0112   meanFreePathTable = 0;
0113   rangeTest = new G4RDRangeNoTest;
0114   generatorName = "geant4.6.2";
0115   ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator");              // default generator
0116
0117
0118   if (verboseLevel > 0)
0119     {
0120       G4cout << GetProcessName() << " is created " << G4endl
0121          << "Energy range: "
0122          << lowEnergyLimit / keV << " keV - "
0123          << highEnergyLimit / GeV << " GeV"
0124          << G4endl;
0125     }
0126 }
0127
0128 G4LowEnergyPhotoElectric::~G4LowEnergyPhotoElectric()
0129 {
0130   delete crossSectionHandler;
0131   delete shellCrossSectionHandler;
0132   delete meanFreePathTable;
0133   delete rangeTest;
0134   delete ElectronAngularGenerator;
0135 }
0136
0137 void G4LowEnergyPhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
0138 {
0139
0140   crossSectionHandler->Clear();
0141   G4String crossSectionFile = "phot/pe-cs-";
0142   crossSectionHandler->LoadData(crossSectionFile);
0143
0144   shellCrossSectionHandler->Clear();
0145   G4String shellCrossSectionFile = "phot/pe-ss-cs-";
0146   shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
0147
0148   delete meanFreePathTable;
0149   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
0150 }
0151
0152 G4VParticleChange* G4LowEnergyPhotoElectric::PostStepDoIt(const G4Track& aTrack,
0153                               const G4Step& aStep)
0154 {
0155   // Fluorescence generated according to:
0156   // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
0157   // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
0158   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
0159
0160   aParticleChange.Initialize(aTrack);
0161
0162   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
0163   G4double photonEnergy = incidentPhoton->GetKineticEnergy();
0164   if (photonEnergy <= lowEnergyLimit)
0165     {
0166       aParticleChange.ProposeTrackStatus(fStopAndKill);
0167       aParticleChange.ProposeEnergy(0.);
0168       aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
0169       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
0170     }
0171
0172   G4ThreeVector photonDirection = incidentPhoton->GetMomentumDirection(); // Returns the normalized direction of the momentum
0173
0174   // Select randomly one element in the current material
0175   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0176   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
0177
0178   // Select the ionised shell in the current atom according to shell cross sections
0179   size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
0180
0181   // Retrieve the corresponding identifier and binding energy of the selected shell
0182   const G4RDAtomicTransitionManager* transitionManager = G4RDAtomicTransitionManager::Instance();
0183   const G4RDAtomicShell* shell = transitionManager->Shell(Z,shellIndex);
0184   G4double bindingEnergy = shell->BindingEnergy();
0185   G4int shellId = shell->ShellId();
0186
0187   // Create lists of pointers to DynamicParticles (photons and electrons)
0188   // (Is the electron vector necessary? To be checked)
0189   std::vector<G4DynamicParticle*>* photonVector = 0;
0190   std::vector<G4DynamicParticle*> electronVector;
0191
0192   G4double energyDeposit = 0.0;
0193
0194   // Primary outcoming electron
0195   G4double eKineticEnergy = photonEnergy - bindingEnergy;
0196
0197   // There may be cases where the binding energy of the selected shell is > photon energy
0198   // In such cases do not generate secondaries
0199   if (eKineticEnergy > 0.)
0200     {
0201       // Generate the electron only if with large enough range w.r.t. cuts and safety
0202       G4double safety = aStep.GetPostStepPoint()->GetSafety();
0203
0204       if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
0205     {
0206
0207       // Calculate direction of the photoelectron
0208       G4ThreeVector gammaPolarization = incidentPhoton->GetPolarization();
0209       G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
0210
0211       // The electron is created ...
0212       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
0213                                    electronDirection,
0214                                    eKineticEnergy);
0215       electronVector.push_back(electron);
0216     }
0217       else
0218     {
0219       energyDeposit += eKineticEnergy;
0220     }
0221     }
0222   else
0223     {
0224       bindingEnergy = photonEnergy;
0225     }
0226
0227   G4int nElectrons = electronVector.size();
0228   size_t nTotPhotons = 0;
0229   G4int nPhotons=0;
0230   const G4ProductionCutsTable* theCoupleTable=
0231         G4ProductionCutsTable::GetProductionCutsTable();
0232
0233   size_t index = couple->GetIndex();
0234   G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
0235   cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
0236
0237   G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
0238   cute = std::min(cutForLowEnergySecondaryPhotons,cute);
0239
0240   G4DynamicParticle* aPhoton;
0241
0242   // Generation of fluorescence
0243   // Data in EADL are available only for Z > 5
0244   // Protection to avoid generating photons in the unphysical case of
0245   // shell binding energy > photon energy
0246   if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
0247     {
0248       photonVector = deexcitationManager.GenerateParticles(Z,shellId);
0249       nTotPhotons = photonVector->size();
0250       for (size_t k=0; k<nTotPhotons; k++)
0251     {
0252       aPhoton = (*photonVector)[k];
0253       if (aPhoton)
0254         {
0255               G4double itsCut = cutg;
0256               if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
0257           G4double itsEnergy = aPhoton->GetKineticEnergy();
0258
0259           if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
0260         {
0261           nPhotons++;
0262           // Local energy deposit is given as the sum of the
0263           // energies of incident photons minus the energies
0264           // of the outcoming fluorescence photons
0265           bindingEnergy -= itsEnergy;
0266
0267         }
0268           else
0269         {
0270                   delete aPhoton;
0271                   (*photonVector)[k] = 0;
0272                 }
0273         }
0274     }
0275     }
0276
0277   energyDeposit += bindingEnergy;
0278
0279   G4int nSecondaries  = nElectrons + nPhotons;
0280   aParticleChange.SetNumberOfSecondaries(nSecondaries);
0281
0282   for (G4int l = 0; l<nElectrons; l++ )
0283     {
0284       aPhoton = electronVector[l];
0285       if(aPhoton) {
0286         aParticleChange.AddSecondary(aPhoton);
0287       }
0288     }
0289   for ( size_t ll = 0; ll < nTotPhotons; ll++)
0290     {
0291       aPhoton = (*photonVector)[ll];
0292       if(aPhoton) {
0293         aParticleChange.AddSecondary(aPhoton);
0294       }
0295     }
0296
0297   delete photonVector;
0298
0299   if (energyDeposit < 0)
0300     {
0301       G4cout << "WARNING - "
0302          << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
0303          << G4endl;
0304       energyDeposit = 0;
0305     }
0306
0307   // Kill the incident photon
0308   aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
0309   aParticleChange.ProposeEnergy( 0. );
0310
0311   aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
0312   aParticleChange.ProposeTrackStatus( fStopAndKill );
0313
0314   // Reset NbOfInteractionLengthLeft and return aParticleChange
0315   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
0316 }
0317
0318 G4bool G4LowEnergyPhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
0319 {
0320   return ( &particle == G4Gamma::Gamma() );
0321 }
0322
0323 G4double G4LowEnergyPhotoElectric::GetMeanFreePath(const G4Track& track,
0324                            G4double, // previousStepSize
0325                            G4ForceCondition*)
0326 {
0327   const G4DynamicParticle* photon = track.GetDynamicParticle();
0328   G4double energy = photon->GetKineticEnergy();
0329   G4Material* material = track.GetMaterial();
0330   //  size_t materialIndex = material->GetIndex();
0331
0332   G4double meanFreePath = DBL_MAX;
0333
0334   //  if (energy > highEnergyLimit)
0335   //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
0336   //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
0337   //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
0338
0339   G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
0340   if(cross > 0.0) meanFreePath = 1.0/cross;
0341
0342   return meanFreePath;
0343 }
0344
0345 void G4LowEnergyPhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
0346 {
0347   cutForLowEnergySecondaryPhotons = cut;
0348   deexcitationManager.SetCutForSecondaryPhotons(cut);
0349 }
0350
0351 void G4LowEnergyPhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
0352 {
0353   cutForLowEnergySecondaryElectrons = cut;
0354   deexcitationManager.SetCutForAugerElectrons(cut);
0355 }
0356
0357 void G4LowEnergyPhotoElectric::ActivateAuger(G4bool val)
0358 {
0359   deexcitationManager.ActivateAugerElectronProduction(val);
0360 }
0361
0362 void G4LowEnergyPhotoElectric::SetAngularGenerator(G4RDVPhotoElectricAngularDistribution* distribution)
0363 {
0364   ElectronAngularGenerator = distribution;
0365   ElectronAngularGenerator->PrintGeneratorInformation();
0366 }
0367
0368 void G4LowEnergyPhotoElectric::SetAngularGenerator(const G4String& name)
0369 {
0370   if (name == "default")
0371     {
0372       delete ElectronAngularGenerator;
0373       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
0374       generatorName = name;
0375     }
0376   else if (name == "standard")
0377     {
0378       delete ElectronAngularGenerator;
0379       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
0380       generatorName = name;
0381     }
0382   else if (name == "polarized")
0383     {
0384       delete ElectronAngularGenerator;
0385       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
0386       generatorName = name;
0387     }
0388   else
0389     {
0390       G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator()",
0391                   "InvalidSetup", FatalException,
0392                   "Generator does not exist!");
0393     }
0394
0395   ElectronAngularGenerator->PrintGeneratorInformation();
0396 }