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