36 #include "G4LivermorePhotoElectricModel.hh" 42 #include "G4ParticleChangeForGamma.hh" 53 std::vector<G4double>* G4LivermorePhotoElectricModel::fParam[] = {
nullptr};
54 G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
55 G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
56 G4ElementData* G4LivermorePhotoElectricModel::fShellCrossSection =
nullptr;
57 G4Material* G4LivermorePhotoElectricModel::fWater =
nullptr;
58 G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
64 G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(
66 :
G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
67 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
68 fAtomDeexcitation(nullptr)
85 G4cout <<
"Livermore PhotoElectric is constructed " 86 <<
" nShellLimit= " << nShellLimit <<
G4endl;
91 fSandiaCof.resize(4,0.0);
97 G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
100 delete fShellCrossSection;
104 delete fCrossSection[i];
105 fCrossSection[i] = 0;
106 delete fCrossSectionLE[i];
107 fCrossSectionLE[i] = 0;
119 G4cout <<
"Calling G4LivermorePhotoElectricModel::Initialise()" <<
G4endl;
126 if(fWater) { fWaterEnergyLimit = 13.6*
eV; }
129 if(!fShellCrossSection) { fShellCrossSection =
new G4ElementData(); }
131 char* path = getenv(
"G4LEDATA");
137 for(
G4int i=0; i<numOfCouples; ++i) {
144 for (
G4int j=0; j<nelm; ++j) {
148 if(!fCrossSection[Z]) {
ReadData(Z, path); }
154 G4cout <<
"Loaded cross section files for LivermorePhotoElectric model" 163 fDeexcitationActive =
false;
164 if(fAtomDeexcitation) {
165 fDeexcitationActive = fAtomDeexcitation->
IsFluoActive();
169 G4cout <<
"LivermorePhotoElectric model is initialized " <<
G4endl 176 G4double G4LivermorePhotoElectricModel::CrossSectionPerVolume(
183 if(fWater && (material == fWater ||
185 if(energy <= fWaterEnergyLimit) {
186 fWater->GetSandiaTable()->GetSandiaCofWater(energy, fSandiaCof);
193 (fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
194 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4);
197 if(0.0 == fCurrSection) {
212 G4cout <<
"G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():" 213 <<
" Z= " << ZZ <<
" R(keV)= " << energy/
keV <<
G4endl;
217 if(Z < 1 || Z >=
maxZ) {
return cs; }
221 if(!fCrossSection[Z]) {
223 if(!fCrossSection[Z]) {
return cs; }
226 G4int idx = fNShells[
Z]*6 - 4;
227 if (energy < (*(fParam[Z]))[idx-1]) { energy = (*(fParam[
Z]))[idx-1]; }
234 if(energy >= (*(fParam[Z]))[0]) {
236 cs = x1*((*(fParam[
Z]))[idx] + x1*(*(fParam[
Z]))[idx+1]
237 + x2*(*(fParam[
Z]))[idx+2] + x3*(*(fParam[
Z]))[idx+3]
238 + x4*(*(fParam[
Z]))[idx+4]);
240 }
else if(energy >= (*(fParam[Z]))[1]) {
241 cs = x3*(fCrossSection[
Z])->
Value(energy);
245 cs = x3*(fCrossSectionLE[
Z])->
Value(energy);
248 G4cout <<
"LivermorePhotoElectricModel: E(keV)= " << energy/
keV 249 <<
" Z= " << Z <<
" cross(barn)= " << cs/
barn <<
G4endl;
257 G4LivermorePhotoElectricModel::SampleSecondaries(
258 std::vector<G4DynamicParticle*>* fvect,
266 G4cout <<
"G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= " 276 if(fWater && (material == fWater ||
278 if(gammaEnergy <= fWaterEnergyLimit) {
299 if(!fCrossSection[Z]) {
306 size_t nn = fNShellsUsed[
Z];
309 if(gammaEnergy >= (*(fParam[Z]))[0]) {
314 G4int idx = nn*6 - 4;
318 + x1*(*(fParam[
Z]))[idx+1]
319 +
x2*(*(fParam[
Z]))[idx+2]
320 + x3*(*(fParam[
Z]))[idx+3]
321 + x4*(*(fParam[
Z]))[idx+4]);
322 for(shellIdx=0; shellIdx<
nn; ++shellIdx) {
323 idx = shellIdx*6 + 2;
324 if(gammaEnergy > (*(fParam[Z]))[idx-1]) {
325 G4double cs = (*(fParam[
Z]))[idx] +
x1*(*(fParam[
Z]))[idx+1]
326 +
x2*(*(fParam[
Z]))[idx+2] + x3*(*(fParam[
Z]))[idx+3]
327 + x4*(*(fParam[
Z]))[idx+4];
328 if(cs >= cs0) {
break; }
331 if(shellIdx >= nn) { shellIdx = nn-1; }
339 if(gammaEnergy >= (*(fParam[Z]))[1]) {
340 cs *= (fCrossSection[
Z])->
Value(gammaEnergy);
342 cs *= (fCrossSectionLE[
Z])->
Value(gammaEnergy);
345 for(
size_t j=0; j<
nn; ++j) {
346 shellIdx = (size_t)fShellCrossSection->GetComponentID(Z, j);
347 if(gammaEnergy > (*(fParam[Z]))[6*shellIdx+1]) {
348 cs -= fShellCrossSection->GetValueForComponent(Z, j, gammaEnergy);
350 if(cs <= 0.0 || j+1 == nn) {
break; }
355 G4double bindingEnergy = (*(fParam[
Z]))[shellIdx*6 + 1];
364 if(fDeexcitationActive && shellIdx + 1 < nn) {
366 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
371 if(gammaEnergy < bindingEnergy) {
377 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
391 fvect->push_back(electron);
396 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
397 G4int nbefore = fvect->size();
399 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
400 G4int nafter = fvect->size();
401 if(nafter > nbefore) {
403 for (
G4int j=nbefore; j<nafter; ++j) {
405 G4double e = ((*fvect)[j])->GetKineticEnergy();
406 if(esec + e > edep) {
409 ((*fvect)[j])->SetKineticEnergy(e);
412 for (
G4int jj=nafter-1; jj>j; --jj) {
433 G4LivermorePhotoElectricModel::ReadData(
G4int Z,
const char* path)
437 G4cout <<
"Calling ReadData() of G4LivermoreGammaConversionModel" 441 if(fCrossSection[Z]) {
return; }
443 const char* datadir = path;
447 datadir = getenv(
"G4LEDATA");
450 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
452 "Environment variable G4LEDATA not defined");
459 fCrossSection[
Z]->SetSpline(
true);
461 std::ostringstream ost;
462 ost << datadir <<
"/livermore/phot/pe-cs-" << Z <<
".dat";
463 std::ifstream
fin(ost.str().c_str());
464 if( !
fin.is_open()) {
466 ed <<
"G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
467 <<
"> is not opened!" <<
G4endl;
468 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
470 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
474 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;}
475 fCrossSection[
Z]->Retrieve(
fin,
true);
476 fCrossSection[
Z]->ScaleVector(
MeV,
barn);
480 fParam[
Z] =
new std::vector<G4double>;
486 std::ostringstream ost1;
487 ost1 << datadir <<
"/livermore/phot/pe-" << Z <<
".dat";
488 std::ifstream fin1(ost1.str().c_str());
489 if( !fin1.is_open()) {
491 ed <<
"G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
492 <<
"> is not opened!" <<
G4endl;
493 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
495 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
499 G4cout <<
"File " << ost1.str().c_str()
500 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
503 if(fin1.fail()) {
return; }
504 if(0 > n1 || n1 >=
INT_MAX) { n1 = 0; }
507 if(fin1.fail()) {
return; }
508 if(0 > n2 || n2 >=
INT_MAX) { n2 = 0; }
511 if(fin1.fail()) {
return; }
514 fParam[
Z]->reserve(6*n1+1);
515 fParam[
Z]->push_back(x*
MeV);
516 for(
G4int i=0; i<n1; ++i) {
517 for(
G4int j=0; j<6; ++j) {
519 if(0 == j) { x *=
MeV; }
521 fParam[
Z]->push_back(x);
527 if(nShellLimit < n2) { n2 = nShellLimit; }
528 fShellCrossSection->InitialiseForComponent(Z, n2);
529 fNShellsUsed[
Z] = n2;
532 std::ostringstream ost2;
533 ost2 << datadir <<
"/livermore/phot/pe-ss-cs-" << Z <<
".dat";
534 std::ifstream fin2(ost2.str().c_str());
535 if( !fin2.is_open()) {
537 ed <<
"G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
538 <<
"> is not opened!" <<
G4endl;
539 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
541 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
545 G4cout <<
"File " << ost2.str().c_str()
546 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
551 for(
G4int i=0; i<n2; ++i) {
552 fin2 >> x >> y >> n3 >> n4;
554 for(
G4int j=0; j<n3; ++j) {
558 fShellCrossSection->AddComponent(Z, n4, v);
565 if(1 < fNShells[Z]) {
568 std::ostringstream ost3;
569 ost3 << datadir <<
"/livermore/phot/pe-le-cs-" << Z <<
".dat";
570 std::ifstream fin3(ost3.str().c_str());
571 if( !fin3.is_open()) {
573 ed <<
"G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
574 <<
"> is not opened!" <<
G4endl;
575 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
577 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
581 G4cout <<
"File " << ost3.str().c_str()
582 <<
" is opened by G4LivermorePhotoElectricModel" <<
G4endl;
584 fCrossSectionLE[
Z]->Retrieve(fin3,
true);
585 fCrossSectionLE[
Z]->ScaleVector(MeV,
barn);
596 void G4LivermorePhotoElectricModel::InitialiseForElement(
599 G4AutoLock l(&LivermorePhotoElectricModelMutex);
602 if(!fCrossSection[Z]) {
ReadData(Z); }
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4LossTableManager * Instance()
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
const G4Material * GetMaterial() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
G4bool IsFluoActive() const
G4VEmAngularDistribution * GetAngularDistribution()
G4double GetDensity() const
#define G4MUTEX_INITIALIZER
G4ParticleChangeForGamma * fParticleChange
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
ComputeCrossSectionPerAtom
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
const G4ThreeVector & GetMomentumDirection() const
void SetAngularDistribution(G4VEmAngularDistribution *)
size_t GetNumberOfElements() const
size_t GetTableSize() const
static G4Electron * Electron()
G4VAtomDeexcitation * AtomDeexcitation()
const G4ElementVector * GetElementVector() const
void ReadData(size_t Z, const char *path=0)
void SetDeexcitationFlag(G4bool val)
const G4Material * GetBaseMaterial() const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)