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 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermorePhotoElectricModel.cc (Version 6.2.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 // Author: Sebastien Incerti                      
 28 //         30 October 2008                        
 29 //         on base of G4LowEnergyPhotoElectric    
 30 //                                                
 31 // 22 Oct 2012   A & V Ivanchenko Migration da    
 32 // 1 June 2017   M Bandieramonte: New model ba    
 33 //               evaluated data - parameteriza    
 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.h    
 46 #include "G4SystemOfUnits.hh"                     
 47 #include "G4VAtomDeexcitation.hh"                 
 48 #include "G4EmParameters.hh"                      
 49 #include <thread>                                 
 50                                                   
 51 //....oooOO0OOooo........oooOO0OOooo........oo    
 52                                                   
 53 G4ElementData* G4LivermorePhotoElectricModel::    
 54 G4ElementData* G4LivermorePhotoElectricModel::    
 55 std::vector<G4double>* G4LivermorePhotoElectri    
 56 std::vector<G4double>* G4LivermorePhotoElectri    
 57 G4int G4LivermorePhotoElectricModel::fNShells[    
 58 G4int G4LivermorePhotoElectricModel::fNShellsU    
 59 G4Material* G4LivermorePhotoElectricModel::fWa    
 60 G4double G4LivermorePhotoElectricModel::fWater    
 61 G4String G4LivermorePhotoElectricModel::fDataD    
 62                                                   
 63 static std::once_flag applyOnce;                  
 64                                                   
 65 namespace                                         
 66 {                                                 
 67   G4Mutex livPhotoeffMutex = G4MUTEX_INITIALIZ    
 68 }                                                 
 69                                                   
 70 //....oooOO0OOooo........oooOO0OOooo........oo    
 71                                                   
 72 G4LivermorePhotoElectricModel::G4LivermorePhot    
 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 o    
 80   // 4 = entering in methods                      
 81                                                   
 82   theGamma = G4Gamma::Gamma();                    
 83   theElectron = G4Electron::Electron();           
 84                                                   
 85   // default generator                            
 86   SetAngularDistribution(new G4SauterGavrilaAn    
 87                                                   
 88   if (verboseLevel > 0) {                         
 89     G4cout << "Livermore PhotoElectric is cons    
 90            << " nShellLimit= " << nShellLimit     
 91   }                                               
 92                                                   
 93   // Mark this model as "applicable" for atomi    
 94   SetDeexcitationFlag(true);                      
 95                                                   
 96   // For water                                    
 97   fSandiaCof.resize(4, 0.0);                      
 98 }                                                 
 99                                                   
100 //....oooOO0OOooo........oooOO0OOooo........oo    
101                                                   
102 G4LivermorePhotoElectricModel::~G4LivermorePho    
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........oo    
119                                                   
120 void G4LivermorePhotoElectricModel::Initialise    
121                  const G4DataVector&)             
122 {                                                 
123   if (verboseLevel > 1) {                         
124     G4cout << "Calling G4LivermorePhotoElectri    
125   }                                               
126                                                   
127   // initialise static tables only once           
128   std::call_once(applyOnce, [this]() { isIniti    
129                                                   
130   if (isInitializer) {                            
131     G4AutoLock l(&livPhotoeffMutex);              
132     FindDirectoryPath();                          
133     if (fWater == nullptr) {                      
134       fWater = G4Material::GetMaterial("G4_WAT    
135       if (fWater == nullptr) {                    
136         fWater = G4Material::GetMaterial("Wate    
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(ZMAX    
147       fCrossSectionLE->SetName("PhotoEffLowXS"    
148     }                                             
149                                                   
150     const G4ElementTable* elemTable = G4Elemen    
151     std::size_t numElems = (*elemTable).size()    
152     for (std::size_t ie = 0; ie < numElems; ++    
153       const G4Element* elem = (*elemTable)[ie]    
154       G4int Z = elem->GetZasInt();                
155       if (Z < ZMAXPE) {                           
156   if (fCrossSection->GetElementData(Z) == null    
157     ReadData(Z);                                  
158   }                                               
159       }                                           
160     }                                             
161     l.unlock();                                   
162   }                                               
163                                                   
164   if (verboseLevel > 1) {                         
165     G4cout << "Loaded cross section files for     
166   }                                               
167   if (nullptr == fParticleChange) {               
168     fParticleChange = GetParticleChangeForGamm    
169     fAtomDeexcitation = G4LossTableManager::In    
170   }                                               
171                                                   
172   fDeexcitationActive = false;                    
173   if (nullptr != fAtomDeexcitation) {             
174     fDeexcitationActive = fAtomDeexcitation->I    
175   }                                               
176                                                   
177   if (verboseLevel > 1) {                         
178     G4cout << "LivermorePhotoElectric model is    
179   }                                               
180 }                                                 
181                                                   
182 //....oooOO0OOooo........oooOO0OOooo........oo    
183                                                   
184 G4double                                          
185 G4LivermorePhotoElectricModel::CrossSectionPer    
186                                                   
187                                                   
188                  G4double, G4double)              
189 {                                                 
190   fCurrSection = 0.0;                             
191   if (fWater && (material == fWater || materia    
192     if (energy <= fWaterEnergyLimit) {            
193       fWater->GetSandiaTable()->GetSandiaCofWa    
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    
201       fSandiaCof[2] / energy3 + fSandiaCof[3]     
202     }                                             
203   }                                               
204   if (0.0 == fCurrSection) {                      
205     fCurrSection = G4VEmModel::CrossSectionPer    
206   }                                               
207   return fCurrSection;                            
208 }                                                 
209                                                   
210 //....oooOO0OOooo........oooOO0OOooo........oo    
211                                                   
212 G4double                                          
213 G4LivermorePhotoElectricModel::ComputeCrossSec    
214                                                   
215                 G4double ZZ, G4double,            
216                 G4double, G4double)               
217 {                                                 
218   if (verboseLevel > 3) {                         
219     G4cout << "\n G4LivermorePhotoElectricMode    
220            << " Z= " << ZZ << "  R(keV)= " <<     
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) == null    
231     InitialiseOnFly(Z);                           
232     if (fCrossSection->GetElementData(Z) == nu    
233   }                                               
234                                                   
235   // 7: rows in the parameterization file; 5:     
236   G4int idx = fNShells[Z] * 7 - 5;                
237                                                   
238   energy = std::max(energy, (*(fParamHigh[Z]))    
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 * (*(fPara    
251        + x2 * (*(fParamHigh[Z]))[idx + 2] + x3    
252        + x4 * (*(fParamHigh[Z]))[idx + 4] + x5    
253   }                                               
254   // low energy parameterisation                  
255   else if (energy >= (*(fParamLow[Z]))[0]) {      
256     G4double x4 = x2 * x2;                        
257     G4double x5 = x4 * x1;  // this variable u    
258     cs = x1 *                                     
259       ((*(fParamLow[Z]))[idx] + x1 * (*(fParam    
260        + x2 * (*(fParamLow[Z]))[idx + 2] + x3     
261        + x4 * (*(fParamLow[Z]))[idx + 4] + x5     
262   }                                               
263   // Tabulated values above k-shell ionization    
264   else if (energy >= (*(fParamHigh[Z]))[1]) {     
265     cs = x3 * (fCrossSection->GetElementData(Z    
266   }                                               
267   // Tabulated values below k-shell ionization    
268   else {                                          
269     cs = x3 * (fCrossSectionLE->GetElementData    
270   }                                               
271   if (verboseLevel > 1) {                         
272     G4cout << "G4LivermorePhotoElectricModel:     
273            << " cross(barn)= " << cs / barn <<    
274   }                                               
275   return cs;                                      
276 }                                                 
277                                                   
278 //....oooOO0OOooo........oooOO0OOooo........oo    
279                                                   
280 void G4LivermorePhotoElectricModel::SampleSeco    
281                                                   
282                                                   
283                                                   
284 {                                                 
285   G4double gammaEnergy = aDynamicGamma->GetKin    
286   if (verboseLevel > 3) {                         
287     G4cout << "G4LivermorePhotoElectricModel::    
288            << gammaEnergy / keV << G4endl;        
289   }                                               
290                                                   
291   // kill incident photon                         
292   fParticleChange->ProposeTrackStatus(fStopAnd    
293   fParticleChange->SetProposedKineticEnergy(0.    
294                                                   
295   // low-energy photo-effect in water - full a    
296   const G4Material* material = couple->GetMate    
297   if (fWater && (material == fWater || materia    
298     if (gammaEnergy <= fWaterEnergyLimit) {       
299       fParticleChange->ProposeLocalEnergyDepos    
300       return;                                     
301     }                                             
302   }                                               
303                                                   
304   // Returns the normalized direction of the m    
305   G4ThreeVector photonDirection = aDynamicGamm    
306                                                   
307   // Select randomly one element in the curren    
308   const G4Element* elm = SelectRandomAtom(mate    
309   G4int Z = elm->GetZasInt();                     
310                                                   
311   // Select the ionised shell in the current a    
312   // cross sections                               
313                                                   
314   // If element was not initialised gamma shou    
315   if (Z >= ZMAXPE || fCrossSection->GetElement    
316     fParticleChange->ProposeLocalEnergyDeposit    
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 n    
332       // so cross section is not real             
333                                                   
334       G4double rand = G4UniformRand();            
335       G4double cs0 = rand *                       
336         ((*(fParamHigh[Z]))[idx] + x1 * (*(fPa    
337         + x2 * (*(fParamHigh[Z]))[idx + 2] + x    
338         + x4 * (*(fParamHigh[Z]))[idx + 4] + x    
339                                                   
340       for (shellIdx = 0; shellIdx < nn; ++shel    
341         idx = shellIdx * 7 + 2;                   
342         if (gammaEnergy > (*(fParamHigh[Z]))[i    
343           G4double cs = (*(fParamHigh[Z]))[idx    
344             + x2 * (*(fParamHigh[Z]))[idx + 2]    
345             + x4 * (*(fParamHigh[Z]))[idx + 4]    
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]))[    
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 n    
364       // so cross section is not real             
365       G4double cs0 = G4UniformRand() *            
366         ((*(fParamLow[Z]))[idx] + x1 * (*(fPar    
367         + x2 * (*(fParamLow[Z]))[idx + 2] + x3    
368         + x4 * (*(fParamLow[Z]))[idx + 4] + x5    
369       for (shellIdx = 0; shellIdx < nn; ++shel    
370         idx = shellIdx * 7 + 2;                   
371         if (gammaEnergy > (*(fParamLow[Z]))[id    
372           G4double cs = (*(fParamLow[Z]))[idx]    
373             + x2 * (*(fParamLow[Z]))[idx + 2]     
374             + x4 * (*(fParamLow[Z]))[idx + 4]     
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 n    
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)    
392       }                                           
393       else {                                      
394         // below K-shell binding energy           
395         cs *= fCrossSectionLE->GetElementData(    
396       }                                           
397                                                   
398       for (G4int j = 0; j < (G4int)nn; ++j) {     
399         shellIdx = (std::size_t)fCrossSection-    
400         if (gammaEnergy > (*(fParamLow[Z]))[7     
401           cs -= fCrossSection->GetValueForComp    
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]))[    
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(    
418   }                                               
419                                                   
420   // If binding energy of the selected shell i    
421   //    do not generate secondaries               
422   if (gammaEnergy < bindingEnergy) {              
423     fParticleChange->ProposeLocalEnergyDeposit    
424     return;                                       
425   }                                               
426                                                   
427   // Primary outcoming electron                   
428   G4double eKineticEnergy = gammaEnergy - bind    
429   G4double edep = bindingEnergy;                  
430                                                   
431   // Calculate direction of the photoelectron     
432   G4ThreeVector electronDirection = GetAngular    
433     aDynamicGamma, eKineticEnergy, (G4int)shel    
434                                                   
435   // The electron is created                      
436   auto electron =                                 
437     new G4DynamicParticle(theElectron, electro    
438   fvect->push_back(electron);                     
439                                                   
440   // Sample deexcitation                          
441   if (nullptr != shell) {                         
442     G4int index = couple->GetIndex();             
443     if (fAtomDeexcitation->CheckDeexcitationAc    
444       std::size_t nbefore = fvect->size();        
445                                                   
446       fAtomDeexcitation->GenerateParticles(fve    
447       std::size_t nafter = fvect->size();         
448       if (nafter > nbefore) {                     
449         G4double esec = 0.0;                      
450         for (std::size_t j = nbefore; j < naft    
451           G4double e = ((*fvect)[j])->GetKinet    
452           if (esec + e > edep) {                  
453             // correct energy in order to have    
454             e = edep - esec;                      
455             ((*fvect)[j])->SetKineticEnergy(e)    
456             esec += e;                            
457             // delete the rest of secondaries     
458             for (std::size_t jj = nafter - 1;     
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    
473   }                                               
474 }                                                 
475                                                   
476 //....oooOO0OOooo........oooOO0OOooo........oo    
477                                                   
478 const G4String& G4LivermorePhotoElectricModel:    
479 {                                                 
480   // no check in this method - environment var    
481   if (fDataDirectory.empty()) {                   
482     auto param = G4EmParameters::Instance();      
483     std::ostringstream ost;                       
484     if (param->LivermoreDataDir() == "livermor    
485       ost << param->GetDirLEDATA() << "/liverm    
486     }                                             
487     else {                                        
488       ost << param->GetDirLEDATA() << "/epics2    
489     }                                             
490     fDataDirectory = ost.str();                   
491   }                                               
492   return fDataDirectory;                          
493 }                                                 
494                                                   
495 //....oooOO0OOooo........oooOO0OOooo........oo    
496                                                   
497 void G4LivermorePhotoElectricModel::ReadData(G    
498 {                                                 
499   if (Z <= 0 || Z >= ZMAXPE) {                    
500     G4cout << "G4LivermorePhotoElectricModel::    
501      << Z << " is out of range - request ignor    
502     return;                                       
503   }                                               
504   if (verboseLevel > 1) {                         
505     G4cout << "G4LivermorePhotoElectricModel::    
506   }                                               
507                                                   
508   if (fCrossSection->GetElementData(Z) != null    
509     return;                                       
510   }                                               
511                                                   
512   // spline for photoeffect total x-section ab    
513   // but below the parameterized ones             
514                                                   
515   G4bool spline = (G4EmParameters::Instance()-    
516   G4int number = G4EmParameters::Instance()->N    
517   auto pv = new G4PhysicsFreeVector(spline);      
518                                                   
519   // fDataDirectory will be defined after thes    
520   std::ostringstream ost;                         
521   ost << FindDirectoryPath() << "pe-cs-" << Z     
522   std::ifstream fin(ost.str().c_str());           
523   if (!fin.is_open()) {                           
524     G4ExceptionDescription ed;                    
525     ed << "G4LivermorePhotoElectricModel data     
526        << ost.str().c_str() << "> is not opene    
527     G4Exception("G4LivermorePhotoElectricModel    
528                 "G4LEDATA version should be G4    
529     return;                                       
530   }                                               
531   if (verboseLevel > 3) {                         
532     G4cout << "File " << ost.str().c_str() <<     
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 <<    
549   std::ifstream fin1(ost1.str().c_str());         
550   if (!fin1.is_open()) {                          
551     G4ExceptionDescription ed;                    
552     ed << "G4LivermorePhotoElectricModel data     
553        << ost1.str().c_str() << "> is not open    
554     G4Exception("G4LivermorePhotoElectricModel    
555                 "G4LEDATA version should be G4    
556     return;                                       
557   }                                               
558   if (verboseLevel > 3) {                         
559     G4cout << "File " << ost1.str().c_str() <<    
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    
607   std::ifstream fin1_low(ost1_low.str().c_str(    
608   if (!fin1_low.is_open()) {                      
609     G4ExceptionDescription ed;                    
610     ed << "G4LivermorePhotoElectricModel data     
611        << "> is not opened!" << G4endl;           
612     G4Exception("G4LivermorePhotoElectricModel    
613                 "G4LEDATA version should be G4    
614     return;                                       
615   }                                               
616   if (verboseLevel > 3) {                         
617     G4cout << "File " << ost1_low.str().c_str(    
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 s    
659   if (nShellLimit < n2) {                         
660     n2 = nShellLimit;                             
661   }                                               
662   fCrossSection->InitialiseForComponent(Z, n2)    
663   fNShellsUsed[Z] = n2;                           
664                                                   
665   if (1 < n2) {                                   
666     std::ostringstream ost2;                      
667     ost2 << fDataDirectory << "pe-ss-cs-" << Z    
668     std::ifstream fin2(ost2.str().c_str());       
669     if (!fin2.is_open()) {                        
670       G4ExceptionDescription ed;                  
671       ed << "G4LivermorePhotoElectricModel dat    
672          << G4endl;                               
673       G4Exception("G4LivermorePhotoElectricMod    
674                   "G4LEDATA version should be     
675       return;                                     
676     }                                             
677     if (verboseLevel > 3) {                       
678       G4cout << "File " << ost2.str().c_str()     
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,     
687       for (G4int j = 0; j < n3; ++j) {            
688         fin2 >> x >> y;                           
689         v->PutValues(j, x * CLHEP::MeV, y * CL    
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    
698   if (1 < fNShells[Z]) {                          
699     auto pv1 = new G4PhysicsFreeVector(false);    
700     std::ostringstream ost3;                      
701     ost3 << fDataDirectory << "pe-le-cs-" << Z    
702     std::ifstream fin3(ost3.str().c_str());       
703     if (!fin3.is_open()) {                        
704       G4ExceptionDescription ed;                  
705       ed << "G4LivermorePhotoElectricModel dat    
706          << G4endl;                               
707       G4Exception("G4LivermorePhotoElectricMod    
708                   "G4LEDATA version should be     
709       return;                                     
710     }                                             
711     if (verboseLevel > 3) {                       
712       G4cout << "File " << ost3.str().c_str()     
713              << G4endl;                           
714     }                                             
715     pv1->Retrieve(fin3, true);                    
716     pv1->ScaleVector(CLHEP::MeV, CLHEP::barn);    
717     pv1->EnableLogBinSearch(number);              
718     fCrossSectionLE->InitialiseForElement(Z, p    
719     fin3.close();                                 
720   }                                               
721 }                                                 
722                                                   
723 //....oooOO0OOooo........oooOO0OOooo........oo    
724                                                   
725 G4double G4LivermorePhotoElectricModel::GetBin    
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) == null    
733       shell < 0 || shell >= fNShellsUsed[Z]) {    
734     return -1;                                    
735   }                                               
736                                                   
737   if (Z > 2) {                                    
738     return fCrossSection->GetComponentDataByIn    
739   }                                               
740   else {                                          
741     return fCrossSection->GetElementData(Z)->E    
742   }                                               
743 }                                                 
744                                                   
745 //....oooOO0OOooo........oooOO0OOooo........oo    
746                                                   
747 void                                              
748 G4LivermorePhotoElectricModel::InitialiseForEl    
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........oo    
760                                                   
761 void G4LivermorePhotoElectricModel::Initialise    
762 {                                                 
763   if (fCrossSection->GetElementData(Z) == null    
764     G4AutoLock l(&livPhotoeffMutex);              
765     if (fCrossSection->GetElementData(Z) == nu    
766       ReadData(Z);                                
767     }                                             
768     l.unlock();                                   
769   }                                               
770 }                                                 
771                                                   
772 //....oooOO0OOooo........oooOO0OOooo........oo    
773