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