65 nShellLimit(100), fDeexcitationActive(false), isInitialised(false), fAtomDeexcitation(0)
89 G4cout <<
"Livermore Polarized PhotoElectric is constructed "
124 G4cout <<
"Calling G4LivermorePolarizedPhotoElectricGDModel::Initialise()" <<
G4endl;
138 char* path = getenv(
"G4LEDATA");
144 for(
G4int i=0; i<numOfCouples; ++i) {
151 for (
G4int j=0; j<nelm; ++j) {
152 G4int Z = (
G4int)(*theElementVector)[j]->GetZ();
161 G4cout <<
"Loaded cross section files for LivermorePhotoElectric model"
176 G4cout <<
"LivermorePolarizedPhotoElectric model is initialized " <<
G4endl
192 G4cout <<
"G4LivermorePolarizedPhotoElectricGDModel::ComputeCrossSectionPerAtom():"
193 <<
" Z= " << ZZ <<
" R(keV)= " << GammaEnergy/
keV <<
G4endl;
197 if(Z < 1 || Z >=
maxZ) {
return cs; }
207 if (GammaEnergy < (*(
fParam[Z]))[idx-1]) { GammaEnergy = (*(
fParam[Z]))[idx-1]; }
214 if(GammaEnergy >= (*(
fParam[Z]))[0]) {
218 + x4*(*(
fParam[Z]))[idx+4]);
220 }
else if (GammaEnergy >= (*(
fParam[Z]))[1]) {
228 G4cout <<
"LivermorePolarizedPhotoElectricGDModel: E(keV)= " << GammaEnergy/
keV
229 <<
" Z= " << Z <<
" cross(barn)= " << cs/
barn <<
G4endl;
238 std::vector<G4DynamicParticle*>* fvect,
245 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricGDModel"
251 G4cout <<
"G4LivermorePolarizedPhotoElectricGDModel::SampleSecondaries() Egamma(keV)= "
280 if(!(gammaPolarization0.isOrthogonal(photonDirection, 1e-6))||(gammaPolarization0.mag()==0))
286 if ( gammaPolarization0.howOrthogonal(photonDirection) != 0)
321 if(photonEnergy >= (*(
fParam[Z]))[0]) {
326 G4int idx = nn*6 - 4;
330 + x1*(*(
fParam[Z]))[idx+1]
331 + x2*(*(
fParam[Z]))[idx+2]
332 + x3*(*(
fParam[Z]))[idx+3]
333 + x4*(*(
fParam[Z]))[idx+4]);
334 for(shellIdx=0; shellIdx<
nn; ++shellIdx) {
335 idx = shellIdx*6 + 2;
336 if(photonEnergy > (*(
fParam[Z]))[idx-1]) {
339 + x4*(*(
fParam[Z]))[idx+4];
340 if(cs >= cs0) {
break; }
343 if(shellIdx >= nn) { shellIdx = nn-1; }
351 if(photonEnergy >= (*(
fParam[Z]))[1]) {
357 for(
size_t j=0; j<
nn; ++j) {
359 if(photonEnergy > (*(
fParam[Z]))[6*shellIdx+1]) {
362 if(cs <= 0.0 || j+1 == nn) {
break; }
383 if(photonEnergy < bindingEnergy) {
395 G4double sintheta = sqrt(1. - costheta*costheta);
405 fvect->push_back(electron);
413 G4int nbefore = fvect->size();
416 G4int nafter = fvect->size();
417 if(nafter > nbefore) {
419 for (
G4int j=nbefore; j<nafter; ++j) {
421 G4double e = ((*fvect)[j])->GetKineticEnergy();
422 if(esec + e > edep) {
425 ((*fvect)[j])->SetKineticEnergy(e);
428 for (
G4int jj=nafter-1; jj>j; --jj) {
453 G4cout <<
"Calling ReadData() of G4LivermorePolarizedPhotoElectricGDModel"
459 const char* datadir = path;
463 datadir = getenv(
"G4LEDATA");
466 G4Exception(
"G4LivermorePhotoElectricModel::ReadData()",
468 "Environment variable G4LEDATA not defined");
477 std::ostringstream ost;
478 ost << datadir <<
"/livermore/phot/pe-cs-" << Z <<
".dat";
479 std::ifstream fin(ost.str().c_str());
480 if( !fin.is_open()) {
482 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost.str().c_str()
483 <<
"> is not opened!" <<
G4endl;
484 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
486 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
490 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;}
496 fParam[Z] =
new std::vector<G4double>;
502 std::ostringstream ost1;
503 ost1 << datadir <<
"/livermore/phot/pe-" << Z <<
".dat";
504 std::ifstream fin1(ost1.str().c_str());
505 if( !fin1.is_open()) {
507 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost1.str().c_str()
508 <<
"> is not opened!" <<
G4endl;
509 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
511 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
515 G4cout <<
"File " << ost1.str().c_str()
516 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
519 if(fin1.fail()) {
return; }
520 if(0 > n1 || n1 >=
INT_MAX) { n1 = 0; }
523 if(fin1.fail()) {
return; }
524 if(0 > n2 || n2 >=
INT_MAX) { n2 = 0; }
527 if(fin1.fail()) {
return; }
530 fParam[Z]->reserve(6*n1+1);
532 for(
G4int i=0; i<n1; ++i) {
533 for(
G4int j=0; j<6; ++j) {
535 if(0 == j) { x *=
MeV; }
548 std::ostringstream ost2;
549 ost2 << datadir <<
"/livermore/phot/pe-ss-cs-" << Z <<
".dat";
550 std::ifstream fin2(ost2.str().c_str());
551 if( !fin2.is_open()) {
553 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost2.str().c_str()
554 <<
"> is not opened!" <<
G4endl;
555 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
557 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
561 G4cout <<
"File " << ost2.str().c_str()
562 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
567 for(
G4int i=0; i<n2; ++i) {
568 fin2 >> x >> y >> n3 >> n4;
570 for(
G4int j=0; j<n3; ++j) {
584 std::ostringstream ost3;
585 ost3 << datadir <<
"/livermore/phot/pe-le-cs-" << Z <<
".dat";
586 std::ifstream fin3(ost3.str().c_str());
587 if( !fin3.is_open()) {
589 ed <<
"G4LivermorePolarizedPhotoElectricGDModel data file <" << ost3.str().c_str()
590 <<
"> is not opened!" <<
G4endl;
591 G4Exception(
"G4LivermorePolarizedPhotoElectricGDModel::ReadData()",
593 ed,
"G4LEDATA version should be G4EMLOW6.32 or later.");
597 G4cout <<
"File " << ost3.str().c_str()
598 <<
" is opened by G4LivermorePolarizedPhotoElectricGDModel" <<
G4endl;
612 G4double rand1,rand2,onemcost,greject;
615 G4double gamma = 1. + energyE/masarep;
618 G4double beta = sqrt((gamma2 - 1.)/gamma2);
622 G4double g1 = 0.5*beta*gamma*(gamma-1.)*(gamma-2.);
626 G4double grejectmax = 2.*(g1+1./alfa);
631 onemcost = 2.*alfa*(2.*rand1 + alfap2 * sqrt(rand1))/
632 (alfap2*alfap2 - 4.*rand1);
633 greject = (2. - onemcost)*(g1+1./(alfa+onemcost));
636 while (rand2*grejectmax > greject);
647 G4double k = Ph_energy/electron_mass_c2;
650 G4double beta = sqrt((gamma2 - 1.)/gamma2);
652 G4double d = (2./(k*gamma*(1-beta*costheta))-1)*(1/k);
665 phiprob = 1 +2*d*cos(phi)*cos(phi);
667 while (rnd2*norm_factor > phiprob);
704 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
705 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
706 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
731 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
749 G4double direction_x = direction1.getX();
750 G4double direction_y = direction1.getY();
751 G4double direction_z = direction1.getZ();
753 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
765 G4AutoLock l(&LivermorePolarizedPhotoElectricGDModelMutex);
G4ParticleChangeForGamma * fParticleChange
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
G4VAtomDeexcitation * fAtomDeexcitation
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
static G4LPhysicsFreeVector * fCrossSectionLE[99]
G4bool IsFluoActive() const
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
static G4int fNShellsUsed[99]
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 fWaterEnergyLimit
G4ParticleDefinition * theElectron
const G4ElementVector * GetElementVector() const
#define G4MUTEX_INITIALIZER
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
std::vector< G4double > fSandiaCof
static G4Material * fWater
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4bool fDeexcitationActive
G4int GetComponentID(G4int Z, size_t idx)
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition * theGamma
size_t GetTableSize() const
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
static G4LPhysicsFreeVector * fCrossSection[99]
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)
G4ThreeVector GetRandomPolarization(G4ThreeVector &direction0)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void ReadData(G4int Z, const char *path=0)
static G4ElementData * fShellCrossSection
static G4int fNShells[99]
static std::vector< G4double > * fParam[99]
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
const G4double x[NPOINTSGL]
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)
size_t GetNumberOfElements() const
G4VAtomDeexcitation * AtomDeexcitation()
G4ThreeVector SetPerpendicularVector(G4ThreeVector &a)
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")
G4double bindingEnergy(G4int A, G4int Z)
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4double SetCosTheta(G4double)
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
G4double SetPhi(G4double, G4double, G4double)