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

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

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