source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectricModel.cc @ 991

Last change on this file since 991 was 991, checked in by garnier, 15 years ago

update
File size: 15.7 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// $Id: G4PenelopePhotoElectricModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 08 Oct 2008   L Pandola    Migration from process to model
34//
35
36#include "G4PenelopePhotoElectricModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4PhysicsTable.hh"
42#include "G4ElementTable.hh"
43#include "G4Element.hh"
44#include "G4CrossSectionHandler.hh"
45#include "G4AtomicTransitionManager.hh"
46#include "G4AtomicShell.hh"
47#include "G4Gamma.hh"
48#include "G4Electron.hh"
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52
53G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition*,
54                                             const G4String& nam)
55  :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
56   shellCrossSectionHandler(0)
57{
58  fIntrinsicLowEnergyLimit = 100.0*eV;
59  fIntrinsicHighEnergyLimit = 100.0*GeV;
60  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
61  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
62  //
63  fUseAtomicDeexcitation = true;
64  verboseLevel= 0;
65  // Verbosity scale:
66  // 0 = nothing
67  // 1 = warning for energy non-conservation
68  // 2 = details of energy budget
69  // 3 = calculation of cross sections, file openings, sampling of atoms
70  // 4 = entering in methods
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel()
76{ 
77  if (crossSectionHandler) delete crossSectionHandler;
78  if (shellCrossSectionHandler) delete shellCrossSectionHandler;
79}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition*,
84                                       const G4DataVector& )
85{
86  if (verboseLevel > 3)
87    G4cout << "Calling  G4PenelopePhotoElectricModel::Initialise()" << G4endl;
88  if (crossSectionHandler)
89    {
90      crossSectionHandler->Clear();
91      delete crossSectionHandler;
92    }
93  if (shellCrossSectionHandler)
94    {
95      shellCrossSectionHandler->Clear();
96      delete shellCrossSectionHandler;
97    }
98
99  //Check energy limits
100  if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
101    {
102      G4cout << "G4PenelopePhotoElectricModel: low energy limit increased from " <<
103        LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << 
104        G4endl;
105      SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
106    }
107  if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
108    {
109      G4cout << "G4PenelopePhotoElectricModel: high energy limit decreased from " <<
110        HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" 
111             << G4endl;
112      SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
113    }
114
115
116  //Re-initialize cross section handlers
117  crossSectionHandler = new G4CrossSectionHandler();
118  crossSectionHandler->Clear();
119  G4String crossSectionFile = "penelope/ph-cs-pen-";
120  crossSectionHandler->LoadData(crossSectionFile);
121  shellCrossSectionHandler = new G4CrossSectionHandler();
122  shellCrossSectionHandler->Clear();
123  crossSectionFile = "penelope/ph-ss-cs-pen-";
124  shellCrossSectionHandler->LoadShellData(crossSectionFile);
125  //This is used to retrieve cross section values later on
126  crossSectionHandler->BuildMeanFreePathForMaterials();
127
128  if (verboseLevel > 2) 
129    G4cout << "Loaded cross section files for PenelopePhotoElectric" << G4endl;
130
131  G4cout << "Penelope Photo-Electric model is initialized " << G4endl
132         << "Energy range: "
133         << LowEnergyLimit() / MeV << " MeV - "
134         << HighEnergyLimit() / GeV << " GeV"
135         << G4endl;
136
137  if(isInitialised) return;
138
139  if(pParticleChange)
140    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
141  else
142    fParticleChange = new G4ParticleChangeForGamma();
143  isInitialised = true;
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom(
149                                       const G4ParticleDefinition*,
150                                             G4double energy,
151                                             G4double Z, G4double,
152                                             G4double, G4double)
153{
154  //
155  // Penelope model. Use data-driven approach for cross section estimate (and
156  // also shell sampling from a given atom). Data are from the Livermore database
157  //  D.E. Cullen et al., Report UCRL-50400 (1989)
158  //
159
160  if (verboseLevel > 3)
161    G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
162
163  G4int iZ = (G4int) Z;
164  if (!crossSectionHandler)
165    {
166      G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
167      G4cout << "The cross section handler is not correctly initialized" << G4endl;
168      G4Exception();
169    }
170  G4double cs = crossSectionHandler->FindValue(iZ,energy);
171 
172  if (verboseLevel > 2)
173    G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
174      " = " << cs/barn << " barn" << G4endl;
175  return cs;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
181                                              const G4MaterialCutsCouple* couple,
182                                              const G4DynamicParticle* aDynamicGamma,
183                                              G4double,
184                                              G4double)
185{
186  //
187  // Photoelectric effect, Penelope model
188  //
189  // The target atom and the target shell are sampled according to the Livermore
190  // database
191  //  D.E. Cullen et al., Report UCRL-50400 (1989)
192  // The angular distribution of the electron in the final state is sampled
193  // according to the Sauter distribution from
194  //  F. Sauter, Ann. Phys. 11 (1931) 454
195  // The energy of the final electron is given by the initial photon energy minus
196  // the binding energy. Fluorescence de-excitation is subsequently produced
197  // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
198  //  J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
199
200  if (verboseLevel > 3)
201    G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
202
203  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
204
205  if (photonEnergy <= LowEnergyLimit())
206  {
207      fParticleChange->ProposeTrackStatus(fStopAndKill);
208      fParticleChange->SetProposedKineticEnergy(0.);
209      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
210      return ;
211  }
212
213  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
214
215  // Select randomly one element in the current material
216  if (verboseLevel > 2)
217    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
218  //use crossSectionHandler instead of G4EmElementSelector because in this case
219  //the dimension of the table is equal to the dimension of the database (less interpolation errors)
220  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
221  if (verboseLevel > 2)
222    G4cout << "Selected Z = " << Z << G4endl;
223
224  // Select the ionised shell in the current atom according to shell cross sections
225  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
226
227  // Retrieve the corresponding identifier and binding energy of the selected shell
228  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
229  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
230  G4double bindingEnergy = shell->BindingEnergy();
231  G4int shellId = shell->ShellId();
232
233  G4double localEnergyDeposit = 0.0;
234
235  // Primary outcoming electron
236  G4double eKineticEnergy = photonEnergy - bindingEnergy;
237 
238  G4double cutForLowEnergySecondaryParticles = 250.0*eV;
239  const G4ProductionCutsTable* theCoupleTable=
240    G4ProductionCutsTable::GetProductionCutsTable();
241  size_t indx = couple->GetIndex();
242  G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
243  cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
244
245  // There may be cases where the binding energy of the selected shell is > photon energy
246  // In such cases do not generate secondaries
247  if (eKineticEnergy > 0.)
248    {
249      //Now check if the electron is above cuts: if so, it is created explicitely
250      if (eKineticEnergy > cutE)
251        {
252          // The electron is created
253          // Direction sampled from the Sauter distribution
254          G4double cosTheta = SampleElectronDirection(eKineticEnergy);
255          G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
256          G4double phi = twopi * G4UniformRand() ;
257          G4double dirx = sinTheta * std::cos(phi);
258          G4double diry = sinTheta * std::sin(phi);
259          G4double dirz = cosTheta ;
260          G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
261          electronDirection.rotateUz(photonDirection);
262          G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
263                                                               electronDirection, 
264                                                               eKineticEnergy);
265          fvect->push_back(electron);
266        } 
267      else 
268        localEnergyDeposit += eKineticEnergy;   
269    }
270  else
271      bindingEnergy = photonEnergy;
272   
273 
274  G4double energyInFluorescence = 0; //testing purposes
275
276  //Now, take care of fluorescence, if required
277  if (fUseAtomicDeexcitation)
278    {
279      G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
280      cutG = std::min(cutForLowEnergySecondaryParticles,cutG);
281     
282      std::vector<G4DynamicParticle*>* photonVector = 0;
283
284      // Protection to avoid generating photons in the unphysical case of
285      // shell binding energy > photon energy
286      if (Z > 5  && (bindingEnergy > cutG || bindingEnergy > cutE))
287        {
288          photonVector = deexcitationManager.GenerateParticles(Z,shellId); 
289          //Check for single photons (if they are above threshold)
290          for (size_t k=0; k< photonVector->size(); k++)
291            {
292              G4DynamicParticle* aPhoton = (*photonVector)[k];
293              if (aPhoton)
294                {
295                  G4double itsCut = cutG;
296                  if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
297                  G4double itsEnergy = aPhoton->GetKineticEnergy();
298                 
299                  if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
300                    {
301                      // Local energy deposit is given as the sum of the
302                      // energies of incident photons minus the energies
303                      // of the outcoming fluorescence photons
304                      bindingEnergy -= itsEnergy;
305                     
306                    }
307                  else
308                    { 
309                      (*photonVector)[k] = 0;
310                    }
311                }
312            }
313        }
314      //Register valid secondaries
315      if (photonVector)
316        {
317          for ( size_t ll = 0; ll <photonVector->size(); ll++) 
318            {
319              G4DynamicParticle* aPhoton = (*photonVector)[ll];
320              if (aPhoton) 
321                {
322                  energyInFluorescence += aPhoton->GetKineticEnergy();
323                  fvect->push_back(aPhoton);
324                }
325            }
326          delete photonVector;
327        }
328    }
329  //Residual energy is deposited locally
330  localEnergyDeposit += bindingEnergy;
331     
332
333  if (localEnergyDeposit < 0)
334    {
335      G4cout << "WARNING - "
336             << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
337             << G4endl;
338      localEnergyDeposit = 0;
339    }
340
341 //Update the status of the primary gamma (kill it)
342  fParticleChange->SetProposedKineticEnergy(0.);
343  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
344  fParticleChange->ProposeTrackStatus(fStopAndKill);
345
346  if (verboseLevel > 1)
347    {
348      G4cout << "-----------------------------------------------------------" << G4endl;
349      G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
350      G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
351      G4cout << "-----------------------------------------------------------" << G4endl;
352      if (eKineticEnergy)
353        G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
354      G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
355      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
356      G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << 
357        " keV" << G4endl;
358      G4cout << "-----------------------------------------------------------" << G4endl;
359    }
360  if (verboseLevel > 0)
361    {
362      G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
363      if (energyDiff > 0.05*keV)
364        G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " << 
365          (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << " keV (final) vs. " << 
366          photonEnergy/keV << " keV (initial)" << G4endl;
367    }
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371
372void G4PenelopePhotoElectricModel::ActivateAuger(G4bool augerbool)
373{
374  if (!fUseAtomicDeexcitation)
375    {
376      G4cout << "WARNING - G4PenelopePhotoElectricModel" << G4endl;
377      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
378      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
379    }
380  deexcitationManager.ActivateAugerElectronProduction(augerbool);
381  if (verboseLevel > 1)
382    G4cout << "Auger production set to " << augerbool << G4endl;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy)
388{
389  G4double costheta = 1.0;
390  if (energy>1*GeV) return costheta;
391 
392  //1) initialize energy-dependent variables
393  // Variable naming according to Eq. (2.24) of Penelope Manual
394  // (pag. 44)
395  G4double gamma = 1.0 + energy/electron_mass_c2;
396  G4double gamma2 = gamma*gamma;
397  G4double beta = std::sqrt((gamma2-1.0)/gamma2);
398   
399  // ac corresponds to "A" of Eq. (2.31)
400  //
401  G4double ac = (1.0/beta) - 1.0;
402  G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
403  G4double a2 = ac + 2.0;
404  G4double gtmax = 2.0*(a1 + 1.0/ac);
405 
406  G4double tsam = 0;
407  G4double gtr = 0;
408
409  //2) sampling. Eq. (2.31) of Penelope Manual
410  // tsam = 1-std::cos(theta)
411  // gtr = rejection function according to Eq. (2.28)
412  do{
413    G4double rand = G4UniformRand();
414    tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
415    gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
416  }while(G4UniformRand()*gtmax > gtr);
417  costheta = 1.0-tsam;
418  return costheta;
419}
420
Note: See TracBrowser for help on using the repository browser.