ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopePhotoElectricModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PenelopePhotoElectricModel.cc
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 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 08 Jan 2010 L Pandola First implementation
32 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger
33 // is active.
34 // Make sure that fluorescence/Auger is generated
35 // only if above threshold
36 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope)
37 // 10 Jun 2011 L Pandola Migrate atomic deexcitation interface
38 // 07 Oct 2011 L Pandola Bug fix (potential violation of energy conservation)
39 // 27 Sep 2013 L Pandola Migrate to MT paradigm, only master model manages
40 // tables.
41 // 02 Oct 2013 L Pandola Rewrite sampling algorithm of SelectRandomShell()
42 // to improve CPU performances
43 //
44 
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ParticleDefinition.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4DynamicParticle.hh"
51 #include "G4PhysicsTable.hh"
52 #include "G4PhysicsFreeVector.hh"
53 #include "G4ElementTable.hh"
54 #include "G4Element.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4AutoLock.hh"
60 #include "G4LossTableManager.hh"
61 #include "G4Exp.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
66  const G4String& nam)
67  :G4VEmModel(nam),fParticleChange(0),fParticle(0),
68  isInitialised(false),fAtomDeexcitation(0),logAtomicShellXS(0),fLocalTable(false)
69 {
72  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
74  //
75 
76  if (part)
77  SetParticle(part);
78 
79  verboseLevel= 0;
80  // Verbosity scale:
81  // 0 = nothing
82  // 1 = warning for energy non-conservation
83  // 2 = details of energy budget
84  // 3 = calculation of cross sections, file openings, sampling of atoms
85  // 4 = entering in methods
86 
87  //Mark this model as "applicable" for atomic deexcitation
88  SetDeexcitationFlag(true);
89 
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
96 {
97  if (IsMaster() || fLocalTable)
98  {
99  if (logAtomicShellXS)
100  {
101  for (auto& item : (*logAtomicShellXS))
102  {
103  //G4PhysicsTable* tab = item.second;
104  //tab->clearAndDestroy();
105  delete item.second;
106  }
107  }
108  delete logAtomicShellXS;
109  }
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
113 
115  const G4DataVector& cuts)
116 {
117  if (verboseLevel > 3)
118  G4cout << "Calling G4PenelopePhotoElectricModel::Initialise()" << G4endl;
119 
121  //Issue warning if the AtomicDeexcitation has not been declared
122  if (!fAtomDeexcitation)
123  {
124  G4cout << G4endl;
125  G4cout << "WARNING from G4PenelopePhotoElectricModel " << G4endl;
126  G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
127  G4cout << "any fluorescence/Auger emission." << G4endl;
128  G4cout << "Please make sure this is intended" << G4endl;
129  }
130 
131  SetParticle(particle);
132 
133  //Only the master model creates/fills/destroys the tables
134  if (IsMaster() && particle == fParticle)
135  {
136 
137  // logAtomicShellXS is created only once, since it is never cleared
138  if (!logAtomicShellXS)
139  logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>;
140 
141  G4ProductionCutsTable* theCoupleTable =
143 
144  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
145  {
146  const G4Material* material =
147  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
148  const G4ElementVector* theElementVector = material->GetElementVector();
149 
150  for (size_t j=0;j<material->GetNumberOfElements();j++)
151  {
152  G4int iZ = (G4int) theElementVector->at(j)->GetZ();
153  //read data files only in the master
154  if (!logAtomicShellXS->count(iZ))
155  ReadDataFile(iZ);
156  }
157  }
158 
159 
160  InitialiseElementSelectors(particle,cuts);
161 
162  if (verboseLevel > 0) {
163  G4cout << "Penelope Photo-Electric model v2008 is initialized " << G4endl
164  << "Energy range: "
165  << LowEnergyLimit() / MeV << " MeV - "
166  << HighEnergyLimit() / GeV << " GeV";
167  }
168  }
169 
170  if(isInitialised) return;
172  isInitialised = true;
173 
174 }
175 
177  G4VEmModel *masterModel)
178 {
179  if (verboseLevel > 3)
180  G4cout << "Calling G4PenelopePhotoElectricModel::InitialiseLocal()" << G4endl;
181 
182  //
183  //Check that particle matches: one might have multiple master models (e.g.
184  //for e+ and e-).
185  //
186  if (part == fParticle)
187  {
189 
190  //Get the const table pointers from the master to the workers
191  const G4PenelopePhotoElectricModel* theModel =
192  static_cast<G4PenelopePhotoElectricModel*> (masterModel);
193 
194  logAtomicShellXS = theModel->logAtomicShellXS;
195 
196  //Same verbosity for all workers, as the master
197  verboseLevel = theModel->verboseLevel;
198  }
199 
200  return;
201 }
202 
203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
204 namespace { G4Mutex PenelopePhotoElectricModelMutex = G4MUTEX_INITIALIZER; }
206  const G4ParticleDefinition*,
210 {
211  //
212  // Penelope model v2008
213  //
214 
215  if (verboseLevel > 3)
216  G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
217 
218  G4int iZ = (G4int) Z;
219 
220  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
221  //not invoked
222  if (!logAtomicShellXS)
223  {
224  //create a **thread-local** version of the table. Used only for G4EmCalculator and
225  //Unit Tests
226  fLocalTable = true;
227  logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>;
228  }
229 
230  //now it should be ok
231  if (!logAtomicShellXS->count(iZ))
232  {
233  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
234  //not filled up. This can happen in a UnitTest or via G4EmCalculator
235  if (verboseLevel > 0)
236  {
237  //Issue a G4Exception (warning) only in verbose mode
239  ed << "Unable to retrieve the shell cross section table for Z=" << iZ << G4endl;
240  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
241  G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
242  "em2038",JustWarning,ed);
243  }
244  //protect file reading via autolock
245  G4AutoLock lock(&PenelopePhotoElectricModelMutex);
246  ReadDataFile(iZ);
247  lock.unlock();
248 
249  }
250 
251  G4double cross = 0;
252 
253  G4PhysicsTable* theTable = logAtomicShellXS->find(iZ)->second;
254  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
255 
256  if (!totalXSLog)
257  {
258  G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
259  "em2039",FatalException,
260  "Unable to retrieve the total cross section table");
261  return 0;
262  }
263  G4double logene = std::log(energy);
264  G4double logXS = totalXSLog->Value(logene);
265  cross = G4Exp(logXS);
266 
267  if (verboseLevel > 2)
268  G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
269  " = " << cross/barn << " barn" << G4endl;
270  return cross;
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 
275 void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
276  const G4MaterialCutsCouple* couple,
277  const G4DynamicParticle* aDynamicGamma,
278  G4double,
279  G4double)
280 {
281  //
282  // Photoelectric effect, Penelope model v2008
283  //
284  // The target atom and the target shell are sampled according to the Livermore
285  // database
286  // D.E. Cullen et al., Report UCRL-50400 (1989)
287  // The angular distribution of the electron in the final state is sampled
288  // according to the Sauter distribution from
289  // F. Sauter, Ann. Phys. 11 (1931) 454
290  // The energy of the final electron is given by the initial photon energy minus
291  // the binding energy. Fluorescence de-excitation is subsequently produced
292  // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
293  // J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
294 
295  if (verboseLevel > 3)
296  G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
297 
298  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
299 
300  // always kill primary
303 
304  if (photonEnergy <= fIntrinsicLowEnergyLimit)
305  {
307  return ;
308  }
309 
310  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
311 
312  // Select randomly one element in the current material
313  if (verboseLevel > 2)
314  G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
315 
316  // atom can be selected efficiently if element selectors are initialised
317  const G4Element* anElement =
318  SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy);
319  G4int Z = (G4int) anElement->GetZ();
320  if (verboseLevel > 2)
321  G4cout << "Selected " << anElement->GetName() << G4endl;
322 
323  // Select the ionised shell in the current atom according to shell cross sections
324  //shellIndex = 0 --> K shell
325  // 1-3 --> L shells
326  // 4-8 --> M shells
327  // 9 --> outer shells cumulatively
328  //
329  size_t shellIndex = SelectRandomShell(Z,photonEnergy);
330 
331  if (verboseLevel > 2)
332  G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl;
333 
334  // Retrieve the corresponding identifier and binding energy of the selected shell
336 
337  //The number of shell cross section possibly reported in the Penelope database
338  //might be different from the number of shells in the G4AtomicTransitionManager
339  //(namely, Penelope may contain more shell, especially for very light elements).
340  //In order to avoid a warning message from the G4AtomicTransitionManager, I
341  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
342  //has a shellID>maxID, it sets the shellID to the last valid shell.
343  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
344  if (shellIndex >= numberOfShells)
345  shellIndex = numberOfShells-1;
346 
347  const G4AtomicShell* shell = fTransitionManager->Shell(Z,shellIndex);
349  //G4int shellId = shell->ShellId();
350 
351  //Penelope considers only K, L and M shells. Cross sections of outer shells are
352  //not included in the Penelope database. If SelectRandomShell() returns
353  //shellIndex = 9, it means that an outer shell was ionized. In this case the
354  //Penelope recipe is to set bindingEnergy = 0 (the energy is entirely assigned
355  //to the electron) and to disregard fluorescence.
356  if (shellIndex == 9)
357  bindingEnergy = 0.*eV;
358 
359 
360  G4double localEnergyDeposit = 0.0;
361  G4double cosTheta = 1.0;
362 
363  // Primary outcoming electron
364  G4double eKineticEnergy = photonEnergy - bindingEnergy;
365 
366  // There may be cases where the binding energy of the selected shell is > photon energy
367  // In such cases do not generate secondaries
368  if (eKineticEnergy > 0.)
369  {
370  // The electron is created
371  // Direction sampled from the Sauter distribution
372  cosTheta = SampleElectronDirection(eKineticEnergy);
373  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
375  G4double dirx = sinTheta * std::cos(phi);
376  G4double diry = sinTheta * std::sin(phi);
377  G4double dirz = cosTheta ;
378  G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
379  electronDirection.rotateUz(photonDirection);
381  electronDirection,
382  eKineticEnergy);
383  fvect->push_back(electron);
384  }
385  else
386  bindingEnergy = photonEnergy;
387 
388 
389  G4double energyInFluorescence = 0; //testing purposes
390  G4double energyInAuger = 0; //testing purposes
391 
392  //Now, take care of fluorescence, if required. According to the Penelope
393  //recipe, I have to skip fluoresence completely if shellIndex == 9
394  //(= sampling of a shell outer than K,L,M)
395  if (fAtomDeexcitation && shellIndex<9)
396  {
397  G4int index = couple->GetIndex();
399  {
400  size_t nBefore = fvect->size();
401  fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
402  size_t nAfter = fvect->size();
403 
404  if (nAfter > nBefore) //actual production of fluorescence
405  {
406  for (size_t j=nBefore;j<nAfter;j++) //loop on products
407  {
408  G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
409  if (itsEnergy < bindingEnergy) // valid secondary, generate it
410  {
411  bindingEnergy -= itsEnergy;
412  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
413  energyInFluorescence += itsEnergy;
414  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
415  energyInAuger += itsEnergy;
416  }
417  else //invalid secondary: takes more than the available energy: delete it
418  {
419  delete (*fvect)[j];
420  (*fvect)[j] = nullptr;
421  }
422  }
423  }
424  }
425  }
426 
427  //Residual energy is deposited locally
428  localEnergyDeposit += bindingEnergy;
429 
430  if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
431  {
432  G4Exception("G4PenelopePhotoElectricModel::SampleSecondaries()",
433  "em2099",JustWarning,"WARNING: Negative local energy deposit");
434  localEnergyDeposit = 0;
435  }
436 
437  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
438 
439  if (verboseLevel > 1)
440  {
441  G4cout << "-----------------------------------------------------------" << G4endl;
442  G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
443  G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " <<
444  anElement->GetName() << G4endl;
445  G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
446  G4cout << "-----------------------------------------------------------" << G4endl;
447  if (eKineticEnergy)
448  G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
449  if (energyInFluorescence)
450  G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
451  if (energyInAuger)
452  G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
453  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
454  G4cout << "Total final state: " <<
455  (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV <<
456  " keV" << G4endl;
457  G4cout << "-----------------------------------------------------------" << G4endl;
458  }
459  if (verboseLevel > 0)
460  {
461  G4double energyDiff =
462  std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy);
463  if (energyDiff > 0.05*keV)
464  {
465  G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
466  (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV
467  << " keV (final) vs. " <<
468  photonEnergy/keV << " keV (initial)" << G4endl;
469  G4cout << "-----------------------------------------------------------" << G4endl;
470  G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
471  G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " <<
472  anElement->GetName() << G4endl;
473  G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
474  G4cout << "-----------------------------------------------------------" << G4endl;
475  if (eKineticEnergy)
476  G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
477  if (energyInFluorescence)
478  G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
479  if (energyInAuger)
480  G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
481  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
482  G4cout << "Total final state: " <<
483  (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV <<
484  " keV" << G4endl;
485  G4cout << "-----------------------------------------------------------" << G4endl;
486  }
487  }
488 }
489 
490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
491 
493 {
494  G4double costheta = 1.0;
495  if (energy>1*GeV) return costheta;
496 
497  //1) initialize energy-dependent variables
498  // Variable naming according to Eq. (2.24) of Penelope Manual
499  // (pag. 44)
500  G4double gamma = 1.0 + energy/electron_mass_c2;
501  G4double gamma2 = gamma*gamma;
502  G4double beta = std::sqrt((gamma2-1.0)/gamma2);
503 
504  // ac corresponds to "A" of Eq. (2.31)
505  //
506  G4double ac = (1.0/beta) - 1.0;
507  G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
508  G4double a2 = ac + 2.0;
509  G4double gtmax = 2.0*(a1 + 1.0/ac);
510 
511  G4double tsam = 0;
512  G4double gtr = 0;
513 
514  //2) sampling. Eq. (2.31) of Penelope Manual
515  // tsam = 1-std::cos(theta)
516  // gtr = rejection function according to Eq. (2.28)
517  do{
518  G4double rand = G4UniformRand();
519  tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
520  gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
521  }while(G4UniformRand()*gtmax > gtr);
522  costheta = 1.0-tsam;
523 
524 
525  return costheta;
526 }
527 
528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
529 
531 {
532  if (!IsMaster())
533  //Should not be here!
534  G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
535  "em0100",FatalException,"Worker thread in this method");
536 
537  if (verboseLevel > 2)
538  {
539  G4cout << "G4PenelopePhotoElectricModel::ReadDataFile()" << G4endl;
540  G4cout << "Going to read PhotoElectric data files for Z=" << Z << G4endl;
541  }
542 
543  char* path = std::getenv("G4LEDATA");
544  if (!path)
545  {
546  G4String excep = "G4PenelopePhotoElectricModel - G4LEDATA environment variable not set!";
547  G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
548  "em0006",FatalException,excep);
549  return;
550  }
551 
552  /*
553  Read the cross section file
554  */
555  std::ostringstream ost;
556  if (Z>9)
557  ost << path << "/penelope/photoelectric/pdgph" << Z << ".p08";
558  else
559  ost << path << "/penelope/photoelectric/pdgph0" << Z << ".p08";
560  std::ifstream file(ost.str().c_str());
561  if (!file.is_open())
562  {
563  G4String excep = "G4PenelopePhotoElectricModel - data file " + G4String(ost.str()) + " not found!";
564  G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
565  "em0003",FatalException,excep);
566  }
567  //I have to know in advance how many points are in the data list
568  //to initialize the G4PhysicsFreeVector()
569  size_t ndata=0;
570  G4String line;
571  while( getline(file, line) )
572  ndata++;
573  ndata -= 1;
574  //G4cout << "Found: " << ndata << " lines" << G4endl;
575 
576  file.clear();
577  file.close();
578  file.open(ost.str().c_str());
579 
580  G4int readZ =0;
581  size_t nShells= 0;
582  file >> readZ >> nShells;
583 
584  if (verboseLevel > 3)
585  G4cout << "Element Z=" << Z << " , nShells = " << nShells << G4endl;
586 
587  //check the right file is opened.
588  if (readZ != Z || nShells <= 0 || nShells > 50) //protect nShell against large values
589  {
591  ed << "Corrupted data file for Z=" << Z << G4endl;
592  G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
593  "em0005",FatalException,ed);
594  return;
595  }
596  G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
597 
598  //the table has to contain nShell+1 G4PhysicsFreeVectors,
599  //(theTable)[0] --> total cross section
600  //(theTable)[ishell] --> cross section for shell (ishell-1)
601 
602  //reserve space for the vectors
603  //everything is log-log
604  for (size_t i=0;i<nShells+1;i++)
605  thePhysicsTable->push_back(new G4PhysicsFreeVector(ndata));
606 
607  size_t k =0;
608  for (k=0;k<ndata && !file.eof();k++)
609  {
610  G4double energy = 0;
611  G4double aValue = 0;
612  file >> energy ;
613  energy *= eV;
614  G4double logene = std::log(energy);
615  //loop on the columns
616  for (size_t i=0;i<nShells+1;i++)
617  {
618  file >> aValue;
619  aValue *= barn;
620  G4PhysicsFreeVector* theVec = (G4PhysicsFreeVector*) ((*thePhysicsTable)[i]);
621  if (aValue < 1e-40*cm2) //protection against log(0)
622  aValue = 1e-40*cm2;
623  theVec->PutValue(k,logene,std::log(aValue));
624  }
625  }
626 
627  if (verboseLevel > 2)
628  {
629  G4cout << "G4PenelopePhotoElectricModel: read " << k << " points for element Z = "
630  << Z << G4endl;
631  }
632 
633  logAtomicShellXS->insert(std::make_pair(Z,thePhysicsTable));
634 
635  file.close();
636  return;
637 }
638 
639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
640 
642 {
643  if (!IsMaster())
644  //Should not be here!
645  G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
646  "em0100",FatalException,"Worker thread in this method");
647 
648  //read data files
649  if (!logAtomicShellXS->count(Z))
650  ReadDataFile(Z);
651  //now it should be ok
652  if (!logAtomicShellXS->count(Z))
653  {
655  ed << "Cannot find shell cross section data for Z=" << Z << G4endl;
656  G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
657  "em2038",FatalException,ed);
658  }
659  //one vector is allocated for the _total_ cross section
660  size_t nEntries = logAtomicShellXS->find(Z)->second->entries();
661  return (nEntries-1);
662 }
663 
664 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
665 
667 {
668  //this forces also the loading of the data
669  size_t entries = GetNumberOfShellXS(Z);
670 
671  if (shellID >= entries)
672  {
673  G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl;
674  G4cout << "so shellID should be from 0 to " << entries-1 << G4endl;
675  return 0;
676  }
677 
678  G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
679  //[0] is the total XS, shellID is in the element [shellID+1]
680  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1];
681 
682  if (!totalXSLog)
683  {
684  G4Exception("G4PenelopePhotoElectricModel::GetShellCrossSection()",
685  "em2039",FatalException,
686  "Unable to retrieve the total cross section table");
687  return 0;
688  }
689  G4double logene = std::log(energy);
690  G4double logXS = totalXSLog->Value(logene);
691  G4double cross = G4Exp(logXS);
692  if (cross < 2e-40*cm2) cross = 0;
693  return cross;
694 }
695 
696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
697 
699 {
700  G4String theShell = "outer shell";
701  if (shellID == 0)
702  theShell = "K";
703  else if (shellID == 1)
704  theShell = "L1";
705  else if (shellID == 2)
706  theShell = "L2";
707  else if (shellID == 3)
708  theShell = "L3";
709  else if (shellID == 4)
710  theShell = "M1";
711  else if (shellID == 5)
712  theShell = "M2";
713  else if (shellID == 6)
714  theShell = "M3";
715  else if (shellID == 7)
716  theShell = "M4";
717  else if (shellID == 8)
718  theShell = "M5";
719 
720  return theShell;
721 }
722 
723 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
724 
726 {
727  if(!fParticle) {
728  fParticle = p;
729  }
730 }
731 
732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
733 
735 {
736  G4double logEnergy = std::log(energy);
737 
738  //Check if data have been read (it should be!)
739  if (!logAtomicShellXS->count(Z))
740  {
742  ed << "Cannot find shell cross section data for Z=" << Z << G4endl;
743  G4Exception("G4PenelopePhotoElectricModel::SelectRandomShell()",
744  "em2038",FatalException,ed);
745  }
746 
747  // size_t shellIndex = 0;
748 
749  G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
750 
751  G4double sum = 0;
752 
753  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
754  G4double logXS = totalXSLog->Value(logEnergy);
755  G4double totalXS = G4Exp(logXS);
756 
757  //Notice: totalXS is the total cross section and it does *not* correspond to
758  //the sum of partialXS's, since these include only K, L and M shells.
759  //
760  // Therefore, here one have to consider the possibility of ionisation of
761  // an outer shell. Conventionally, it is indicated with id=10 in Penelope
762  //
763  G4double random = G4UniformRand()*totalXS;
764 
765  for (size_t k=1;k<theTable->entries();k++)
766  {
767  //Add one shell
768  G4PhysicsFreeVector* partialXSLog = (G4PhysicsFreeVector*) (*theTable)[k];
769  G4double logXSLocal = partialXSLog->Value(logEnergy);
770  G4double partialXS = G4Exp(logXSLocal);
771  sum += partialXS;
772  if (random <= sum)
773  return k-1;
774  }
775  //none of the shells K, L, M: return outer shell
776  return 9;
777 }