Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectricModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectricModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopePhotoElectricModel.cc (Version 9.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4PenelopePhotoElectricModel.cc,v 1.13 2010-11-26 11:51:11 pandola Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
 26 //                                                 28 //
 27 // Author: Luciano Pandola                         29 // Author: Luciano Pandola
 28 //                                                 30 //
 29 // History:                                        31 // History:
 30 // --------                                        32 // --------
 31 // 08 Jan 2010   L Pandola  First implementati <<  33 // 08 Oct 2008   L Pandola  Migration from process to model 
 32 // 01 Feb 2011   L Pandola  Suppress fake ener <<  34 // 08 Jan 2009   L Pandola  Check shell index to avoid mismatch between 
 33 //                          is active.         <<  35 //                          the Penelope cross section database and the 
 34 //                          Make sure that flu <<  36 //                          G4AtomicTransitionManager database. It suppresses 
 35 //                          only if above thre <<  37 //                          a warning from G4AtomicTransitionManager only. 
 36 // 25 May 2011   L Pandola  Renamed (make v200 <<  38 //                          Results are unchanged.
 37 // 10 Jun 2011   L Pandola  Migrate atomic dee <<  39 // 25 Mar 2009   L Pandola  Small fix to avoid wrong energy-violation warnings
 38 // 07 Oct 2011   L Pandola  Bug fix (potential <<  40 // 17 Apr 2009   V Ivanchenko Cleanup initialisation and generation of secondaries:
 39 // 27 Sep 2013   L Pandola  Migrate to MT para <<  41 //                  - apply internal high-energy limit only in constructor 
 40 //                          tables.            <<  42 //                  - do not apply low-energy limit (default is 0)
 41 // 02 Oct 2013   L Pandola  Rewrite sampling a <<  43 //                  - do not apply production threshold on secondaries
 42 //                          to improve CPU per <<  44 // 19 May 2009   L Pandola    Explicitely set to zero pointers deleted in 
                                                   >>  45 //                            Initialise(), since they might be checked later on
                                                   >>  46 // 21 Oct 2009   L Pandola    Remove un-necessary fUseAtomicDeexcitation flag - now managed by
                                                   >>  47 //                            G4VEmModel::DeexcitationFlag()
                                                   >>  48 // 15 Mar 2010   L Pandola    Explicitely initialize Auger to false
 43 //                                                 49 //
 44                                                    50 
 45 #include "G4PenelopePhotoElectricModel.hh"         51 #include "G4PenelopePhotoElectricModel.hh"
 46 #include "G4PhysicalConstants.hh"              << 
 47 #include "G4SystemOfUnits.hh"                  << 
 48 #include "G4ParticleDefinition.hh"                 52 #include "G4ParticleDefinition.hh"
 49 #include "G4MaterialCutsCouple.hh"                 53 #include "G4MaterialCutsCouple.hh"
                                                   >>  54 #include "G4ProductionCutsTable.hh"
 50 #include "G4DynamicParticle.hh"                    55 #include "G4DynamicParticle.hh"
 51 #include "G4PhysicsTable.hh"                       56 #include "G4PhysicsTable.hh"
 52 #include "G4PhysicsFreeVector.hh"              << 
 53 #include "G4ElementTable.hh"                       57 #include "G4ElementTable.hh"
 54 #include "G4Element.hh"                            58 #include "G4Element.hh"
                                                   >>  59 #include "G4CrossSectionHandler.hh"
 55 #include "G4AtomicTransitionManager.hh"            60 #include "G4AtomicTransitionManager.hh"
 56 #include "G4AtomicShell.hh"                        61 #include "G4AtomicShell.hh"
 57 #include "G4Gamma.hh"                              62 #include "G4Gamma.hh"
 58 #include "G4Electron.hh"                           63 #include "G4Electron.hh"
 59 #include "G4AutoLock.hh"                       <<  64 #include "G4VEMDataSet.hh"
 60 #include "G4LossTableManager.hh"               << 
 61 #include "G4Exp.hh"                            << 
 62                                                    65 
 63 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 64                                                    67 
 65 const G4int G4PenelopePhotoElectricModel::fMax << 
 66 G4PhysicsTable* G4PenelopePhotoElectricModel:: << 
 67                                                    68 
 68 //....oooOO0OOooo........oooOO0OOooo........oo <<  69 G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition*,
 69                                                <<  70                                              const G4String& nam)
 70 G4PenelopePhotoElectricModel::G4PenelopePhotoE <<  71   :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
 71                  const G4String& nam)          <<  72    shellCrossSectionHandler(0)
 72   :G4VEmModel(nam),fParticleChange(nullptr),fP << 
 73    fAtomDeexcitation(nullptr),fIsInitialised(f << 
 74 {                                                  73 {
 75   fIntrinsicLowEnergyLimit = 100.0*eV;             74   fIntrinsicLowEnergyLimit = 100.0*eV;
 76   fIntrinsicHighEnergyLimit = 100.0*GeV;           75   fIntrinsicHighEnergyLimit = 100.0*GeV;
 77   //  SetLowEnergyLimit(fIntrinsicLowEnergyLim     76   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
 78   SetHighEnergyLimit(fIntrinsicHighEnergyLimit     77   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
 79   //                                               78   //
 80                                                <<  79   verboseLevel= 0;
 81   if (part)                                    << 
 82     SetParticle(part);                         << 
 83                                                << 
 84   fVerboseLevel= 0;                            << 
 85   // Verbosity scale:                              80   // Verbosity scale:
 86   // 0 = nothing                               <<  81   // 0 = nothing 
 87   // 1 = warning for energy non-conservation   <<  82   // 1 = warning for energy non-conservation 
 88   // 2 = details of energy budget                  83   // 2 = details of energy budget
 89   // 3 = calculation of cross sections, file o     84   // 3 = calculation of cross sections, file openings, sampling of atoms
 90   // 4 = entering in methods                       85   // 4 = entering in methods
 91                                                    86 
 92   //Mark this model as "applicable" for atomic <<  87   //by default the model will inkove the atomic deexcitation
 93   SetDeexcitationFlag(true);                   <<  88   SetDeexcitationFlag(true);  
 94                                                <<  89   ActivateAuger(false);
 95   fTransitionManager = G4AtomicTransitionManag << 
 96 }                                                  90 }
 97                                                    91 
 98 //....oooOO0OOooo........oooOO0OOooo........oo     92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 99                                                    93 
100 G4PenelopePhotoElectricModel::~G4PenelopePhoto     94 G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel()
101 {                                              <<  95 {  
102   if (IsMaster() || fLocalTable)               <<  96   if (crossSectionHandler) delete crossSectionHandler;
103     {                                          <<  97   if (shellCrossSectionHandler) delete shellCrossSectionHandler;
104       for(G4int i=0; i<=fMaxZ; ++i)            << 
105   {                                            << 
106     if(fLogAtomicShellXS[i]) {                 << 
107       fLogAtomicShellXS[i]->clearAndDestroy(); << 
108       delete fLogAtomicShellXS[i];             << 
109       fLogAtomicShellXS[i] = nullptr;          << 
110     }                                          << 
111   }                                            << 
112     }                                          << 
113 }                                                  98 }
114                                                    99 
115 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116                                                   101 
117 void G4PenelopePhotoElectricModel::Initialise( << 102 void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition*,
118                 const G4DataVector& cuts)      << 103                                        const G4DataVector& )
119 {                                                 104 {
120   if (fVerboseLevel > 3)                       << 105   if (verboseLevel > 3)
121     G4cout << "Calling  G4PenelopePhotoElectri    106     G4cout << "Calling  G4PenelopePhotoElectricModel::Initialise()" << G4endl;
122                                                << 107   if (crossSectionHandler)
123   fAtomDeexcitation = G4LossTableManager::Inst << 
124   //Issue warning if the AtomicDeexcitation ha << 
125   if (!fAtomDeexcitation)                      << 
126     {                                             108     {
127       G4cout << G4endl;                        << 109       crossSectionHandler->Clear();
128       G4cout << "WARNING from G4PenelopePhotoE << 110       delete crossSectionHandler;
129       G4cout << "Atomic de-excitation module i << 111       crossSectionHandler = 0;
130       G4cout << "any fluorescence/Auger emissi << 112     }
131       G4cout << "Please make sure this is inte << 113   if (shellCrossSectionHandler)
132     }                                          << 114     {
133                                                << 115       shellCrossSectionHandler->Clear();
134   SetParticle(particle);                       << 116       delete shellCrossSectionHandler;
135                                                << 117       shellCrossSectionHandler =0;
136   //Only the master model creates/fills/destro << 118     }
137   if (IsMaster() && particle == fParticle)     << 119 
138     {                                          << 120   //Re-initialize cross section handlers
139       G4ProductionCutsTable* theCoupleTable =  << 121   crossSectionHandler = new G4CrossSectionHandler();
140   G4ProductionCutsTable::GetProductionCutsTabl << 122   crossSectionHandler->Clear();
141                                                << 123   G4String crossSectionFile = "penelope/ph-cs-pen-";
142       for (G4int i=0;i<(G4int)theCoupleTable-> << 124   crossSectionHandler->LoadData(crossSectionFile);
143   {                                            << 125   shellCrossSectionHandler = new G4CrossSectionHandler();
144     const G4Material* material =               << 126   shellCrossSectionHandler->Clear();
145       theCoupleTable->GetMaterialCutsCouple(i) << 127   crossSectionFile = "penelope/ph-ss-cs-pen-";
146     const G4ElementVector* theElementVector =  << 128   shellCrossSectionHandler->LoadShellData(crossSectionFile);
147                                                << 129   //This is used to retrieve cross section values later on
148     for (std::size_t j=0;j<material->GetNumber << 130   G4VEMDataSet* emdata = 
149       {                                        << 131     crossSectionHandler->BuildMeanFreePathForMaterials();
150         G4int iZ = theElementVector->at(j)->Ge << 132   //The method BuildMeanFreePathForMaterials() is required here only to force 
151         //read data files only in the master   << 133   //the building of an internal table: the output pointer can be deleted
152         if (!fLogAtomicShellXS[iZ])            << 134   delete emdata;
153     ReadDataFile(iZ);                          << 135 
154       }                                        << 136   if (verboseLevel > 2) 
155   }                                            << 137     G4cout << "Loaded cross section files for PenelopePhotoElectric" << G4endl;
156                                                << 138 
157       InitialiseElementSelectors(particle,cuts << 139   if (verboseLevel > 0) { 
158                                                << 140     G4cout << "Penelope Photo-Electric model is initialized " << G4endl
159       if (fVerboseLevel > 0) {                 << 141      << "Energy range: "
160   G4cout << "Penelope Photo-Electric model v20 << 142      << LowEnergyLimit() / MeV << " MeV - "
161          << "Energy range: "                   << 143      << HighEnergyLimit() / GeV << " GeV"
162          << LowEnergyLimit() / MeV << " MeV -  << 144      << G4endl;
163          << HighEnergyLimit() / GeV << " GeV"; << 145   }
164       }                                        << 
165     }                                          << 
166                                                   146 
167   if(fIsInitialised) return;                   << 147   if(isInitialised) return;
168   fParticleChange = GetParticleChangeForGamma(    148   fParticleChange = GetParticleChangeForGamma();
169   fIsInitialised = true;                       << 149   isInitialised = true;
170                                                << 
171 }                                              << 
172                                                << 
173 void G4PenelopePhotoElectricModel::InitialiseL << 
174                  G4VEmModel *masterModel)      << 
175 {                                              << 
176   if (fVerboseLevel > 3)                       << 
177     G4cout << "Calling  G4PenelopePhotoElectri << 
178   //                                           << 
179   //Check that particle matches: one might hav << 
180   //for e+ and e-).                            << 
181   //                                           << 
182   if (part == fParticle)                       << 
183     {                                          << 
184       SetElementSelectors(masterModel->GetElem << 
185                                                << 
186       //Get the const table pointers from the  << 
187       const G4PenelopePhotoElectricModel* theM << 
188   static_cast<G4PenelopePhotoElectricModel*> ( << 
189       for(G4int i=0; i<=fMaxZ; ++i)            << 
190   fLogAtomicShellXS[i] = theModel->fLogAtomicS << 
191       //Same verbosity for all workers, as the << 
192       fVerboseLevel = theModel->fVerboseLevel; << 
193     }                                          << 
194                                                << 
195  return;                                       << 
196 }                                                 150 }
197                                                   151 
198 //....oooOO0OOooo........oooOO0OOooo........oo    152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 namespace { G4Mutex  PenelopePhotoElectricMode << 153 
200 G4double G4PenelopePhotoElectricModel::Compute    154 G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom(
201                   const G4ParticleDefinition*, << 155                                        const G4ParticleDefinition*,
202                   G4double energy,             << 156                                              G4double energy,
203                   G4double Z, G4double,        << 157                                              G4double Z, G4double,
204                   G4double, G4double)          << 158                                              G4double, G4double)
205 {                                                 159 {
206   //                                              160   //
207   // Penelope model v2008                      << 161   // Penelope model. Use data-driven approach for cross section estimate (and 
                                                   >> 162   // also shell sampling from a given atom). Data are from the Livermore database
                                                   >> 163   //  D.E. Cullen et al., Report UCRL-50400 (1989)
208   //                                              164   //
209   if (fVerboseLevel > 3)                       << 
210     G4cout << "Calling ComputeCrossSectionPerA << 
211                                                << 
212   G4int iZ = G4int(Z);                         << 
213                                                << 
214   if (!fLogAtomicShellXS[iZ])                  << 
215     {                                          << 
216       //If we are here, it means that Initiali << 
217       //not filled up. This can happen in a Un << 
218       if (fVerboseLevel > 0)                   << 
219   {                                            << 
220     //Issue a G4Exception (warning) only in ve << 
221     G4ExceptionDescription ed;                 << 
222     ed << "Unable to retrieve the shell cross  << 
223     ed << "This can happen only in Unit Tests  << 
224     G4Exception("G4PenelopePhotoElectricModel: << 
225           "em2038",JustWarning,ed);            << 
226   }                                            << 
227       //protect file reading via autolock      << 
228       G4AutoLock lock(&PenelopePhotoElectricMo << 
229       ReadDataFile(iZ);                        << 
230       lock.unlock();                           << 
231     }                                          << 
232                                                   165 
233   G4double cross = 0;                          << 166   if (verboseLevel > 3)
234   G4PhysicsTable* theTable =  fLogAtomicShellX << 167     G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
235   G4PhysicsFreeVector* totalXSLog = (G4Physics << 
236                                                << 
237    if (!totalXSLog)                            << 
238      {                                         << 
239        G4Exception("G4PenelopePhotoElectricMod << 
240        "em2039",FatalException,                << 
241        "Unable to retrieve the total cross sec << 
242        return 0;                               << 
243      }                                         << 
244    G4double logene = G4Log(energy);            << 
245    G4double logXS = totalXSLog->Value(logene); << 
246    cross = G4Exp(logXS);                       << 
247                                                   168 
248   if (fVerboseLevel > 2)                       << 169   G4int iZ = (G4int) Z;
                                                   >> 170   //  if (!crossSectionHandler) // VI: should not be 
                                                   >> 171   //  {
                                                   >> 172   //    G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
                                                   >> 173   //    G4cout << "The cross section handler is not correctly initialized" << G4endl;
                                                   >> 174   //    G4Exception();
                                                   >> 175   //  }
                                                   >> 176   G4double cs = crossSectionHandler->FindValue(iZ,energy);
                                                   >> 177  
                                                   >> 178   if (verboseLevel > 2)
249     G4cout << "Photoelectric cross section at     179     G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
250       " = " << cross/barn << " barn" << G4endl << 180       " = " << cs/barn << " barn" << G4endl;
251   return cross;                                << 181   return cs;
252 }                                                 182 }
253                                                   183 
254 //....oooOO0OOooo........oooOO0OOooo........oo    184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255                                                   185 
256 void G4PenelopePhotoElectricModel::SampleSecon    186 void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
257                  const G4MaterialCutsCouple* c << 187                 const G4MaterialCutsCouple* couple,
258                  const G4DynamicParticle* aDyn << 188                 const G4DynamicParticle* aDynamicGamma,
259                  G4double,                     << 189                 G4double,
260                  G4double)                     << 190                 G4double)
261 {                                                 191 {
262   //                                              192   //
263   // Photoelectric effect, Penelope model v200 << 193   // Photoelectric effect, Penelope model
264   //                                              194   //
265   // The target atom and the target shell are  << 195   // The target atom and the target shell are sampled according to the Livermore 
266   // database                                  << 196   // database 
267   //  D.E. Cullen et al., Report UCRL-50400 (1    197   //  D.E. Cullen et al., Report UCRL-50400 (1989)
268   // The angular distribution of the electron  << 198   // The angular distribution of the electron in the final state is sampled 
269   // according to the Sauter distribution from << 199   // according to the Sauter distribution from 
270   //  F. Sauter, Ann. Phys. 11 (1931) 454         200   //  F. Sauter, Ann. Phys. 11 (1931) 454
271   // The energy of the final electron is given << 201   // The energy of the final electron is given by the initial photon energy minus 
272   // the binding energy. Fluorescence de-excit << 202   // the binding energy. Fluorescence de-excitation is subsequently produced 
273   // (to fill the vacancy) according to the ge    203   // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
274   //  J. Stepanek, Comp. Phys. Comm. 1206 pp 1    204   //  J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
275                                                   205 
276   if (fVerboseLevel > 3)                       << 206   if (verboseLevel > 3)
277     G4cout << "Calling SamplingSecondaries() o    207     G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
278                                                   208 
279   G4double photonEnergy = aDynamicGamma->GetKi    209   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
280                                                   210 
281   // always kill primary                          211   // always kill primary
282   fParticleChange->ProposeTrackStatus(fStopAnd    212   fParticleChange->ProposeTrackStatus(fStopAndKill);
283   fParticleChange->SetProposedKineticEnergy(0.    213   fParticleChange->SetProposedKineticEnergy(0.);
284                                                   214 
285   if (photonEnergy <= fIntrinsicLowEnergyLimit    215   if (photonEnergy <= fIntrinsicLowEnergyLimit)
286     {                                             216     {
287       fParticleChange->ProposeLocalEnergyDepos    217       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
288       return ;                                    218       return ;
289     }                                             219     }
290                                                   220 
291   G4ParticleMomentum photonDirection = aDynami    221   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
292                                                   222 
293   // Select randomly one element in the curren    223   // Select randomly one element in the current material
294   if (fVerboseLevel > 2)                       << 224   if (verboseLevel > 2)
295     G4cout << "Going to select element in " <<    225     G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
296                                                << 226   //use crossSectionHandler instead of G4EmElementSelector because in this case 
297   // atom can be selected efficiently if eleme << 227   //the dimension of the table is equal to the dimension of the database 
298   const G4Element* anElement =                 << 228   //(less interpolation errors)
299     SelectRandomAtom(couple,G4Gamma::GammaDefi << 229   G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
300   G4int Z = anElement->GetZasInt();            << 230   if (verboseLevel > 2)
301   if (fVerboseLevel > 2)                       << 231     G4cout << "Selected Z = " << Z << G4endl;
302     G4cout << "Selected " << anElement->GetNam << 
303                                                   232 
304   // Select the ionised shell in the current a    233   // Select the ionised shell in the current atom according to shell cross sections
305   //shellIndex = 0 --> K shell                 << 234   size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
306   //             1-3 --> L shells              << 
307   //             4-8 --> M shells              << 
308   //             9 --> outer shells cumulative << 
309   //                                           << 
310   std::size_t shellIndex = SelectRandomShell(Z << 
311                                                << 
312   if (fVerboseLevel > 2)                       << 
313     G4cout << "Selected shell " << shellIndex  << 
314                                                   235 
315   // Retrieve the corresponding identifier and    236   // Retrieve the corresponding identifier and binding energy of the selected shell
316   const G4AtomicTransitionManager* transitionM    237   const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
317                                                   238 
318   //The number of shell cross section possibly << 239   //The number of shell cross section possibly reported in the Penelope database 
319   //might be different from the number of shel    240   //might be different from the number of shells in the G4AtomicTransitionManager
320   //(namely, Penelope may contain more shell,     241   //(namely, Penelope may contain more shell, especially for very light elements).
321   //In order to avoid a warning message from t << 242   //In order to avoid a warning message from the G4AtomicTransitionManager, I 
322   //add this protection. Results are anyway ch    243   //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
323   //has a shellID>maxID, it sets the shellID t << 244   //has a shellID>maxID, it sets the shellID to the last valid shell. 
324   std::size_t numberOfShells = (std::size_t) t << 245   size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
325   if (shellIndex >= numberOfShells)               246   if (shellIndex >= numberOfShells)
326     shellIndex = numberOfShells-1;                247     shellIndex = numberOfShells-1;
327                                                   248 
328   const G4AtomicShell* shell = fTransitionMana << 249   const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
329   G4double bindingEnergy = shell->BindingEnerg    250   G4double bindingEnergy = shell->BindingEnergy();
330                                                << 251   G4int shellId = shell->ShellId();
331   //Penelope considers only K, L and M shells. << 
332   //not included in the Penelope database. If  << 
333   //shellIndex = 9, it means that an outer she << 
334   //Penelope recipe is to set bindingEnergy =  << 
335   //to the electron) and to disregard fluoresc << 
336   if (shellIndex == 9)                         << 
337     bindingEnergy = 0.*eV;                     << 
338                                                   252 
339   G4double localEnergyDeposit = 0.0;              253   G4double localEnergyDeposit = 0.0;
340   G4double cosTheta = 1.0;                     << 
341                                                   254 
342   // Primary outcoming electron                   255   // Primary outcoming electron
343   G4double eKineticEnergy = photonEnergy - bin    256   G4double eKineticEnergy = photonEnergy - bindingEnergy;
344                                                << 257   
345   // There may be cases where the binding ener    258   // There may be cases where the binding energy of the selected shell is > photon energy
346   // In such cases do not generate secondaries    259   // In such cases do not generate secondaries
347   if (eKineticEnergy > 0.)                        260   if (eKineticEnergy > 0.)
348     {                                          << 261     {    
349       // The electron is created                  262       // The electron is created
350       // Direction sampled from the Sauter dis    263       // Direction sampled from the Sauter distribution
351       cosTheta = SampleElectronDirection(eKine << 264       G4double cosTheta = SampleElectronDirection(eKineticEnergy);
352       G4double sinTheta = std::sqrt(1-cosTheta    265       G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
353       G4double phi = twopi * G4UniformRand() ;    266       G4double phi = twopi * G4UniformRand() ;
354       G4double dirx = sinTheta * std::cos(phi)    267       G4double dirx = sinTheta * std::cos(phi);
355       G4double diry = sinTheta * std::sin(phi)    268       G4double diry = sinTheta * std::sin(phi);
356       G4double dirz = cosTheta ;                  269       G4double dirz = cosTheta ;
357       G4ThreeVector electronDirection(dirx,dir    270       G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
358       electronDirection.rotateUz(photonDirecti    271       electronDirection.rotateUz(photonDirection);
359       G4DynamicParticle* electron = new G4Dyna << 272       G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
360                  electronDirection,            << 273                  electronDirection, 
361                  eKineticEnergy);                 274                  eKineticEnergy);
362       fvect->push_back(electron);                 275       fvect->push_back(electron);
363     }                                          << 276     } 
364   else                                            277   else
365     bindingEnergy = photonEnergy;              << 278     {
366                                                << 279       bindingEnergy = photonEnergy;
                                                   >> 280     }
367   G4double energyInFluorescence = 0; //testing    281   G4double energyInFluorescence = 0; //testing purposes
368   G4double energyInAuger = 0; //testing purpos << 
369                                                   282 
370   //Now, take care of fluorescence, if require << 283   //Now, take care of fluorescence, if required
371   //recipe, I have to skip fluoresence complet << 284   if(DeexcitationFlag() && Z > 5) 
372   //(= sampling of a shell outer than K,L,M)   << 
373   if (fAtomDeexcitation && shellIndex<9)       << 
374     {                                             285     {
375       G4int index = couple->GetIndex();        << 286       const G4ProductionCutsTable* theCoupleTable=
376       if (fAtomDeexcitation->CheckDeexcitation << 287   G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 288       size_t indx = couple->GetIndex();
                                                   >> 289       G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
                                                   >> 290       G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
                                                   >> 291       
                                                   >> 292       // Protection to avoid generating photons in the unphysical case of 
                                                   >> 293       // shell binding energy > photon energy
                                                   >> 294       if (bindingEnergy > cutG || bindingEnergy > cutE)
377   {                                               295   {
378     std::size_t nBefore = fvect->size();       << 296     deexcitationManager.SetCutForSecondaryPhotons(cutG);
379     fAtomDeexcitation->GenerateParticles(fvect << 297     deexcitationManager.SetCutForAugerElectrons(cutE);
380     std::size_t nAfter = fvect->size();        << 298     std::vector<G4DynamicParticle*>* photonVector = 
381                                                << 299       deexcitationManager.GenerateParticles(Z,shellId); 
382     if (nAfter > nBefore) //actual production  << 300     //Check for secondaries
                                                   >> 301           if(photonVector) 
383       {                                           302       {
384         for (std::size_t j=nBefore;j<nAfter;++ << 303         for (size_t k=0; k< photonVector->size(); k++)
385     {                                             304     {
386       G4double itsEnergy = ((*fvect)[j])->GetK << 305       G4DynamicParticle* aPhoton = (*photonVector)[k];
387       if (itsEnergy < bindingEnergy) // valid  << 306       if (aPhoton)
388         {                                         307         {
389           bindingEnergy -= itsEnergy;          << 308           G4double itsEnergy = aPhoton->GetKineticEnergy();
390           if (((*fvect)[j])->GetParticleDefini << 309           if (itsEnergy <= bindingEnergy)
391       energyInFluorescence += itsEnergy;       << 310       {
392           else if (((*fvect)[j])->GetParticleD << 311         if(aPhoton->GetDefinition() == G4Gamma::Gamma())
393       energyInAuger += itsEnergy;              << 312           energyInFluorescence += itsEnergy;
                                                   >> 313         bindingEnergy -= itsEnergy;
                                                   >> 314         fvect->push_back(aPhoton);
                                                   >> 315       }
                                                   >> 316           else
                                                   >> 317       {
                                                   >> 318         delete aPhoton;
                                                   >> 319         (*photonVector)[k] = 0;
                                                   >> 320       }
394         }                                         321         }
395       else //invalid secondary: takes more tha << 
396         {                                      << 
397           delete (*fvect)[j];                  << 
398           (*fvect)[j] = nullptr;               << 
399         }                                      << 
400     }                                             322     }
                                                   >> 323         delete photonVector;
401       }                                           324       }
402   }                                               325   }
403     }                                             326     }
404                                                << 
405   //Residual energy is deposited locally          327   //Residual energy is deposited locally
406   localEnergyDeposit += bindingEnergy;            328   localEnergyDeposit += bindingEnergy;
407                                                << 329       
408   if (localEnergyDeposit < 0) //Should not be: << 330   if (localEnergyDeposit < 0)
409     {                                             331     {
410       G4Exception("G4PenelopePhotoElectricMode << 332       G4cout << "WARNING - "
411       "em2099",JustWarning,"WARNING: Negative  << 333        << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
                                                   >> 334        << G4endl;
412       localEnergyDeposit = 0;                     335       localEnergyDeposit = 0;
413     }                                             336     }
414                                                   337 
415   fParticleChange->ProposeLocalEnergyDeposit(l    338   fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
416                                                   339 
417   if (fVerboseLevel > 1)                       << 340   if (verboseLevel > 1)
418     {                                             341     {
419       G4cout << "-----------------------------    342       G4cout << "-----------------------------------------------------------" << G4endl;
420       G4cout << "Energy balance from G4Penelop    343       G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
421       G4cout << "Selected shell: " << WriteTar << 
422   anElement->GetName() << G4endl;              << 
423       G4cout << "Incoming photon energy: " <<     344       G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
424       G4cout << "-----------------------------    345       G4cout << "-----------------------------------------------------------" << G4endl;
425       if (eKineticEnergy)                         346       if (eKineticEnergy)
426   G4cout << "Outgoing electron " << eKineticEn    347   G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
427       if (energyInFluorescence)                << 348       G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
428   G4cout << "Fluorescence x-rays: " << energyI << 
429       if (energyInAuger)                       << 
430   G4cout << "Auger electrons: " << energyInAug << 
431       G4cout << "Local energy deposit " << loc    349       G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
432       G4cout << "Total final state: " <<       << 350       G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << 
433   (eKineticEnergy+energyInFluorescence+localEn << 
434   " keV" << G4endl;                               351   " keV" << G4endl;
435       G4cout << "-----------------------------    352       G4cout << "-----------------------------------------------------------" << G4endl;
436     }                                             353     }
437   if (fVerboseLevel > 0)                       << 354   if (verboseLevel > 0)
438     {                                             355     {
439       G4double energyDiff =                    << 356       G4double energyDiff = 
440   std::fabs(eKineticEnergy+energyInFluorescenc << 357   std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
441       if (energyDiff > 0.05*keV)                  358       if (energyDiff > 0.05*keV)
442   {                                            << 359   G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " << 
443     G4cout << "Warning from G4PenelopePhotoEle << 360     (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV 
444       (eKineticEnergy+energyInFluorescence+loc << 361          << " keV (final) vs. " << 
445      << " keV (final) vs. " <<                 << 362     photonEnergy/keV << " keV (initial)" << G4endl;
446       photonEnergy/keV << " keV (initial)" <<  << 363     }
447     G4cout << "------------------------------- << 364 }
448     G4cout << "Energy balance from G4PenelopeP << 365 
449     G4cout << "Selected shell: " << WriteTarge << 366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450       anElement->GetName() << G4endl;          << 367 
451     G4cout << "Incoming photon energy: " << ph << 368 void G4PenelopePhotoElectricModel::ActivateAuger(G4bool augerbool)
452     G4cout << "------------------------------- << 369 {
453     if (eKineticEnergy)                        << 370   if (!DeexcitationFlag() && augerbool)
454       G4cout << "Outgoing electron " << eKinet << 371     {
455     if (energyInFluorescence)                  << 372       G4cout << "WARNING - G4PenelopePhotoElectricModel" << G4endl;
456       G4cout << "Fluorescence x-rays: " << ene << 373       G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
457     if (energyInAuger)                         << 374       G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
458       G4cout << "Auger electrons: " << energyI << 
459     G4cout << "Local energy deposit " << local << 
460     G4cout << "Total final state: " <<         << 
461       (eKineticEnergy+energyInFluorescence+loc << 
462       " keV" << G4endl;                        << 
463     G4cout << "------------------------------- << 
464   }                                            << 
465     }                                             375     }
                                                   >> 376   deexcitationManager.ActivateAugerElectronProduction(augerbool);
                                                   >> 377   if (verboseLevel > 1)
                                                   >> 378     G4cout << "Auger production set to " << augerbool << G4endl;
466 }                                                 379 }
467                                                   380 
468 //....oooOO0OOooo........oooOO0OOooo........oo    381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469                                                   382 
470 G4double G4PenelopePhotoElectricModel::SampleE    383 G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy)
471 {                                                 384 {
472   G4double costheta = 1.0;                        385   G4double costheta = 1.0;
473   if (energy>1*GeV) return costheta;              386   if (energy>1*GeV) return costheta;
474                                                << 387  
475   //1) initialize energy-dependent variables      388   //1) initialize energy-dependent variables
476   // Variable naming according to Eq. (2.24) o    389   // Variable naming according to Eq. (2.24) of Penelope Manual
477   // (pag. 44)                                    390   // (pag. 44)
478   G4double gamma = 1.0 + energy/electron_mass_    391   G4double gamma = 1.0 + energy/electron_mass_c2;
479   G4double gamma2 = gamma*gamma;                  392   G4double gamma2 = gamma*gamma;
480   G4double beta = std::sqrt((gamma2-1.0)/gamma    393   G4double beta = std::sqrt((gamma2-1.0)/gamma2);
481                                                << 394    
482   // ac corresponds to "A" of Eq. (2.31)          395   // ac corresponds to "A" of Eq. (2.31)
483   //                                              396   //
484   G4double ac = (1.0/beta) - 1.0;                 397   G4double ac = (1.0/beta) - 1.0;
485   G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(ga    398   G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
486   G4double a2 = ac + 2.0;                         399   G4double a2 = ac + 2.0;
487   G4double gtmax = 2.0*(a1 + 1.0/ac);             400   G4double gtmax = 2.0*(a1 + 1.0/ac);
488                                                << 401  
489   G4double tsam = 0;                              402   G4double tsam = 0;
490   G4double gtr = 0;                               403   G4double gtr = 0;
491                                                   404 
492   //2) sampling. Eq. (2.31) of Penelope Manual    405   //2) sampling. Eq. (2.31) of Penelope Manual
493   // tsam = 1-std::cos(theta)                     406   // tsam = 1-std::cos(theta)
494   // gtr = rejection function according to Eq.    407   // gtr = rejection function according to Eq. (2.28)
495   do{                                             408   do{
496     G4double rand = G4UniformRand();              409     G4double rand = G4UniformRand();
497     tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(r    410     tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
498     gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));    411     gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
499   }while(G4UniformRand()*gtmax > gtr);            412   }while(G4UniformRand()*gtmax > gtr);
500   costheta = 1.0-tsam;                            413   costheta = 1.0-tsam;
501                                                << 
502   return costheta;                                414   return costheta;
503 }                                                 415 }
504                                                   416 
505 //....oooOO0OOooo........oooOO0OOooo........oo << 
506                                                << 
507 void G4PenelopePhotoElectricModel::ReadDataFil << 
508 {                                              << 
509   if (!IsMaster())                             << 
510       //Should not be here!                    << 
511     G4Exception("G4PenelopePhotoElectricModel: << 
512     "em0100",FatalException,"Worker thread in  << 
513                                                << 
514   if (fVerboseLevel > 2)                       << 
515     {                                          << 
516       G4cout << "G4PenelopePhotoElectricModel: << 
517       G4cout << "Going to read PhotoElectric d << 
518     }                                          << 
519                                                << 
520     const char* path = G4FindDataDir("G4LEDATA << 
521     if(!path)                                  << 
522     {                                          << 
523       G4String excep = "G4PenelopePhotoElectri << 
524       G4Exception("G4PenelopePhotoElectricMode << 
525       "em0006",FatalException,excep);          << 
526       return;                                  << 
527     }                                          << 
528                                                << 
529   /*                                           << 
530     Read the cross section file                << 
531   */                                           << 
532   std::ostringstream ost;                      << 
533   if (Z>9)                                     << 
534     ost << path << "/penelope/photoelectric/pd << 
535   else                                         << 
536     ost << path << "/penelope/photoelectric/pd << 
537   std::ifstream file(ost.str().c_str());       << 
538   if (!file.is_open())                         << 
539     {                                          << 
540       G4String excep = "G4PenelopePhotoElectri << 
541       G4Exception("G4PenelopePhotoElectricMode << 
542       "em0003",FatalException,excep);          << 
543     }                                          << 
544   //I have to know in advance how many points  << 
545   //to initialize the G4PhysicsFreeVector()    << 
546   std::size_t ndata=0;                         << 
547   G4String line;                               << 
548   while( getline(file, line) )                 << 
549     ndata++;                                   << 
550   ndata -= 1;                                  << 
551   //G4cout << "Found: " << ndata << " lines" < << 
552                                                << 
553   file.clear();                                << 
554   file.close();                                << 
555   file.open(ost.str().c_str());                << 
556                                                << 
557   G4int readZ =0;                              << 
558   std::size_t nShells= 0;                      << 
559   file >> readZ >> nShells;                    << 
560                                                << 
561   if (fVerboseLevel > 3)                       << 
562     G4cout << "Element Z=" << Z << " , nShells << 
563                                                << 
564   //check the right file is opened.            << 
565   if (readZ != Z || nShells <= 0 || nShells >  << 
566     {                                          << 
567       G4ExceptionDescription ed;               << 
568       ed << "Corrupted data file for Z=" << Z  << 
569       G4Exception("G4PenelopePhotoElectricMode << 
570       "em0005",FatalException,ed);             << 
571       return;                                  << 
572     }                                          << 
573   G4PhysicsTable* thePhysicsTable = new G4Phys << 
574                                                << 
575   //the table has to contain nShell+1 G4Physic << 
576   //(theTable)[0] --> total cross section      << 
577   //(theTable)[ishell] --> cross section for s << 
578                                                << 
579   //reserve space for the vectors              << 
580   //everything is log-log                      << 
581   for (std::size_t i=0;i<nShells+1;++i)        << 
582     thePhysicsTable->push_back(new G4PhysicsFr << 
583                                                << 
584   std::size_t k =0;                            << 
585   for (k=0;k<ndata && !file.eof();++k)         << 
586     {                                          << 
587       G4double energy = 0;                     << 
588       G4double aValue = 0;                     << 
589       file >> energy ;                         << 
590       energy *= eV;                            << 
591       G4double logene = G4Log(energy);         << 
592       //loop on the columns                    << 
593       for (std::size_t i=0;i<nShells+1;++i)    << 
594   {                                            << 
595     file >> aValue;                            << 
596     aValue *= barn;                            << 
597     G4PhysicsFreeVector* theVec = (G4PhysicsFr << 
598     if (aValue < 1e-40*cm2) //protection again << 
599       aValue = 1e-40*cm2;                      << 
600     theVec->PutValue(k,logene,G4Log(aValue));  << 
601   }                                            << 
602     }                                          << 
603                                                << 
604   if (fVerboseLevel > 2)                       << 
605     {                                          << 
606       G4cout << "G4PenelopePhotoElectricModel: << 
607        << Z << G4endl;                         << 
608     }                                          << 
609                                                << 
610   fLogAtomicShellXS[Z] = thePhysicsTable;      << 
611                                                << 
612   file.close();                                << 
613   return;                                      << 
614 }                                              << 
615                                                << 
616 //....oooOO0OOooo........oooOO0OOooo........oo << 
617                                                << 
618 std::size_t G4PenelopePhotoElectricModel::GetN << 
619 {                                              << 
620   if (!IsMaster())                             << 
621     //Should not be here!                      << 
622     G4Exception("G4PenelopePhotoElectricModel: << 
623     "em0100",FatalException,"Worker thread in  << 
624                                                << 
625   //read data files                            << 
626   if (!fLogAtomicShellXS[Z])                   << 
627     ReadDataFile(Z);                           << 
628   //now it should be ok                        << 
629   if (!fLogAtomicShellXS[Z])                   << 
630      {                                         << 
631        G4ExceptionDescription ed;              << 
632        ed << "Cannot find shell cross section  << 
633        G4Exception("G4PenelopePhotoElectricMod << 
634        "em2038",FatalException,ed);            << 
635      }                                         << 
636   //one vector is allocated for the _total_ cr << 
637   std::size_t nEntries = fLogAtomicShellXS[Z]- << 
638   return  (nEntries-1);                        << 
639 }                                              << 
640                                                << 
641 //....oooOO0OOooo........oooOO0OOooo........oo << 
642                                                << 
643 G4double G4PenelopePhotoElectricModel::GetShel << 
644 {                                              << 
645   //this forces also the loading of the data   << 
646   std::size_t entries = GetNumberOfShellXS(Z); << 
647                                                << 
648   if (shellID >= entries)                      << 
649     {                                          << 
650       G4cout << "Element Z=" << Z << " has dat << 
651       G4cout << "so shellID should be from 0 t << 
652       return 0;                                << 
653     }                                          << 
654                                                << 
655   G4PhysicsTable* theTable =  fLogAtomicShellX << 
656   //[0] is the total XS, shellID is in the ele << 
657   G4PhysicsFreeVector* totalXSLog = (G4Physics << 
658                                                << 
659   if (!totalXSLog)                             << 
660      {                                         << 
661        G4Exception("G4PenelopePhotoElectricMod << 
662        "em2039",FatalException,                << 
663        "Unable to retrieve the total cross sec << 
664        return 0;                               << 
665      }                                         << 
666    G4double logene = G4Log(energy);            << 
667    G4double logXS = totalXSLog->Value(logene); << 
668    G4double cross = G4Exp(logXS);              << 
669    if (cross < 2e-40*cm2) cross = 0;           << 
670    return cross;                               << 
671 }                                              << 
672                                                << 
673 //....oooOO0OOooo........oooOO0OOooo........oo << 
674                                                << 
675 G4String G4PenelopePhotoElectricModel::WriteTa << 
676 {                                              << 
677   G4String theShell = "outer shell";           << 
678   if (shellID == 0)                            << 
679     theShell = "K";                            << 
680   else if (shellID == 1)                       << 
681     theShell = "L1";                           << 
682   else if (shellID == 2)                       << 
683     theShell = "L2";                           << 
684   else if (shellID == 3)                       << 
685     theShell = "L3";                           << 
686   else if (shellID == 4)                       << 
687     theShell = "M1";                           << 
688   else if (shellID == 5)                       << 
689     theShell = "M2";                           << 
690   else if (shellID == 6)                       << 
691     theShell = "M3";                           << 
692   else if (shellID == 7)                       << 
693     theShell = "M4";                           << 
694   else if (shellID == 8)                       << 
695     theShell = "M5";                           << 
696                                                << 
697   return theShell;                             << 
698 }                                              << 
699                                                << 
700 //....oooOO0OOooo........oooOO0OOooo........oo << 
701                                                << 
702 void G4PenelopePhotoElectricModel::SetParticle << 
703 {                                              << 
704   if(!fParticle) {                             << 
705     fParticle = p;                             << 
706   }                                            << 
707 }                                              << 
708                                                << 
709 //....oooOO0OOooo........oooOO0OOooo........oo << 
710                                                << 
711 std::size_t G4PenelopePhotoElectricModel::Sele << 
712 {                                              << 
713   G4double logEnergy = G4Log(energy);          << 
714                                                << 
715   //Check if data have been read (it should be << 
716   if (!fLogAtomicShellXS[Z])                   << 
717      {                                         << 
718        G4ExceptionDescription ed;              << 
719        ed << "Cannot find shell cross section  << 
720        G4Exception("G4PenelopePhotoElectricMod << 
721        "em2038",FatalException,ed);            << 
722      }                                         << 
723                                                << 
724   G4PhysicsTable* theTable =  fLogAtomicShellX << 
725                                                << 
726   G4double sum = 0;                            << 
727   G4PhysicsFreeVector* totalXSLog = (G4Physics << 
728   G4double logXS = totalXSLog->Value(logEnergy << 
729   G4double totalXS = G4Exp(logXS);             << 
730                                                << 
731   //Notice: totalXS is the total cross section << 
732   //the sum of partialXS's, since these includ << 
733   //                                           << 
734   // Therefore, here one have to consider the  << 
735   // an outer shell. Conventionally, it is ind << 
736   //                                           << 
737   G4double random = G4UniformRand()*totalXS;   << 
738                                                << 
739   for (std::size_t k=1;k<theTable->entries();+ << 
740     {                                          << 
741       //Add one shell                          << 
742       G4PhysicsFreeVector* partialXSLog = (G4P << 
743       G4double logXSLocal = partialXSLog->Valu << 
744       G4double partialXS = G4Exp(logXSLocal);  << 
745       sum += partialXS;                        << 
746       if (random <= sum)                       << 
747   return k-1;                                  << 
748     }                                          << 
749   //none of the shells K, L, M: return outer s << 
750   return 9;                                    << 
751 }                                              << 
752                                                   417