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

Last change on this file since 989 was 968, checked in by garnier, 15 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.