source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.cc @ 1228

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

update geant4.9.3 tag
File size: 12.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: G4LivermorePhotoElectricModel.cc,v 1.9 2009/10/23 09:31:03 pandola Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29//
30// Author: Sebastien Inserti
31//         30 October 2008
32//
33// History:
34// --------
35// 15 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
36//                  - apply internal high-energy limit only in constructor
37//                  - do not apply low-energy limit (default is 0)
38//                  - remove GetMeanFreePath method and table
39//                  - simplify sampling of deexcitation by using cut in energy
40//                  - added protection against numerical problem in energy sampling
41//                  - use G4ElementSelector
42// 23 Oct 2009   L Pandola
43//                  - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
44//                    set as "true" (default would be false)
45//
46
47#include "G4LivermorePhotoElectricModel.hh"
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51using namespace std;
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4ParticleDefinition*,
56                                             const G4String& nam)
57:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
58 crossSectionHandler(0),shellCrossSectionHandler(0),ElectronAngularGenerator(0)
59{
60  lowEnergyLimit = 250 * eV; 
61  highEnergyLimit = 100 * GeV;
62  //  SetLowEnergyLimit(lowEnergyLimit);
63  SetHighEnergyLimit(highEnergyLimit);
64
65  //Set atomic deexcitation by default
66  SetDeexcitationFlag(true);
67  ActivateAuger(false);
68   
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  if(verboseLevel>0) {
77    G4cout << "Livermore PhotoElectric is constructed " << G4endl
78           << "Energy range: "
79           << lowEnergyLimit / eV << " eV - "
80           << highEnergyLimit / GeV << " GeV"
81           << G4endl;
82  }
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
88{ 
89  if (crossSectionHandler) delete crossSectionHandler;
90  if (shellCrossSectionHandler) delete shellCrossSectionHandler;
91  if (ElectronAngularGenerator) delete ElectronAngularGenerator;
92}
93
94//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95
96void 
97G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
98                                          const G4DataVector&)
99{
100  if (verboseLevel > 3)
101    G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
102
103  if (crossSectionHandler)
104  {
105    crossSectionHandler->Clear();
106    delete crossSectionHandler;
107  }
108 
109  if (shellCrossSectionHandler)
110  {
111    shellCrossSectionHandler->Clear();
112    delete shellCrossSectionHandler;
113  }
114 
115  // Read data tables for all materials
116 
117  crossSectionHandler = new G4CrossSectionHandler();
118  crossSectionHandler->Clear();
119  G4String crossSectionFile = "phot/pe-cs-";
120  crossSectionHandler->LoadData(crossSectionFile);
121
122  shellCrossSectionHandler = new G4CrossSectionHandler();
123  shellCrossSectionHandler->Clear();
124  G4String shellCrossSectionFile = "phot/pe-ss-cs-";
125  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
126 
127  // SI - Simple generator is buggy
128  //generatorName = "geant4.6.2";
129  //ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator");              // default generator
130  ElectronAngularGenerator = 
131    new G4PhotoElectricAngularGeneratorSauterGavrila("GEANTSauterGavrilaGenerator");       
132
133  // 
134  if (verboseLevel > 2) 
135    G4cout << "Loaded cross section files for Livermore PhotoElectric model" << G4endl;
136
137  //  InitialiseElementSelectors(particle,cuts);
138
139  if (verboseLevel > 0) { 
140    G4cout << "Livermore PhotoElectric model is initialized " << G4endl
141           << "Energy range: "
142           << LowEnergyLimit() / eV << " eV - "
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 G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(
155                                       const G4ParticleDefinition*,
156                                             G4double GammaEnergy,
157                                             G4double Z, G4double,
158                                             G4double, G4double)
159{
160  if (verboseLevel > 3)
161    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel" 
162           << G4endl;
163
164  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
165    return 0;
166
167  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
168  return cs;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
173void 
174G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
175                                                 const G4MaterialCutsCouple* couple,
176                                                 const G4DynamicParticle* aDynamicGamma,
177                                                 G4double,
178                                                 G4double)
179{
180
181  // Fluorescence generated according to:
182  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
183  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
184  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
185 
186  if (verboseLevel > 3)
187    G4cout << "Calling SampleSecondaries() of G4LivermorePhotoElectricModel" << G4endl;
188
189  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
190 
191  // kill incident photon
192  fParticleChange->SetProposedKineticEnergy(0.);
193  fParticleChange->ProposeTrackStatus(fStopAndKill);   
194
195  // low-energy gamma is absorpted by this process
196  if (photonEnergy <= lowEnergyLimit)
197    {
198      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
199      return;
200    }
201 
202  // Returns the normalized direction of the momentum
203  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection(); 
204
205  // Select randomly one element in the current material
206  //  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
207  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
208  const G4Element* elm = SelectRandomAtom(couple->GetMaterial(),particle,photonEnergy);
209  G4int Z = (G4int)elm->GetZ();
210
211  // Select the ionised shell in the current atom according to shell cross sections
212  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
213
214  // Retrieve the corresponding identifier and binding energy of the selected shell
215  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
216  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
217  G4double bindingEnergy = shell->BindingEnergy();
218  G4int shellId = shell->ShellId();
219
220  // Primary outcoming electron
221  G4double eKineticEnergy = photonEnergy - bindingEnergy;
222
223  // There may be cases where the binding energy of the selected shell is > photon energy
224  // In such cases do not generate secondaries
225  if (eKineticEnergy > 0.)
226    {
227      // Calculate direction of the photoelectron
228      G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
229      G4ThreeVector electronDirection = 
230        ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,
231                                                            eKineticEnergy,
232                                                            gammaPolarization,
233                                                            shellId);
234
235      // The electron is created ...
236      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
237                                                           electronDirection,
238                                                           eKineticEnergy);
239      fvect->push_back(electron);
240    }
241  else
242    {
243      bindingEnergy = photonEnergy;
244    }
245
246  // deexcitation
247  if(DeexcitationFlag() && Z > 5) {
248
249    const G4ProductionCutsTable* theCoupleTable=
250      G4ProductionCutsTable::GetProductionCutsTable();
251
252    size_t index = couple->GetIndex();
253    G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
254    G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
255
256    // Generation of fluorescence
257    // Data in EADL are available only for Z > 5
258    // Protection to avoid generating photons in the unphysical case of
259    // shell binding energy > photon energy
260    if (bindingEnergy > cutg || bindingEnergy > cute)
261      {
262        G4DynamicParticle* aPhoton;
263        deexcitationManager.SetCutForSecondaryPhotons(cutg);
264        deexcitationManager.SetCutForAugerElectrons(cute);
265 
266        std::vector<G4DynamicParticle*>* photonVector = 
267          deexcitationManager.GenerateParticles(Z,shellId);
268        size_t nTotPhotons = photonVector->size();
269        for (size_t k=0; k<nTotPhotons; k++)
270          {
271            aPhoton = (*photonVector)[k];
272            if (aPhoton)
273              {
274                G4double itsEnergy = aPhoton->GetKineticEnergy();
275                if (itsEnergy <= bindingEnergy)
276                  {
277                    // Local energy deposit is given as the sum of the
278                    // energies of incident photons minus the energies
279                    // of the outcoming fluorescence photons
280                    bindingEnergy -= itsEnergy;
281                    fvect->push_back(aPhoton);
282                  }
283                else
284                  {
285                    // abnormal case of energy non-conservation
286                    delete aPhoton;
287                    (*photonVector)[k] = 0;
288                  }
289              }
290          }
291        delete photonVector;
292      }
293  }
294  // excitation energy left
295  fParticleChange->ProposeLocalEnergyDeposit(bindingEnergy);
296}
297
298//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299
300void G4LivermorePhotoElectricModel::SetCutForLowEnSecPhotons(G4double cut)
301{
302  cutForLowEnergySecondaryPhotons = cut;
303  deexcitationManager.SetCutForSecondaryPhotons(cut);
304}
305
306//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307
308void G4LivermorePhotoElectricModel::SetCutForLowEnSecElectrons(G4double cut)
309{
310  cutForLowEnergySecondaryElectrons = cut;
311  deexcitationManager.SetCutForAugerElectrons(cut);
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316void G4LivermorePhotoElectricModel::ActivateAuger(G4bool val)
317{
318  deexcitationManager.ActivateAugerElectronProduction(val);
319}
320
321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322
323void 
324G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* dist)
325{
326  ElectronAngularGenerator = dist;
327  ElectronAngularGenerator->PrintGeneratorInformation();
328}
329
330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
331
332void G4LivermorePhotoElectricModel::SetAngularGenerator(const G4String& name)
333{
334  if (name == "default") 
335    {
336      delete ElectronAngularGenerator;
337      ElectronAngularGenerator = 
338        new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
339      generatorName = name;
340    }
341  else if (name == "standard")
342    {
343      delete ElectronAngularGenerator;
344      ElectronAngularGenerator = 
345        new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
346      generatorName = name;
347    }
348  else if (name == "polarized")
349    {
350      delete ElectronAngularGenerator;
351      ElectronAngularGenerator = 
352        new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
353      generatorName = name;
354    }
355  else
356    {
357      G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator - generator does not exist");
358    }
359
360  ElectronAngularGenerator->PrintGeneratorInformation();
361}
362
Note: See TracBrowser for help on using the repository browser.