Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4LowEnergyPhotoElectric.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/eRosita/physics/src/G4LowEnergyPhotoElectric.cc (Version 11.3.0) and /examples/advanced/eRosita/physics/src/G4LowEnergyPhotoElectric.cc (Version 9.5.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
                                                   >>  27 // $Id: G4LowEnergyPhotoElectric.cc,v 1.56 2006/06/29 19:40:23 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
 27 //                                                 29 //
 28 // Author: A. Forti                                30 // Author: A. Forti
 29 //         Maria Grazia Pia (Maria.Grazia.Pia@     31 //         Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
 30 //                                                 32 //
 31 // History:                                        33 // History:
 32 // --------                                        34 // --------
 33 // October 1998 - low energy modifications by      35 // October 1998 - low energy modifications by Alessandra Forti
 34 // Added Livermore data table construction met     36 // Added Livermore data table construction methods A. Forti
 35 // Modified BuildMeanFreePath to read new data     37 // Modified BuildMeanFreePath to read new data tables A. Forti
 36 // Added EnergySampling method A. Forti            38 // Added EnergySampling method A. Forti
 37 // Modified PostStepDoIt to insert sampling wi     39 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
 38 // Added SelectRandomAtom A. Forti                 40 // Added SelectRandomAtom A. Forti
 39 // Added map of the elements A. Forti              41 // Added map of the elements A. Forti
 40 //   10.04.2000 VL                                 42 //   10.04.2000 VL
 41 // - Correcting Fluorescence transition probab     43 // - Correcting Fluorescence transition probabilities in order to take into account 
 42 //   non-radiative transitions. No Auger elect     44 //   non-radiative transitions. No Auger electron simulated yet: energy is locally deposited.
 43 // 17.02.2000 Veronique Lefebure                   45 // 17.02.2000 Veronique Lefebure
 44 // - bugs corrected in fluorescence simulation     46 // - bugs corrected in fluorescence simulation: 
 45 //   . when final use of binding energy: no ph     47 //   . when final use of binding energy: no photon was ever created
 46 //   . no Fluorescence was simulated when the      48 //   . no Fluorescence was simulated when the photo-electron energy
 47 //     was below production threshold.             49 //     was below production threshold.
 48 //                                                 50 //
 49 // 07-09-99,  if no e- emitted: edep=photon en     51 // 07-09-99,  if no e- emitted: edep=photon energy, mma
 50 // 24.04.01   V.Ivanchenko     remove RogueWav     52 // 24.04.01   V.Ivanchenko     remove RogueWave 
 51 // 12.08.2001 MGP              Revised accordi     53 // 12.08.2001 MGP              Revised according to a design iteration
 52 // 16.09.2001 E. Guardincerri  Added fluoresce     54 // 16.09.2001 E. Guardincerri  Added fluorescence generation
 53 // 06.10.2001 MGP              Added protectio     55 // 06.10.2001 MGP              Added protection to avoid negative electron energies
 54 //                             when binding en     56 //                             when binding energy of selected shell > photon energy
 55 // 18.04.2001 V.Ivanchenko     Fix problem wit     57 // 18.04.2001 V.Ivanchenko     Fix problem with low energy gammas from fluorescence
 56 //                             MeanFreePath is     58 //                             MeanFreePath is calculated by crosSectionHandler directly
 57 // 31.05.2002 V.Ivanchenko     Add path of Flu     59 // 31.05.2002 V.Ivanchenko     Add path of Fluo + Auger cuts to AtomicDeexcitation
 58 // 14.06.2002 V.Ivanchenko     By default do n     60 // 14.06.2002 V.Ivanchenko     By default do not cheak range of e-
 59 // 21.01.2003 V.Ivanchenko     Cut per region      61 // 21.01.2003 V.Ivanchenko     Cut per region
 60 // 10.05.2004 P.Rodrigues      Changes to acco     62 // 10.05.2004 P.Rodrigues      Changes to accommodate new angular generators
 61 // 20.01.2006 A.Trindade       Changes to acco     63 // 20.01.2006 A.Trindade       Changes to accommodate polarized angular generator
 62 //                                                 64 //
 63 // -------------------------------------------     65 // --------------------------------------------------------------
 64                                                    66 
 65 #include "G4LowEnergyPhotoElectric.hh"             67 #include "G4LowEnergyPhotoElectric.hh"
 66                                                    68 
 67 #include "G4RDVPhotoElectricAngularDistributio     69 #include "G4RDVPhotoElectricAngularDistribution.hh"
 68 #include "G4RDPhotoElectricAngularGeneratorSim     70 #include "G4RDPhotoElectricAngularGeneratorSimple.hh"
 69 #include "G4RDPhotoElectricAngularGeneratorSau     71 #include "G4RDPhotoElectricAngularGeneratorSauterGavrila.hh"
 70 #include "G4RDPhotoElectricAngularGeneratorPol     72 #include "G4RDPhotoElectricAngularGeneratorPolarized.hh"
 71                                                    73 
 72 #include "G4PhysicalConstants.hh"              << 
 73 #include "G4SystemOfUnits.hh"                  << 
 74 #include "G4ParticleDefinition.hh"                 74 #include "G4ParticleDefinition.hh"
 75 #include "G4Track.hh"                              75 #include "G4Track.hh"
 76 #include "G4Step.hh"                               76 #include "G4Step.hh"
 77 #include "G4ForceCondition.hh"                     77 #include "G4ForceCondition.hh"
 78 #include "G4Gamma.hh"                              78 #include "G4Gamma.hh"
 79 #include "G4Electron.hh"                           79 #include "G4Electron.hh"
 80 #include "G4DynamicParticle.hh"                    80 #include "G4DynamicParticle.hh"
 81 #include "G4VParticleChange.hh"                    81 #include "G4VParticleChange.hh"
 82 #include "G4ThreeVector.hh"                        82 #include "G4ThreeVector.hh"
 83 #include "G4RDVCrossSectionHandler.hh"             83 #include "G4RDVCrossSectionHandler.hh"
 84 #include "G4RDCrossSectionHandler.hh"              84 #include "G4RDCrossSectionHandler.hh"
 85 #include "G4RDVEMDataSet.hh"                       85 #include "G4RDVEMDataSet.hh"
 86 #include "G4RDCompositeEMDataSet.hh"               86 #include "G4RDCompositeEMDataSet.hh"
 87 #include "G4RDVDataSetAlgorithm.hh"                87 #include "G4RDVDataSetAlgorithm.hh"
 88 #include "G4RDLogLogInterpolation.hh"              88 #include "G4RDLogLogInterpolation.hh"
 89 #include "G4RDVRangeTest.hh"                       89 #include "G4RDVRangeTest.hh"
 90 #include "G4RDRangeNoTest.hh"                      90 #include "G4RDRangeNoTest.hh"
 91 #include "G4RDAtomicTransitionManager.hh"          91 #include "G4RDAtomicTransitionManager.hh"
 92 #include "G4RDAtomicShell.hh"                      92 #include "G4RDAtomicShell.hh"
 93 #include "G4ProductionCutsTable.hh"                93 #include "G4ProductionCutsTable.hh"
 94                                                    94 
 95 G4LowEnergyPhotoElectric::G4LowEnergyPhotoElec     95 G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric(const G4String& processName)
 96   : G4VDiscreteProcess(processName), lowEnergy     96   : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
 97     intrinsicLowEnergyLimit(10*eV),                97     intrinsicLowEnergyLimit(10*eV),
 98     intrinsicHighEnergyLimit(100*GeV),             98     intrinsicHighEnergyLimit(100*GeV),
 99     cutForLowEnergySecondaryPhotons(250.*eV),      99     cutForLowEnergySecondaryPhotons(250.*eV),
100     cutForLowEnergySecondaryElectrons(250.*eV)    100     cutForLowEnergySecondaryElectrons(250.*eV)
101 {                                                 101 {
102   if (lowEnergyLimit < intrinsicLowEnergyLimit    102   if (lowEnergyLimit < intrinsicLowEnergyLimit || 
103       highEnergyLimit > intrinsicHighEnergyLim    103       highEnergyLimit > intrinsicHighEnergyLimit)
104     {                                             104     {
105       G4Exception("G4LowEnergyPhotoElectric::G    105       G4Exception("G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric()",
106                   "OutOfRange", FatalException    106                   "OutOfRange", FatalException,
107                   "Energy limit outside intrin    107                   "Energy limit outside intrinsic process validity range!");
108     }                                             108     }
109                                                   109 
110   crossSectionHandler = new G4RDCrossSectionHa    110   crossSectionHandler = new G4RDCrossSectionHandler();
111   shellCrossSectionHandler = new G4RDCrossSect    111   shellCrossSectionHandler = new G4RDCrossSectionHandler();
112   meanFreePathTable = 0;                          112   meanFreePathTable = 0;
113   rangeTest = new G4RDRangeNoTest;                113   rangeTest = new G4RDRangeNoTest;
114   generatorName = "geant4.6.2";                   114   generatorName = "geant4.6.2";
115   ElectronAngularGenerator = new G4RDPhotoElec    115   ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator");              // default generator
116                                                   116 
117                                                   117 
118   if (verboseLevel > 0)                           118   if (verboseLevel > 0)
119     {                                             119     {
120       G4cout << GetProcessName() << " is creat    120       G4cout << GetProcessName() << " is created " << G4endl
121        << "Energy range: "                        121        << "Energy range: "
122        << lowEnergyLimit / keV << " keV - "       122        << lowEnergyLimit / keV << " keV - "
123        << highEnergyLimit / GeV << " GeV"         123        << highEnergyLimit / GeV << " GeV"
124        << G4endl;                                 124        << G4endl;
125     }                                             125     }
126 }                                                 126 }
127                                                   127 
128 G4LowEnergyPhotoElectric::~G4LowEnergyPhotoEle    128 G4LowEnergyPhotoElectric::~G4LowEnergyPhotoElectric()
129 {                                                 129 {
130   delete crossSectionHandler;                     130   delete crossSectionHandler;
131   delete shellCrossSectionHandler;                131   delete shellCrossSectionHandler;
132   delete meanFreePathTable;                       132   delete meanFreePathTable;
133   delete rangeTest;                               133   delete rangeTest;
134   delete ElectronAngularGenerator;                134   delete ElectronAngularGenerator;
135 }                                                 135 }
136                                                   136 
137 void G4LowEnergyPhotoElectric::BuildPhysicsTab    137 void G4LowEnergyPhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
138 {                                                 138 {
139                                                   139 
140   crossSectionHandler->Clear();                   140   crossSectionHandler->Clear();
141   G4String crossSectionFile = "phot/pe-cs-";      141   G4String crossSectionFile = "phot/pe-cs-";
142   crossSectionHandler->LoadData(crossSectionFi    142   crossSectionHandler->LoadData(crossSectionFile);
143                                                   143 
144   shellCrossSectionHandler->Clear();              144   shellCrossSectionHandler->Clear();
145   G4String shellCrossSectionFile = "phot/pe-ss    145   G4String shellCrossSectionFile = "phot/pe-ss-cs-";
146   shellCrossSectionHandler->LoadShellData(shel    146   shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
147                                                   147 
148   delete meanFreePathTable;                       148   delete meanFreePathTable;
149   meanFreePathTable = crossSectionHandler->Bui    149   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
150 }                                                 150 }
151                                                   151 
152 G4VParticleChange* G4LowEnergyPhotoElectric::P    152 G4VParticleChange* G4LowEnergyPhotoElectric::PostStepDoIt(const G4Track& aTrack,
153                 const G4Step& aStep)              153                 const G4Step& aStep)
154 {                                                 154 {
155   // Fluorescence generated according to:         155   // Fluorescence generated according to:
156   // J. Stepanek ,"A program to determine the     156   // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
157   // subshell ionisation by a particle or due     157   // subshell ionisation by a particle or due to deexcitation or decay of radionuclides", 
158   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)       158   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
159                                                   159  
160   aParticleChange.Initialize(aTrack);             160   aParticleChange.Initialize(aTrack);
161                                                   161   
162   const G4DynamicParticle* incidentPhoton = aT    162   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
163   G4double photonEnergy = incidentPhoton->GetK    163   G4double photonEnergy = incidentPhoton->GetKineticEnergy();
164   if (photonEnergy <= lowEnergyLimit)             164   if (photonEnergy <= lowEnergyLimit)
165     {                                             165     {
166       aParticleChange.ProposeTrackStatus(fStop    166       aParticleChange.ProposeTrackStatus(fStopAndKill);
167       aParticleChange.ProposeEnergy(0.);          167       aParticleChange.ProposeEnergy(0.);
168       aParticleChange.ProposeLocalEnergyDeposi    168       aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
169       return G4VDiscreteProcess::PostStepDoIt(    169       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
170     }                                             170     }
171                                                   171  
172   G4ThreeVector photonDirection = incidentPhot    172   G4ThreeVector photonDirection = incidentPhoton->GetMomentumDirection(); // Returns the normalized direction of the momentum
173                                                   173 
174   // Select randomly one element in the curren    174   // Select randomly one element in the current material
175   const G4MaterialCutsCouple* couple = aTrack.    175   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
176   G4int Z = crossSectionHandler->SelectRandomA    176   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
177                                                   177 
178   // Select the ionised shell in the current a    178   // Select the ionised shell in the current atom according to shell cross sections
179   size_t shellIndex = shellCrossSectionHandler    179   size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
180                                                   180 
181   // Retrieve the corresponding identifier and    181   // Retrieve the corresponding identifier and binding energy of the selected shell
182   const G4RDAtomicTransitionManager* transitio    182   const G4RDAtomicTransitionManager* transitionManager = G4RDAtomicTransitionManager::Instance();
183   const G4RDAtomicShell* shell = transitionMan    183   const G4RDAtomicShell* shell = transitionManager->Shell(Z,shellIndex);
184   G4double bindingEnergy = shell->BindingEnerg    184   G4double bindingEnergy = shell->BindingEnergy();
185   G4int shellId = shell->ShellId();               185   G4int shellId = shell->ShellId();
186                                                   186 
187   // Create lists of pointers to DynamicPartic    187   // Create lists of pointers to DynamicParticles (photons and electrons)
188   // (Is the electron vector necessary? To be     188   // (Is the electron vector necessary? To be checked)
189   std::vector<G4DynamicParticle*>* photonVecto    189   std::vector<G4DynamicParticle*>* photonVector = 0;
190   std::vector<G4DynamicParticle*> electronVect    190   std::vector<G4DynamicParticle*> electronVector;
191                                                   191 
192   G4double energyDeposit = 0.0;                   192   G4double energyDeposit = 0.0;
193                                                   193 
194   // Primary outcoming electron                   194   // Primary outcoming electron
195   G4double eKineticEnergy = photonEnergy - bin    195   G4double eKineticEnergy = photonEnergy - bindingEnergy;
196                                                   196 
197   // There may be cases where the binding ener    197   // There may be cases where the binding energy of the selected shell is > photon energy
198   // In such cases do not generate secondaries    198   // In such cases do not generate secondaries
199   if (eKineticEnergy > 0.)                        199   if (eKineticEnergy > 0.)
200     {                                             200     {
201       // Generate the electron only if with la    201       // Generate the electron only if with large enough range w.r.t. cuts and safety
202       G4double safety = aStep.GetPostStepPoint    202       G4double safety = aStep.GetPostStepPoint()->GetSafety();
203                                                   203 
204       if (rangeTest->Escape(G4Electron::Electr    204       if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
205   {                                               205   {
206                                                   206 
207     // Calculate direction of the photoelectro    207     // Calculate direction of the photoelectron
208     G4ThreeVector gammaPolarization = incident    208     G4ThreeVector gammaPolarization = incidentPhoton->GetPolarization();
209     G4ThreeVector electronDirection = Electron    209     G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
210                                                   210 
211     // The electron is created ...                211     // The electron is created ...
212     G4DynamicParticle* electron = new G4Dynami    212     G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
213                      electronDirection,           213                      electronDirection,
214                      eKineticEnergy);             214                      eKineticEnergy);
215     electronVector.push_back(electron);           215     electronVector.push_back(electron);
216   }                                               216   }
217       else                                        217       else
218   {                                               218   {
219     energyDeposit += eKineticEnergy;              219     energyDeposit += eKineticEnergy;
220   }                                               220   }
221     }                                             221     }
222   else                                            222   else
223     {                                             223     {
224       bindingEnergy = photonEnergy;               224       bindingEnergy = photonEnergy;
225     }                                             225     }
226                                                   226 
227   G4int nElectrons = electronVector.size();       227   G4int nElectrons = electronVector.size();
228   size_t nTotPhotons = 0;                         228   size_t nTotPhotons = 0;
229   G4int nPhotons=0;                               229   G4int nPhotons=0;
230   const G4ProductionCutsTable* theCoupleTable=    230   const G4ProductionCutsTable* theCoupleTable=
231         G4ProductionCutsTable::GetProductionCu    231         G4ProductionCutsTable::GetProductionCutsTable();
232                                                   232 
233   size_t index = couple->GetIndex();              233   size_t index = couple->GetIndex();
234   G4double cutg = (*(theCoupleTable->GetEnergy    234   G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
235   cutg = std::min(cutForLowEnergySecondaryPhot    235   cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
236                                                   236   
237   G4double cute = (*(theCoupleTable->GetEnergy    237   G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
238   cute = std::min(cutForLowEnergySecondaryPhot    238   cute = std::min(cutForLowEnergySecondaryPhotons,cute);
239                                                   239 
240   G4DynamicParticle* aPhoton;                     240   G4DynamicParticle* aPhoton;
241                                                   241 
242   // Generation of fluorescence                   242   // Generation of fluorescence
243   // Data in EADL are available only for Z > 5    243   // Data in EADL are available only for Z > 5
244   // Protection to avoid generating photons in    244   // Protection to avoid generating photons in the unphysical case of
245   // shell binding energy > photon energy         245   // shell binding energy > photon energy
246   if (Z > 5  && (bindingEnergy > cutg || bindi    246   if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
247     {                                             247     {
248       photonVector = deexcitationManager.Gener    248       photonVector = deexcitationManager.GenerateParticles(Z,shellId);
249       nTotPhotons = photonVector->size();         249       nTotPhotons = photonVector->size();
250       for (size_t k=0; k<nTotPhotons; k++)        250       for (size_t k=0; k<nTotPhotons; k++)
251   {                                               251   {
252     aPhoton = (*photonVector)[k];                 252     aPhoton = (*photonVector)[k];
253     if (aPhoton)                                  253     if (aPhoton)
254       {                                           254       {
255               G4double itsCut = cutg;             255               G4double itsCut = cutg;
256               if(aPhoton->GetDefinition() == G    256               if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
257         G4double itsEnergy = aPhoton->GetKinet    257         G4double itsEnergy = aPhoton->GetKineticEnergy();
258                                                   258 
259         if (itsEnergy > itsCut && itsEnergy <=    259         if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
260     {                                             260     {
261       nPhotons++;                                 261       nPhotons++;
262       // Local energy deposit is given as the     262       // Local energy deposit is given as the sum of the
263       // energies of incident photons minus th    263       // energies of incident photons minus the energies
264       // of the outcoming fluorescence photons    264       // of the outcoming fluorescence photons
265       bindingEnergy -= itsEnergy;                 265       bindingEnergy -= itsEnergy;
266                                                   266 
267     }                                             267     }
268         else                                      268         else
269     {                                             269     {
270                   delete aPhoton;                 270                   delete aPhoton;
271                   (*photonVector)[k] = 0;         271                   (*photonVector)[k] = 0;
272                 }                                 272                 }
273       }                                           273       }
274   }                                               274   }
275     }                                             275     }
276                                                   276 
277   energyDeposit += bindingEnergy;                 277   energyDeposit += bindingEnergy;
278                                                   278 
279   G4int nSecondaries  = nElectrons + nPhotons;    279   G4int nSecondaries  = nElectrons + nPhotons;
280   aParticleChange.SetNumberOfSecondaries(nSeco    280   aParticleChange.SetNumberOfSecondaries(nSecondaries);
281                                                   281 
282   for (G4int l = 0; l<nElectrons; l++ )           282   for (G4int l = 0; l<nElectrons; l++ )
283     {                                             283     {
284       aPhoton = electronVector[l];                284       aPhoton = electronVector[l];
285       if(aPhoton) {                               285       if(aPhoton) {
286         aParticleChange.AddSecondary(aPhoton);    286         aParticleChange.AddSecondary(aPhoton);
287       }                                           287       }
288     }                                             288     }
289   for ( size_t ll = 0; ll < nTotPhotons; ll++)    289   for ( size_t ll = 0; ll < nTotPhotons; ll++)
290     {                                             290     {
291       aPhoton = (*photonVector)[ll];              291       aPhoton = (*photonVector)[ll];
292       if(aPhoton) {                               292       if(aPhoton) {
293         aParticleChange.AddSecondary(aPhoton);    293         aParticleChange.AddSecondary(aPhoton);
294       }                                           294       }
295     }                                             295     }
296                                                   296 
297   delete photonVector;                            297   delete photonVector;
298                                                   298 
299   if (energyDeposit < 0)                          299   if (energyDeposit < 0)
300     {                                             300     {
301       G4cout << "WARNING - "                      301       G4cout << "WARNING - "
302        << "G4LowEnergyPhotoElectric::PostStepD    302        << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
303        << G4endl;                                 303        << G4endl;
304       energyDeposit = 0;                          304       energyDeposit = 0;
305     }                                             305     }
306                                                   306 
307   // Kill the incident photon                     307   // Kill the incident photon
308   aParticleChange.ProposeMomentumDirection( 0.    308   aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
309   aParticleChange.ProposeEnergy( 0. );            309   aParticleChange.ProposeEnergy( 0. );
310                                                   310 
311   aParticleChange.ProposeLocalEnergyDeposit(en    311   aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
312   aParticleChange.ProposeTrackStatus( fStopAnd    312   aParticleChange.ProposeTrackStatus( fStopAndKill );
313                                                   313 
314   // Reset NbOfInteractionLengthLeft and retur    314   // Reset NbOfInteractionLengthLeft and return aParticleChange
315   return G4VDiscreteProcess::PostStepDoIt( aTr    315   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
316 }                                                 316 }
317                                                   317 
318 G4bool G4LowEnergyPhotoElectric::IsApplicable(    318 G4bool G4LowEnergyPhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
319 {                                                 319 {
320   return ( &particle == G4Gamma::Gamma() );       320   return ( &particle == G4Gamma::Gamma() );
321 }                                                 321 }
322                                                   322 
323 G4double G4LowEnergyPhotoElectric::GetMeanFree    323 G4double G4LowEnergyPhotoElectric::GetMeanFreePath(const G4Track& track,
324                G4double, // previousStepSize      324                G4double, // previousStepSize
325                  G4ForceCondition*)               325                  G4ForceCondition*)
326 {                                                 326 {
327   const G4DynamicParticle* photon = track.GetD    327   const G4DynamicParticle* photon = track.GetDynamicParticle();
328   G4double energy = photon->GetKineticEnergy()    328   G4double energy = photon->GetKineticEnergy();
329   G4Material* material = track.GetMaterial();     329   G4Material* material = track.GetMaterial();
330   //  size_t materialIndex = material->GetInde    330   //  size_t materialIndex = material->GetIndex();
331                                                   331 
332   G4double meanFreePath = DBL_MAX;                332   G4double meanFreePath = DBL_MAX;
333                                                   333 
334   //  if (energy > highEnergyLimit)               334   //  if (energy > highEnergyLimit) 
335   //    meanFreePath = meanFreePathTable->Find    335   //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
336   //  else if (energy < lowEnergyLimit) meanFr    336   //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
337   //  else meanFreePath = meanFreePathTable->F    337   //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
338                                                   338 
339   G4double cross = shellCrossSectionHandler->V    339   G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
340   if(cross > 0.0) meanFreePath = 1.0/cross;       340   if(cross > 0.0) meanFreePath = 1.0/cross;
341                                                   341 
342   return meanFreePath;                            342   return meanFreePath;
343 }                                                 343 }
344                                                   344 
345 void G4LowEnergyPhotoElectric::SetCutForLowEnS    345 void G4LowEnergyPhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
346 {                                                 346 {
347   cutForLowEnergySecondaryPhotons = cut;          347   cutForLowEnergySecondaryPhotons = cut;
348   deexcitationManager.SetCutForSecondaryPhoton    348   deexcitationManager.SetCutForSecondaryPhotons(cut);
349 }                                                 349 }
350                                                   350 
351 void G4LowEnergyPhotoElectric::SetCutForLowEnS    351 void G4LowEnergyPhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
352 {                                                 352 {
353   cutForLowEnergySecondaryElectrons = cut;        353   cutForLowEnergySecondaryElectrons = cut;
354   deexcitationManager.SetCutForAugerElectrons(    354   deexcitationManager.SetCutForAugerElectrons(cut);
355 }                                                 355 }
356                                                   356 
357 void G4LowEnergyPhotoElectric::ActivateAuger(G    357 void G4LowEnergyPhotoElectric::ActivateAuger(G4bool val)
358 {                                                 358 {
359   deexcitationManager.ActivateAugerElectronPro    359   deexcitationManager.ActivateAugerElectronProduction(val);
360 }                                                 360 }
361                                                   361 
362 void G4LowEnergyPhotoElectric::SetAngularGener    362 void G4LowEnergyPhotoElectric::SetAngularGenerator(G4RDVPhotoElectricAngularDistribution* distribution)
363 {                                                 363 {
364   ElectronAngularGenerator = distribution;        364   ElectronAngularGenerator = distribution;
365   ElectronAngularGenerator->PrintGeneratorInfo    365   ElectronAngularGenerator->PrintGeneratorInformation();
366 }                                                 366 }
367                                                   367 
368 void G4LowEnergyPhotoElectric::SetAngularGener    368 void G4LowEnergyPhotoElectric::SetAngularGenerator(const G4String& name)
369 {                                                 369 {
370   if (name == "default")                          370   if (name == "default") 
371     {                                             371     {
372       delete ElectronAngularGenerator;            372       delete ElectronAngularGenerator;
373       ElectronAngularGenerator = new G4RDPhoto    373       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
374       generatorName = name;                       374       generatorName = name;
375     }                                             375     }
376   else if (name == "standard")                    376   else if (name == "standard")
377     {                                             377     {
378       delete ElectronAngularGenerator;            378       delete ElectronAngularGenerator;
379       ElectronAngularGenerator = new G4RDPhoto    379       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
380       generatorName = name;                       380       generatorName = name;
381     }                                             381     }
382   else if (name == "polarized")                   382   else if (name == "polarized")
383     {                                             383     {
384       delete ElectronAngularGenerator;            384       delete ElectronAngularGenerator;
385       ElectronAngularGenerator = new G4RDPhoto    385       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
386       generatorName = name;                       386       generatorName = name;
387     }                                             387     }
388   else                                            388   else
389     {                                             389     {
390       G4Exception("G4LowEnergyPhotoElectric::S    390       G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator()",
391                   "InvalidSetup", FatalExcepti    391                   "InvalidSetup", FatalException,
392                   "Generator does not exist!")    392                   "Generator does not exist!");
393     }                                             393     }
394                                                   394 
395   ElectronAngularGenerator->PrintGeneratorInfo    395   ElectronAngularGenerator->PrintGeneratorInformation();
396 }                                                 396 }
397                                                   397