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

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

update
File size: 14.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// $Id: G4LivermorePhotoElectricModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-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)
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  delete meanFreePathTable;
67  delete crossSectionHandler;
68  delete shellCrossSectionHandler;
69  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  InitialiseElementSelectors(particle,cuts);
81
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);
179      // SI - IS THE FOLLOWING RETURN NECESSARY ?
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.