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

Last change on this file since 973 was 961, checked in by garnier, 17 years ago

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