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

Last change on this file since 992 was 991, checked in by garnier, 17 years ago

update
File size: 15.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//
[991]26// $Id: G4PenelopePhotoElectricModel.cc,v 1.2 2008/12/04 14:09:36 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
[991]33// 08 Oct 2008 L Pandola Migration from process to model
[968]34//
35
36#include "G4PenelopePhotoElectricModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4PhysicsTable.hh"
42#include "G4ElementTable.hh"
43#include "G4Element.hh"
44#include "G4CrossSectionHandler.hh"
45#include "G4AtomicTransitionManager.hh"
46#include "G4AtomicShell.hh"
47#include "G4Gamma.hh"
48#include "G4Electron.hh"
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52
53G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition*,
54 const G4String& nam)
55 :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
56 shellCrossSectionHandler(0)
57{
58 fIntrinsicLowEnergyLimit = 100.0*eV;
59 fIntrinsicHighEnergyLimit = 100.0*GeV;
60 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
61 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
62 //
63 fUseAtomicDeexcitation = true;
64 verboseLevel= 0;
65 // Verbosity scale:
66 // 0 = nothing
67 // 1 = warning for energy non-conservation
68 // 2 = details of energy budget
69 // 3 = calculation of cross sections, file openings, sampling of atoms
70 // 4 = entering in methods
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74
75G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel()
76{
77 if (crossSectionHandler) delete crossSectionHandler;
78 if (shellCrossSectionHandler) delete shellCrossSectionHandler;
79}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition*,
84 const G4DataVector& )
85{
86 if (verboseLevel > 3)
87 G4cout << "Calling G4PenelopePhotoElectricModel::Initialise()" << G4endl;
88 if (crossSectionHandler)
89 {
90 crossSectionHandler->Clear();
91 delete crossSectionHandler;
92 }
93 if (shellCrossSectionHandler)
94 {
95 shellCrossSectionHandler->Clear();
96 delete shellCrossSectionHandler;
97 }
98
99 //Check energy limits
100 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
101 {
102 G4cout << "G4PenelopePhotoElectricModel: low energy limit increased from " <<
103 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" <<
104 G4endl;
105 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
106 }
107 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
108 {
109 G4cout << "G4PenelopePhotoElectricModel: high energy limit decreased from " <<
110 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV"
111 << G4endl;
112 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
113 }
114
115
116 //Re-initialize cross section handlers
117 crossSectionHandler = new G4CrossSectionHandler();
118 crossSectionHandler->Clear();
119 G4String crossSectionFile = "penelope/ph-cs-pen-";
120 crossSectionHandler->LoadData(crossSectionFile);
121 shellCrossSectionHandler = new G4CrossSectionHandler();
122 shellCrossSectionHandler->Clear();
123 crossSectionFile = "penelope/ph-ss-cs-pen-";
124 shellCrossSectionHandler->LoadShellData(crossSectionFile);
125 //This is used to retrieve cross section values later on
126 crossSectionHandler->BuildMeanFreePathForMaterials();
127
128 if (verboseLevel > 2)
129 G4cout << "Loaded cross section files for PenelopePhotoElectric" << G4endl;
130
131 G4cout << "Penelope Photo-Electric model is initialized " << G4endl
132 << "Energy range: "
133 << LowEnergyLimit() / MeV << " MeV - "
134 << HighEnergyLimit() / GeV << " GeV"
135 << G4endl;
136
137 if(isInitialised) return;
138
139 if(pParticleChange)
140 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
141 else
142 fParticleChange = new G4ParticleChangeForGamma();
143 isInitialised = true;
144}
145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom(
149 const G4ParticleDefinition*,
150 G4double energy,
151 G4double Z, G4double,
152 G4double, G4double)
153{
154 //
155 // Penelope model. Use data-driven approach for cross section estimate (and
156 // also shell sampling from a given atom). Data are from the Livermore database
157 // D.E. Cullen et al., Report UCRL-50400 (1989)
158 //
159
160 if (verboseLevel > 3)
161 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
162
163 G4int iZ = (G4int) Z;
164 if (!crossSectionHandler)
165 {
166 G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
167 G4cout << "The cross section handler is not correctly initialized" << G4endl;
168 G4Exception();
169 }
170 G4double cs = crossSectionHandler->FindValue(iZ,energy);
171
172 if (verboseLevel > 2)
173 G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
174 " = " << cs/barn << " barn" << G4endl;
175 return cs;
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
181 const G4MaterialCutsCouple* couple,
182 const G4DynamicParticle* aDynamicGamma,
183 G4double,
184 G4double)
185{
186 //
187 // Photoelectric effect, Penelope model
188 //
189 // The target atom and the target shell are sampled according to the Livermore
190 // database
191 // D.E. Cullen et al., Report UCRL-50400 (1989)
192 // The angular distribution of the electron in the final state is sampled
193 // according to the Sauter distribution from
194 // F. Sauter, Ann. Phys. 11 (1931) 454
195 // The energy of the final electron is given by the initial photon energy minus
196 // the binding energy. Fluorescence de-excitation is subsequently produced
197 // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
198 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
199
200 if (verboseLevel > 3)
201 G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
202
203 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
204
205 if (photonEnergy <= LowEnergyLimit())
206 {
207 fParticleChange->ProposeTrackStatus(fStopAndKill);
208 fParticleChange->SetProposedKineticEnergy(0.);
209 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
210 return ;
211 }
212
213 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
214
215 // Select randomly one element in the current material
216 if (verboseLevel > 2)
217 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
218 //use crossSectionHandler instead of G4EmElementSelector because in this case
219 //the dimension of the table is equal to the dimension of the database (less interpolation errors)
220 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
221 if (verboseLevel > 2)
222 G4cout << "Selected Z = " << Z << G4endl;
223
224 // Select the ionised shell in the current atom according to shell cross sections
225 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
226
227 // Retrieve the corresponding identifier and binding energy of the selected shell
228 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
229 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
230 G4double bindingEnergy = shell->BindingEnergy();
231 G4int shellId = shell->ShellId();
232
233 G4double localEnergyDeposit = 0.0;
234
235 // Primary outcoming electron
236 G4double eKineticEnergy = photonEnergy - bindingEnergy;
237
238 G4double cutForLowEnergySecondaryParticles = 250.0*eV;
239 const G4ProductionCutsTable* theCoupleTable=
240 G4ProductionCutsTable::GetProductionCutsTable();
241 size_t indx = couple->GetIndex();
242 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
243 cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
244
245 // There may be cases where the binding energy of the selected shell is > photon energy
246 // In such cases do not generate secondaries
247 if (eKineticEnergy > 0.)
248 {
249 //Now check if the electron is above cuts: if so, it is created explicitely
250 if (eKineticEnergy > cutE)
251 {
252 // The electron is created
253 // Direction sampled from the Sauter distribution
254 G4double cosTheta = SampleElectronDirection(eKineticEnergy);
255 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
256 G4double phi = twopi * G4UniformRand() ;
257 G4double dirx = sinTheta * std::cos(phi);
258 G4double diry = sinTheta * std::sin(phi);
259 G4double dirz = cosTheta ;
260 G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
261 electronDirection.rotateUz(photonDirection);
262 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
263 electronDirection,
264 eKineticEnergy);
265 fvect->push_back(electron);
266 }
267 else
268 localEnergyDeposit += eKineticEnergy;
269 }
270 else
271 bindingEnergy = photonEnergy;
272
273
274 G4double energyInFluorescence = 0; //testing purposes
275
276 //Now, take care of fluorescence, if required
277 if (fUseAtomicDeexcitation)
278 {
279 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
280 cutG = std::min(cutForLowEnergySecondaryParticles,cutG);
281
282 std::vector<G4DynamicParticle*>* photonVector = 0;
283
284 // Protection to avoid generating photons in the unphysical case of
285 // shell binding energy > photon energy
286 if (Z > 5 && (bindingEnergy > cutG || bindingEnergy > cutE))
287 {
288 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
289 //Check for single photons (if they are above threshold)
290 for (size_t k=0; k< photonVector->size(); k++)
291 {
292 G4DynamicParticle* aPhoton = (*photonVector)[k];
293 if (aPhoton)
294 {
295 G4double itsCut = cutG;
296 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
297 G4double itsEnergy = aPhoton->GetKineticEnergy();
298
299 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
300 {
301 // Local energy deposit is given as the sum of the
302 // energies of incident photons minus the energies
303 // of the outcoming fluorescence photons
304 bindingEnergy -= itsEnergy;
305
306 }
307 else
308 {
309 (*photonVector)[k] = 0;
310 }
311 }
312 }
313 }
314 //Register valid secondaries
315 if (photonVector)
316 {
317 for ( size_t ll = 0; ll <photonVector->size(); ll++)
318 {
319 G4DynamicParticle* aPhoton = (*photonVector)[ll];
320 if (aPhoton)
321 {
322 energyInFluorescence += aPhoton->GetKineticEnergy();
323 fvect->push_back(aPhoton);
324 }
325 }
326 delete photonVector;
327 }
328 }
329 //Residual energy is deposited locally
330 localEnergyDeposit += bindingEnergy;
331
332
333 if (localEnergyDeposit < 0)
334 {
335 G4cout << "WARNING - "
336 << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
337 << G4endl;
338 localEnergyDeposit = 0;
339 }
340
341 //Update the status of the primary gamma (kill it)
342 fParticleChange->SetProposedKineticEnergy(0.);
343 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
344 fParticleChange->ProposeTrackStatus(fStopAndKill);
345
346 if (verboseLevel > 1)
347 {
348 G4cout << "-----------------------------------------------------------" << G4endl;
349 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
350 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
351 G4cout << "-----------------------------------------------------------" << G4endl;
352 if (eKineticEnergy)
353 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
354 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
355 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
356 G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV <<
357 " keV" << G4endl;
358 G4cout << "-----------------------------------------------------------" << G4endl;
359 }
360 if (verboseLevel > 0)
361 {
362 G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
363 if (energyDiff > 0.05*keV)
364 G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
365 (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << " keV (final) vs. " <<
366 photonEnergy/keV << " keV (initial)" << G4endl;
367 }
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371
372void G4PenelopePhotoElectricModel::ActivateAuger(G4bool augerbool)
373{
374 if (!fUseAtomicDeexcitation)
375 {
376 G4cout << "WARNING - G4PenelopePhotoElectricModel" << G4endl;
377 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
378 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
379 }
380 deexcitationManager.ActivateAugerElectronProduction(augerbool);
381 if (verboseLevel > 1)
382 G4cout << "Auger production set to " << augerbool << G4endl;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy)
388{
389 G4double costheta = 1.0;
390 if (energy>1*GeV) return costheta;
391
392 //1) initialize energy-dependent variables
393 // Variable naming according to Eq. (2.24) of Penelope Manual
394 // (pag. 44)
395 G4double gamma = 1.0 + energy/electron_mass_c2;
396 G4double gamma2 = gamma*gamma;
397 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
398
399 // ac corresponds to "A" of Eq. (2.31)
400 //
401 G4double ac = (1.0/beta) - 1.0;
402 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
403 G4double a2 = ac + 2.0;
404 G4double gtmax = 2.0*(a1 + 1.0/ac);
405
406 G4double tsam = 0;
407 G4double gtr = 0;
408
409 //2) sampling. Eq. (2.31) of Penelope Manual
410 // tsam = 1-std::cos(theta)
411 // gtr = rejection function according to Eq. (2.28)
412 do{
413 G4double rand = G4UniformRand();
414 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
415 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
416 }while(G4UniformRand()*gtmax > gtr);
417 costheta = 1.0-tsam;
418 return costheta;
419}
420
Note: See TracBrowser for help on using the repository browser.