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

Last change on this file since 1005 was 1005, checked in by garnier, 17 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.