53G4ElementData* G4LivermorePhotoElectricModel::fCrossSection =
nullptr;
54G4ElementData* G4LivermorePhotoElectricModel::fCrossSectionLE =
nullptr;
55std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {
nullptr};
56std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {
nullptr};
57G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
58G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
59G4Material* G4LivermorePhotoElectricModel::fWater =
nullptr;
60G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
61G4String G4LivermorePhotoElectricModel::fDataDirectory =
"";
63static std::once_flag applyOnce;
88 if (verboseLevel > 0) {
89 G4cout <<
"Livermore PhotoElectric is constructed "
90 <<
" nShellLimit= " << nShellLimit <<
G4endl;
97 fSandiaCof.resize(4, 0.0);
105 for (
G4int i = 0; i < ZMAXPE; ++i) {
107 delete fParamHigh[i];
108 fParamHigh[i] =
nullptr;
112 fParamLow[i] =
nullptr;
123 if (verboseLevel > 1) {
124 G4cout <<
"Calling G4LivermorePhotoElectricModel::Initialise() " <<
G4endl;
128 std::call_once(applyOnce, [
this]() { isInitializer =
true; });
133 if (fWater ==
nullptr) {
135 if (fWater ==
nullptr) {
138 if (fWater !=
nullptr) {
139 fWaterEnergyLimit = 13.6 * CLHEP::eV;
143 if (fCrossSection ==
nullptr) {
145 fCrossSection->
SetName(
"PhotoEffXS");
147 fCrossSectionLE->
SetName(
"PhotoEffLowXS");
151 std::size_t numElems = (*elemTable).size();
152 for (std::size_t ie = 0; ie < numElems; ++ie) {
164 if (verboseLevel > 1) {
165 G4cout <<
"Loaded cross section files for new LivermorePhotoElectric model" <<
G4endl;
172 fDeexcitationActive =
false;
173 if (
nullptr != fAtomDeexcitation) {
174 fDeexcitationActive = fAtomDeexcitation->
IsFluoActive();
177 if (verboseLevel > 1) {
178 G4cout <<
"LivermorePhotoElectric model is initialized " <<
G4endl;
191 if (fWater && (material == fWater || material->
GetBaseMaterial() == fWater)) {
192 if (energy <= fWaterEnergyLimit) {
196 G4double energy3 = energy * energy2;
197 G4double energy4 = energy2 * energy2;
200 * (fSandiaCof[0] / energy + fSandiaCof[1] / energy2 +
201 fSandiaCof[2] / energy3 + fSandiaCof[3] / energy4);
204 if (0.0 == fCurrSection) {
218 if (verboseLevel > 3) {
219 G4cout <<
"\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
220 <<
" Z= " << ZZ <<
" R(keV)= " << energy / keV <<
G4endl;
224 if (Z >= ZMAXPE || Z <= 0) {
236 G4int idx = fNShells[Z] * 7 - 5;
238 energy = std::max(energy, (*(fParamHigh[Z]))[idx - 1]);
245 if (energy >= (*(fParamHigh[Z]))[0]) {
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]);
255 else if (energy >= (*(fParamLow[Z]))[0]) {
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]);
264 else if (energy >= (*(fParamHigh[Z]))[1]) {
271 if (verboseLevel > 1) {
272 G4cout <<
"G4LivermorePhotoElectricModel: E(keV)= " << energy / keV <<
" Z= " << Z
273 <<
" cross(barn)= " << cs / barn <<
G4endl;
286 if (verboseLevel > 3) {
287 G4cout <<
"G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
288 << gammaEnergy / keV <<
G4endl;
297 if (fWater && (material == fWater || material->
GetBaseMaterial() == fWater)) {
298 if (gammaEnergy <= fWaterEnergyLimit) {
321 std::size_t shellIdx = 0;
322 std::size_t nn = fNShellsUsed[Z];
324 if (gammaEnergy >= (*(fParamHigh[Z]))[0]) {
330 std::size_t idx = nn * 7 - 5;
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]);
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];
352 if (shellIdx >= nn) {
356 else if (gammaEnergy >= (*(fParamLow[Z]))[0]) {
362 std::size_t idx = nn * 7 - 5;
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];
380 if (shellIdx >= nn) {
389 if (gammaEnergy >= (*(fParamHigh[Z]))[1]) {
400 if (gammaEnergy > (*(fParamLow[Z]))[7 * shellIdx + 1]) {
403 if (cs <= 0.0 || j + 1 == (
G4int)nn) {
411 G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx * 7 + 1];
415 if (fDeexcitationActive && shellIdx + 1 < nn) {
422 if (gammaEnergy < bindingEnergy) {
428 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
438 fvect->push_back(electron);
441 if (
nullptr != shell) {
444 std::size_t nbefore = fvect->size();
447 std::size_t nafter = fvect->size();
448 if (nafter > nbefore) {
450 for (std::size_t j = nbefore; j < nafter; ++j) {
451 G4double e = ((*fvect)[j])->GetKineticEnergy();
452 if (esec + e > edep) {
455 ((*fvect)[j])->SetKineticEnergy(e);
458 for (std::size_t jj = nafter - 1; jj > j; --jj) {
478const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
481 if (fDataDirectory.empty()) {
483 std::ostringstream ost;
484 if (param->LivermoreDataDir() ==
"livermore") {
485 ost << param->GetDirLEDATA() <<
"/livermore/phot_epics2014/";
488 ost << param->GetDirLEDATA() <<
"/epics2017/phot/";
490 fDataDirectory = ost.str();
492 return fDataDirectory;
497void G4LivermorePhotoElectricModel::ReadData(
G4int Z)
499 if (Z <= 0 || Z >= ZMAXPE) {
500 G4cout <<
"G4LivermorePhotoElectricModel::ReadData() Warning: Z="
501 << Z <<
" is out of range - request ignored" <<
G4endl;
504 if (verboseLevel > 1) {
505 G4cout <<
"G4LivermorePhotoElectricModel::ReadData() for Z=" << Z <<
G4endl;
520 std::ostringstream ost;
521 ost << FindDirectoryPath() <<
"pe-cs-" << Z <<
".dat";
522 std::ifstream fin(ost.str().c_str());
523 if (!fin.is_open()) {
525 ed <<
"G4LivermorePhotoElectricModel data file <"
526 << ost.str().c_str() <<
"> is not opened!" <<
G4endl;
528 "G4LEDATA version should be G4EMLOW8.0 or later.");
531 if (verboseLevel > 3) {
532 G4cout <<
"File " << ost.str().c_str() <<
" is opened by G4LivermorePhotoElectricModel"
535 pv->Retrieve(fin,
true);
536 pv->ScaleVector(MeV, barn);
537 pv->FillSecondDerivatives();
538 pv->EnableLogBinSearch(number);
543 fParamHigh[Z] =
new std::vector<G4double>;
547 std::ostringstream ost1;
548 ost1 << fDataDirectory <<
"pe-high-" << Z <<
".dat";
549 std::ifstream fin1(ost1.str().c_str());
550 if (!fin1.is_open()) {
552 ed <<
"G4LivermorePhotoElectricModel data file <"
553 << ost1.str().c_str() <<
"> is not opened!" <<
G4endl;
555 "G4LEDATA version should be G4EMLOW7.2 or later.");
558 if (verboseLevel > 3) {
559 G4cout <<
"File " << ost1.str().c_str() <<
" is opened by G4LivermorePhotoElectricModel"
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) {
595 fParamHigh[Z]->push_back(x);
601 fParamLow[Z] =
new std::vector<G4double>;
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()) {
610 ed <<
"G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
611 <<
"> is not opened!" <<
G4endl;
613 "G4LEDATA version should be G4EMLOW8.0 or later.");
616 if (verboseLevel > 3) {
617 G4cout <<
"File " << ost1_low.str().c_str() <<
" is opened by G4LivermorePhotoElectricModel"
621 if (fin1_low.fail()) {
624 if (0 > n1_low || n1_low >=
INT_MAX) {
629 if (fin1_low.fail()) {
632 if (0 > n2_low || n2_low >=
INT_MAX) {
637 if (fin1_low.fail()) {
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) {
653 fParamLow[Z]->push_back(x_low);
659 if (nShellLimit < n2) {
663 fNShellsUsed[Z] = 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()) {
671 ed <<
"G4LivermorePhotoElectricModel data file <" << ost2.str().c_str() <<
"> is not opened!"
674 "G4LEDATA version should be G4EMLOW7.2 or later.");
677 if (verboseLevel > 3) {
678 G4cout <<
"File " << ost2.str().c_str() <<
" is opened by G4LivermorePhotoElectricModel"
684 for (
G4int i = 0; i < n2; ++i) {
685 fin2 >> x >> y >> n3 >> n4;
687 for (
G4int j = 0; j < n3; ++j) {
689 v->PutValues(j, x * CLHEP::MeV, y * CLHEP::barn);
691 v->EnableLogBinSearch(number);
698 if (1 < fNShells[Z]) {
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()) {
705 ed <<
"G4LivermorePhotoElectricModel data file <" << ost3.str().c_str() <<
"> is not opened!"
708 "G4LEDATA version should be G4EMLOW8.0 or later.");
711 if (verboseLevel > 3) {
712 G4cout <<
"File " << ost3.str().c_str() <<
" is opened by G4LivermorePhotoElectricModel"
715 pv1->Retrieve(fin3,
true);
716 pv1->ScaleVector(CLHEP::MeV, CLHEP::barn);
717 pv1->EnableLogBinSearch(number);
727 if (Z < 1 || Z >= ZMAXPE) {
733 shell < 0 || shell >= fNShellsUsed[Z]) {
750 if (fCrossSection ==
nullptr) {
752 fCrossSection->
SetName(
"PhotoEffXS");
754 fCrossSectionLE->
SetName(
"PhotoEffLowXS");
761void G4LivermorePhotoElectricModel::InitialiseOnFly(
G4int Z)
763 if (fCrossSection->
GetElementData(Z) ==
nullptr && Z > 0 && Z < ZMAXPE) {
std::vector< G4Element * > G4ElementTable
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
#define G4MUTEX_INITIALIZER
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
G4PhysicsVector * GetElementData(G4int Z) const
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4PhysicsVector * GetComponentDataByIndex(G4int Z, std::size_t idx) const
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
G4int GetComponentID(G4int Z, std::size_t idx) const
G4double GetValueForComponent(G4int Z, std::size_t idx, G4double kinEnergy) const
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
void SetName(const G4String &nam)
static G4ElementTable * GetElementTable()
static G4EmParameters * Instance()
const G4String & LivermoreDataDir()
G4int NumberForFreeVector() const
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double energy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX) override
G4ParticleChangeForGamma * fParticleChange
G4double GetBindingEnergy(G4int Z, G4int shell)
void InitialiseForElement(const G4ParticleDefinition *, G4int Z) override
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
~G4LivermorePhotoElectricModel() override
G4LivermorePhotoElectricModel(const G4String &nam="LivermorePhElectric")
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double energy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) override
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
G4double GetDensity() const
const G4Material * GetBaseMaterial() const
G4SandiaTable * GetSandiaTable() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double Energy(const std::size_t index) const
G4double Value(const G4double energy, std::size_t &lastidx) const
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)
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)