46 std::vector<G4double>* G4LivermorePolarizedPhotoElectricGDModel::fParam[] = {0};
47 G4int G4LivermorePolarizedPhotoElectricGDModel::fNShells[] = {0};
48 G4int G4LivermorePolarizedPhotoElectricGDModel::fNShellsUsed[] = {0};
49 G4ElementData* G4LivermorePolarizedPhotoElectricGDModel::fShellCrossSection =
nullptr;
50 G4Material* G4LivermorePolarizedPhotoElectricGDModel::fWater =
nullptr;
51 G4double G4LivermorePolarizedPhotoElectricGDModel::fWaterEnergyLimit = 0.0;
62 :
G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
63 nShellLimit(100), fDeexcitationActive(false), isInitialised(false),
64 fAtomDeexcitation(nullptr)
78 fSandiaCof.resize(4,0.0);
81 if (verboseLevel > 0) {
82 G4cout <<
"Livermore Polarized PhotoElectric is constructed "
93 delete fShellCrossSection;
94 for(
G4int i=0; i<maxZ; ++i) {
97 delete fCrossSection[i];
99 delete fCrossSectionLE[i];
100 fCrossSectionLE[i] = 0;
112 if (verboseLevel > 2) {
113 G4cout <<
"Calling G4LivermorePolarizedPhotoElectricGDModel::Initialise()" <<
G4endl;
120 if(fWater) { fWaterEnergyLimit = 13.6*
eV; }
123 if(!fShellCrossSection) { fShellCrossSection =
new G4ElementData(); }
125 char* path = getenv(
"G4LEDATA");
131 for(
G4int i=0; i<numOfCouples; ++i) {
138 for (
G4int j=0; j<nelm; ++j) {
141 else if(Z > maxZ) { Z = maxZ; }
142 if(!fCrossSection[Z]) { ReadData(Z, path); }
147 if (verboseLevel > 2) {
148 G4cout <<
"Loaded cross section files for LivermorePhotoElectric model"
152 isInitialised =
true;
157 fDeexcitationActive =
false;
158 if(fAtomDeexcitation) {
159 fDeexcitationActive = fAtomDeexcitation->
IsFluoActive();
162 if (verboseLevel > 0) {
163 G4cout <<
"LivermorePolarizedPhotoElectric model is initialized " <<
G4endl
176 if (verboseLevel > 3) {
177 G4cout <<
"G4LivermorePolarizedPhotoElectricGDModel::ComputeCrossSectionPerAtom():"
178 <<
" Z= " << ZZ <<
" R(keV)= " << GammaEnergy/
keV <<
G4endl;
182 if(Z < 1 || Z >= maxZ) {
return cs; }
186 if(!fCrossSection[Z]) {
188 if(!fCrossSection[Z]) {
return cs; }
191 G4int idx = fNShells[
Z]*6 - 4;
192 if (GammaEnergy < (*(fParam[Z]))[idx-1]) { GammaEnergy = (*(fParam[
Z]))[idx-1]; }
199 if(GammaEnergy >= (*(fParam[Z]))[0]) {
201 cs = x1*((*(fParam[
Z]))[idx] + x1*(*(fParam[
Z]))[idx+1]
202 + x2*(*(fParam[
Z]))[idx+2] + x3*(*(fParam[
Z]))[idx+3]
203 + x4*(*(fParam[
Z]))[idx+4]);
205 }
else if (GammaEnergy >= (*(fParam[Z]))[1]) {
206 cs = x3*(fCrossSection[
Z])->
Value(GammaEnergy);
210 cs = x3*(fCrossSectionLE[
Z])->
Value(GammaEnergy);
212 if (verboseLevel > 1) {
213 G4cout <<
"LivermorePolarizedPhotoElectricGDModel: E(keV)= " << GammaEnergy/
keV
214 <<
" Z= " << Z <<
" cross(barn)= " << cs/
barn <<
G4endl;
223 std::vector<G4DynamicParticle*>* fvect,
229 if (verboseLevel > 3) {
230 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricGDModel"
235 if (verboseLevel > 3) {
236 G4cout <<
"G4LivermorePolarizedPhotoElectricGDModel::SampleSecondaries() Egamma(keV)= "
250 if(fWater && (material == fWater ||
252 if(photonEnergy <= fWaterEnergyLimit) {
265 if(!(gammaPolarization0.
isOrthogonal(photonDirection, 1e-6))||(gammaPolarization0.
mag()==0))
267 gammaPolarization0 = GetRandomPolarization(photonDirection);
273 gammaPolarization0 = GetPerpendicularPolarization(photonDirection, gammaPolarization0);
293 if(Z >= maxZ) { Z = maxZ-1; }
296 if(!fCrossSection[Z]) {
303 size_t nn = fNShellsUsed[
Z];
306 if(photonEnergy >= (*(fParam[Z]))[0]) {
311 G4int idx = nn*6 - 4;
315 + x1*(*(fParam[
Z]))[idx+1]
316 + x2*(*(fParam[
Z]))[idx+2]
317 + x3*(*(fParam[
Z]))[idx+3]
318 + x4*(*(fParam[
Z]))[idx+4]);
319 for(shellIdx=0; shellIdx<
nn; ++shellIdx) {
320 idx = shellIdx*6 + 2;
321 if(photonEnergy > (*(fParam[Z]))[idx-1]) {
322 G4double cs = (*(fParam[
Z]))[idx] + x1*(*(fParam[
Z]))[idx+1]
323 + x2*(*(fParam[
Z]))[idx+2] + x3*(*(fParam[
Z]))[idx+3]
324 + x4*(*(fParam[
Z]))[idx+4];
325 if(cs >= cs0) {
break; }
328 if(shellIdx >= nn) { shellIdx = nn-1; }
336 if(photonEnergy >= (*(fParam[Z]))[1]) {
337 cs *= (fCrossSection[
Z])->
Value(photonEnergy);
339 cs *= (fCrossSectionLE[
Z])->
Value(photonEnergy);
342 for(
size_t j=0; j<
nn; ++j) {
344 if(photonEnergy > (*(fParam[Z]))[6*shellIdx+1]) {
347 if(cs <= 0.0 || j+1 == nn) {
break; }
361 if(fDeexcitationActive && shellIdx + 1 < nn) {
368 if(photonEnergy < bindingEnergy) {
379 G4double costheta = SetCosTheta(eKineticEnergy);
380 G4double sintheta = sqrt(1. - costheta*costheta);
381 G4double phi = SetPhi(photonEnergy,eKineticEnergy,costheta);
386 SystemOfRefChange(photonDirection, electronDirection, gammaPolarization0);
390 fvect->push_back(electron);
397 G4int nbefore = fvect->size();
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) {
432 void G4LivermorePolarizedPhotoElectricGDModel::ReadData(
G4int Z,
const char* path)
434 if (verboseLevel > 1)
436 G4cout <<
"Calling ReadData() of G4LivermorePolarizedPhotoElectricGDModel"
440 if(fCrossSection[Z]) {
return; }
442 const char* datadir = path;
446 datadir = getenv(
"G4LEDATA");
449 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
451 "Environment variable G4LEDATA not defined");
460 std::ostringstream ost;
461 ost << datadir <<
"/livermore/phot/pe-cs-" << Z <<
".dat";
462 std::ifstream fin(ost.str().c_str());
463 if( !fin.is_open()) {
465 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost.str().c_str()
466 <<
"> is not opened!" <<
G4endl;
467 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
469 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
472 if(verboseLevel > 3) {
G4cout <<
"File " << ost.str().c_str()
473 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;}
479 fParam[
Z] =
new std::vector<G4double>;
485 std::ostringstream ost1;
486 ost1 << datadir <<
"/livermore/phot/pe-" << Z <<
".dat";
487 std::ifstream fin1(ost1.str().c_str());
488 if( !fin1.is_open()) {
490 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost1.str().c_str()
491 <<
"> is not opened!" <<
G4endl;
492 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
494 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
497 if(verboseLevel > 3) {
498 G4cout <<
"File " << ost1.str().c_str()
499 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
502 if(fin1.fail()) {
return; }
503 if(0 > n1 || n1 >=
INT_MAX) { n1 = 0; }
506 if(fin1.fail()) {
return; }
507 if(0 > n2 || n2 >=
INT_MAX) { n2 = 0; }
510 if(fin1.fail()) {
return; }
513 fParam[
Z]->reserve(6*n1+1);
514 fParam[
Z]->push_back(x*
MeV);
515 for(
G4int i=0; i<n1; ++i) {
516 for(
G4int j=0; j<6; ++j) {
518 if(0 == j) { x *=
MeV; }
520 fParam[
Z]->push_back(x);
526 if(nShellLimit < n2) { n2 = nShellLimit; }
528 fNShellsUsed[
Z] = n2;
531 std::ostringstream ost2;
532 ost2 << datadir <<
"/livermore/phot/pe-ss-cs-" << Z <<
".dat";
533 std::ifstream fin2(ost2.str().c_str());
534 if( !fin2.is_open()) {
536 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost2.str().c_str()
537 <<
"> is not opened!" <<
G4endl;
538 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
540 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
543 if(verboseLevel > 3) {
544 G4cout <<
"File " << ost2.str().c_str()
545 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
550 for(
G4int i=0; i<n2; ++i) {
551 fin2 >> x >> y >> n3 >> n4;
553 for(
G4int j=0; j<n3; ++j) {
564 if(1 < fNShells[Z]) {
567 std::ostringstream ost3;
568 ost3 << datadir <<
"/livermore/phot/pe-le-cs-" << Z <<
".dat";
569 std::ifstream fin3(ost3.str().c_str());
570 if( !fin3.is_open()) {
572 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost3.str().c_str()
573 <<
"> is not opened!" <<
G4endl;
574 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
576 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
579 if(verboseLevel > 3) {
580 G4cout <<
"File " << ost3.str().c_str()
581 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
583 fCrossSectionLE[
Z]->
Retrieve(fin3,
true);
593 G4double G4LivermorePolarizedPhotoElectricGDModel::SetCosTheta(
G4double energyE)
595 G4double rand1,rand2,onemcost,greject;
598 G4double gamma = 1. + energyE/masarep;
601 G4double beta = sqrt((gamma2 - 1.)/gamma2);
605 G4double g1 = 0.5*beta*gamma*(gamma-1.)*(gamma-2.);
609 G4double grejectmax = 2.*(g1+1./alfa);
614 onemcost = 2.*alfa*(2.*rand1 + alfap2 * sqrt(rand1))/
615 (alfap2*alfap2 - 4.*rand1);
616 greject = (2. - onemcost)*(g1+1./(alfa+onemcost));
619 while (rand2*grejectmax > greject);
625 G4double G4LivermorePolarizedPhotoElectricGDModel::SetPhi(
G4double Ph_energy,
633 G4double beta = sqrt((gamma2 - 1.)/gamma2);
635 G4double d = (2./(k*gamma*(1-beta*costheta))-1)*(1/k);
648 phiprob = 1 +2*d*cos(phi)*cos(phi);
650 while (rnd2*norm_factor > phiprob);
666 return x < z ?
G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
668 return y < z ?
G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
687 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
688 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
689 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
699 G4ThreeVector G4LivermorePolarizedPhotoElectricGDModel::GetPerpendicularPolarization
714 return gammaPolarization - gammaPolarization.
dot(gammaDirection)/gammaDirection.
dot(gammaDirection) * gammaDirection;
721 void G4LivermorePolarizedPhotoElectricGDModel::SystemOfRefChange
736 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
748 G4AutoLock l(&LivermorePolarizedPhotoElectricGDModelMutex);
751 if(!fCrossSection[Z]) { ReadData(Z); }
G4ParticleChangeForGamma * fParticleChange
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4bool IsFluoActive() const
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
double dot(const Hep3Vector &) const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
static G4double angle[DIM]
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
const G4ElementVector * GetElementVector() const
#define G4MUTEX_INITIALIZER
double howOrthogonal(const Hep3Vector &v) const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
static constexpr double electron_mass_c2
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4int GetComponentID(G4int Z, size_t idx)
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
const G4ThreeVector & GetMomentumDirection() const
static constexpr double eV
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual ~G4LivermorePolarizedPhotoElectricGDModel()
virtual void ScaleVector(G4double factorE, G4double factorV)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4ProductionCutsTable * GetProductionCutsTable()
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4ThreeVector & GetPolarization() const
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
G4double GetValueForComponent(G4int Z, size_t idx, G4double kinEnergy)
const G4Material * GetBaseMaterial() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void PutValues(size_t index, G4double e, G4double dataValue)
static constexpr double MeV
size_t GetNumberOfElements() const
Hep3Vector cross(const Hep3Vector &) const
G4VAtomDeexcitation * AtomDeexcitation()
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4LivermorePolarizedPhotoElectricGDModel(const G4String &nam="LivermorePolarizedPhotoElectric")
static constexpr double barn
G4double bindingEnergy(G4int A, G4int Z)
void InitialiseForComponent(G4int Z, G4int nComponents=0)
static constexpr double keV
double epsilon(double density, double temperature)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
const G4Material * GetMaterial() const