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

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

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