68 G4cout <<
"Livermore Polarized GammaConversion is constructed "
95 G4cout <<
"Calling1 G4LivermorePolarizedGammaConversionModel::Initialise()"
113 char* path = getenv(
"G4LEDATA");
120 for(
G4int i=0; i<numOfCouples; ++i)
127 for (
G4int j=0; j<nelm; ++j)
165 G4cout <<
"Calling ReadData() of G4LivermorePolarizedGammaConversionModel"
169 if(
data[Z]) {
return; }
171 const char* datadir = path;
175 datadir = getenv(
"G4LEDATA");
178 G4Exception(
"G4LivermorePolarizedGammaConversionModel::ReadData()",
180 "Environment variable G4LEDATA not defined");
191 std::ostringstream ost;
192 ost << datadir <<
"/livermore/pair/pp-cs-" << Z <<
".dat";
193 std::ifstream
fin(ost.str().c_str());
198 ed <<
"G4LivermorePolarizedGammaConversionModel data file <" << ost.str().c_str()
199 <<
"> is not opened!" <<
G4endl;
200 G4Exception(
"G4LivermorePolarizedGammaConversionModel::ReadData()",
202 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
209 <<
" is opened by G4LivermorePolarizedGammaConversionModel" <<
G4endl;}
215 data[
Z] ->SetSpline(
true);
228 G4cout <<
"G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()"
237 if(intZ < 1 || intZ >
maxZ) {
return xs; }
247 if(!pv) {
return xs; }
250 xs = pv->
Value(GammaEnergy);
255 G4cout <<
"****** DEBUG: tcs value for Z=" << Z <<
" at energy (MeV)="
257 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
258 G4cout <<
" -> first cs value in EADL data file (iu) =" << (*pv)[0] <<
G4endl;
259 G4cout <<
" -> last cs value in EADL data file (iu) =" << (*pv)[
n] <<
G4endl;
260 G4cout <<
"*********************************************************" <<
G4endl;
282 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" <<
G4endl;
301 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1
e-6))||(gammaPolarization0.
mag()==0))
323 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
336 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - element = 0" <<
G4endl;
345 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - ionisation = 0" <<
G4endl;
352 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
360 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
362 G4double epsilonRange = 0.5 - epsilonMin ;
376 epsilon = 0.5 - epsilonRange * pow(
G4UniformRand(), 0.3333) ;
377 screen = screenFactor / (epsilon * (1. -
epsilon));
383 screen = screenFactor / (epsilon * (1 -
epsilon));
401 electronTotEnergy = (1. -
epsilon) * photonEnergy;
402 positronTotEnergy = epsilon * photonEnergy;
406 positronTotEnergy = (1. -
epsilon) * photonEnergy;
407 electronTotEnergy = epsilon * photonEnergy;
445 phi =
SetPhi(photonEnergy);
446 psi =
SetPsi(photonEnergy,phi);
521 dirX = sinTheta*cos(phip);
522 dirY = sinTheta*sin(phip);
532 positronDirection, positronKineEnergy);
535 fvect->push_back(particle1);
536 fvect->push_back(particle2);
552 if (screenVariable > 1.)
553 value = 42.24 - 8.368 * log(screenVariable + 0.952);
555 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
568 if (screenVariable > 1.)
569 value = 42.24 - 8.368 * log(screenVariable + 0.952);
571 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
583 G4double Momentum = sqrt(Energy*Energy -1);
586 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
587 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
611 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
612 const G4double aw = 0.0151, bw = 10.7, cw = -410.;
614 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
616 pl[0] =
Fln(ay0,by0,Ene);
617 pl[1] = aa0 + ba0*(Ene);
618 pl[2] =
Poli(aw,bw,cw,Ene);
619 pl[3] =
Poli(axc,bxc,cxc,Ene);
621 const G4double abf = 3.1216, bbf = 2.68;
623 pt[1] = abf + bbf/Ene;
639 const G4double aw = 0.21, bw = 10.8, cw = -58.;
640 const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
642 pl[0] =
Fln(ay0, by0, Ene);
643 pl[1] =
Fln(aa0, ba0, Ene);
644 pl[2] =
Poli(aw,bw,cw,Ene);
645 pl[3] =
Poli(axc,bxc,cxc,Ene);
717 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
718 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
719 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
720 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
721 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
722 const G4double axcp = 2.8E-3,bxcp = -3.133;
723 const G4double abf0 = 3.1213, bbf0 = 2.61;
724 const G4double abfpm = 3.1231, bbfpm = 2.84;
726 p0l[0] =
Fln(ay00, by00, Ene);
727 p0l[1] =
Fln(aa00, ba00, Ene);
728 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
729 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
731 ppml[0] =
Fln(ay0p, by0p, Ene);
732 ppml[1] = aap + bap*(Ene);
733 ppml[2] =
Poli(awp, bwp, cwp, Ene);
734 ppml[3] =
Fln(axcp,bxcp,Ene);
737 p0t[1] = abf0 + bbf0/Ene;
739 ppmt[1] = abfpm + bbfpm/Ene;
745 xe0 =
Encu(p0l, p0t, xi);
747 xepm =
Encu(ppml, ppmt, xi);
752 const G4double ay00 = 2.82, by00 = 6.35;
753 const G4double aa00 = -1.75, ba00 = 0.25;
755 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
756 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
757 const G4double ay0p = 1.56, by0p = 3.6;
758 const G4double aap = 0.86, bap = 8.3E-3;
759 const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
762 p0l[0] =
Fln(ay00, by00, Ene);
763 p0l[1] = aa00+pow(Ene, ba00);
764 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
765 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
766 ppml[0] =
Fln(ay0p, by0p, Ene);
767 ppml[1] = aap + bap*(Ene);
768 ppml[2] =
Poli(awp, bwp, cwp, Ene);
779 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
783 b =
Ftan(ppmt,PhiLocal);
787 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
791 a =
Ftan(p0t,PhiLocal);
796 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
797 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
818 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
831 value =(a + b/x + c/(x*x*
x));
865 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
866 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
868 if(x > xmax) {
return xmax; }
873 if(std::fabs(fx) <= x*1.0e-6) {
break; }
876 if(x < 0.0) { x = 0.0; }
889 value = 1./(
pi*(w*w + 4.*(x-xc)*(x-xc)));
902 value = (y0 *
pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
915 value = (-16.*A*w*(x-xc))/
916 (
pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
929 value = y0*x + A*atan( 2*(x-xc)/w) /
pi;
943 nor = atan(2.*(
pi-xc)/w)/(2.*
pi*
w) - atan(2.*(x-xc)/
w)/(2.*
pi*w);
944 value = xc - (w/2.)*tan(-2.*r*nor*
pi*w+atan(2*(xc-x)/w));
967 value = -1.*a / ((x-b)*(x-b));
990 value = b*(1-
G4Exp(r*cnor/a));
1030 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
1031 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
1032 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
1056 return gammaPolarization - gammaPolarization.
dot(gammaDirection)/gammaDirection.
dot(gammaDirection) * gammaDirection;
1077 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
1091 G4AutoLock l(&LivermorePolarizedGammaConversionModelMutex);
G4double Flor(G4double *, G4double)
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreeVector GetRandomPolarization(G4ThreeVector &direction0)
CLHEP::Hep3Vector G4ThreeVector
void ReadData(size_t Z, const char *path=0)
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
G4double Fdlor(G4double *, G4double)
static constexpr double MeV
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
const G4ThreeVector & GetMomentumDirection() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ElementVector * GetElementVector() const
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
double dot(const Hep3Vector &) const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double Finvlor(G4double *, G4double, G4double)
G4double Value(G4double theEnergy, size_t &lastidx) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
G4double Glor(G4double *, G4double)
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
G4double LowEnergyLimit() const
const XML_Char const XML_Char * data
#define G4MUTEX_INITIALIZER
G4ParticleDefinition * GetDefinition() const
G4double Finvtan(G4double *, G4double, G4double)
const XML_Char int const XML_Char * value
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
static constexpr double electron_mass_c2
virtual ~G4LivermorePolarizedGammaConversionModel()
static constexpr double twopi
double A(double temperature)
G4double Poli(G4double, G4double, G4double, G4double)
static G4Positron * Positron()
double howOrthogonal(const Hep3Vector &v) const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4LivermorePolarizedGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedGammaConversion")
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
static G4Electron * Electron()
size_t GetTableSize() const
Hep3Vector cross(const Hep3Vector &) const
G4double SetPhi(G4double)
double epsilon(double density, double temperature)
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double Finttan(G4double *, G4double)
G4ThreeVector SetPerpendicularVector(G4ThreeVector &a)
const G4ThreeVector & GetPolarization() const
std::vector< G4Element * > G4ElementVector
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
void SetTheta(G4double *, G4double *, G4double)
G4double ScreenFunction2(G4double screenVariable)
G4double Fdtan(G4double *, G4double)
G4IonisParamElm * GetIonisation() const
G4double SetPsi(G4double, G4double)
G4ParticleChangeForGamma * fParticleChange
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4Material * GetMaterial() const
static constexpr double pi
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double Encu(G4double *, G4double *, G4double)
G4double GetfCoulomb() const
G4double GetlogZ3() const
G4double HighEnergyLimit() const
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
G4double ScreenFunction1(G4double screenVariable)
G4double Ftan(G4double *, G4double)
G4double Fintlor(G4double *, G4double)
size_t GetNumberOfElements() const
size_t GetVectorLength() const
G4double Fln(G4double, G4double, G4double)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4LPhysicsFreeVector * data[100]