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

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

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