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

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

update geant4.9.3 tag
File size: 12.7 KB
RevLine 
[968]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//
[1192]26// $Id: G4LivermorePhotoElectricModel.cc,v 1.9 2009/10/23 09:31:03 pandola Exp $
[1228]27// GEANT4 tag $Name: geant4-09-03 $
[968]28//
[1055]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
[1192]42// 23 Oct 2009 L Pandola
43// - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
44// set as "true" (default would be false)
45//
[968]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)
[1055]57:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),
58 crossSectionHandler(0),shellCrossSectionHandler(0),ElectronAngularGenerator(0)
[968]59{
[1055]60 lowEnergyLimit = 250 * eV;
[968]61 highEnergyLimit = 100 * GeV;
[1055]62 // SetLowEnergyLimit(lowEnergyLimit);
[968]63 SetHighEnergyLimit(highEnergyLimit);
[1055]64
[1192]65 //Set atomic deexcitation by default
66 SetDeexcitationFlag(true);
[1055]67 ActivateAuger(false);
68
[968]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
[1055]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 }
[968]83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
88{
[1055]89 if (crossSectionHandler) delete crossSectionHandler;
90 if (shellCrossSectionHandler) delete shellCrossSectionHandler;
91 if (ElectronAngularGenerator) delete ElectronAngularGenerator;
[968]92}
93
94//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95
[1055]96void
97G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
98 const G4DataVector&)
[968]99{
100 if (verboseLevel > 3)
101 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
102
[1055]103 if (crossSectionHandler)
[968]104 {
[1055]105 crossSectionHandler->Clear();
106 delete crossSectionHandler;
[968]107 }
[1055]108
109 if (shellCrossSectionHandler)
[968]110 {
[1055]111 shellCrossSectionHandler->Clear();
112 delete shellCrossSectionHandler;
[968]113 }
[1055]114
[968]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
[1055]127 // SI - Simple generator is buggy
[968]128 //generatorName = "geant4.6.2";
129 //ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator"); // default generator
[1055]130 ElectronAngularGenerator =
131 new G4PhotoElectricAngularGeneratorSauterGavrila("GEANTSauterGavrilaGenerator");
[968]132
[1055]133 //
[968]134 if (verboseLevel > 2)
135 G4cout << "Loaded cross section files for Livermore PhotoElectric model" << G4endl;
136
[1055]137 // InitialiseElementSelectors(particle,cuts);
[968]138
[1055]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
[968]147 if(isInitialised) return;
[1055]148 fParticleChange = GetParticleChangeForGamma();
[968]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)
[1055]161 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel"
162 << G4endl;
[968]163
[1055]164 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
165 return 0;
166
[968]167 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
168 return cs;
169}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172
[1055]173void
174G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
175 const G4MaterialCutsCouple* couple,
176 const G4DynamicParticle* aDynamicGamma,
177 G4double,
178 G4double)
[968]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
[1055]191 // kill incident photon
192 fParticleChange->SetProposedKineticEnergy(0.);
193 fParticleChange->ProposeTrackStatus(fStopAndKill);
194
195 // low-energy gamma is absorpted by this process
[968]196 if (photonEnergy <= lowEnergyLimit)
197 {
198 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
[1055]199 return;
[968]200 }
201
[1055]202 // Returns the normalized direction of the momentum
203 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
[968]204
205 // Select randomly one element in the current material
[1055]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();
[968]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 {
[1055]227 // Calculate direction of the photoelectron
228 G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
229 G4ThreeVector electronDirection =
230 ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,
231 eKineticEnergy,
232 gammaPolarization,
233 shellId);
[968]234
[1055]235 // The electron is created ...
236 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
237 electronDirection,
238 eKineticEnergy);
239 fvect->push_back(electron);
[968]240 }
241 else
242 {
243 bindingEnergy = photonEnergy;
244 }
245
[1055]246 // deexcitation
247 if(DeexcitationFlag() && Z > 5) {
[968]248
[1055]249 const G4ProductionCutsTable* theCoupleTable=
250 G4ProductionCutsTable::GetProductionCutsTable();
[968]251
[1055]252 size_t index = couple->GetIndex();
253 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
254 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
[968]255
[1055]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;
[968]292 }
[1055]293 }
294 // excitation energy left
295 fParticleChange->ProposeLocalEnergyDeposit(bindingEnergy);
[968]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
[1055]323void
324G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* dist)
[968]325{
[1055]326 ElectronAngularGenerator = dist;
[968]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;
[1055]337 ElectronAngularGenerator =
338 new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
[968]339 generatorName = name;
340 }
341 else if (name == "standard")
342 {
343 delete ElectronAngularGenerator;
[1055]344 ElectronAngularGenerator =
345 new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
[968]346 generatorName = name;
347 }
348 else if (name == "polarized")
349 {
350 delete ElectronAngularGenerator;
[1055]351 ElectronAngularGenerator =
352 new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
[968]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.