source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectric.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 13.8 KB
Line 
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: G4PenelopePhotoElectric.cc,v 1.16 2009/06/11 15:47:08 mantero Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Author: L. Pandola
31//
32// History:
33// --------
34// January 2003 - Created
35// 12 Feb 2003   MG Pia     Migration to "cuts per region"
36// 10 Mar 2003 V.Ivanchenko   Remome CutPerMaterial warning
37// 31 May 2005  L. Pandola  Added Sauter formula for the sampling of
38//                          the electron direction
39// 08 Jan 2009  L. Pandola  Check shell index to avoid mismatch between
40//                          the Penelope cross section database and the
41//                          G4AtomicTransitionManager database. It suppresses
42//                          a warning from G4AtomicTransitionManager only.
43//                          Results are unchanged.
44// --------------------------------------------------------------
45
46#include "G4PenelopePhotoElectric.hh"
47
48#include "G4ParticleDefinition.hh"
49#include "G4Track.hh"
50#include "G4Step.hh"
51#include "G4ForceCondition.hh"
52#include "G4Gamma.hh"
53#include "G4Electron.hh"
54#include "G4DynamicParticle.hh"
55#include "G4VParticleChange.hh"
56#include "G4ThreeVector.hh"
57#include "G4VCrossSectionHandler.hh"
58#include "G4CrossSectionHandler.hh"
59#include "G4VEMDataSet.hh"
60#include "G4CompositeEMDataSet.hh"
61#include "G4VDataSetAlgorithm.hh"
62#include "G4LogLogInterpolation.hh"
63#include "G4VRangeTest.hh"
64#include "G4RangeNoTest.hh"
65#include "G4AtomicTransitionManager.hh"
66#include "G4AtomicShell.hh"
67#include "G4MaterialCutsCouple.hh"
68#include "G4ProductionCutsTable.hh"
69
70G4PenelopePhotoElectric::G4PenelopePhotoElectric(const G4String& processName)
71  : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
72    intrinsicLowEnergyLimit(10*eV),
73    intrinsicHighEnergyLimit(100*GeV),
74    cutForLowEnergySecondaryPhotons(250.*eV),
75    cutForLowEnergySecondaryElectrons(250.*eV)
76{
77  if (lowEnergyLimit < intrinsicLowEnergyLimit || 
78      highEnergyLimit > intrinsicHighEnergyLimit)
79    {
80      G4Exception("G4PenelopePhotoElectric::G4PenelopePhotoElectric - energy limit outside intrinsic process validity range");
81    }
82
83  crossSectionHandler = new G4CrossSectionHandler();
84  shellCrossSectionHandler = new G4CrossSectionHandler();
85  meanFreePathTable = 0;
86  rangeTest = new G4RangeNoTest;
87
88  if (verboseLevel > 0) 
89    {
90      G4cout << GetProcessName() << " is created " << G4endl
91             << "Energy range: " 
92             << lowEnergyLimit / keV << " keV - "
93             << highEnergyLimit / GeV << " GeV" 
94             << G4endl;
95    }
96
97   G4cout << G4endl;
98   G4cout << "*******************************************************************************" << G4endl;
99   G4cout << "*******************************************************************************" << G4endl;
100   G4cout << "   The class G4PenelopePhotoElectric is NOT SUPPORTED ANYMORE. " << G4endl;
101   G4cout << "   It will be REMOVED with the next major release of Geant4. " << G4endl;
102   G4cout << "   Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
103   G4cout << "*******************************************************************************" << G4endl;
104   G4cout << "*******************************************************************************" << G4endl;
105   G4cout << G4endl;
106
107}
108
109G4PenelopePhotoElectric::~G4PenelopePhotoElectric()
110{
111  delete crossSectionHandler;
112  delete shellCrossSectionHandler;
113  delete meanFreePathTable;
114  delete rangeTest;
115}
116
117void G4PenelopePhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
118{
119 
120  crossSectionHandler->Clear();
121  G4String crossSectionFile = "penelope/ph-cs-pen-";
122  crossSectionHandler->LoadData(crossSectionFile);
123
124  shellCrossSectionHandler->Clear();
125  G4String shellCrossSectionFile = "penelope/ph-ss-cs-pen-";
126  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
127
128  delete meanFreePathTable;
129  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
130}
131
132G4VParticleChange* G4PenelopePhotoElectric::PostStepDoIt(const G4Track& aTrack,
133                                                         const G4Step& aStep)
134{
135  // Fluorescence generated according to:
136  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
137  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
138  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
139
140  aParticleChange.Initialize(aTrack);
141
142  const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
143  G4double photonEnergy = incidentPhoton->GetKineticEnergy();
144  if (photonEnergy <= lowEnergyLimit)
145    {
146      aParticleChange.ProposeTrackStatus(fStopAndKill);
147      aParticleChange.ProposeEnergy(0.);
148      aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
149      return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
150    }
151 
152  G4ParticleMomentum photonDirection = incidentPhoton->GetMomentumDirection();
153   
154  // Select randomly one element in the current material
155  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
156 
157  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
158
159  // Select the ionised shell in the current atom according to shell cross sections
160  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
161
162  // Retrieve the corresponding identifier and binding energy of the selected shell
163  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
164
165  //The number of shell cross section possibly reported in the Penelope database
166  //might be different from the number of shells in the G4AtomicTransitionManager
167  //(namely, Penelope may contain more shell, especially for very light elements).
168  //In order to avoid a warning message from the G4AtomicTransitionManager, I
169  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
170  //has a shellID>maxID, it sets the shellID to the last valid shell.
171  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
172  if (shellIndex >= numberOfShells)
173    shellIndex = numberOfShells-1;
174
175  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
176  G4double bindingEnergy = shell->BindingEnergy();
177  G4int shellId = shell->ShellId();
178
179  // Create lists of pointers to DynamicParticles (photons and electrons)
180  // (Is the electron vector necessary? To be checked)
181  std::vector<G4DynamicParticle*>* photonVector = 0;
182  std::vector<G4DynamicParticle*> electronVector;
183
184  G4double energyDeposit = 0.0;
185
186  // Primary outcoming electron
187  G4double eKineticEnergy = photonEnergy - bindingEnergy;
188
189  // There may be cases where the binding energy of the selected shell is > photon energy
190  // In such cases do not generate secondaries
191  if (eKineticEnergy > 0.)
192    {
193      // Generate the electron only if with large enough range w.r.t. cuts and safety
194      G4double safety = aStep.GetPostStepPoint()->GetSafety();
195
196      if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
197        {
198          // The electron is created
199          // Direction sampled from the Sauter distribution
200          G4double cosTheta = SampleElectronDirection(eKineticEnergy);
201          G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
202          G4double phi = twopi * G4UniformRand() ;
203          G4double dirx = sinTheta * std::cos(phi);
204          G4double diry = sinTheta * std::sin(phi);
205          G4double dirz = cosTheta ;
206          G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
207          electronDirection.rotateUz(photonDirection);
208
209          G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
210                                                               electronDirection, 
211                                                               eKineticEnergy);
212          electronVector.push_back(electron);
213        } 
214      else 
215        {
216          energyDeposit += eKineticEnergy;   
217        }
218    }
219  else
220    {
221      bindingEnergy = photonEnergy;
222    }
223
224  G4int nElectrons = electronVector.size();
225  size_t nTotPhotons = 0;
226  G4int nPhotons=0;
227 
228  const G4ProductionCutsTable* theCoupleTable=
229        G4ProductionCutsTable::GetProductionCutsTable();
230
231  size_t index = couple->GetIndex();
232  G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
233  cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
234 
235  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
236  cute = std::min(cutForLowEnergySecondaryPhotons,cute);
237
238  G4DynamicParticle* aPhoton; 
239
240  // Generation of fluorescence
241  // Data in EADL are available only for Z > 5
242  // Protection to avoid generating photons in the unphysical case of
243  // shell binding energy > photon energy
244  if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
245    {
246      photonVector = deexcitationManager.GenerateParticles(Z,shellId); 
247      nTotPhotons = photonVector->size();
248      for (size_t k=0; k<nTotPhotons; k++)
249        {
250          aPhoton = (*photonVector)[k];
251          if (aPhoton)
252            {
253              G4double itsCut = cutg;
254              if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
255              G4double itsEnergy = aPhoton->GetKineticEnergy();
256
257              if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
258                {
259                  nPhotons++;
260                  // Local energy deposit is given as the sum of the
261                  // energies of incident photons minus the energies
262                  // of the outcoming fluorescence photons
263                  bindingEnergy -= itsEnergy;
264                 
265                }
266              else
267                { 
268                  delete aPhoton; 
269                  (*photonVector)[k] = 0;
270                }
271            }
272        }
273    }
274
275  energyDeposit += bindingEnergy;
276
277  G4int nSecondaries  = nElectrons + nPhotons;
278  aParticleChange.SetNumberOfSecondaries(nSecondaries);
279
280  for (G4int l = 0; l<nElectrons; l++ )
281    {
282      aPhoton = electronVector[l];
283      if(aPhoton) {
284        aParticleChange.AddSecondary(aPhoton);
285      }
286    }
287  for ( size_t ll = 0; ll < nTotPhotons; ll++) 
288    {
289      aPhoton = (*photonVector)[ll];
290      if(aPhoton) {
291        aParticleChange.AddSecondary(aPhoton);
292      } 
293    }
294 
295  delete photonVector;
296 
297  if (energyDeposit < 0)
298    {
299      G4cout << "WARNING - "
300             << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
301             << G4endl;
302      energyDeposit = 0;
303    }
304 
305  // Kill the incident photon
306  aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
307  aParticleChange.ProposeEnergy( 0. );
308 
309  aParticleChange.ProposeLocalEnergyDeposit(energyDeposit); 
310  aParticleChange.ProposeTrackStatus( fStopAndKill ); 
311 
312  // Reset NbOfInteractionLengthLeft and return aParticleChange
313  return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
314}
315
316G4bool G4PenelopePhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
317{
318  return ( &particle == G4Gamma::Gamma() ); 
319}
320
321G4double G4PenelopePhotoElectric::GetMeanFreePath(const G4Track& track, 
322                                                  G4double, // previousStepSize
323                                                  G4ForceCondition*)
324{
325  const G4DynamicParticle* photon = track.GetDynamicParticle();
326  G4double energy = photon->GetKineticEnergy();
327  G4Material* material = track.GetMaterial();
328
329  G4double meanFreePath = DBL_MAX;
330
331  G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
332  if(cross > 0.0) meanFreePath = 1.0/cross;
333
334  return meanFreePath;
335}
336
337void G4PenelopePhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
338{
339  cutForLowEnergySecondaryPhotons = cut;
340  deexcitationManager.SetCutForSecondaryPhotons(cut);
341}
342
343void G4PenelopePhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
344{
345  cutForLowEnergySecondaryElectrons = cut;
346  deexcitationManager.SetCutForAugerElectrons(cut);
347}
348
349void G4PenelopePhotoElectric::ActivateAuger(G4bool val)
350{
351  deexcitationManager.ActivateAugerElectronProduction(val);
352}
353
354
355G4double G4PenelopePhotoElectric::SampleElectronDirection(G4double energy)
356{
357  G4double costheta = 1.0;
358  if (energy>1*GeV) return costheta;
359
360  //1) initialize energy-dependent variables
361  // Variable naming according to Eq. (2.24) of Penelope Manual
362  // (pag. 44)
363  G4double gamma = 1.0 + energy/electron_mass_c2;
364  G4double gamma2 = gamma*gamma;
365  G4double beta = std::sqrt((gamma2-1.0)/gamma2);
366 
367  // ac corresponds to "A" of Eq. (2.31)
368  //
369  G4double ac = (1.0/beta) - 1.0;
370  G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
371  G4double a2 = ac + 2.0;
372  G4double gtmax = 2.0*(a1 + 1.0/ac);
373
374  G4double tsam = 0;
375  G4double gtr = 0;
376 
377  //2) sampling. Eq. (2.31) of Penelope Manual
378  // tsam = 1-std::cos(theta)
379  // gtr = rejection function according to Eq. (2.28)
380  do{
381    G4double rand = G4UniformRand();
382    tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
383    gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
384  }while(G4UniformRand()*gtmax > gtr);
385  costheta = 1.0-tsam;
386  return costheta;
387}
388
389
390
391
392
Note: See TracBrowser for help on using the repository browser.