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