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

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

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