Geant4  10.02.p01
G4LowEnergyPhotoElectric.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
29 //
30 // Author: A. Forti
31 // Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
32 //
33 // History:
34 // --------
35 // October 1998 - low energy modifications by Alessandra Forti
36 // Added Livermore data table construction methods A. Forti
37 // Modified BuildMeanFreePath to read new data tables A. Forti
38 // Added EnergySampling method A. Forti
39 // Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
40 // Added SelectRandomAtom A. Forti
41 // Added map of the elements A. Forti
42 // 10.04.2000 VL
43 // - Correcting Fluorescence transition probabilities in order to take into account
44 // non-radiative transitions. No Auger electron simulated yet: energy is locally deposited.
45 // 17.02.2000 Veronique Lefebure
46 // - bugs corrected in fluorescence simulation:
47 // . when final use of binding energy: no photon was ever created
48 // . no Fluorescence was simulated when the photo-electron energy
49 // was below production threshold.
50 //
51 // 07-09-99, if no e- emitted: edep=photon energy, mma
52 // 24.04.01 V.Ivanchenko remove RogueWave
53 // 12.08.2001 MGP Revised according to a design iteration
54 // 16.09.2001 E. Guardincerri Added fluorescence generation
55 // 06.10.2001 MGP Added protection to avoid negative electron energies
56 // when binding energy of selected shell > photon energy
57 // 18.04.2001 V.Ivanchenko Fix problem with low energy gammas from fluorescence
58 // MeanFreePath is calculated by crosSectionHandler directly
59 // 31.05.2002 V.Ivanchenko Add path of Fluo + Auger cuts to AtomicDeexcitation
60 // 14.06.2002 V.Ivanchenko By default do not cheak range of e-
61 // 21.01.2003 V.Ivanchenko Cut per region
62 // 10.05.2004 P.Rodrigues Changes to accommodate new angular generators
63 // 20.01.2006 A.Trindade Changes to accommodate polarized angular generator
64 //
65 // --------------------------------------------------------------
66 
68 
73 
74 #include "G4PhysicalConstants.hh"
75 #include "G4SystemOfUnits.hh"
76 #include "G4ParticleDefinition.hh"
77 #include "G4Track.hh"
78 #include "G4Step.hh"
79 #include "G4ForceCondition.hh"
80 #include "G4Gamma.hh"
81 #include "G4Electron.hh"
82 #include "G4DynamicParticle.hh"
83 #include "G4VParticleChange.hh"
84 #include "G4ThreeVector.hh"
87 #include "G4RDVEMDataSet.hh"
89 #include "G4RDVDataSetAlgorithm.hh"
91 #include "G4RDVRangeTest.hh"
92 #include "G4RDRangeNoTest.hh"
94 #include "G4RDAtomicShell.hh"
95 #include "G4ProductionCutsTable.hh"
96 
98  : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
99  intrinsicLowEnergyLimit(10*eV),
100  intrinsicHighEnergyLimit(100*GeV),
101  cutForLowEnergySecondaryPhotons(250.*eV),
102  cutForLowEnergySecondaryElectrons(250.*eV)
103 {
106  {
107  G4Exception("G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric()",
108  "OutOfRange", FatalException,
109  "Energy limit outside intrinsic process validity range!");
110  }
111 
114  meanFreePathTable = 0;
116  generatorName = "geant4.6.2";
117  ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator"); // default generator
118 
119 
120  if (verboseLevel > 0)
121  {
122  G4cout << GetProcessName() << " is created " << G4endl
123  << "Energy range: "
124  << lowEnergyLimit / keV << " keV - "
125  << highEnergyLimit / GeV << " GeV"
126  << G4endl;
127  }
128 }
129 
131 {
132  delete crossSectionHandler;
134  delete meanFreePathTable;
135  delete rangeTest;
137 }
138 
140 {
141 
143  G4String crossSectionFile = "phot/pe-cs-";
144  crossSectionHandler->LoadData(crossSectionFile);
145 
147  G4String shellCrossSectionFile = "phot/pe-ss-cs-";
148  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
149 
150  delete meanFreePathTable;
152 }
153 
155  const G4Step& aStep)
156 {
157  // Fluorescence generated according to:
158  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
159  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
160  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
161 
162  aParticleChange.Initialize(aTrack);
163 
164  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
165  G4double photonEnergy = incidentPhoton->GetKineticEnergy();
166  if (photonEnergy <= lowEnergyLimit)
167  {
171  return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
172  }
173 
174  G4ThreeVector photonDirection = incidentPhoton->GetMomentumDirection(); // Returns the normalized direction of the momentum
175 
176  // Select randomly one element in the current material
177  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
178  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
179 
180  // Select the ionised shell in the current atom according to shell cross sections
181  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
182 
183  // Retrieve the corresponding identifier and binding energy of the selected shell
185  const G4RDAtomicShell* shell = transitionManager->Shell(Z,shellIndex);
187  G4int shellId = shell->ShellId();
188 
189  // Create lists of pointers to DynamicParticles (photons and electrons)
190  // (Is the electron vector necessary? To be checked)
191  std::vector<G4DynamicParticle*>* photonVector = 0;
192  std::vector<G4DynamicParticle*> electronVector;
193 
194  G4double energyDeposit = 0.0;
195 
196  // Primary outcoming electron
197  G4double eKineticEnergy = photonEnergy - bindingEnergy;
198 
199  // There may be cases where the binding energy of the selected shell is > photon energy
200  // In such cases do not generate secondaries
201  if (eKineticEnergy > 0.)
202  {
203  // Generate the electron only if with large enough range w.r.t. cuts and safety
204  G4double safety = aStep.GetPostStepPoint()->GetSafety();
205 
206  if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
207  {
208 
209  // Calculate direction of the photoelectron
210  G4ThreeVector gammaPolarization = incidentPhoton->GetPolarization();
211  G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
212 
213  // The electron is created ...
215  electronDirection,
216  eKineticEnergy);
217  electronVector.push_back(electron);
218  }
219  else
220  {
221  energyDeposit += eKineticEnergy;
222  }
223  }
224  else
225  {
226  bindingEnergy = photonEnergy;
227  }
228 
229  G4int nElectrons = electronVector.size();
230  size_t nTotPhotons = 0;
231  G4int nPhotons=0;
232  const G4ProductionCutsTable* theCoupleTable=
234 
235  size_t index = couple->GetIndex();
236  G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
238 
239  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
241 
242  G4DynamicParticle* aPhoton;
243 
244  // Generation of fluorescence
245  // Data in EADL are available only for Z > 5
246  // Protection to avoid generating photons in the unphysical case of
247  // shell binding energy > photon energy
248  if (Z > 5 && (bindingEnergy > cutg || bindingEnergy > cute))
249  {
250  photonVector = deexcitationManager.GenerateParticles(Z,shellId);
251  nTotPhotons = photonVector->size();
252  for (size_t k=0; k<nTotPhotons; k++)
253  {
254  aPhoton = (*photonVector)[k];
255  if (aPhoton)
256  {
257  G4double itsCut = cutg;
258  if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
259  G4double itsEnergy = aPhoton->GetKineticEnergy();
260 
261  if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
262  {
263  nPhotons++;
264  // Local energy deposit is given as the sum of the
265  // energies of incident photons minus the energies
266  // of the outcoming fluorescence photons
267  bindingEnergy -= itsEnergy;
268 
269  }
270  else
271  {
272  delete aPhoton;
273  (*photonVector)[k] = 0;
274  }
275  }
276  }
277  }
278 
279  energyDeposit += bindingEnergy;
280 
281  G4int nSecondaries = nElectrons + nPhotons;
283 
284  for (G4int l = 0; l<nElectrons; l++ )
285  {
286  aPhoton = electronVector[l];
287  if(aPhoton) {
288  aParticleChange.AddSecondary(aPhoton);
289  }
290  }
291  for ( size_t ll = 0; ll < nTotPhotons; ll++)
292  {
293  aPhoton = (*photonVector)[ll];
294  if(aPhoton) {
295  aParticleChange.AddSecondary(aPhoton);
296  }
297  }
298 
299  delete photonVector;
300 
301  if (energyDeposit < 0)
302  {
303  G4cout << "WARNING - "
304  << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
305  << G4endl;
306  energyDeposit = 0;
307  }
308 
309  // Kill the incident photon
312 
315 
316  // Reset NbOfInteractionLengthLeft and return aParticleChange
317  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
318 }
319 
321 {
322  return ( &particle == G4Gamma::Gamma() );
323 }
324 
326  G4double, // previousStepSize
328 {
330  G4double energy = photon->GetKineticEnergy();
331  G4Material* material = track.GetMaterial();
332  // size_t materialIndex = material->GetIndex();
333 
334  G4double meanFreePath = DBL_MAX;
335 
336  // if (energy > highEnergyLimit)
337  // meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
338  // else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
339  // else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
340 
341  G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
342  if(cross > 0.0) meanFreePath = 1.0/cross;
343 
344  return meanFreePath;
345 }
346 
348 {
351 }
352 
354 {
357 }
358 
360 {
362 }
363 
365 {
366  ElectronAngularGenerator = distribution;
368 }
369 
371 {
372  if (name == "default")
373  {
375  ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
377  }
378  else if (name == "standard")
379  {
381  ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
383  }
384  else if (name == "polarized")
385  {
387  ElectronAngularGenerator = new G4RDPhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
389  }
390  else
391  {
392  G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator()",
393  "InvalidSetup", FatalException,
394  "Generator does not exist!");
395  }
396 
398 }
G4double ValueForMaterial(const G4Material *material, G4double e) const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
void SetCutForSecondaryPhotons(G4double cut)
G4int verboseLevel
Definition: G4VProcess.hh:368
std::vector< G4DynamicParticle * > * GenerateParticles(G4int Z, G4int shellId)
void ActivateAugerElectronProduction(G4bool val)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsApplicable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
G4String name
Definition: TRTMaterials.hh:40
void BuildPhysicsTable(const G4ParticleDefinition &photon)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
void LoadShellData(const G4String &dataFile)
G4RDVPhotoElectricAngularDistribution * ElectronAngularGenerator
G4ParticleDefinition * GetDefinition() const
G4RDVEMDataSet * BuildMeanFreePathForMaterials(const G4DataVector *energyCuts=0)
int G4int
Definition: G4Types.hh:78
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetAngularGenerator(G4RDVPhotoElectricAngularDistribution *distribution)
static G4RDAtomicTransitionManager * Instance()
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double BindingEnergy() const
G4int SelectRandomShell(G4int Z, G4double e) const
G4LowEnergyPhotoElectric(const G4String &processName="LowEnPhotoElec")
G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
G4RDVCrossSectionHandler * crossSectionHandler
static const double GeV
Definition: G4SIunits.hh:214
Definition: G4Step.hh:76
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual G4bool Escape(const G4ParticleDefinition *particle, const G4MaterialCutsCouple *couple, G4double energy, G4double safety) const =0
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:212
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
G4StepPoint * GetPostStepPoint() const
void SetCutForAugerElectrons(G4double cut)
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const
const G4ThreeVector & GetPolarization() const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4RDAtomicShell * Shell(G4int Z, size_t shellIndex) const
G4double GetSafety() const
void AddSecondary(G4Track *aSecondary)
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4RDVCrossSectionHandler * shellCrossSectionHandler
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
virtual G4ThreeVector GetPhotoElectronDirection(const G4ThreeVector &direction, const G4double kineticEnergy, const G4ThreeVector &polarization, const G4int shellID) const =0
void LoadData(const G4String &dataFile)
G4RDAtomicDeexcitation deexcitationManager
static const double keV
Definition: G4SIunits.hh:213
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4double bindingEnergy(G4int A, G4int Z)
#define DBL_MAX
Definition: templates.hh:83
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
G4int ShellId() const