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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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