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

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

update
File size: 14.1 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: G4LivermorePhotoElectricModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29
30#include "G4LivermorePhotoElectricModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4ParticleDefinition*,
39 const G4String& nam)
[991]40:G4VEmModel(nam),isInitialised(false)
[968]41{
42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
43 highEnergyLimit = 100 * GeV;
44 SetLowEnergyLimit(lowEnergyLimit);
45 SetHighEnergyLimit(highEnergyLimit);
46
47 G4cout << "Livermore Compton is constructed " << G4endl
48 << "Energy range: "
49 << lowEnergyLimit / keV << " keV - "
50 << highEnergyLimit / GeV << " GeV"
51 << G4endl;
52
53 verboseLevel= 0;
54 // Verbosity scale:
55 // 0 = nothing
56 // 1 = warning for energy non-conservation
57 // 2 = details of energy budget
58 // 3 = calculation of cross sections, file openings, sampling of atoms
59 // 4 = entering in methods
60}
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
65{
[991]66 delete meanFreePathTable;
67 delete crossSectionHandler;
68 delete shellCrossSectionHandler;
69 delete ElectronAngularGenerator;
[968]70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74void G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
75 const G4DataVector& cuts)
76{
77 if (verboseLevel > 3)
78 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
79
[991]80 InitialiseElementSelectors(particle,cuts);
81
[968]82 // Energy limits
83
84 if (LowEnergyLimit() < lowEnergyLimit)
85 {
86 G4cout << "G4LivermorePhotoElectricModel: low energy limit increased from " <<
87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" <<
88 G4endl;
89 SetLowEnergyLimit(lowEnergyLimit);
90 }
91
92 if (HighEnergyLimit() > highEnergyLimit)
93 {
94 G4cout << "G4LivermorePhotoElectricModel: high energy limit decreased from " <<
95 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV"
96 << G4endl;
97 SetHighEnergyLimit(highEnergyLimit);
98 }
99
100 // Read data tables for all materials
101
102 crossSectionHandler = new G4CrossSectionHandler();
103 crossSectionHandler->Clear();
104 G4String crossSectionFile = "phot/pe-cs-";
105 crossSectionHandler->LoadData(crossSectionFile);
106
107 meanFreePathTable = 0;
108 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
109
110 shellCrossSectionHandler = new G4CrossSectionHandler();
111 shellCrossSectionHandler->Clear();
112 G4String shellCrossSectionFile = "phot/pe-ss-cs-";
113 shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
114
115 // SI - Buggy default ?
116 //generatorName = "geant4.6.2";
117 //ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator"); // default generator
118
119 //
120
121 if (verboseLevel > 2)
122 G4cout << "Loaded cross section files for Livermore PhotoElectric model" << G4endl;
123
124 G4cout << "Livermore PhotoElectric model is initialized " << G4endl
125 << "Energy range: "
126 << LowEnergyLimit() / keV << " keV - "
127 << HighEnergyLimit() / GeV << " GeV"
128 << G4endl;
129
130 if(isInitialised) return;
131
132 if(pParticleChange)
133 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
134 else
135 fParticleChange = new G4ParticleChangeForGamma();
136
137 isInitialised = true;
138}
139
140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141
142G4double G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(
143 const G4ParticleDefinition*,
144 G4double GammaEnergy,
145 G4double Z, G4double,
146 G4double, G4double)
147{
148 if (verboseLevel > 3)
149 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel" << G4endl;
150
151 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
152 return cs;
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157void G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
158 const G4MaterialCutsCouple* couple,
159 const G4DynamicParticle* aDynamicGamma,
160 G4double,
161 G4double)
162{
163
164 // Fluorescence generated according to:
165 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
166 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
167 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
168
169 if (verboseLevel > 3)
170 G4cout << "Calling SampleSecondaries() of G4LivermorePhotoElectricModel" << G4endl;
171
172 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
173
174 if (photonEnergy <= lowEnergyLimit)
175 {
176 fParticleChange->ProposeTrackStatus(fStopAndKill);
177 fParticleChange->SetProposedKineticEnergy(0.);
178 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
[991]179 // SI - IS THE FOLLOWING RETURN NECESSARY ?
[968]180 return ;
181 }
182
183 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection(); // Returns the normalized direction of the momentum
184
185 // Select randomly one element in the current material
186 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
187
188 // Select the ionised shell in the current atom according to shell cross sections
189 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
190
191 // Retrieve the corresponding identifier and binding energy of the selected shell
192 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
193 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
194 G4double bindingEnergy = shell->BindingEnergy();
195 G4int shellId = shell->ShellId();
196
197 // Create lists of pointers to DynamicParticles (photons and electrons)
198 // (Is the electron vector necessary? To be checked)
199 std::vector<G4DynamicParticle*>* photonVector = 0;
200 std::vector<G4DynamicParticle*> electronVector;
201
202 G4double energyDeposit = 0.0;
203
204 // Primary outcoming electron
205 G4double eKineticEnergy = photonEnergy - bindingEnergy;
206
207 // There may be cases where the binding energy of the selected shell is > photon energy
208 // In such cases do not generate secondaries
209 if (eKineticEnergy > 0.)
210 {
211 // SI - Removed safety
212
213 // Generate the electron only if with large enough range w.r.t. cuts and safety
214 //G4double safety = aStep.GetPostStepPoint()->GetSafety();
215
216 //if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
217 {
218
219 // Calculate direction of the photoelectron
220 G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
221 G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
222
223 // The electron is created ...
224 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
225 electronDirection,
226 eKineticEnergy);
227 electronVector.push_back(electron);
228 }
229 /*else
230 {
231 energyDeposit += eKineticEnergy;
232 }*/
233 }
234 else
235 {
236 bindingEnergy = photonEnergy;
237 }
238
239 G4int nElectrons = electronVector.size();
240 size_t nTotPhotons = 0;
241 G4int nPhotons=0;
242 const G4ProductionCutsTable* theCoupleTable=
243 G4ProductionCutsTable::GetProductionCutsTable();
244
245 size_t index = couple->GetIndex();
246 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
247 cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
248
249 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
250 cute = std::min(cutForLowEnergySecondaryPhotons,cute);
251
252 G4DynamicParticle* aPhoton;
253
254 // Generation of fluorescence
255 // Data in EADL are available only for Z > 5
256 // Protection to avoid generating photons in the unphysical case of
257 // shell binding energy > photon energy
258 if (Z > 5 && (bindingEnergy > cutg || bindingEnergy > cute))
259 {
260 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
261 nTotPhotons = photonVector->size();
262 for (size_t k=0; k<nTotPhotons; k++)
263 {
264 aPhoton = (*photonVector)[k];
265 if (aPhoton)
266 {
267 G4double itsCut = cutg;
268 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
269 G4double itsEnergy = aPhoton->GetKineticEnergy();
270
271 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
272 {
273 nPhotons++;
274 // Local energy deposit is given as the sum of the
275 // energies of incident photons minus the energies
276 // of the outcoming fluorescence photons
277 bindingEnergy -= itsEnergy;
278
279 }
280 else
281 {
282 delete aPhoton;
283 (*photonVector)[k] = 0;
284 }
285 }
286 }
287 }
288
289 energyDeposit += bindingEnergy;
290
291 // Final state
292
293 for (G4int l = 0; l<nElectrons; l++ )
294 {
295 aPhoton = electronVector[l];
296 if(aPhoton) {
297 fvect->push_back(aPhoton);
298 }
299 }
300 for ( size_t ll = 0; ll < nTotPhotons; ll++)
301 {
302 aPhoton = (*photonVector)[ll];
303 if(aPhoton) {
304 fvect->push_back(aPhoton);
305 }
306 }
307
308 delete photonVector;
309
310 if (energyDeposit < 0)
311 {
312 G4cout << "WARNING - "
313 << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
314 << G4endl;
315 energyDeposit = 0;
316 }
317
318 // kill incident photon
319 fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
320 fParticleChange->SetProposedKineticEnergy(0.);
321 fParticleChange->ProposeTrackStatus(fStopAndKill);
322 fParticleChange->ProposeLocalEnergyDeposit(energyDeposit);
323}
324
325//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
326
327void G4LivermorePhotoElectricModel::SetCutForLowEnSecPhotons(G4double cut)
328{
329 cutForLowEnergySecondaryPhotons = cut;
330 deexcitationManager.SetCutForSecondaryPhotons(cut);
331}
332
333//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
334
335void G4LivermorePhotoElectricModel::SetCutForLowEnSecElectrons(G4double cut)
336{
337 cutForLowEnergySecondaryElectrons = cut;
338 deexcitationManager.SetCutForAugerElectrons(cut);
339}
340
341//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342
343void G4LivermorePhotoElectricModel::ActivateAuger(G4bool val)
344{
345 deexcitationManager.ActivateAugerElectronProduction(val);
346}
347
348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349
350void G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* distribution)
351{
352 ElectronAngularGenerator = distribution;
353 ElectronAngularGenerator->PrintGeneratorInformation();
354}
355
356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357
358void G4LivermorePhotoElectricModel::SetAngularGenerator(const G4String& name)
359{
360 if (name == "default")
361 {
362 delete ElectronAngularGenerator;
363 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
364 generatorName = name;
365 }
366 else if (name == "standard")
367 {
368 delete ElectronAngularGenerator;
369 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
370 generatorName = name;
371 }
372 else if (name == "polarized")
373 {
374 delete ElectronAngularGenerator;
375 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
376 generatorName = name;
377 }
378 else
379 {
380 G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator - generator does not exist");
381 }
382
383 ElectronAngularGenerator->PrintGeneratorInformation();
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
387
388G4double G4LivermorePhotoElectricModel::GetMeanFreePath(const G4Track& track,
389 G4double, // previousStepSize
390 G4ForceCondition*)
391{
392 const G4DynamicParticle* photon = track.GetDynamicParticle();
393 G4double energy = photon->GetKineticEnergy();
394 G4Material* material = track.GetMaterial();
395 // size_t materialIndex = material->GetIndex();
396
397 G4double meanFreePath = DBL_MAX;
398
399 // if (energy > highEnergyLimit)
400 // meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
401 // else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
402 // else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
403
404 G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
405 if(cross > 0.0) meanFreePath = 1.0/cross;
406
407 return meanFreePath;
408}
409
Note: See TracBrowser for help on using the repository browser.