52std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {
nullptr};
53std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {
nullptr};
54G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
55G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
56G4ElementData* G4LivermorePhotoElectricModel::fShellCrossSection =
nullptr;
57G4Material* G4LivermorePhotoElectricModel::fWater =
nullptr;
58G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
59G4String G4LivermorePhotoElectricModel::fDataDirectory =
"";
70 :
G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
71 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
72 fAtomDeexcitation(nullptr)
89 G4cout <<
"Livermore PhotoElectric is constructed "
90 <<
" nShellLimit= " << nShellLimit <<
G4endl;
95 fSandiaCof.resize(4,0.0);
104 delete fShellCrossSection;
105 fShellCrossSection =
nullptr;
106 for(
G4int i=0; i<maxZ; ++i) {
107 delete fParamHigh[i];
108 fParamHigh[i] =
nullptr;
110 fParamLow[i] =
nullptr;
111 delete fCrossSection[i];
112 fCrossSection[i] =
nullptr;
113 delete fCrossSectionLE[i];
114 fCrossSectionLE[i] =
nullptr;
125 if (verboseLevel > 2) {
126 G4cout <<
"Calling G4LivermorePhotoElectricModel::Initialise() " <<
G4endl;
134 if(fWater) { fWaterEnergyLimit = 13.6*eV; }
137 if(!fShellCrossSection) { fShellCrossSection =
new G4ElementData(); }
143 for(
G4int i=0; i<numOfCouples; ++i) {
150 for (
G4int j=0; j<nelm; ++j) {
151 G4int Z = std::min(maxZ, (*theElementVector)[j]->GetZasInt());
152 if(!fCrossSection[Z]) { ReadData(Z); }
157 if (verboseLevel > 2) {
158 G4cout <<
"Loaded cross section files for new LivermorePhotoElectric model"
162 isInitialised =
true;
167 fDeexcitationActive =
false;
168 if(fAtomDeexcitation) {
169 fDeexcitationActive = fAtomDeexcitation->
IsFluoActive();
172 if (verboseLevel > 0) {
173 G4cout <<
"LivermorePhotoElectric model is initialized " <<
G4endl
187 if(fWater && (material == fWater ||
189 if(energy <= fWaterEnergyLimit) {
197 (fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
198 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4);
201 if(0.0 == fCurrSection) {
215 if (verboseLevel > 3) {
216 G4cout <<
"\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
217 <<
" Z= " << ZZ <<
" R(keV)= " << energy/keV <<
G4endl;
221 if(Z >= maxZ) {
return cs; }
228 G4int idx = fNShells[Z]*7 - 5;
230 energy = std::max(energy, (*(fParamHigh[Z]))[idx-1]);
237 if(energy >= (*(fParamHigh[Z]))[0]) {
242 cs = x1*((*(fParamHigh[Z]))[idx] + x1*(*(fParamHigh[Z]))[idx+1]
243 + x2*(*(fParamHigh[Z]))[idx+2] + x3*(*(fParamHigh[Z]))[idx+3]
244 + x4*(*(fParamHigh[Z]))[idx+4]+ x5*(*(fParamHigh[Z]))[idx+5]);
248 else if(energy >= (*(fParamLow[Z]))[0]) {
252 cs = x1*((*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
253 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
254 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5]);
258 else if(energy >= (*(fParamHigh[Z]))[1]) {
259 cs = x3*(fCrossSection[Z])->
Value(energy);
264 cs = x3*(fCrossSectionLE[Z])->
Value(energy);
266 if (verboseLevel > 1) {
267 G4cout <<
"G4LivermorePhotoElectricModel: E(keV)= " << energy/keV
268 <<
" Z= " << Z <<
" cross(barn)= " << cs/barn <<
G4endl;
277 std::vector<G4DynamicParticle*>* fvect,
283 if (verboseLevel > 3) {
284 G4cout <<
"G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
285 << gammaEnergy/keV <<
G4endl;
294 if(fWater && (material == fWater ||
296 if(gammaEnergy <= fWaterEnergyLimit) {
314 if(Z >= maxZ) { Z = maxZ-1; }
317 if(!fCrossSection[Z]) {
324 size_t nn = fNShellsUsed[Z];
327 if(gammaEnergy >= (*(fParamHigh[Z]))[0])
334 G4int idx = nn*7 - 5;
339 G4double cs0 = rand*( (*(fParamHigh[Z]))[idx]
340 + x1*(*(fParamHigh[Z]))[idx+1]
341 + x2*(*(fParamHigh[Z]))[idx+2]
342 + x3*(*(fParamHigh[Z]))[idx+3]
343 + x4*(*(fParamHigh[Z]))[idx+4]
344 + x5*(*(fParamHigh[Z]))[idx+5]);
346 for(shellIdx=0; shellIdx<nn; ++shellIdx)
348 idx = shellIdx*7 + 2;
349 if(gammaEnergy > (*(fParamHigh[Z]))[idx-1])
352 (*(fParamHigh[Z]))[idx]
353 + x1*(*(fParamHigh[Z]))[idx+1]
354 + x2*(*(fParamHigh[Z]))[idx+2]
355 + x3*(*(fParamHigh[Z]))[idx+3]
356 + x4*(*(fParamHigh[Z]))[idx+4]
357 + x5*(*(fParamHigh[Z]))[idx+5];
359 if(cs >= cs0) {
break; }
362 if(shellIdx >= nn) { shellIdx = nn-1; }
364 else if(gammaEnergy >= (*(fParamLow[Z]))[0])
371 G4int idx = nn*7 - 5;
375 + x1*(*(fParamLow[Z]))[idx+1]
376 + x2*(*(fParamLow[Z]))[idx+2]
377 + x3*(*(fParamLow[Z]))[idx+3]
378 + x4*(*(fParamLow[Z]))[idx+4]
379 + x5*(*(fParamLow[Z]))[idx+5]);
380 for(shellIdx=0; shellIdx<nn; ++shellIdx)
382 idx = shellIdx*7 + 2;
383 if(gammaEnergy > (*(fParamLow[Z]))[idx-1])
385 G4double cs = (*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
386 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
387 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5];
388 if(cs >= cs0) {
break; }
391 if(shellIdx >= nn) {shellIdx = nn-1;}
399 if(gammaEnergy >= (*(fParamHigh[Z]))[1]) {
401 cs*= (fCrossSection[Z])->
Value(gammaEnergy);
406 cs *= (fCrossSectionLE[Z])->
Value(gammaEnergy);
409 for(
size_t j=0; j<nn; ++j) {
412 if(gammaEnergy > (*(fParamLow[Z]))[7*shellIdx+1]) {
415 if(cs <= 0.0 || j+1 == nn) {
break;}
421 G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx*7 + 1];
430 if(fDeexcitationActive && shellIdx + 1 < nn) {
437 if(gammaEnergy < bindingEnergy) {
443 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
457 fvect->push_back(electron);
463 G4int nbefore = fvect->size();
466 G4int nafter = fvect->size();
467 if(nafter > nbefore) {
469 for (
G4int j=nbefore; j<nafter; ++j) {
471 G4double e = ((*fvect)[j])->GetKineticEnergy();
472 if(esec + e > edep) {
475 ((*fvect)[j])->SetKineticEnergy(e);
478 for (
G4int jj=nafter-1; jj>j; --jj) {
498const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
502 if(fDataDirectory.empty()) {
503 const char* path = std::getenv(
"G4LEDATA");
505 std::ostringstream ost;
506 ost << path <<
"/livermore/phot_epics2014/";
507 fDataDirectory = ost.str();
509 G4Exception(
"G4SeltzerBergerModel::FindDirectoryPath()",
"em0006",
511 "Environment variable G4LEDATA not defined");
514 return fDataDirectory;
519void G4LivermorePhotoElectricModel::ReadData(
G4int Z)
521 if (verboseLevel > 1)
523 G4cout <<
"Calling ReadData() of G4LivermorePhotoElectricModel"
527 if(fCrossSection[Z]) {
return; }
533 std::ostringstream ost;
534 ost << FindDirectoryPath() <<
"pe-cs-" << Z <<
".dat";
535 std::ifstream fin(ost.str().c_str());
536 if( !fin.is_open()) {
538 ed <<
"G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
539 <<
"> is not opened!" <<
G4endl;
540 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
542 ed,
"G4LEDATA version should be G4EMLOW7.2 or later.");
545 if(verboseLevel > 3) {
546 G4cout <<
"File " << ost.str().c_str()
547 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
549 fCrossSection[Z]->
Retrieve(fin,
true);
554 fParamHigh[Z] =
new std::vector<G4double>;
558 std::ostringstream ost1;
559 ost1 << fDataDirectory <<
"pe-high-" << Z <<
".dat";
560 std::ifstream fin1(ost1.str().c_str());
561 if( !fin1.is_open()) {
563 ed <<
"G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
564 <<
"> is not opened!" <<
G4endl;
565 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
567 ed,
"G4LEDATA version should be G4EMLOW7.2 or later.");
570 if(verboseLevel > 3) {
571 G4cout <<
"File " << ost1.str().c_str()
572 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
575 if(fin1.fail()) {
return; }
576 if(0 > n1 || n1 >=
INT_MAX) { n1 = 0; }
579 if(fin1.fail()) {
return; }
580 if(0 > n2 || n2 >=
INT_MAX) { n2 = 0; }
583 if(fin1.fail()) {
return; }
586 fParamHigh[Z]->reserve(7*n1+1);
587 fParamHigh[Z]->push_back(x*MeV);
588 for(
G4int i=0; i<n1; ++i) {
589 for(
G4int j=0; j<7; ++j) {
591 if(0 == j) { x *= MeV; }
593 fParamHigh[Z]->push_back(x);
599 fParamLow[Z] =
new std::vector<G4double>;
603 std::ostringstream ost1_low;
604 ost1_low << fDataDirectory <<
"pe-low-" << Z <<
".dat";
605 std::ifstream fin1_low(ost1_low.str().c_str());
606 if( !fin1_low.is_open()) {
608 ed <<
"G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
609 <<
"> is not opened!" <<
G4endl;
610 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
612 ed,
"G4LEDATA version should be G4EMLOW7.2 or later.");
615 if(verboseLevel > 3) {
616 G4cout <<
"File " << ost1_low.str().c_str()
617 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
620 if(fin1_low.fail()) {
return; }
621 if(0 > n1_low || n1_low >=
INT_MAX) { n1_low = 0; }
624 if(fin1_low.fail()) {
return; }
625 if(0 > n2_low || n2_low >=
INT_MAX) { n2_low = 0; }
628 if(fin1_low.fail()) {
return; }
630 fNShells[Z] = n1_low;
631 fParamLow[Z]->reserve(7*n1_low+1);
632 fParamLow[Z]->push_back(x_low*MeV);
633 for(
G4int i=0; i<n1_low; ++i) {
634 for(
G4int j=0; j<7; ++j) {
636 if(0 == j) { x_low *= MeV; }
637 else { x_low *= barn; }
638 fParamLow[Z]->push_back(x_low);
644 if(nShellLimit < n2) { n2 = nShellLimit; }
646 fNShellsUsed[Z] = n2;
649 std::ostringstream ost2;
650 ost2 << fDataDirectory <<
"pe-ss-cs-" << Z <<
".dat";
651 std::ifstream fin2(ost2.str().c_str());
652 if( !fin2.is_open()) {
654 ed <<
"G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
655 <<
"> is not opened!" <<
G4endl;
656 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
658 ed,
"G4LEDATA version should be G4EMLOW7.2 or later.");
661 if(verboseLevel > 3) {
662 G4cout <<
"File " << ost2.str().c_str()
663 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
668 for(
G4int i=0; i<n2; ++i) {
669 fin2 >> x >> y >> n3 >> n4;
671 for(
G4int j=0; j<n3; ++j) {
681 if(1 < fNShells[Z]) {
683 std::ostringstream ost3;
684 ost3 << fDataDirectory <<
"pe-le-cs-" << Z <<
".dat";
685 std::ifstream fin3(ost3.str().c_str());
686 if( !fin3.is_open()) {
688 ed <<
"G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
689 <<
"> is not opened!" <<
G4endl;
690 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
692 ed,
"G4LEDATA version should be G4EMLOW7.2 or later.");
695 if(verboseLevel > 3) {
696 G4cout <<
"File " << ost3.str().c_str()
697 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
699 fCrossSectionLE[Z]->
Retrieve(fin3,
true);
709 if(Z < 1 || Z >= maxZ) {
return -1;}
712 if(!fCrossSection[Z] || shell < 0 || shell >= fNShellsUsed[Z]) {
return -1; }
717 return fCrossSection[Z]->
Energy(0);
725 if (!fCrossSection[Z]) {
726#ifdef G4MULTITHREADED
728 if (!fCrossSection[Z]) {
731#ifdef G4MULTITHREADED
std::vector< G4Element * > G4ElementVector
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
#define G4MUTEX_INITIALIZER
#define G4MUTEXLOCK(mutex)
#define G4MUTEXUNLOCK(mutex)
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4int GetComponentID(G4int Z, size_t idx)
G4PhysicsVector * GetComponentDataByIndex(G4int Z, size_t idx)
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
G4double GetValueForComponent(G4int Z, size_t idx, G4double kinEnergy)
void PutValues(std::size_t index, G4double e, G4double dataValue)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double energy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4ParticleChangeForGamma * fParticleChange
G4double GetBindingEnergy(G4int Z, G4int shell)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual ~G4LivermorePhotoElectricModel()
G4LivermorePhotoElectricModel(const G4String &nam="LivermorePhElectric")
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double energy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
G4double GetDensity() const
const G4ElementVector * GetElementVector() const
const G4Material * GetBaseMaterial() const
size_t GetNumberOfElements() const
G4SandiaTable * GetSandiaTable() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Energy(std::size_t index) const
virtual void ScaleVector(G4double factorE, G4double factorV)
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
static G4ProductionCutsTable * GetProductionCutsTable()
void GetSandiaCofWater(G4double energy, std::vector< G4double > &coeff) const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
G4bool IsFluoActive() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4VEmAngularDistribution * GetAngularDistribution()
G4ParticleChangeForGamma * GetParticleChangeForGamma()
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void SetDeexcitationFlag(G4bool val)
void SetAngularDistribution(G4VEmAngularDistribution *)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)