source: trunk/source/processes/electromagnetic/lowenergy/src/G4Penelope08PhotoElectricModel.cc @ 1317

Last change on this file since 1317 was 1316, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 22.6 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: G4Penelope08PhotoElectricModel.cc,v 1.3 2010/03/31 11:12:08 pandola Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 08 Jan 2010   L Pandola  First implementation
34
35#include "G4Penelope08PhotoElectricModel.hh"
36#include "G4ParticleDefinition.hh"
37#include "G4MaterialCutsCouple.hh"
38#include "G4ProductionCutsTable.hh"
39#include "G4DynamicParticle.hh"
40#include "G4PhysicsTable.hh"
41#include "G4PhysicsFreeVector.hh"
42#include "G4ElementTable.hh"
43#include "G4Element.hh"
44#include "G4CrossSectionHandler.hh"
45#include "G4AtomicTransitionManager.hh"
46#include "G4AtomicShell.hh"
47#include "G4Gamma.hh"
48#include "G4Electron.hh"
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52
53G4Penelope08PhotoElectricModel::G4Penelope08PhotoElectricModel(const G4ParticleDefinition*,
54                                             const G4String& nam)
55  :G4VEmModel(nam),isInitialised(false),logAtomicShellXS(0)
56{
57  fIntrinsicLowEnergyLimit = 100.0*eV;
58  fIntrinsicHighEnergyLimit = 100.0*GeV;
59  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
60  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
61  //
62  verboseLevel= 0;
63  // Verbosity scale:
64  // 0 = nothing
65  // 1 = warning for energy non-conservation
66  // 2 = details of energy budget
67  // 3 = calculation of cross sections, file openings, sampling of atoms
68  // 4 = entering in methods
69
70  //by default the model will inkove the atomic deexcitation
71  SetDeexcitationFlag(true); 
72  ActivateAuger(false);
73}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77G4Penelope08PhotoElectricModel::~G4Penelope08PhotoElectricModel()
78{ 
79  std::map <const G4int,G4PhysicsTable*>::iterator i;
80  if (logAtomicShellXS)
81    {
82      for (i=logAtomicShellXS->begin();i != logAtomicShellXS->end();i++)
83        {
84          G4PhysicsTable* tab = i->second;
85          tab->clearAndDestroy();
86          delete tab;
87        }
88    }
89  delete logAtomicShellXS;
90}
91
92//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93
94void G4Penelope08PhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
95                                       const G4DataVector& cuts)
96{
97  if (verboseLevel > 3)
98    G4cout << "Calling  G4Penelope08PhotoElectricModel::Initialise()" << G4endl;
99
100  // logAtomicShellXS is created only once, since it is  never cleared
101  if (!logAtomicShellXS)
102    logAtomicShellXS = new std::map<const G4int,G4PhysicsTable*>;
103
104  InitialiseElementSelectors(particle,cuts);
105
106  if (verboseLevel > 0) { 
107    G4cout << "Penelope Photo-Electric model is initialized " << G4endl
108           << "Energy range: "
109           << LowEnergyLimit() / MeV << " MeV - "
110           << HighEnergyLimit() / GeV << " GeV"
111           << G4endl;
112  }
113
114  if(isInitialised) return;
115  fParticleChange = GetParticleChangeForGamma();
116  isInitialised = true;
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121G4double G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom(
122                                       const G4ParticleDefinition*,
123                                             G4double energy,
124                                             G4double Z, G4double,
125                                             G4double, G4double)
126{
127  //
128  // Penelope model.
129  //
130
131  if (verboseLevel > 3)
132    G4cout << "Calling ComputeCrossSectionPerAtom() of G4Penelope08PhotoElectricModel" << G4endl;
133
134  G4int iZ = (G4int) Z;
135
136  //read data files
137  if (!logAtomicShellXS->count(iZ))
138    ReadDataFile(iZ);
139  //now it should be ok
140  if (!logAtomicShellXS->count(iZ))
141     {
142       G4cout << "Problem in G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom"
143              << G4endl;
144       G4Exception();
145     }
146
147  G4double cross = 0;
148
149  G4PhysicsTable* theTable =  logAtomicShellXS->find(iZ)->second;
150  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
151
152   if (!totalXSLog)
153     {
154       G4cout << "Problem in G4Penelope08PhotoElectricModel::ComputeCrossSectionPerAtom"
155         << G4endl;
156       G4Exception();
157     }
158   G4double logene = log(energy);
159   G4double logXS = totalXSLog->Value(logene);
160   cross = exp(logXS);
161 
162  if (verboseLevel > 2)
163    G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
164      " = " << cross/barn << " barn" << G4endl;
165  return cross;
166}
167
168//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169
170void G4Penelope08PhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
171                                                       const G4MaterialCutsCouple* couple,
172                                                       const G4DynamicParticle* aDynamicGamma,
173                                                       G4double,
174                                                       G4double)
175{
176  //
177  // Photoelectric effect, Penelope model
178  //
179  // The target atom and the target shell are sampled according to the Livermore
180  // database
181  //  D.E. Cullen et al., Report UCRL-50400 (1989)
182  // The angular distribution of the electron in the final state is sampled
183  // according to the Sauter distribution from
184  //  F. Sauter, Ann. Phys. 11 (1931) 454
185  // The energy of the final electron is given by the initial photon energy minus
186  // the binding energy. Fluorescence de-excitation is subsequently produced
187  // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
188  //  J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
189
190  if (verboseLevel > 3)
191    G4cout << "Calling SamplingSecondaries() of G4Penelope08PhotoElectricModel" << G4endl;
192
193  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
194
195  // always kill primary
196  fParticleChange->ProposeTrackStatus(fStopAndKill);
197  fParticleChange->SetProposedKineticEnergy(0.);
198
199  if (photonEnergy <= fIntrinsicLowEnergyLimit)
200    {
201      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
202      return ;
203    }
204
205  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
206
207  // Select randomly one element in the current material
208  if (verboseLevel > 2)
209    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
210
211  // atom can be selected efficiently if element selectors are initialised
212  const G4Element* anElement =
213    SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy);
214  G4int Z = (G4int) anElement->GetZ();
215  if (verboseLevel > 2)
216    G4cout << "Selected " << anElement->GetName() << G4endl;
217 
218  // Select the ionised shell in the current atom according to shell cross sections
219  //shellIndex = 0 --> K shell
220  //             1-3 --> L shells
221  //             4-8 --> M shells
222  //             9 --> outer shells cumulatively
223  //
224  size_t shellIndex = SelectRandomShell(Z,photonEnergy);
225
226  if (verboseLevel > 2)
227    G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl;
228
229  // Retrieve the corresponding identifier and binding energy of the selected shell
230  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
231
232  //The number of shell cross section possibly reported in the Penelope database
233  //might be different from the number of shells in the G4AtomicTransitionManager
234  //(namely, Penelope may contain more shell, especially for very light elements).
235  //In order to avoid a warning message from the G4AtomicTransitionManager, I
236  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
237  //has a shellID>maxID, it sets the shellID to the last valid shell.
238  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
239  if (shellIndex >= numberOfShells)
240    shellIndex = numberOfShells-1;
241
242  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
243  G4double bindingEnergy = shell->BindingEnergy();
244  G4int shellId = shell->ShellId();
245
246  //Penelope considers only K, L and M shells. Cross sections of outer shells are
247  //not included in the Penelope database. If SelectRandomShell() returns
248  //shellIndex = 9, it means that an outer shell was ionized. In this case the
249  //Penelope recipe is to set bindingEnergy = 0 (the energy is entirely assigned
250  //to the electron) and to disregard fluorescence.
251  if (shellIndex == 9)
252    bindingEnergy = 0.*eV;
253
254
255  G4double localEnergyDeposit = 0.0;
256  G4double cosTheta = 1.0;
257
258  // Primary outcoming electron
259  G4double eKineticEnergy = photonEnergy - bindingEnergy;
260 
261  // There may be cases where the binding energy of the selected shell is > photon energy
262  // In such cases do not generate secondaries
263  if (eKineticEnergy > 0.)
264    {
265      // The electron is created
266      // Direction sampled from the Sauter distribution
267      cosTheta = SampleElectronDirection(eKineticEnergy);
268      G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
269      G4double phi = twopi * G4UniformRand() ;
270      G4double dirx = sinTheta * std::cos(phi);
271      G4double diry = sinTheta * std::sin(phi);
272      G4double dirz = cosTheta ;
273      G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
274      electronDirection.rotateUz(photonDirection);
275      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
276                                                           electronDirection, 
277                                                           eKineticEnergy);
278      fvect->push_back(electron);
279    } 
280  else
281    {
282      bindingEnergy = photonEnergy;
283    }
284
285  G4double energyInFluorescence = 0; //testing purposes
286
287  //Now, take care of fluorescence, if required
288  if(DeexcitationFlag() && Z > 5) 
289    {
290      const G4ProductionCutsTable* theCoupleTable=
291        G4ProductionCutsTable::GetProductionCutsTable();
292      size_t indx = couple->GetIndex();
293      G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
294      G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
295     
296      // Protection to avoid generating photons in the unphysical case of
297      // shell binding energy > photon energy
298      if (bindingEnergy > cutG || bindingEnergy > cutE)
299        {
300          deexcitationManager.SetCutForSecondaryPhotons(cutG);
301          deexcitationManager.SetCutForAugerElectrons(cutE);
302          std::vector<G4DynamicParticle*>* photonVector = 
303            deexcitationManager.GenerateParticles(Z,shellId); 
304          //Check for secondaries
305          if(photonVector) 
306            {
307              for (size_t k=0; k< photonVector->size(); k++)
308                {
309                  G4DynamicParticle* aPhoton = (*photonVector)[k];
310                  if (aPhoton)
311                    {
312                      G4double itsEnergy = aPhoton->GetKineticEnergy();
313                      if (itsEnergy <= bindingEnergy)
314                        {
315                          if(aPhoton->GetDefinition() == G4Gamma::Gamma())
316                            energyInFluorescence += itsEnergy;
317                          bindingEnergy -= itsEnergy;
318                          fvect->push_back(aPhoton);
319                        }
320                      else
321                        {
322                          delete aPhoton;
323                          (*photonVector)[k] = 0;
324                        }
325                    }
326                }
327              delete photonVector;
328            }
329        }
330    }
331  //Residual energy is deposited locally
332  localEnergyDeposit += bindingEnergy;
333     
334  if (localEnergyDeposit < 0)
335    {
336      G4cout << "WARNING - "
337             << "G4Penelope08PhotoElectric::PostStepDoIt - Negative energy deposit"
338             << G4endl;
339      localEnergyDeposit = 0;
340    }
341
342  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
343
344  if (verboseLevel > 1)
345    {
346      G4cout << "-----------------------------------------------------------" << G4endl;
347      G4cout << "Energy balance from G4Penelope08PhotoElectric" << G4endl;
348      G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " << 
349        anElement->GetName() << G4endl;
350      G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
351      G4cout << "-----------------------------------------------------------" << G4endl;
352      if (eKineticEnergy)
353        G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
354      G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
355      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
356      G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << 
357        " keV" << G4endl;
358      G4cout << "-----------------------------------------------------------" << G4endl;
359    }
360  if (verboseLevel > 0)
361    {
362      G4double energyDiff = 
363        std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
364      if (energyDiff > 0.05*keV)
365        G4cout << "Warning from G4Penelope08PhotoElectric: problem with energy conservation: " << 
366          (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV
367               << " keV (final) vs. " << 
368          photonEnergy/keV << " keV (initial)" << G4endl;
369    }
370}
371
372//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373
374void G4Penelope08PhotoElectricModel::ActivateAuger(G4bool augerbool)
375{
376  if (!DeexcitationFlag() && augerbool)
377    {
378      G4cout << "WARNING - G4Penelope08PhotoElectricModel" << G4endl;
379      G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
380      G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
381    }
382  deexcitationManager.ActivateAugerElectronProduction(augerbool);
383  if (verboseLevel > 1)
384    G4cout << "Auger production set to " << augerbool << G4endl;
385}
386
387//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388
389G4double G4Penelope08PhotoElectricModel::SampleElectronDirection(G4double energy)
390{
391  G4double costheta = 1.0;
392  if (energy>1*GeV) return costheta;
393 
394  //1) initialize energy-dependent variables
395  // Variable naming according to Eq. (2.24) of Penelope Manual
396  // (pag. 44)
397  G4double gamma = 1.0 + energy/electron_mass_c2;
398  G4double gamma2 = gamma*gamma;
399  G4double beta = std::sqrt((gamma2-1.0)/gamma2);
400   
401  // ac corresponds to "A" of Eq. (2.31)
402  //
403  G4double ac = (1.0/beta) - 1.0;
404  G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
405  G4double a2 = ac + 2.0;
406  G4double gtmax = 2.0*(a1 + 1.0/ac);
407 
408  G4double tsam = 0;
409  G4double gtr = 0;
410
411  //2) sampling. Eq. (2.31) of Penelope Manual
412  // tsam = 1-std::cos(theta)
413  // gtr = rejection function according to Eq. (2.28)
414  do{
415    G4double rand = G4UniformRand();
416    tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
417    gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
418  }while(G4UniformRand()*gtmax > gtr);
419  costheta = 1.0-tsam;
420 
421
422  return costheta;
423}
424
425//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426
427void G4Penelope08PhotoElectricModel::ReadDataFile(G4int Z)
428{
429  if (verboseLevel > 2)
430    {
431      G4cout << "G4Penelope08PhotoElectricModel::ReadDataFile()" << G4endl;
432      G4cout << "Going to read PhotoElectric data files for Z=" << Z << G4endl;
433    }
434 
435  char* path = getenv("G4LEDATA");
436  if (!path)
437    {
438      G4String excep = "G4Penelope08PhotoElectricModel - G4LEDATA environment variable not set!";
439      G4Exception(excep);
440    }
441 
442  /*
443    Read the cross section file
444  */
445  std::ostringstream ost;
446  if (Z>9)
447    ost << path << "/penelope/photoelectric/pdgph" << Z << ".p08";
448  else
449    ost << path << "/penelope/photoelectric/pdgph0" << Z << ".p08";
450  std::ifstream file(ost.str().c_str());
451  if (!file.is_open())
452    {
453      G4String excep = "G4Penelope08PhotoElectricModel - data file " + G4String(ost.str()) + " not found!";
454      G4Exception(excep);
455    }
456  //I have to know in advance how many points are in the data list
457  //to initialize the G4PhysicsFreeVector()
458  size_t ndata=0;
459  G4String line;
460  while( getline(file, line) )
461    ndata++;
462  ndata -= 1;
463  //G4cout << "Found: " << ndata << " lines" << G4endl;
464
465  file.clear();
466  file.close();
467  file.open(ost.str().c_str());
468
469  G4int readZ =0;
470  size_t nShells= 0;
471  file >> readZ >> nShells;
472
473  if (verboseLevel > 3)
474    G4cout << "Element Z=" << Z << " , nShells = " << nShells << G4endl;
475
476  //check the right file is opened.
477  if (readZ != Z || nShells <= 0)
478    {
479      G4cout << "G4Penelope08PhotoElectricModel::ReadDataFile()" << G4endl;
480      G4cout << "Corrupted data file for Z=" << Z << G4endl;
481      G4Exception();
482    }
483  G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
484
485  //the table has to contain nShell+1 G4PhysicsFreeVectors,
486  //(theTable)[0] --> total cross section
487  //(theTable)[ishell] --> cross section for shell (ishell-1)
488
489
490
491  //reserve space for the vectors
492  //everything is log-log
493  for (size_t i=0;i<nShells+1;i++)
494    thePhysicsTable->push_back(new G4PhysicsFreeVector(ndata));
495
496  size_t k =0;
497  for (k=0;k<ndata && !file.eof();k++)
498    {
499      G4double energy = 0;
500      G4double aValue = 0;
501      file >> energy ;
502      energy *= eV;
503      G4double logene = log(energy);
504      //loop on the columns
505      for (size_t i=0;i<nShells+1;i++)
506        {
507          file >> aValue;
508          aValue *= barn;
509          G4PhysicsFreeVector* theVec = (G4PhysicsFreeVector*) ((*thePhysicsTable)[i]); 
510          if (aValue < 1e-40*cm2) //protection against log(0)
511            aValue = 1e-40*cm2;
512          theVec->PutValue(k,logene,log(aValue));
513        }
514    }
515
516  if (verboseLevel > 2)
517    {
518      G4cout << "G4Penelope08PhotoElectricModel: read " << k << " points for element Z = " 
519             << Z << G4endl;
520    }
521
522  logAtomicShellXS->insert(std::make_pair(Z,thePhysicsTable));
523 
524  file.close();
525  return;
526}
527
528//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
529
530size_t G4Penelope08PhotoElectricModel::SelectRandomShell(G4int Z,G4double energy)
531{
532  G4double logEnergy = log(energy);
533
534  //Check if data have been read (it should be!)
535  if (!logAtomicShellXS->count(Z))
536     {
537       G4cout << "Problem in G4Penelope08PhotoElectricModel::SelectRandomShell" << G4endl;
538       G4cout << "Cannot find data for Z=" << Z << G4endl;
539       G4Exception();
540     }
541
542  size_t shellIndex = 0;
543 
544  G4PhysicsTable* theTable =  logAtomicShellXS->find(Z)->second;
545
546  G4DataVector* tempVector = new G4DataVector();
547
548  G4double sum = 0;
549  //loop on shell partial XS, retrieve the value for the given energy and store on
550  //a temporary vector
551  tempVector->push_back(sum); //first element is zero
552
553  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
554  G4double logXS = totalXSLog->Value(logEnergy);
555  G4double totalXS = exp(logXS);
556                                           
557  //Notice: totalXS is the total cross section and it does *not* correspond to
558  //the sum of partialXS's, since these include only K, L and M shells.
559  //
560  // Therefore, here one have to consider the possibility of ionisation of
561  // an outer shell. Conventionally, it is indicated with id=10 in Penelope
562  //
563 
564  for (size_t k=1;k<theTable->entries();k++)
565    {
566      G4PhysicsFreeVector* partialXSLog = (G4PhysicsFreeVector*) (*theTable)[k];
567      G4double logXS = partialXSLog->Value(logEnergy);
568      G4double partialXS = exp(logXS);
569      sum += partialXS;
570      tempVector->push_back(sum);     
571    }
572
573  tempVector->push_back(totalXS); //last element
574
575  G4double random = G4UniformRand()*totalXS; 
576
577  /*
578  for (size_t i=0;i<tempVector->size(); i++)
579    G4cout << i << " " << (*tempVector)[i]/totalXS << G4endl;
580  */
581 
582  //locate bin of tempVector
583  //Now one has to sample according to the elements in tempVector
584  //This gives the left edge of the interval...
585  size_t lowerBound = 0;
586  size_t upperBound = tempVector->size()-1; 
587  while (lowerBound <= upperBound)
588   {
589     size_t midBin = (lowerBound + upperBound)/2;
590     if( random < (*tempVector)[midBin])
591       upperBound = midBin-1; 
592     else
593       lowerBound = midBin+1; 
594   }
595 
596  shellIndex = upperBound;
597
598  delete tempVector;
599  return shellIndex;
600}
601
602//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
603
604size_t G4Penelope08PhotoElectricModel::GetNumberOfShellXS(G4int Z)
605{
606  //read data files
607  if (!logAtomicShellXS->count(Z))
608    ReadDataFile(Z);
609  //now it should be ok
610  if (!logAtomicShellXS->count(Z))
611     {
612       G4cout << "Problem in G4Penelope08PhotoElectricModel::GetNumberOfShellXS()"
613              << G4endl;
614       G4Exception();
615     }
616  //one vector is allocated for the _total_ cross section
617  size_t nEntries = logAtomicShellXS->find(Z)->second->entries();
618  return  (nEntries-1);
619}
620
621//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
622
623G4double G4Penelope08PhotoElectricModel::GetShellCrossSection(G4int Z,size_t shellID,G4double energy)
624{
625  //this forces also the loading of the data
626  size_t entries = GetNumberOfShellXS(Z);
627
628  if (shellID >= entries)
629    {
630      G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl;
631      G4cout << "so shellID should be from 0 to " << entries-1 << G4endl;
632      return 0;
633    }
634 
635  G4PhysicsTable* theTable =  logAtomicShellXS->find(Z)->second;
636  //[0] is the total XS, shellID is in the element [shellID+1]
637  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1];
638 
639  if (!totalXSLog)
640     {
641       G4cout << "Problem in G4Penelope08PhotoElectricModel::GetShellCrossSection()"
642         << G4endl;
643       G4Exception();
644     }
645   G4double logene = log(energy);
646   G4double logXS = totalXSLog->Value(logene);
647   G4double cross = exp(logXS);
648   if (cross < 2e-40*cm2) cross = 0;
649   return cross;
650}
651
652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
653
654G4String G4Penelope08PhotoElectricModel::WriteTargetShell(size_t shellID)
655{
656  G4String theShell = "outer shell";
657  if (shellID == 0)
658    theShell = "K";
659  else if (shellID == 1)
660    theShell = "L1";
661  else if (shellID == 2)
662    theShell = "L2";
663  else if (shellID == 3)
664    theShell = "L3";
665  else if (shellID == 4)
666    theShell = "M1";
667  else if (shellID == 5)
668    theShell = "M2";
669  else if (shellID == 6)
670    theShell = "M3";
671  else if (shellID == 7)
672    theShell = "M4";
673  else if (shellID == 8)
674    theShell = "M5";
675     
676  return theShell;
677}
Note: See TracBrowser for help on using the repository browser.