Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 // Author: Sebastien Incerti
 28 //         30 October 2008
 29 //         on base of G4LowEnergyPhotoElectric developed by A.Forti and M.G.Pia
 30 //
 31 // 22 Oct 2012   A & V Ivanchenko Migration data structure to G4PhysicsVector
 32 // 1 June 2017   M Bandieramonte: New model based on livermore/epics2014
 33 //               evaluated data - parameterization fits in two ranges
 34 
 35 #include "G4LivermorePhotoElectricModel.hh"
 36 
 37 #include "G4AtomicShell.hh"
 38 #include "G4AutoLock.hh"
 39 #include "G4CrossSectionHandler.hh"
 40 #include "G4Electron.hh"
 41 #include "G4Gamma.hh"
 42 #include "G4LossTableManager.hh"
 43 #include "G4ParticleChangeForGamma.hh"
 44 #include "G4PhysicsFreeVector.hh"
 45 #include "G4SauterGavrilaAngularDistribution.hh"
 46 #include "G4SystemOfUnits.hh"
 47 #include "G4VAtomDeexcitation.hh"
 48 #include "G4EmParameters.hh"
 49 #include <thread>
 50 
 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 52 
 53 G4ElementData* G4LivermorePhotoElectricModel::fCrossSection = nullptr;
 54 G4ElementData* G4LivermorePhotoElectricModel::fCrossSectionLE = nullptr;
 55 std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {nullptr};
 56 std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {nullptr};
 57 G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
 58 G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
 59 G4Material* G4LivermorePhotoElectricModel::fWater = nullptr;
 60 G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
 61 G4String G4LivermorePhotoElectricModel::fDataDirectory = "";
 62 
 63 static std::once_flag applyOnce;
 64 
 65 namespace
 66 {
 67   G4Mutex livPhotoeffMutex = G4MUTEX_INITIALIZER;
 68 }
 69 
 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 71 
 72 G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4String& nam) : G4VEmModel(nam)
 73 {
 74   verboseLevel = 0;
 75   // Verbosity scale:
 76   // 0 = nothing
 77   // 1 = warning for energy non-conservation
 78   // 2 = details of energy budget
 79   // 3 = calculation of cross sections, file openings, sampling of atoms
 80   // 4 = entering in methods
 81 
 82   theGamma = G4Gamma::Gamma();
 83   theElectron = G4Electron::Electron();
 84 
 85   // default generator
 86   SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
 87 
 88   if (verboseLevel > 0) {
 89     G4cout << "Livermore PhotoElectric is constructed "
 90            << " nShellLimit= " << nShellLimit << G4endl;
 91   }
 92 
 93   // Mark this model as "applicable" for atomic deexcitation
 94   SetDeexcitationFlag(true);
 95 
 96   // For water
 97   fSandiaCof.resize(4, 0.0);
 98 }
 99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
103 {
104   if (isInitializer) {
105     for (G4int i = 0; i < ZMAXPE; ++i) {
106       if (fParamHigh[i]) {
107         delete fParamHigh[i];
108         fParamHigh[i] = nullptr;
109       }
110       if (fParamLow[i]) {
111         delete fParamLow[i];
112         fParamLow[i] = nullptr;
113       }
114     }
115   }
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
120 void G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
121                  const G4DataVector&)
122 {
123   if (verboseLevel > 1) {
124     G4cout << "Calling G4LivermorePhotoElectricModel::Initialise() " << G4endl;
125   }
126 
127   // initialise static tables only once
128   std::call_once(applyOnce, [this]() { isInitializer = true; });
129 
130   if (isInitializer) {
131     G4AutoLock l(&livPhotoeffMutex);
132     FindDirectoryPath();
133     if (fWater == nullptr) {
134       fWater = G4Material::GetMaterial("G4_WATER", false);
135       if (fWater == nullptr) {
136         fWater = G4Material::GetMaterial("Water", false);
137       }
138       if (fWater != nullptr) {
139         fWaterEnergyLimit = 13.6 * CLHEP::eV;
140       }
141     }
142 
143     if (fCrossSection == nullptr) {
144       fCrossSection = new G4ElementData(ZMAXPE);
145       fCrossSection->SetName("PhotoEffXS");
146       fCrossSectionLE = new G4ElementData(ZMAXPE);
147       fCrossSectionLE->SetName("PhotoEffLowXS");
148     }
149 
150     const G4ElementTable* elemTable = G4Element::GetElementTable();
151     std::size_t numElems = (*elemTable).size();
152     for (std::size_t ie = 0; ie < numElems; ++ie) {
153       const G4Element* elem = (*elemTable)[ie];
154       G4int Z = elem->GetZasInt();
155       if (Z < ZMAXPE) {
156   if (fCrossSection->GetElementData(Z) == nullptr) {
157     ReadData(Z);
158   }
159       }
160     }
161     l.unlock();
162   }
163 
164   if (verboseLevel > 1) {
165     G4cout << "Loaded cross section files for new LivermorePhotoElectric model" << G4endl;
166   }
167   if (nullptr == fParticleChange) {
168     fParticleChange = GetParticleChangeForGamma();
169     fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
170   }
171 
172   fDeexcitationActive = false;
173   if (nullptr != fAtomDeexcitation) {
174     fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
175   }
176 
177   if (verboseLevel > 1) {
178     G4cout << "LivermorePhotoElectric model is initialized " << G4endl;
179   }
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 
184 G4double
185 G4LivermorePhotoElectricModel::CrossSectionPerVolume(const G4Material* material,
186                                                      const G4ParticleDefinition* p,
187                                                      G4double energy,
188                  G4double, G4double)
189 {
190   fCurrSection = 0.0;
191   if (fWater && (material == fWater || material->GetBaseMaterial() == fWater)) {
192     if (energy <= fWaterEnergyLimit) {
193       fWater->GetSandiaTable()->GetSandiaCofWater(energy, fSandiaCof);
194 
195       G4double energy2 = energy * energy;
196       G4double energy3 = energy * energy2;
197       G4double energy4 = energy2 * energy2;
198 
199       fCurrSection = material->GetDensity()
200                      * (fSandiaCof[0] / energy + fSandiaCof[1] / energy2 +
201       fSandiaCof[2] / energy3 + fSandiaCof[3] / energy4);
202     }
203   }
204   if (0.0 == fCurrSection) {
205     fCurrSection = G4VEmModel::CrossSectionPerVolume(material, p, energy);
206   }
207   return fCurrSection;
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 
212 G4double
213 G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
214                                                           G4double energy,
215                 G4double ZZ, G4double,
216                 G4double, G4double)
217 {
218   if (verboseLevel > 3) {
219     G4cout << "\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
220            << " Z= " << ZZ << "  R(keV)= " << energy / keV << G4endl;
221   }
222   G4double cs = 0.0;
223   G4int Z = G4lrint(ZZ);
224   if (Z >= ZMAXPE || Z <= 0) {
225     return cs;
226   }
227   // if element was not initialised
228 
229   // do initialisation safely for MT mode
230   if (fCrossSection->GetElementData(Z) == nullptr) {
231     InitialiseOnFly(Z);
232     if (fCrossSection->GetElementData(Z) == nullptr) { return cs; }
233   }
234 
235   // 7: rows in the parameterization file; 5: number of parameters
236   G4int idx = fNShells[Z] * 7 - 5;
237 
238   energy = std::max(energy, (*(fParamHigh[Z]))[idx - 1]);
239 
240   G4double x1 = 1.0 / energy;
241   G4double x2 = x1 * x1;
242   G4double x3 = x2 * x1;
243 
244   // high energy parameterisation
245   if (energy >= (*(fParamHigh[Z]))[0]) {
246     G4double x4 = x2 * x2;
247     G4double x5 = x4 * x1;
248 
249     cs = x1 *
250       ((*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
251        + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
252        + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5]);
253   }
254   // low energy parameterisation
255   else if (energy >= (*(fParamLow[Z]))[0]) {
256     G4double x4 = x2 * x2;
257     G4double x5 = x4 * x1;  // this variable usage can probably be optimized
258     cs = x1 *
259       ((*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
260        + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
261        + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5]);
262   }
263   // Tabulated values above k-shell ionization energy
264   else if (energy >= (*(fParamHigh[Z]))[1]) {
265     cs = x3 * (fCrossSection->GetElementData(Z))->Value(energy);
266   }
267   // Tabulated values below k-shell ionization energy
268   else {
269     cs = x3 * (fCrossSectionLE->GetElementData(Z))->Value(energy);
270   }
271   if (verboseLevel > 1) {
272     G4cout << "G4LivermorePhotoElectricModel: E(keV)= " << energy / keV << " Z= " << Z
273            << " cross(barn)= " << cs / barn << G4endl;
274   }
275   return cs;
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
280 void G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
281                                                       const G4MaterialCutsCouple* couple,
282                                                       const G4DynamicParticle* aDynamicGamma,
283                                                       G4double, G4double)
284 {
285   G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
286   if (verboseLevel > 3) {
287     G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
288            << gammaEnergy / keV << G4endl;
289   }
290 
291   // kill incident photon
292   fParticleChange->ProposeTrackStatus(fStopAndKill);
293   fParticleChange->SetProposedKineticEnergy(0.);
294 
295   // low-energy photo-effect in water - full absorption
296   const G4Material* material = couple->GetMaterial();
297   if (fWater && (material == fWater || material->GetBaseMaterial() == fWater)) {
298     if (gammaEnergy <= fWaterEnergyLimit) {
299       fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
300       return;
301     }
302   }
303 
304   // Returns the normalized direction of the momentum
305   G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
306 
307   // Select randomly one element in the current material
308   const G4Element* elm = SelectRandomAtom(material, theGamma, gammaEnergy);
309   G4int Z = elm->GetZasInt();
310 
311   // Select the ionised shell in the current atom according to shell
312   // cross sections
313 
314   // If element was not initialised gamma should be absorbed
315   if (Z >= ZMAXPE || fCrossSection->GetElementData(Z) == nullptr) {
316     fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
317     return;
318   }
319 
320   // SAMPLING OF THE SHELL INDEX
321   std::size_t shellIdx = 0;
322   std::size_t nn = fNShellsUsed[Z];
323   if (nn > 1) {
324     if (gammaEnergy >= (*(fParamHigh[Z]))[0]) {
325       G4double x1 = 1.0 / gammaEnergy;
326       G4double x2 = x1 * x1;
327       G4double x3 = x2 * x1;
328       G4double x4 = x3 * x1;
329       G4double x5 = x4 * x1;
330       std::size_t idx = nn * 7 - 5;
331       // when do sampling common factors are not taken into account
332       // so cross section is not real
333 
334       G4double rand = G4UniformRand();
335       G4double cs0 = rand *
336         ((*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
337         + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
338         + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5]);
339 
340       for (shellIdx = 0; shellIdx < nn; ++shellIdx) {
341         idx = shellIdx * 7 + 2;
342         if (gammaEnergy > (*(fParamHigh[Z]))[idx - 1]) {
343           G4double cs = (*(fParamHigh[Z]))[idx] + x1 * (*(fParamHigh[Z]))[idx + 1]
344             + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 * (*(fParamHigh[Z]))[idx + 3]
345             + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 * (*(fParamHigh[Z]))[idx + 5];
346 
347           if (cs >= cs0) {
348             break;
349           }
350         }
351       }
352       if (shellIdx >= nn) {
353         shellIdx = nn - 1;
354       }
355     }
356     else if (gammaEnergy >= (*(fParamLow[Z]))[0]) {
357       G4double x1 = 1.0 / gammaEnergy;
358       G4double x2 = x1 * x1;
359       G4double x3 = x2 * x1;
360       G4double x4 = x3 * x1;
361       G4double x5 = x4 * x1;
362       std::size_t idx = nn * 7 - 5;
363       // when do sampling common factors are not taken into account
364       // so cross section is not real
365       G4double cs0 = G4UniformRand() *
366         ((*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
367         + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
368         + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5]);
369       for (shellIdx = 0; shellIdx < nn; ++shellIdx) {
370         idx = shellIdx * 7 + 2;
371         if (gammaEnergy > (*(fParamLow[Z]))[idx - 1]) {
372           G4double cs = (*(fParamLow[Z]))[idx] + x1 * (*(fParamLow[Z]))[idx + 1]
373             + x2 * (*(fParamLow[Z]))[idx + 2] + x3 * (*(fParamLow[Z]))[idx + 3]
374             + x4 * (*(fParamLow[Z]))[idx + 4] + x5 * (*(fParamLow[Z]))[idx + 5];
375           if (cs >= cs0) {
376             break;
377           }
378         }
379       }
380       if (shellIdx >= nn) {
381         shellIdx = nn - 1;
382       }
383     }
384     else {
385       // when do sampling common factors are not taken into account
386       // so cross section is not real
387       G4double cs = G4UniformRand();
388 
389       if (gammaEnergy >= (*(fParamHigh[Z]))[1]) {
390         // above K-shell binding energy
391         cs *= fCrossSection->GetElementData(Z)->Value(gammaEnergy);
392       }
393       else {
394         // below K-shell binding energy
395         cs *= fCrossSectionLE->GetElementData(Z)->Value(gammaEnergy);
396       }
397 
398       for (G4int j = 0; j < (G4int)nn; ++j) {
399         shellIdx = (std::size_t)fCrossSection->GetComponentID(Z, j);
400         if (gammaEnergy > (*(fParamLow[Z]))[7 * shellIdx + 1]) {
401           cs -= fCrossSection->GetValueForComponent(Z, j, gammaEnergy);
402         }
403         if (cs <= 0.0 || j + 1 == (G4int)nn) {
404           break;
405         }
406       }
407     }
408   }
409   // END: SAMPLING OF THE SHELL
410 
411   G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx * 7 + 1];
412   const G4AtomicShell* shell = nullptr;
413 
414   // no de-excitation from the last shell
415   if (fDeexcitationActive && shellIdx + 1 < nn) {
416     auto as = G4AtomicShellEnumerator(shellIdx);
417     shell = fAtomDeexcitation->GetAtomicShell(Z, as);
418   }
419 
420   // If binding energy of the selected shell is larger than photon energy
421   //    do not generate secondaries
422   if (gammaEnergy < bindingEnergy) {
423     fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
424     return;
425   }
426 
427   // Primary outcoming electron
428   G4double eKineticEnergy = gammaEnergy - bindingEnergy;
429   G4double edep = bindingEnergy;
430 
431   // Calculate direction of the photoelectron
432   G4ThreeVector electronDirection = GetAngularDistribution()->SampleDirection(
433     aDynamicGamma, eKineticEnergy, (G4int)shellIdx, couple->GetMaterial());
434 
435   // The electron is created
436   auto electron =
437     new G4DynamicParticle(theElectron, electronDirection, eKineticEnergy);
438   fvect->push_back(electron);
439 
440   // Sample deexcitation
441   if (nullptr != shell) {
442     G4int index = couple->GetIndex();
443     if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
444       std::size_t nbefore = fvect->size();
445 
446       fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
447       std::size_t nafter = fvect->size();
448       if (nafter > nbefore) {
449         G4double esec = 0.0;
450         for (std::size_t j = nbefore; j < nafter; ++j) {
451           G4double e = ((*fvect)[j])->GetKineticEnergy();
452           if (esec + e > edep) {
453             // correct energy in order to have energy balance
454             e = edep - esec;
455             ((*fvect)[j])->SetKineticEnergy(e);
456             esec += e;
457             // delete the rest of secondaries (should not happens)
458             for (std::size_t jj = nafter - 1; jj > j; --jj) {
459               delete (*fvect)[jj];
460               fvect->pop_back();
461             }
462             break;
463           }
464           esec += e;
465         }
466         edep -= esec;
467       }
468     }
469   }
470   // energy balance - excitation energy left
471   if (edep > 0.0) {
472     fParticleChange->ProposeLocalEnergyDeposit(edep);
473   }
474 }
475 
476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
477 
478 const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
479 {
480   // no check in this method - environment variable is check by utility
481   if (fDataDirectory.empty()) {
482     auto param = G4EmParameters::Instance();
483     std::ostringstream ost;
484     if (param->LivermoreDataDir() == "livermore") {
485       ost << param->GetDirLEDATA() << "/livermore/phot_epics2014/";
486     }
487     else {
488       ost << param->GetDirLEDATA() << "/epics2017/phot/";
489     }
490     fDataDirectory = ost.str();
491   }
492   return fDataDirectory;
493 }
494 
495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
496 
497 void G4LivermorePhotoElectricModel::ReadData(G4int Z)
498 {
499   if (Z <= 0 || Z >= ZMAXPE) {
500     G4cout << "G4LivermorePhotoElectricModel::ReadData() Warning: Z="
501      << Z << " is out of range - request ignored" << G4endl;
502     return;
503   }
504   if (verboseLevel > 1) {
505     G4cout << "G4LivermorePhotoElectricModel::ReadData() for Z=" << Z << G4endl;
506   }
507 
508   if (fCrossSection->GetElementData(Z) != nullptr) {
509     return;
510   }
511 
512   // spline for photoeffect total x-section above K-shell when using EPDL97
513   // but below the parameterized ones
514 
515   G4bool spline = (G4EmParameters::Instance()->LivermoreDataDir() == "livermore");
516   G4int number = G4EmParameters::Instance()->NumberForFreeVector();
517   auto pv = new G4PhysicsFreeVector(spline);
518 
519   // fDataDirectory will be defined after these lines
520   std::ostringstream ost;
521   ost << FindDirectoryPath() << "pe-cs-" << Z << ".dat";
522   std::ifstream fin(ost.str().c_str());
523   if (!fin.is_open()) {
524     G4ExceptionDescription ed;
525     ed << "G4LivermorePhotoElectricModel data file <"
526        << ost.str().c_str() << "> is not opened!" << G4endl;
527     G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
528                 "G4LEDATA version should be G4EMLOW8.0 or later.");
529     return;
530   }
531   if (verboseLevel > 3) {
532     G4cout << "File " << ost.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
533            << G4endl;
534   }
535   pv->Retrieve(fin, true);
536   pv->ScaleVector(MeV, barn);
537   pv->FillSecondDerivatives();
538   pv->EnableLogBinSearch(number);
539   fCrossSection->InitialiseForElement(Z, pv);
540   fin.close();
541 
542   // read high-energy fit parameters
543   fParamHigh[Z] = new std::vector<G4double>;
544   G4int n1 = 0;
545   G4int n2 = 0;
546   G4double x;
547   std::ostringstream ost1;
548   ost1 << fDataDirectory << "pe-high-" << Z << ".dat";
549   std::ifstream fin1(ost1.str().c_str());
550   if (!fin1.is_open()) {
551     G4ExceptionDescription ed;
552     ed << "G4LivermorePhotoElectricModel data file <"
553        << ost1.str().c_str() << "> is not opened!" << G4endl;
554     G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
555                 "G4LEDATA version should be G4EMLOW7.2 or later.");
556     return;
557   }
558   if (verboseLevel > 3) {
559     G4cout << "File " << ost1.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
560            << G4endl;
561   }
562   fin1 >> n1;
563   if (fin1.fail()) {
564     return;
565   }
566   if (0 > n1 || n1 >= INT_MAX) {
567     n1 = 0;
568   }
569 
570   fin1 >> n2;
571   if (fin1.fail()) {
572     return;
573   }
574   if (0 > n2 || n2 >= INT_MAX) {
575     n2 = 0;
576   }
577 
578   fin1 >> x;
579   if (fin1.fail()) {
580     return;
581   }
582 
583   fNShells[Z] = n1;
584   fParamHigh[Z]->reserve(7 * n1 + 1);
585   fParamHigh[Z]->push_back(x * MeV);
586   for (G4int i = 0; i < n1; ++i) {
587     for (G4int j = 0; j < 7; ++j) {
588       fin1 >> x;
589       if (0 == j) {
590         x *= MeV;
591       }
592       else {
593         x *= barn;
594       }
595       fParamHigh[Z]->push_back(x);
596     }
597   }
598   fin1.close();
599 
600   // read low-energy fit parameters
601   fParamLow[Z] = new std::vector<G4double>;
602   G4int n1_low = 0;
603   G4int n2_low = 0;
604   G4double x_low;
605   std::ostringstream ost1_low;
606   ost1_low << fDataDirectory << "pe-low-" << Z << ".dat";
607   std::ifstream fin1_low(ost1_low.str().c_str());
608   if (!fin1_low.is_open()) {
609     G4ExceptionDescription ed;
610     ed << "G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
611        << "> is not opened!" << G4endl;
612     G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
613                 "G4LEDATA version should be G4EMLOW8.0 or later.");
614     return;
615   }
616   if (verboseLevel > 3) {
617     G4cout << "File " << ost1_low.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
618            << G4endl;
619   }
620   fin1_low >> n1_low;
621   if (fin1_low.fail()) {
622     return;
623   }
624   if (0 > n1_low || n1_low >= INT_MAX) {
625     n1_low = 0;
626   }
627 
628   fin1_low >> n2_low;
629   if (fin1_low.fail()) {
630     return;
631   }
632   if (0 > n2_low || n2_low >= INT_MAX) {
633     n2_low = 0;
634   }
635 
636   fin1_low >> x_low;
637   if (fin1_low.fail()) {
638     return;
639   }
640 
641   fNShells[Z] = n1_low;
642   fParamLow[Z]->reserve(7 * n1_low + 1);
643   fParamLow[Z]->push_back(x_low * MeV);
644   for (G4int i = 0; i < n1_low; ++i) {
645     for (G4int j = 0; j < 7; ++j) {
646       fin1_low >> x_low;
647       if (0 == j) {
648         x_low *= MeV;
649       }
650       else {
651         x_low *= barn;
652       }
653       fParamLow[Z]->push_back(x_low);
654     }
655   }
656   fin1_low.close();
657 
658   // there is a possibility to use only main shells
659   if (nShellLimit < n2) {
660     n2 = nShellLimit;
661   }
662   fCrossSection->InitialiseForComponent(Z, n2);  // number of shells
663   fNShellsUsed[Z] = n2;
664 
665   if (1 < n2) {
666     std::ostringstream ost2;
667     ost2 << fDataDirectory << "pe-ss-cs-" << Z << ".dat";
668     std::ifstream fin2(ost2.str().c_str());
669     if (!fin2.is_open()) {
670       G4ExceptionDescription ed;
671       ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str() << "> is not opened!"
672          << G4endl;
673       G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
674                   "G4LEDATA version should be G4EMLOW7.2 or later.");
675       return;
676     }
677     if (verboseLevel > 3) {
678       G4cout << "File " << ost2.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
679              << G4endl;
680     }
681 
682     G4int n3, n4;
683     G4double y;
684     for (G4int i = 0; i < n2; ++i) {
685       fin2 >> x >> y >> n3 >> n4;
686       auto v = new G4PhysicsFreeVector(n3, x, y);
687       for (G4int j = 0; j < n3; ++j) {
688         fin2 >> x >> y;
689         v->PutValues(j, x * CLHEP::MeV, y * CLHEP::barn);
690       }
691       v->EnableLogBinSearch(number);
692       fCrossSection->AddComponent(Z, n4, v);
693     }
694     fin2.close();
695   }
696 
697   // no spline for photoeffect total x-section below K-shell
698   if (1 < fNShells[Z]) {
699     auto pv1 = new G4PhysicsFreeVector(false);
700     std::ostringstream ost3;
701     ost3 << fDataDirectory << "pe-le-cs-" << Z << ".dat";
702     std::ifstream fin3(ost3.str().c_str());
703     if (!fin3.is_open()) {
704       G4ExceptionDescription ed;
705       ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str() << "> is not opened!"
706          << G4endl;
707       G4Exception("G4LivermorePhotoElectricModel::ReadData()", "em0003", FatalException, ed,
708                   "G4LEDATA version should be G4EMLOW8.0 or later.");
709       return;
710     }
711     if (verboseLevel > 3) {
712       G4cout << "File " << ost3.str().c_str() << " is opened by G4LivermorePhotoElectricModel"
713              << G4endl;
714     }
715     pv1->Retrieve(fin3, true);
716     pv1->ScaleVector(CLHEP::MeV, CLHEP::barn);
717     pv1->EnableLogBinSearch(number);
718     fCrossSectionLE->InitialiseForElement(Z, pv1);
719     fin3.close();
720   }
721 }
722 
723 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
724 
725 G4double G4LivermorePhotoElectricModel::GetBindingEnergy(G4int Z, G4int shell)
726 {
727   if (Z < 1 || Z >= ZMAXPE) {
728     return -1;
729   }  // If Z is out of the supported return 0
730 
731   InitialiseOnFly(Z);
732   if (fCrossSection->GetElementData(Z) == nullptr ||
733       shell < 0 || shell >= fNShellsUsed[Z]) {
734     return -1;
735   }
736 
737   if (Z > 2) {
738     return fCrossSection->GetComponentDataByIndex(Z, shell)->Energy(0);
739   }
740   else {
741     return fCrossSection->GetElementData(Z)->Energy(0);
742   }
743 }
744 
745 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
746 
747 void
748 G4LivermorePhotoElectricModel::InitialiseForElement(const G4ParticleDefinition*, G4int Z)
749 {
750   if (fCrossSection == nullptr) {
751     fCrossSection = new G4ElementData(ZMAXPE);
752     fCrossSection->SetName("PhotoEffXS");
753     fCrossSectionLE = new G4ElementData(ZMAXPE);
754     fCrossSectionLE->SetName("PhotoEffLowXS");
755   }
756   ReadData(Z);
757 }
758 
759 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
760 
761 void G4LivermorePhotoElectricModel::InitialiseOnFly(G4int Z)
762 {
763   if (fCrossSection->GetElementData(Z) == nullptr && Z > 0 && Z < ZMAXPE) {
764     G4AutoLock l(&livPhotoeffMutex);
765     if (fCrossSection->GetElementData(Z) == nullptr) {
766       ReadData(Z);
767     }
768     l.unlock();
769   }
770 }
771 
772 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
773