source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectric.cc@ 830

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 12.1 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: G4PenelopePhotoElectric.cc,v 1.12 2006/06/29 19:40:51 gunter Exp $
28// GEANT4 tag $Name: $
29//
30// Author: L. Pandola
31//
32// History:
33// --------
34// January 2003 - Created
35// 12 Feb 2003 MG Pia Migration to "cuts per region"
36// 10 Mar 2003 V.Ivanchenko Remome CutPerMaterial warning
37// 31 May 2005 L. Pandola Added Sauter formula for the sampling of
38// the electron direction
39// --------------------------------------------------------------
40
41#include "G4PenelopePhotoElectric.hh"
42
43#include "G4ParticleDefinition.hh"
44#include "G4Track.hh"
45#include "G4Step.hh"
46#include "G4ForceCondition.hh"
47#include "G4Gamma.hh"
48#include "G4Electron.hh"
49#include "G4DynamicParticle.hh"
50#include "G4VParticleChange.hh"
51#include "G4ThreeVector.hh"
52#include "G4VCrossSectionHandler.hh"
53#include "G4CrossSectionHandler.hh"
54#include "G4VEMDataSet.hh"
55#include "G4CompositeEMDataSet.hh"
56#include "G4VDataSetAlgorithm.hh"
57#include "G4LogLogInterpolation.hh"
58#include "G4VRangeTest.hh"
59#include "G4RangeNoTest.hh"
60#include "G4AtomicTransitionManager.hh"
61#include "G4AtomicShell.hh"
62#include "G4MaterialCutsCouple.hh"
63#include "G4ProductionCutsTable.hh"
64
65G4PenelopePhotoElectric::G4PenelopePhotoElectric(const G4String& processName)
66 : G4VDiscreteProcess(processName), lowEnergyLimit(250*eV), highEnergyLimit(100*GeV),
67 intrinsicLowEnergyLimit(10*eV),
68 intrinsicHighEnergyLimit(100*GeV),
69 cutForLowEnergySecondaryPhotons(250.*eV),
70 cutForLowEnergySecondaryElectrons(250.*eV)
71{
72 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
73 highEnergyLimit > intrinsicHighEnergyLimit)
74 {
75 G4Exception("G4PenelopePhotoElectric::G4PenelopePhotoElectric - energy limit outside intrinsic process validity range");
76 }
77
78 crossSectionHandler = new G4CrossSectionHandler();
79 shellCrossSectionHandler = new G4CrossSectionHandler();
80 meanFreePathTable = 0;
81 rangeTest = new G4RangeNoTest;
82
83 if (verboseLevel > 0)
84 {
85 G4cout << GetProcessName() << " is created " << G4endl
86 << "Energy range: "
87 << lowEnergyLimit / keV << " keV - "
88 << highEnergyLimit / GeV << " GeV"
89 << G4endl;
90 }
91}
92
93G4PenelopePhotoElectric::~G4PenelopePhotoElectric()
94{
95 delete crossSectionHandler;
96 delete shellCrossSectionHandler;
97 delete meanFreePathTable;
98 delete rangeTest;
99}
100
101void G4PenelopePhotoElectric::BuildPhysicsTable(const G4ParticleDefinition& )
102{
103
104 crossSectionHandler->Clear();
105 G4String crossSectionFile = "penelope/ph-cs-pen-";
106 crossSectionHandler->LoadData(crossSectionFile);
107
108 shellCrossSectionHandler->Clear();
109 G4String shellCrossSectionFile = "penelope/ph-ss-cs-pen-";
110 shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
111
112 delete meanFreePathTable;
113 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
114}
115
116G4VParticleChange* G4PenelopePhotoElectric::PostStepDoIt(const G4Track& aTrack,
117 const G4Step& aStep)
118{
119 // Fluorescence generated according to:
120 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
121 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
122 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
123
124 aParticleChange.Initialize(aTrack);
125
126 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
127 G4double photonEnergy = incidentPhoton->GetKineticEnergy();
128 if (photonEnergy <= lowEnergyLimit)
129 {
130 aParticleChange.ProposeTrackStatus(fStopAndKill);
131 aParticleChange.ProposeEnergy(0.);
132 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy);
133 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
134 }
135
136 G4ParticleMomentum photonDirection = incidentPhoton->GetMomentumDirection();
137
138 // Select randomly one element in the current material
139 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
140
141 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
142
143 // Select the ionised shell in the current atom according to shell cross sections
144 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
145
146 // Retrieve the corresponding identifier and binding energy of the selected shell
147 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
148 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
149 G4double bindingEnergy = shell->BindingEnergy();
150 G4int shellId = shell->ShellId();
151
152 // Create lists of pointers to DynamicParticles (photons and electrons)
153 // (Is the electron vector necessary? To be checked)
154 std::vector<G4DynamicParticle*>* photonVector = 0;
155 std::vector<G4DynamicParticle*> electronVector;
156
157 G4double energyDeposit = 0.0;
158
159 // Primary outcoming electron
160 G4double eKineticEnergy = photonEnergy - bindingEnergy;
161
162 // There may be cases where the binding energy of the selected shell is > photon energy
163 // In such cases do not generate secondaries
164 if (eKineticEnergy > 0.)
165 {
166 // Generate the electron only if with large enough range w.r.t. cuts and safety
167 G4double safety = aStep.GetPostStepPoint()->GetSafety();
168
169 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
170 {
171 // The electron is created
172 // Direction sampled from the Sauter distribution
173 G4double cosTheta = SampleElectronDirection(eKineticEnergy);
174 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
175 G4double phi = twopi * G4UniformRand() ;
176 G4double dirx = sinTheta * std::cos(phi);
177 G4double diry = sinTheta * std::sin(phi);
178 G4double dirz = cosTheta ;
179 G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
180 electronDirection.rotateUz(photonDirection);
181
182 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
183 electronDirection,
184 eKineticEnergy);
185 electronVector.push_back(electron);
186 }
187 else
188 {
189 energyDeposit += eKineticEnergy;
190 }
191 }
192 else
193 {
194 bindingEnergy = photonEnergy;
195 }
196
197 G4int nElectrons = electronVector.size();
198 size_t nTotPhotons = 0;
199 G4int nPhotons=0;
200
201 const G4ProductionCutsTable* theCoupleTable=
202 G4ProductionCutsTable::GetProductionCutsTable();
203
204 size_t index = couple->GetIndex();
205 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
206 cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
207
208 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
209 cute = std::min(cutForLowEnergySecondaryPhotons,cute);
210
211 G4DynamicParticle* aPhoton;
212
213 // Generation of fluorescence
214 // Data in EADL are available only for Z > 5
215 // Protection to avoid generating photons in the unphysical case of
216 // shell binding energy > photon energy
217 if (Z > 5 && (bindingEnergy > cutg || bindingEnergy > cute))
218 {
219 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
220 nTotPhotons = photonVector->size();
221 for (size_t k=0; k<nTotPhotons; k++)
222 {
223 aPhoton = (*photonVector)[k];
224 if (aPhoton)
225 {
226 G4double itsCut = cutg;
227 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
228 G4double itsEnergy = aPhoton->GetKineticEnergy();
229
230 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
231 {
232 nPhotons++;
233 // Local energy deposit is given as the sum of the
234 // energies of incident photons minus the energies
235 // of the outcoming fluorescence photons
236 bindingEnergy -= itsEnergy;
237
238 }
239 else
240 {
241 delete aPhoton;
242 (*photonVector)[k] = 0;
243 }
244 }
245 }
246 }
247
248 energyDeposit += bindingEnergy;
249
250 G4int nSecondaries = nElectrons + nPhotons;
251 aParticleChange.SetNumberOfSecondaries(nSecondaries);
252
253 for (G4int l = 0; l<nElectrons; l++ )
254 {
255 aPhoton = electronVector[l];
256 if(aPhoton) {
257 aParticleChange.AddSecondary(aPhoton);
258 }
259 }
260 for ( size_t ll = 0; ll < nTotPhotons; ll++)
261 {
262 aPhoton = (*photonVector)[ll];
263 if(aPhoton) {
264 aParticleChange.AddSecondary(aPhoton);
265 }
266 }
267
268 delete photonVector;
269
270 if (energyDeposit < 0)
271 {
272 G4cout << "WARNING - "
273 << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
274 << G4endl;
275 energyDeposit = 0;
276 }
277
278 // Kill the incident photon
279 aParticleChange.ProposeMomentumDirection( 0., 0., 0. );
280 aParticleChange.ProposeEnergy( 0. );
281
282 aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
283 aParticleChange.ProposeTrackStatus( fStopAndKill );
284
285 // Reset NbOfInteractionLengthLeft and return aParticleChange
286 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep );
287}
288
289G4bool G4PenelopePhotoElectric::IsApplicable(const G4ParticleDefinition& particle)
290{
291 return ( &particle == G4Gamma::Gamma() );
292}
293
294G4double G4PenelopePhotoElectric::GetMeanFreePath(const G4Track& track,
295 G4double, // previousStepSize
296 G4ForceCondition*)
297{
298 const G4DynamicParticle* photon = track.GetDynamicParticle();
299 G4double energy = photon->GetKineticEnergy();
300 G4Material* material = track.GetMaterial();
301
302 G4double meanFreePath = DBL_MAX;
303
304 G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
305 if(cross > 0.0) meanFreePath = 1.0/cross;
306
307 return meanFreePath;
308}
309
310void G4PenelopePhotoElectric::SetCutForLowEnSecPhotons(G4double cut)
311{
312 cutForLowEnergySecondaryPhotons = cut;
313 deexcitationManager.SetCutForSecondaryPhotons(cut);
314}
315
316void G4PenelopePhotoElectric::SetCutForLowEnSecElectrons(G4double cut)
317{
318 cutForLowEnergySecondaryElectrons = cut;
319 deexcitationManager.SetCutForAugerElectrons(cut);
320}
321
322void G4PenelopePhotoElectric::ActivateAuger(G4bool val)
323{
324 deexcitationManager.ActivateAugerElectronProduction(val);
325}
326
327
328G4double G4PenelopePhotoElectric::SampleElectronDirection(G4double energy)
329{
330 G4double costheta = 1.0;
331 if (energy>1*GeV) return costheta;
332
333 //1) initialize energy-dependent variables
334 // Variable naming according to Eq. (2.24) of Penelope Manual
335 // (pag. 44)
336 G4double gamma = 1.0 + energy/electron_mass_c2;
337 G4double gamma2 = gamma*gamma;
338 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
339
340 // ac corresponds to "A" of Eq. (2.31)
341 //
342 G4double ac = (1.0/beta) - 1.0;
343 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
344 G4double a2 = ac + 2.0;
345 G4double gtmax = 2.0*(a1 + 1.0/ac);
346
347 G4double tsam = 0;
348 G4double gtr = 0;
349
350 //2) sampling. Eq. (2.31) of Penelope Manual
351 // tsam = 1-std::cos(theta)
352 // gtr = rejection function according to Eq. (2.28)
353 do{
354 G4double rand = G4UniformRand();
355 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
356 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
357 }while(G4UniformRand()*gtmax > gtr);
358 costheta = 1.0-tsam;
359 return costheta;
360}
361
362
363
364
365
Note: See TracBrowser for help on using the repository browser.