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

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

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