45 isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
65 G4cout <<
"Livermore Polarized GammaConversion is constructed " <<
G4endl
91 G4cout <<
"Calling G4LivermorePolarizedGammaConversionModel::Initialise()" <<
G4endl;
111 G4cout <<
"G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " <<
121 G4String crossSectionFile =
"pair/pp-cs-";
126 G4cout <<
"Loaded cross section files for Livermore Polarized GammaConversion model"
132 G4cout <<
"Livermore Polarized GammaConversion model is initialized " << G4endl
155 G4cout <<
"G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()"
179 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" <<
G4endl;
198 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
204 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
214 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
220 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
253 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
257 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
261 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
263 G4double epsilonRange = 0.5 - epsilonMin ;
277 epsilon = 0.5 - epsilonRange * pow(
G4UniformRand(), 0.3333) ;
278 screen = screenFactor / (epsilon * (1. - epsilon));
284 screen = screenFactor / (epsilon * (1 - epsilon));
301 electronTotEnergy = (1. - epsilon) * photonEnergy;
302 positronTotEnergy = epsilon * photonEnergy;
306 positronTotEnergy = (1. - epsilon) * photonEnergy;
307 electronTotEnergy = epsilon * photonEnergy;
329 G4double Ene = electronTotEnergy/electron_mass_c2;
345 phi =
SetPhi(photonEnergy);
346 psi =
SetPsi(photonEnergy,phi);
389 G4double electronKineEnergy =
std::max(0.,electronTotEnergy - electron_mass_c2) ;
400 Ene = positronTotEnergy/electron_mass_c2;
409 dirX = sinTheta*cos(phip);
410 dirY = sinTheta*sin(phip);
414 G4double positronKineEnergy =
std::max(0.,positronTotEnergy - electron_mass_c2) ;
420 positronDirection, positronKineEnergy);
423 fvect->push_back(particle1);
424 fvect->push_back(particle2);
451 if (screenVariable > 1.)
452 value = 42.24 - 8.368 * log(screenVariable + 0.952);
454 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
467 if (screenVariable > 1.)
468 value = 42.24 - 8.368 * log(screenVariable + 0.952);
470 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
482 G4double Momentum = sqrt(Energy*Energy -1);
485 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
486 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
510 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
511 const G4double aw = 0.0151, bw = 10.7, cw = -410.;
513 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
515 pl[0] =
Fln(ay0,by0,Ene);
516 pl[1] = aa0 + ba0*(Ene);
517 pl[2] =
Poli(aw,bw,cw,Ene);
518 pl[3] =
Poli(axc,bxc,cxc,Ene);
520 const G4double abf = 3.1216, bbf = 2.68;
522 pt[1] = abf + bbf/Ene;
538 const G4double aw = 0.21, bw = 10.8, cw = -58.;
539 const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
541 pl[0] =
Fln(ay0, by0, Ene);
542 pl[1] =
Fln(aa0, ba0, Ene);
543 pl[2] =
Poli(aw,bw,cw,Ene);
544 pl[3] =
Poli(axc,bxc,cxc,Ene);
616 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
617 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
618 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
619 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
620 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
621 const G4double axcp = 2.8E-3,bxcp = -3.133;
622 const G4double abf0 = 3.1213, bbf0 = 2.61;
623 const G4double abfpm = 3.1231, bbfpm = 2.84;
625 p0l[0] =
Fln(ay00, by00, Ene);
626 p0l[1] =
Fln(aa00, ba00, Ene);
627 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
628 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
630 ppml[0] =
Fln(ay0p, by0p, Ene);
631 ppml[1] = aap + bap*(Ene);
632 ppml[2] =
Poli(awp, bwp, cwp, Ene);
633 ppml[3] =
Fln(axcp,bxcp,Ene);
636 p0t[1] = abf0 + bbf0/Ene;
638 ppmt[1] = abfpm + bbfpm/Ene;
644 xe0 =
Encu(p0l, p0t, xi);
646 xepm =
Encu(ppml, ppmt, xi);
651 const G4double ay00 = 2.82, by00 = 6.35;
652 const G4double aa00 = -1.75, ba00 = 0.25;
654 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
655 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
656 const G4double ay0p = 1.56, by0p = 3.6;
657 const G4double aap = 0.86, bap = 8.3E-3;
658 const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
661 p0l[0] =
Fln(ay00, by00, Ene);
662 p0l[1] = aa00+pow(Ene, ba00);
663 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
664 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
665 ppml[0] =
Fln(ay0p, by0p, Ene);
666 ppml[1] = aap + bap*(Ene);
667 ppml[2] =
Poli(awp, bwp, cwp, Ene);
678 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
682 b =
Ftan(ppmt,PhiLocal);
686 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
690 a =
Ftan(p0t,PhiLocal);
695 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
696 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
716 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
729 value =(a + b/x + c/(x*x*x));
763 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
764 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
766 if(x > xmax) {
return xmax; }
771 if(std::fabs(fx) <= x*1.0e-6) {
break; }
774 if(x < 0.0) { x = 0.0; }
787 value = 1./(
pi*(w*w + 4.*(x-xc)*(x-xc)));
800 value = (y0 *
pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
813 value = (-16.*A*w*(x-xc))/
814 (
pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
827 value = y0*x + A*atan( 2*(x-xc)/w) /
pi;
841 nor = atan(2.*(
pi-xc)/w)/(2.*
pi*w) - atan(2.*(x-xc)/w)/(2.*
pi*w);
842 value = xc - (w/2.)*tan(-2.*r*nor*
pi*w+atan(2*(xc-x)/w));
865 value = -1.*a / ((x-b)*(x-b));
888 value = b*(1-exp(r*cnor/a));
928 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
929 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
930 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
954 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
971 G4double direction_x = direction1.getX();
972 G4double direction_y = direction1.getY();
973 G4double direction_z = direction1.getZ();
975 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
G4double Poli(G4double, G4double, G4double, G4double)
G4double LowEnergyLimit() const
virtual ~G4LivermorePolarizedGammaConversionModel()
G4LivermorePolarizedGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedGammaConversion")
G4double Finvtan(G4double *, G4double, G4double)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double HighEnergyLimit() const
void SetTheta(G4double *, G4double *, G4double)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double GetfCoulomb() const
G4ParticleDefinition * GetDefinition() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double Flor(G4double *, G4double)
G4double SetPsi(G4double, G4double)
void SetHighEnergyLimit(G4double)
G4double Finvlor(G4double *, G4double, G4double)
G4double Ftan(G4double *, G4double)
G4double FindValue(G4int Z, G4double e) const
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
G4GLOB_DLL std::ostream G4cout
G4double Fdlor(G4double *, G4double)
const G4ThreeVector & GetMomentumDirection() const
G4double Fintlor(G4double *, G4double)
static const G4double A[nN]
G4double GetlogZ3() const
G4double Fln(G4double, G4double, G4double)
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
void Initialise(G4VDataSetAlgorithm *interpolation=0, G4double minE=250 *CLHEP::eV, G4double maxE=100 *CLHEP::GeV, G4int numberOfBins=200, G4double unitE=CLHEP::MeV, G4double unitData=CLHEP::barn, G4int minZ=1, G4int maxZ=99)
static G4Positron * Positron()
G4double ScreenFunction2(G4double screenVariable)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double ScreenFunction1(G4double screenVariable)
G4IonisParamElm * GetIonisation() const
G4double Glor(G4double *, G4double)
G4double SetPhi(G4double)
const G4ThreeVector & GetPolarization() const
void LoadData(const G4String &dataFile)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4ThreeVector SetPerpendicularVector(G4ThreeVector &a)
G4ParticleChangeForGamma * fParticleChange
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
G4VCrossSectionHandler * crossSectionHandler
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4ThreeVector GetRandomPolarization(G4ThreeVector &direction0)
G4double Encu(G4double *, G4double *, G4double)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4double Finttan(G4double *, G4double)
G4double Fdtan(G4double *, G4double)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()