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

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

tag geant4.9.4 beta 1 + modifs locales
File size: 16.1 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//
[1315]26// $Id: G4PenelopePhotoElectricModel.cc,v 1.12 2010/03/26 09:32:50 pandola Exp $
[1337]27// GEANT4 tag $Name: geant4-09-04-beta-01 $
[968]28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
[1055]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
[1192]46// 21 Oct 2009 L Pandola Remove un-necessary fUseAtomicDeexcitation flag - now managed by
47// G4VEmModel::DeexcitationFlag()
[1315]48// 15 Mar 2010 L Pandola Explicitely initialize Auger to false
[968]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;
[1055]75 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
[968]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
[1315]85
86 //by default the model will inkove the atomic deexcitation
87 SetDeexcitationFlag(true);
88 ActivateAuger(false);
[968]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;
[1055]110 crossSectionHandler = 0;
[968]111 }
112 if (shellCrossSectionHandler)
113 {
114 shellCrossSectionHandler->Clear();
115 delete shellCrossSectionHandler;
[1055]116 shellCrossSectionHandler =0;
[968]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
[1055]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 }
[968]141
142 if(isInitialised) return;
[1055]143 fParticleChange = GetParticleChangeForGamma();
[968]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;
[1055]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 // }
[968]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
[1055]206 // always kill primary
207 fParticleChange->ProposeTrackStatus(fStopAndKill);
208 fParticleChange->SetProposedKineticEnergy(0.);
209
210 if (photonEnergy <= fIntrinsicLowEnergyLimit)
211 {
[968]212 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
213 return ;
[1055]214 }
[968]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
[1055]222 //the dimension of the table is equal to the dimension of the database
223 //(less interpolation errors)
[968]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();
[1055]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
[968]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.)
[1315]256 {
[1055]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 }
[968]272 else
[1055]273 {
[968]274 bindingEnergy = photonEnergy;
[1055]275 }
[968]276 G4double energyInFluorescence = 0; //testing purposes
277
278 //Now, take care of fluorescence, if required
[1055]279 if(DeexcitationFlag() && Z > 5)
[968]280 {
[1055]281 const G4ProductionCutsTable* theCoupleTable=
282 G4ProductionCutsTable::GetProductionCutsTable();
283 size_t indx = couple->GetIndex();
[968]284 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
[1055]285 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
[968]286
287 // Protection to avoid generating photons in the unphysical case of
288 // shell binding energy > photon energy
[1055]289 if (bindingEnergy > cutG || bindingEnergy > cutE)
[968]290 {
[1055]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)
[968]297 {
[1055]298 for (size_t k=0; k< photonVector->size(); k++)
[968]299 {
[1055]300 G4DynamicParticle* aPhoton = (*photonVector)[k];
301 if (aPhoton)
[968]302 {
[1055]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 }
[968]316 }
317 }
[1055]318 delete photonVector;
[968]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 {
[1055]351 G4double energyDiff =
352 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
[968]353 if (energyDiff > 0.05*keV)
354 G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
[1055]355 (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV
356 << " keV (final) vs. " <<
[968]357 photonEnergy/keV << " keV (initial)" << G4endl;
358 }
359}
360
361//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
362
363void G4PenelopePhotoElectricModel::ActivateAuger(G4bool augerbool)
364{
[1192]365 if (!DeexcitationFlag() && augerbool)
[968]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.