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.6.p1)


  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$
                                                   >>  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"                  74 #include "G4PhysicalConstants.hh"
 73 #include "G4SystemOfUnits.hh"                      75 #include "G4SystemOfUnits.hh"
 74 #include "G4ParticleDefinition.hh"                 76 #include "G4ParticleDefinition.hh"
 75 #include "G4Track.hh"                              77 #include "G4Track.hh"
 76 #include "G4Step.hh"                               78 #include "G4Step.hh"
 77 #include "G4ForceCondition.hh"                     79 #include "G4ForceCondition.hh"
 78 #include "G4Gamma.hh"                              80 #include "G4Gamma.hh"
 79 #include "G4Electron.hh"                           81 #include "G4Electron.hh"
 80 #include "G4DynamicParticle.hh"                    82 #include "G4DynamicParticle.hh"
 81 #include "G4VParticleChange.hh"                    83 #include "G4VParticleChange.hh"
 82 #include "G4ThreeVector.hh"                        84 #include "G4ThreeVector.hh"
 83 #include "G4RDVCrossSectionHandler.hh"             85 #include "G4RDVCrossSectionHandler.hh"
 84 #include "G4RDCrossSectionHandler.hh"              86 #include "G4RDCrossSectionHandler.hh"
 85 #include "G4RDVEMDataSet.hh"                       87 #include "G4RDVEMDataSet.hh"
 86 #include "G4RDCompositeEMDataSet.hh"               88 #include "G4RDCompositeEMDataSet.hh"
 87 #include "G4RDVDataSetAlgorithm.hh"                89 #include "G4RDVDataSetAlgorithm.hh"
 88 #include "G4RDLogLogInterpolation.hh"              90 #include "G4RDLogLogInterpolation.hh"
 89 #include "G4RDVRangeTest.hh"                       91 #include "G4RDVRangeTest.hh"
 90 #include "G4RDRangeNoTest.hh"                      92 #include "G4RDRangeNoTest.hh"
 91 #include "G4RDAtomicTransitionManager.hh"          93 #include "G4RDAtomicTransitionManager.hh"
 92 #include "G4RDAtomicShell.hh"                      94 #include "G4RDAtomicShell.hh"
 93 #include "G4ProductionCutsTable.hh"                95 #include "G4ProductionCutsTable.hh"
 94                                                    96 
 95 G4LowEnergyPhotoElectric::G4LowEnergyPhotoElec     97 G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric(const G4String& processName)
 96   : G4VDiscreteProcess(processName), lowEnergy     98   : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
 97     intrinsicLowEnergyLimit(10*eV),                99     intrinsicLowEnergyLimit(10*eV),
 98     intrinsicHighEnergyLimit(100*GeV),            100     intrinsicHighEnergyLimit(100*GeV),
 99     cutForLowEnergySecondaryPhotons(250.*eV),     101     cutForLowEnergySecondaryPhotons(250.*eV),
100     cutForLowEnergySecondaryElectrons(250.*eV)    102     cutForLowEnergySecondaryElectrons(250.*eV)
101 {                                                 103 {
102   if (lowEnergyLimit < intrinsicLowEnergyLimit    104   if (lowEnergyLimit < intrinsicLowEnergyLimit || 
103       highEnergyLimit > intrinsicHighEnergyLim    105       highEnergyLimit > intrinsicHighEnergyLimit)
104     {                                             106     {
105       G4Exception("G4LowEnergyPhotoElectric::G    107       G4Exception("G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric()",
106                   "OutOfRange", FatalException    108                   "OutOfRange", FatalException,
107                   "Energy limit outside intrin    109                   "Energy limit outside intrinsic process validity range!");
108     }                                             110     }
109                                                   111 
110   crossSectionHandler = new G4RDCrossSectionHa    112   crossSectionHandler = new G4RDCrossSectionHandler();
111   shellCrossSectionHandler = new G4RDCrossSect    113   shellCrossSectionHandler = new G4RDCrossSectionHandler();
112   meanFreePathTable = 0;                          114   meanFreePathTable = 0;
113   rangeTest = new G4RDRangeNoTest;                115   rangeTest = new G4RDRangeNoTest;
114   generatorName = "geant4.6.2";                   116   generatorName = "geant4.6.2";
115   ElectronAngularGenerator = new G4RDPhotoElec    117   ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator");              // default generator
116                                                   118 
117                                                   119 
118   if (verboseLevel > 0)                           120   if (verboseLevel > 0)
119     {                                             121     {
120       G4cout << GetProcessName() << " is creat    122       G4cout << GetProcessName() << " is created " << G4endl
121        << "Energy range: "                        123        << "Energy range: "
122        << lowEnergyLimit / keV << " keV - "       124        << lowEnergyLimit / keV << " keV - "
123        << highEnergyLimit / GeV << " GeV"         125        << highEnergyLimit / GeV << " GeV"
124        << G4endl;                                 126        << G4endl;
125     }                                             127     }
126 }                                                 128 }
127                                                   129 
128 G4LowEnergyPhotoElectric::~G4LowEnergyPhotoEle    130 G4LowEnergyPhotoElectric::~G4LowEnergyPhotoElectric()
129 {                                                 131 {
130   delete crossSectionHandler;                     132   delete crossSectionHandler;
131   delete shellCrossSectionHandler;                133   delete shellCrossSectionHandler;
132   delete meanFreePathTable;                       134   delete meanFreePathTable;
133   delete rangeTest;                               135   delete rangeTest;
134   delete ElectronAngularGenerator;                136   delete ElectronAngularGenerator;
135 }                                                 137 }
136                                                   138 
137 void G4LowEnergyPhotoElectric::BuildPhysicsTab    139 void G4LowEnergyPhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
138 {                                                 140 {
139                                                   141 
140   crossSectionHandler->Clear();                   142   crossSectionHandler->Clear();
141   G4String crossSectionFile = "phot/pe-cs-";      143   G4String crossSectionFile = "phot/pe-cs-";
142   crossSectionHandler->LoadData(crossSectionFi    144   crossSectionHandler->LoadData(crossSectionFile);
143                                                   145 
144   shellCrossSectionHandler->Clear();              146   shellCrossSectionHandler->Clear();
145   G4String shellCrossSectionFile = "phot/pe-ss    147   G4String shellCrossSectionFile = "phot/pe-ss-cs-";
146   shellCrossSectionHandler->LoadShellData(shel    148   shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
147                                                   149 
148   delete meanFreePathTable;                       150   delete meanFreePathTable;
149   meanFreePathTable = crossSectionHandler->Bui    151   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
150 }                                                 152 }
151                                                   153 
152 G4VParticleChange* G4LowEnergyPhotoElectric::P    154 G4VParticleChange* G4LowEnergyPhotoElectric::PostStepDoIt(const G4Track& aTrack,
153                 const G4Step& aStep)              155                 const G4Step& aStep)
154 {                                                 156 {
155   // Fluorescence generated according to:         157   // Fluorescence generated according to:
156   // J. Stepanek ,"A program to determine the     158   // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
157   // subshell ionisation by a particle or due     159   // subshell ionisation by a particle or due to deexcitation or decay of radionuclides", 
158   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)       160   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
159                                                   161  
160   aParticleChange.Initialize(aTrack);             162   aParticleChange.Initialize(aTrack);
161                                                   163   
162   const G4DynamicParticle* incidentPhoton = aT    164   const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
163   G4double photonEnergy = incidentPhoton->GetK    165   G4double photonEnergy = incidentPhoton->GetKineticEnergy();
164   if (photonEnergy <= lowEnergyLimit)             166   if (photonEnergy <= lowEnergyLimit)
165     {                                             167     {
166       aParticleChange.ProposeTrackStatus(fStop    168       aParticleChange.ProposeTrackStatus(fStopAndKill);
167       aParticleChange.ProposeEnergy(0.);          169       aParticleChange.ProposeEnergy(0.);
168       aParticleChange.ProposeLocalEnergyDeposi    170       aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
169       return G4VDiscreteProcess::PostStepDoIt(    171       return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
170     }                                             172     }
171                                                   173  
172   G4ThreeVector photonDirection = incidentPhot    174   G4ThreeVector photonDirection = incidentPhoton->GetMomentumDirection(); // Returns the normalized direction of the momentum
173                                                   175 
174   // Select randomly one element in the curren    176   // Select randomly one element in the current material
175   const G4MaterialCutsCouple* couple = aTrack.    177   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
176   G4int Z = crossSectionHandler->SelectRandomA    178   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
177                                                   179 
178   // Select the ionised shell in the current a    180   // Select the ionised shell in the current atom according to shell cross sections
179   size_t shellIndex = shellCrossSectionHandler    181   size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
180                                                   182 
181   // Retrieve the corresponding identifier and    183   // Retrieve the corresponding identifier and binding energy of the selected shell
182   const G4RDAtomicTransitionManager* transitio    184   const G4RDAtomicTransitionManager* transitionManager = G4RDAtomicTransitionManager::Instance();
183   const G4RDAtomicShell* shell = transitionMan    185   const G4RDAtomicShell* shell = transitionManager->Shell(Z,shellIndex);
184   G4double bindingEnergy = shell->BindingEnerg    186   G4double bindingEnergy = shell->BindingEnergy();
185   G4int shellId = shell->ShellId();               187   G4int shellId = shell->ShellId();
186                                                   188 
187   // Create lists of pointers to DynamicPartic    189   // Create lists of pointers to DynamicParticles (photons and electrons)
188   // (Is the electron vector necessary? To be     190   // (Is the electron vector necessary? To be checked)
189   std::vector<G4DynamicParticle*>* photonVecto    191   std::vector<G4DynamicParticle*>* photonVector = 0;
190   std::vector<G4DynamicParticle*> electronVect    192   std::vector<G4DynamicParticle*> electronVector;
191                                                   193 
192   G4double energyDeposit = 0.0;                   194   G4double energyDeposit = 0.0;
193                                                   195 
194   // Primary outcoming electron                   196   // Primary outcoming electron
195   G4double eKineticEnergy = photonEnergy - bin    197   G4double eKineticEnergy = photonEnergy - bindingEnergy;
196                                                   198 
197   // There may be cases where the binding ener    199   // There may be cases where the binding energy of the selected shell is > photon energy
198   // In such cases do not generate secondaries    200   // In such cases do not generate secondaries
199   if (eKineticEnergy > 0.)                        201   if (eKineticEnergy > 0.)
200     {                                             202     {
201       // Generate the electron only if with la    203       // Generate the electron only if with large enough range w.r.t. cuts and safety
202       G4double safety = aStep.GetPostStepPoint    204       G4double safety = aStep.GetPostStepPoint()->GetSafety();
203                                                   205 
204       if (rangeTest->Escape(G4Electron::Electr    206       if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
205   {                                               207   {
206                                                   208 
207     // Calculate direction of the photoelectro    209     // Calculate direction of the photoelectron
208     G4ThreeVector gammaPolarization = incident    210     G4ThreeVector gammaPolarization = incidentPhoton->GetPolarization();
209     G4ThreeVector electronDirection = Electron    211     G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
210                                                   212 
211     // The electron is created ...                213     // The electron is created ...
212     G4DynamicParticle* electron = new G4Dynami    214     G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
213                      electronDirection,           215                      electronDirection,
214                      eKineticEnergy);             216                      eKineticEnergy);
215     electronVector.push_back(electron);           217     electronVector.push_back(electron);
216   }                                               218   }
217       else                                        219       else
218   {                                               220   {
219     energyDeposit += eKineticEnergy;              221     energyDeposit += eKineticEnergy;
220   }                                               222   }
221     }                                             223     }
222   else                                            224   else
223     {                                             225     {
224       bindingEnergy = photonEnergy;               226       bindingEnergy = photonEnergy;
225     }                                             227     }
226                                                   228 
227   G4int nElectrons = electronVector.size();       229   G4int nElectrons = electronVector.size();
228   size_t nTotPhotons = 0;                         230   size_t nTotPhotons = 0;
229   G4int nPhotons=0;                               231   G4int nPhotons=0;
230   const G4ProductionCutsTable* theCoupleTable=    232   const G4ProductionCutsTable* theCoupleTable=
231         G4ProductionCutsTable::GetProductionCu    233         G4ProductionCutsTable::GetProductionCutsTable();
232                                                   234 
233   size_t index = couple->GetIndex();              235   size_t index = couple->GetIndex();
234   G4double cutg = (*(theCoupleTable->GetEnergy    236   G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
235   cutg = std::min(cutForLowEnergySecondaryPhot    237   cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
236                                                   238   
237   G4double cute = (*(theCoupleTable->GetEnergy    239   G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
238   cute = std::min(cutForLowEnergySecondaryPhot    240   cute = std::min(cutForLowEnergySecondaryPhotons,cute);
239                                                   241 
240   G4DynamicParticle* aPhoton;                     242   G4DynamicParticle* aPhoton;
241                                                   243 
242   // Generation of fluorescence                   244   // Generation of fluorescence
243   // Data in EADL are available only for Z > 5    245   // Data in EADL are available only for Z > 5
244   // Protection to avoid generating photons in    246   // Protection to avoid generating photons in the unphysical case of
245   // shell binding energy > photon energy         247   // shell binding energy > photon energy
246   if (Z > 5  && (bindingEnergy > cutg || bindi    248   if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
247     {                                             249     {
248       photonVector = deexcitationManager.Gener    250       photonVector = deexcitationManager.GenerateParticles(Z,shellId);
249       nTotPhotons = photonVector->size();         251       nTotPhotons = photonVector->size();
250       for (size_t k=0; k<nTotPhotons; k++)        252       for (size_t k=0; k<nTotPhotons; k++)
251   {                                               253   {
252     aPhoton = (*photonVector)[k];                 254     aPhoton = (*photonVector)[k];
253     if (aPhoton)                                  255     if (aPhoton)
254       {                                           256       {
255               G4double itsCut = cutg;             257               G4double itsCut = cutg;
256               if(aPhoton->GetDefinition() == G    258               if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
257         G4double itsEnergy = aPhoton->GetKinet    259         G4double itsEnergy = aPhoton->GetKineticEnergy();
258                                                   260 
259         if (itsEnergy > itsCut && itsEnergy <=    261         if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
260     {                                             262     {
261       nPhotons++;                                 263       nPhotons++;
262       // Local energy deposit is given as the     264       // Local energy deposit is given as the sum of the
263       // energies of incident photons minus th    265       // energies of incident photons minus the energies
264       // of the outcoming fluorescence photons    266       // of the outcoming fluorescence photons
265       bindingEnergy -= itsEnergy;                 267       bindingEnergy -= itsEnergy;
266                                                   268 
267     }                                             269     }
268         else                                      270         else
269     {                                             271     {
270                   delete aPhoton;                 272                   delete aPhoton;
271                   (*photonVector)[k] = 0;         273                   (*photonVector)[k] = 0;
272                 }                                 274                 }
273       }                                           275       }
274   }                                               276   }
275     }                                             277     }
276                                                   278 
277   energyDeposit += bindingEnergy;                 279   energyDeposit += bindingEnergy;
278                                                   280 
279   G4int nSecondaries  = nElectrons + nPhotons;    281   G4int nSecondaries  = nElectrons + nPhotons;
280   aParticleChange.SetNumberOfSecondaries(nSeco    282   aParticleChange.SetNumberOfSecondaries(nSecondaries);
281                                                   283 
282   for (G4int l = 0; l<nElectrons; l++ )           284   for (G4int l = 0; l<nElectrons; l++ )
283     {                                             285     {
284       aPhoton = electronVector[l];                286       aPhoton = electronVector[l];
285       if(aPhoton) {                               287       if(aPhoton) {
286         aParticleChange.AddSecondary(aPhoton);    288         aParticleChange.AddSecondary(aPhoton);
287       }                                           289       }
288     }                                             290     }
289   for ( size_t ll = 0; ll < nTotPhotons; ll++)    291   for ( size_t ll = 0; ll < nTotPhotons; ll++)
290     {                                             292     {
291       aPhoton = (*photonVector)[ll];              293       aPhoton = (*photonVector)[ll];
292       if(aPhoton) {                               294       if(aPhoton) {
293         aParticleChange.AddSecondary(aPhoton);    295         aParticleChange.AddSecondary(aPhoton);
294       }                                           296       }
295     }                                             297     }
296                                                   298 
297   delete photonVector;                            299   delete photonVector;
298                                                   300 
299   if (energyDeposit < 0)                          301   if (energyDeposit < 0)
300     {                                             302     {
301       G4cout << "WARNING - "                      303       G4cout << "WARNING - "
302        << "G4LowEnergyPhotoElectric::PostStepD    304        << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
303        << G4endl;                                 305        << G4endl;
304       energyDeposit = 0;                          306       energyDeposit = 0;
305     }                                             307     }
306                                                   308 
307   // Kill the incident photon                     309   // Kill the incident photon
308   aParticleChange.ProposeMomentumDirection( 0.    310   aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
309   aParticleChange.ProposeEnergy( 0. );            311   aParticleChange.ProposeEnergy( 0. );
310                                                   312 
311   aParticleChange.ProposeLocalEnergyDeposit(en    313   aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
312   aParticleChange.ProposeTrackStatus( fStopAnd    314   aParticleChange.ProposeTrackStatus( fStopAndKill );
313                                                   315 
314   // Reset NbOfInteractionLengthLeft and retur    316   // Reset NbOfInteractionLengthLeft and return aParticleChange
315   return G4VDiscreteProcess::PostStepDoIt( aTr    317   return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
316 }                                                 318 }
317                                                   319 
318 G4bool G4LowEnergyPhotoElectric::IsApplicable(    320 G4bool G4LowEnergyPhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
319 {                                                 321 {
320   return ( &particle == G4Gamma::Gamma() );       322   return ( &particle == G4Gamma::Gamma() );
321 }                                                 323 }
322                                                   324 
323 G4double G4LowEnergyPhotoElectric::GetMeanFree    325 G4double G4LowEnergyPhotoElectric::GetMeanFreePath(const G4Track& track,
324                G4double, // previousStepSize      326                G4double, // previousStepSize
325                  G4ForceCondition*)               327                  G4ForceCondition*)
326 {                                                 328 {
327   const G4DynamicParticle* photon = track.GetD    329   const G4DynamicParticle* photon = track.GetDynamicParticle();
328   G4double energy = photon->GetKineticEnergy()    330   G4double energy = photon->GetKineticEnergy();
329   G4Material* material = track.GetMaterial();     331   G4Material* material = track.GetMaterial();
330   //  size_t materialIndex = material->GetInde    332   //  size_t materialIndex = material->GetIndex();
331                                                   333 
332   G4double meanFreePath = DBL_MAX;                334   G4double meanFreePath = DBL_MAX;
333                                                   335 
334   //  if (energy > highEnergyLimit)               336   //  if (energy > highEnergyLimit) 
335   //    meanFreePath = meanFreePathTable->Find    337   //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
336   //  else if (energy < lowEnergyLimit) meanFr    338   //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
337   //  else meanFreePath = meanFreePathTable->F    339   //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
338                                                   340 
339   G4double cross = shellCrossSectionHandler->V    341   G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
340   if(cross > 0.0) meanFreePath = 1.0/cross;       342   if(cross > 0.0) meanFreePath = 1.0/cross;
341                                                   343 
342   return meanFreePath;                            344   return meanFreePath;
343 }                                                 345 }
344                                                   346 
345 void G4LowEnergyPhotoElectric::SetCutForLowEnS    347 void G4LowEnergyPhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
346 {                                                 348 {
347   cutForLowEnergySecondaryPhotons = cut;          349   cutForLowEnergySecondaryPhotons = cut;
348   deexcitationManager.SetCutForSecondaryPhoton    350   deexcitationManager.SetCutForSecondaryPhotons(cut);
349 }                                                 351 }
350                                                   352 
351 void G4LowEnergyPhotoElectric::SetCutForLowEnS    353 void G4LowEnergyPhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
352 {                                                 354 {
353   cutForLowEnergySecondaryElectrons = cut;        355   cutForLowEnergySecondaryElectrons = cut;
354   deexcitationManager.SetCutForAugerElectrons(    356   deexcitationManager.SetCutForAugerElectrons(cut);
355 }                                                 357 }
356                                                   358 
357 void G4LowEnergyPhotoElectric::ActivateAuger(G    359 void G4LowEnergyPhotoElectric::ActivateAuger(G4bool val)
358 {                                                 360 {
359   deexcitationManager.ActivateAugerElectronPro    361   deexcitationManager.ActivateAugerElectronProduction(val);
360 }                                                 362 }
361                                                   363 
362 void G4LowEnergyPhotoElectric::SetAngularGener    364 void G4LowEnergyPhotoElectric::SetAngularGenerator(G4RDVPhotoElectricAngularDistribution* distribution)
363 {                                                 365 {
364   ElectronAngularGenerator = distribution;        366   ElectronAngularGenerator = distribution;
365   ElectronAngularGenerator->PrintGeneratorInfo    367   ElectronAngularGenerator->PrintGeneratorInformation();
366 }                                                 368 }
367                                                   369 
368 void G4LowEnergyPhotoElectric::SetAngularGener    370 void G4LowEnergyPhotoElectric::SetAngularGenerator(const G4String& name)
369 {                                                 371 {
370   if (name == "default")                          372   if (name == "default") 
371     {                                             373     {
372       delete ElectronAngularGenerator;            374       delete ElectronAngularGenerator;
373       ElectronAngularGenerator = new G4RDPhoto    375       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
374       generatorName = name;                       376       generatorName = name;
375     }                                             377     }
376   else if (name == "standard")                    378   else if (name == "standard")
377     {                                             379     {
378       delete ElectronAngularGenerator;            380       delete ElectronAngularGenerator;
379       ElectronAngularGenerator = new G4RDPhoto    381       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
380       generatorName = name;                       382       generatorName = name;
381     }                                             383     }
382   else if (name == "polarized")                   384   else if (name == "polarized")
383     {                                             385     {
384       delete ElectronAngularGenerator;            386       delete ElectronAngularGenerator;
385       ElectronAngularGenerator = new G4RDPhoto    387       ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
386       generatorName = name;                       388       generatorName = name;
387     }                                             389     }
388   else                                            390   else
389     {                                             391     {
390       G4Exception("G4LowEnergyPhotoElectric::S    392       G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator()",
391                   "InvalidSetup", FatalExcepti    393                   "InvalidSetup", FatalException,
392                   "Generator does not exist!")    394                   "Generator does not exist!");
393     }                                             395     }
394                                                   396 
395   ElectronAngularGenerator->PrintGeneratorInfo    397   ElectronAngularGenerator->PrintGeneratorInformation();
396 }                                                 398 }
397                                                   399