source: trunk/source/processes/electromagnetic/lowenergy/src/G4LowEnergyPhotoElectric.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 15.0 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//
27// $Id: G4LowEnergyPhotoElectric.cc,v 1.59 2009/06/11 15:47:08 mantero Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Author: A. Forti
31// Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
32//
33// History:
34// --------
35// October 1998 - low energy modifications by Alessandra Forti
36// Added Livermore data table construction methods A. Forti
37// Modified BuildMeanFreePath to read new data tables A. Forti
38// Added EnergySampling method A. Forti
39// Modified PostStepDoIt to insert sampling with EPDL97 data A. Forti
40// Added SelectRandomAtom A. Forti
41// Added map of the elements A. Forti
42// 10.04.2000 VL
43// - Correcting Fluorescence transition probabilities in order to take into account
44// non-radiative transitions. No Auger electron simulated yet: energy is locally deposited.
45// 17.02.2000 Veronique Lefebure
46// - bugs corrected in fluorescence simulation:
47// . when final use of binding energy: no photon was ever created
48// . no Fluorescence was simulated when the photo-electron energy
49// was below production threshold.
50//
51// 07-09-99, if no e- emitted: edep=photon energy, mma
52// 24.04.01 V.Ivanchenko remove RogueWave
53// 12.08.2001 MGP Revised according to a design iteration
54// 16.09.2001 E. Guardincerri Added fluorescence generation
55// 06.10.2001 MGP Added protection to avoid negative electron energies
56// when binding energy of selected shell > photon energy
57// 18.04.2001 V.Ivanchenko Fix problem with low energy gammas from fluorescence
58// MeanFreePath is calculated by crosSectionHandler directly
59// 31.05.2002 V.Ivanchenko Add path of Fluo + Auger cuts to AtomicDeexcitation
60// 14.06.2002 V.Ivanchenko By default do not cheak range of e-
61// 21.01.2003 V.Ivanchenko Cut per region
62// 10.05.2004 P.Rodrigues Changes to accommodate new angular generators
63// 20.01.2006 A.Trindade Changes to accommodate polarized angular generator
64//
65// --------------------------------------------------------------
66
67#include "G4LowEnergyPhotoElectric.hh"
68
69#include "G4VPhotoElectricAngularDistribution.hh"
70#include "G4PhotoElectricAngularGeneratorSimple.hh"
71#include "G4PhotoElectricAngularGeneratorSauterGavrila.hh"
72#include "G4PhotoElectricAngularGeneratorPolarized.hh"
73
74#include "G4ParticleDefinition.hh"
75#include "G4Track.hh"
76#include "G4Step.hh"
77#include "G4ForceCondition.hh"
78#include "G4Gamma.hh"
79#include "G4Electron.hh"
80#include "G4DynamicParticle.hh"
81#include "G4VParticleChange.hh"
82#include "G4ThreeVector.hh"
83#include "G4VCrossSectionHandler.hh"
84#include "G4CrossSectionHandler.hh"
85#include "G4VEMDataSet.hh"
86#include "G4CompositeEMDataSet.hh"
87#include "G4VDataSetAlgorithm.hh"
88#include "G4LogLogInterpolation.hh"
89#include "G4VRangeTest.hh"
90#include "G4RangeNoTest.hh"
91#include "G4AtomicTransitionManager.hh"
92#include "G4AtomicShell.hh"
93#include "G4ProductionCutsTable.hh"
94
95G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric(const G4String& processName)
96 : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
97 intrinsicLowEnergyLimit(10*eV),
98 intrinsicHighEnergyLimit(100*GeV),
99 cutForLowEnergySecondaryPhotons(250.*eV),
100 cutForLowEnergySecondaryElectrons(250.*eV)
101{
102 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
103 highEnergyLimit > intrinsicHighEnergyLimit)
104 {
105 G4Exception("G4LowEnergyPhotoElectric::G4LowEnergyPhotoElectric - energy limit outside intrinsic process validity range");
106 }
107
108 crossSectionHandler = new G4CrossSectionHandler();
109 shellCrossSectionHandler = new G4CrossSectionHandler();
110 meanFreePathTable = 0;
111 rangeTest = new G4RangeNoTest;
112 generatorName = "geant4.6.2";
113 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator"); // default generator
114
115
116 if (verboseLevel > 0)
117 {
118 G4cout << GetProcessName() << " is created " << G4endl
119 << "Energy range: "
120 << lowEnergyLimit / keV << " keV - "
121 << highEnergyLimit / GeV << " GeV"
122 << G4endl;
123 }
124
125 G4cout << G4endl;
126 G4cout << "*******************************************************************************" << G4endl;
127 G4cout << "*******************************************************************************" << G4endl;
128 G4cout << " The class G4LowEnergyPhotoElectric is NOT SUPPORTED ANYMORE. " << G4endl;
129 G4cout << " It will be REMOVED with the next major release of Geant4. " << G4endl;
130 G4cout << " Please consult: https://twiki.cern.ch/twiki/bin/view/Geant4/LoweProcesses" << G4endl;
131 G4cout << "*******************************************************************************" << G4endl;
132 G4cout << "*******************************************************************************" << G4endl;
133 G4cout << G4endl;
134}
135
136G4LowEnergyPhotoElectric::~G4LowEnergyPhotoElectric()
137{
138 delete crossSectionHandler;
139 delete shellCrossSectionHandler;
140 delete meanFreePathTable;
141 delete rangeTest;
142 delete ElectronAngularGenerator;
143}
144
145void G4LowEnergyPhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
146{
147
148 crossSectionHandler->Clear();
149 G4String crossSectionFile = "phot/pe-cs-";
150 crossSectionHandler->LoadData(crossSectionFile);
151
152 shellCrossSectionHandler->Clear();
153 G4String shellCrossSectionFile = "phot/pe-ss-cs-";
154 shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
155
156 delete meanFreePathTable;
157 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
158}
159
160G4VParticleChange* G4LowEnergyPhotoElectric::PostStepDoIt(const G4Track& aTrack,
161 const G4Step& aStep)
162{
163 // Fluorescence generated according to:
164 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
165 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
166 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
167
168 aParticleChange.Initialize(aTrack);
169
170 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
171 G4double photonEnergy = incidentPhoton->GetKineticEnergy();
172 if (photonEnergy <= lowEnergyLimit)
173 {
174 aParticleChange.ProposeTrackStatus(fStopAndKill);
175 aParticleChange.ProposeEnergy(0.);
176 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
177 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
178 }
179
180 G4ThreeVector photonDirection = incidentPhoton->GetMomentumDirection(); // Returns the normalized direction of the momentum
181
182 // Select randomly one element in the current material
183 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
184 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
185
186 // Select the ionised shell in the current atom according to shell cross sections
187 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
188
189 // Retrieve the corresponding identifier and binding energy of the selected shell
190 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
191 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
192 G4double bindingEnergy = shell->BindingEnergy();
193 G4int shellId = shell->ShellId();
194
195 // Create lists of pointers to DynamicParticles (photons and electrons)
196 // (Is the electron vector necessary? To be checked)
197 std::vector<G4DynamicParticle*>* photonVector = 0;
198 std::vector<G4DynamicParticle*> electronVector;
199
200 G4double energyDeposit = 0.0;
201
202 // Primary outcoming electron
203 G4double eKineticEnergy = photonEnergy - bindingEnergy;
204
205 // There may be cases where the binding energy of the selected shell is > photon energy
206 // In such cases do not generate secondaries
207 if (eKineticEnergy > 0.)
208 {
209 // Generate the electron only if with large enough range w.r.t. cuts and safety
210 G4double safety = aStep.GetPostStepPoint()->GetSafety();
211
212 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
213 {
214
215 // Calculate direction of the photoelectron
216 G4ThreeVector gammaPolarization = incidentPhoton->GetPolarization();
217 G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
218
219 // The electron is created ...
220 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
221 electronDirection,
222 eKineticEnergy);
223 electronVector.push_back(electron);
224 }
225 else
226 {
227 energyDeposit += eKineticEnergy;
228 }
229 }
230 else
231 {
232 bindingEnergy = photonEnergy;
233 }
234
235 G4int nElectrons = electronVector.size();
236 size_t nTotPhotons = 0;
237 G4int nPhotons=0;
238 const G4ProductionCutsTable* theCoupleTable=
239 G4ProductionCutsTable::GetProductionCutsTable();
240
241 size_t index = couple->GetIndex();
242 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
243 cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
244
245 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
246 cute = std::min(cutForLowEnergySecondaryPhotons,cute);
247
248 G4DynamicParticle* aPhoton;
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 (Z > 5 && (bindingEnergy > cutg || bindingEnergy > cute))
255 {
256 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
257 nTotPhotons = photonVector->size();
258 for (size_t k=0; k<nTotPhotons; k++)
259 {
260 aPhoton = (*photonVector)[k];
261 if (aPhoton)
262 {
263 G4double itsCut = cutg;
264 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
265 G4double itsEnergy = aPhoton->GetKineticEnergy();
266
267 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
268 {
269 nPhotons++;
270 // Local energy deposit is given as the sum of the
271 // energies of incident photons minus the energies
272 // of the outcoming fluorescence photons
273 bindingEnergy -= itsEnergy;
274
275 }
276 else
277 {
278 delete aPhoton;
279 (*photonVector)[k] = 0;
280 }
281 }
282 }
283 }
284
285 energyDeposit += bindingEnergy;
286
287 G4int nSecondaries = nElectrons + nPhotons;
288 aParticleChange.SetNumberOfSecondaries(nSecondaries);
289
290 for (G4int l = 0; l<nElectrons; l++ )
291 {
292 aPhoton = electronVector[l];
293 if(aPhoton) {
294 aParticleChange.AddSecondary(aPhoton);
295 }
296 }
297 for ( size_t ll = 0; ll < nTotPhotons; ll++)
298 {
299 aPhoton = (*photonVector)[ll];
300 if(aPhoton) {
301 aParticleChange.AddSecondary(aPhoton);
302 }
303 }
304
305 delete photonVector;
306
307 if (energyDeposit < 0)
308 {
309 G4cout << "WARNING - "
310 << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
311 << G4endl;
312 energyDeposit = 0;
313 }
314
315 // Kill the incident photon
316 aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
317 aParticleChange.ProposeEnergy( 0. );
318
319 aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
320 aParticleChange.ProposeTrackStatus( fStopAndKill );
321
322 // Reset NbOfInteractionLengthLeft and return aParticleChange
323 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
324}
325
326G4bool G4LowEnergyPhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
327{
328 return ( &particle == G4Gamma::Gamma() );
329}
330
331G4double G4LowEnergyPhotoElectric::GetMeanFreePath(const G4Track& track,
332 G4double, // previousStepSize
333 G4ForceCondition*)
334{
335 const G4DynamicParticle* photon = track.GetDynamicParticle();
336 G4double energy = photon->GetKineticEnergy();
337 G4Material* material = track.GetMaterial();
338 // size_t materialIndex = material->GetIndex();
339
340 G4double meanFreePath = DBL_MAX;
341
342 // if (energy > highEnergyLimit)
343 // meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
344 // else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
345 // else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
346
347 G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
348 if(cross > 0.0) meanFreePath = 1.0/cross;
349
350 return meanFreePath;
351}
352
353void G4LowEnergyPhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
354{
355 cutForLowEnergySecondaryPhotons = cut;
356 deexcitationManager.SetCutForSecondaryPhotons(cut);
357}
358
359void G4LowEnergyPhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
360{
361 cutForLowEnergySecondaryElectrons = cut;
362 deexcitationManager.SetCutForAugerElectrons(cut);
363}
364
365void G4LowEnergyPhotoElectric::ActivateAuger(G4bool val)
366{
367 deexcitationManager.ActivateAugerElectronProduction(val);
368}
369
370void G4LowEnergyPhotoElectric::SetAngularGenerator(G4VPhotoElectricAngularDistribution* distribution)
371{
372 ElectronAngularGenerator = distribution;
373 ElectronAngularGenerator->PrintGeneratorInformation();
374}
375
376void G4LowEnergyPhotoElectric::SetAngularGenerator(const G4String& name)
377{
378 if (name == "default")
379 {
380 delete ElectronAngularGenerator;
381 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
382 generatorName = name;
383 }
384 else if (name == "standard")
385 {
386 delete ElectronAngularGenerator;
387 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
388 generatorName = name;
389 }
390 else if (name == "polarized")
391 {
392 delete ElectronAngularGenerator;
393 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
394 generatorName = name;
395 }
396 else
397 {
398 G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator - generator does not exist");
399 }
400
401 ElectronAngularGenerator->PrintGeneratorInformation();
402}
Note: See TracBrowser for help on using the repository browser.