72 G4cout <<
"Livermore Polarized GammaConversion is constructed "
101 G4cout <<
"Calling1 G4LivermorePolarizedGammaConversionModel::Initialise()"
119 char* path = getenv(
"G4LEDATA");
126 for(
G4int i=0; i<numOfCouples; ++i)
133 for (
G4int j=0; j<nelm; ++j)
135 G4int Z = (
G4int)(*theElementVector)[j]->GetZ();
172 G4cout <<
"Calling ReadData() of G4LivermorePolarizedGammaConversionModel"
176 if(
data[Z]) {
return; }
178 const char* datadir = path;
182 datadir = getenv(
"G4LEDATA");
185 G4Exception(
"G4LivermorePolarizedGammaConversionModel::ReadData()",
187 "Environment variable G4LEDATA not defined");
198 std::ostringstream ost;
199 ost << datadir <<
"/livermore/pair/pp-cs-" << Z <<
".dat";
200 std::ifstream fin(ost.str().c_str());
205 ed <<
"G4LivermorePolarizedGammaConversionModel data file <" << ost.str().c_str()
206 <<
"> is not opened!" <<
G4endl;
207 G4Exception(
"G4LivermorePolarizedGammaConversionModel::ReadData()",
209 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
216 <<
" is opened by G4LivermorePolarizedGammaConversionModel" <<
G4endl;}
239 G4cout <<
"G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()"
248 if(intZ < 1 || intZ >
maxZ) {
return xs; }
258 if(!pv) {
return xs; }
261 xs = pv->
Value(GammaEnergy);
266 G4cout <<
"****** DEBUG: tcs value for Z=" << Z <<
" at energy (MeV)="
268 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
269 G4cout <<
" -> first cs value in EADL data file (iu) =" << (*pv)[0] <<
G4endl;
270 G4cout <<
" -> last cs value in EADL data file (iu) =" << (*pv)[
n] <<
G4endl;
271 G4cout <<
"*********************************************************" <<
G4endl;
293 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" <<
G4endl;
312 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
318 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
328 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
334 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
347 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - element = 0" <<
G4endl;
356 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - ionisation = 0" <<
G4endl;
363 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
367 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
371 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
373 G4double epsilonRange = 0.5 - epsilonMin ;
387 epsilon = 0.5 - epsilonRange * pow(
G4UniformRand(), 0.3333) ;
388 screen = screenFactor / (epsilon * (1. -
epsilon));
394 screen = screenFactor / (epsilon * (1 -
epsilon));
412 electronTotEnergy = (1. -
epsilon) * photonEnergy;
413 positronTotEnergy = epsilon * photonEnergy;
417 positronTotEnergy = (1. -
epsilon) * photonEnergy;
418 electronTotEnergy = epsilon * photonEnergy;
440 G4double Ene = electronTotEnergy/electron_mass_c2;
456 phi =
SetPhi(photonEnergy);
457 psi =
SetPsi(photonEnergy,phi);
512 G4double electronKineEnergy =
std::max(0.,electronTotEnergy - electron_mass_c2) ;
523 Ene = positronTotEnergy/electron_mass_c2;
532 dirX = sinTheta*cos(phip);
533 dirY = sinTheta*sin(phip);
537 G4double positronKineEnergy =
std::max(0.,positronTotEnergy - electron_mass_c2) ;
543 positronDirection, positronKineEnergy);
546 fvect->push_back(particle1);
547 fvect->push_back(particle2);
574 if (screenVariable > 1.)
575 value = 42.24 - 8.368 * log(screenVariable + 0.952);
577 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
590 if (screenVariable > 1.)
591 value = 42.24 - 8.368 * log(screenVariable + 0.952);
593 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
605 G4double Momentum = sqrt(Energy*Energy -1);
608 *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
609 *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
633 const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
634 const G4double aw = 0.0151, bw = 10.7, cw = -410.;
636 const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
638 pl[0] =
Fln(ay0,by0,Ene);
639 pl[1] = aa0 + ba0*(Ene);
640 pl[2] =
Poli(aw,bw,cw,Ene);
641 pl[3] =
Poli(axc,bxc,cxc,Ene);
643 const G4double abf = 3.1216, bbf = 2.68;
645 pt[1] = abf + bbf/Ene;
661 const G4double aw = 0.21, bw = 10.8, cw = -58.;
662 const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
664 pl[0] =
Fln(ay0, by0, Ene);
665 pl[1] =
Fln(aa0, ba0, Ene);
666 pl[2] =
Poli(aw,bw,cw,Ene);
667 pl[3] =
Poli(axc,bxc,cxc,Ene);
739 const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
740 const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
741 const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
742 const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
743 const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
744 const G4double axcp = 2.8E-3,bxcp = -3.133;
745 const G4double abf0 = 3.1213, bbf0 = 2.61;
746 const G4double abfpm = 3.1231, bbfpm = 2.84;
748 p0l[0] =
Fln(ay00, by00, Ene);
749 p0l[1] =
Fln(aa00, ba00, Ene);
750 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
751 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
753 ppml[0] =
Fln(ay0p, by0p, Ene);
754 ppml[1] = aap + bap*(Ene);
755 ppml[2] =
Poli(awp, bwp, cwp, Ene);
756 ppml[3] =
Fln(axcp,bxcp,Ene);
759 p0t[1] = abf0 + bbf0/Ene;
761 ppmt[1] = abfpm + bbfpm/Ene;
767 xe0 =
Encu(p0l, p0t, xi);
769 xepm =
Encu(ppml, ppmt, xi);
774 const G4double ay00 = 2.82, by00 = 6.35;
775 const G4double aa00 = -1.75, ba00 = 0.25;
777 const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
778 const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
779 const G4double ay0p = 1.56, by0p = 3.6;
780 const G4double aap = 0.86, bap = 8.3E-3;
781 const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
784 p0l[0] =
Fln(ay00, by00, Ene);
785 p0l[1] = aa00+pow(Ene, ba00);
786 p0l[2] =
Poli(aw0, bw0, cw0, Ene);
787 p0l[3] =
Poli(axc0, bxc0, cxc0, Ene);
788 ppml[0] =
Fln(ay0p, by0p, Ene);
789 ppml[1] = aap + bap*(Ene);
790 ppml[2] =
Poli(awp, bwp, cwp, Ene);
801 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
805 b =
Ftan(ppmt,PhiLocal);
809 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
813 a =
Ftan(p0t,PhiLocal);
818 b = (ppml[0]+2*ppml[1]*ppml[2]*
Flor(ppml,PhiLocal));
819 a = (p0l[0]+2*p0l[1]*p0l[2]*
Flor(p0l,PhiLocal));
840 r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
853 value =(a + b/x + c/(x*x*
x));
887 fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
888 (Fdlor(p_p1,x) - Fdtan(p_p2,x));
890 if(x > xmax) {
return xmax; }
895 if(std::fabs(fx) <= x*1.0e-6) {
break; }
898 if(x < 0.0) { x = 0.0; }
911 value = 1./(
pi*(w*w + 4.*(x-xc)*(x-xc)));
924 value = (y0 *
pi*(w*w + 4.*(x-xc)*(x-xc)) + 2.*A*w);
937 value = (-16.*A*w*(x-xc))/
938 (
pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
951 value = y0*x + A*atan( 2*(x-xc)/w) /
pi;
965 nor = atan(2.*(
pi-xc)/w)/(2.*
pi*
w) - atan(2.*(x-xc)/
w)/(2.*
pi*w);
966 value = xc - (w/2.)*tan(-2.*r*nor*
pi*w+atan(2*(xc-x)/w));
989 value = -1.*a / ((x-b)*(x-b));
1012 value = b*(1-exp(r*cnor/a));
1052 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
1053 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
1054 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
1078 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
1095 G4double direction_x = direction1.getX();
1096 G4double direction_y = direction1.getY();
1097 G4double direction_z = direction1.getZ();
1099 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
1113 G4AutoLock l(&LivermorePolarizedGammaConversionModelMutex);
G4double Poli(G4double, G4double, G4double, G4double)
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
G4double LowEnergyLimit() const
virtual ~G4LivermorePolarizedGammaConversionModel()
G4LivermorePolarizedGammaConversionModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedGammaConversion")
G4double Finvtan(G4double *, G4double, G4double)
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
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)
const G4double w[NPOINTSGL]
G4double GetfCoulomb() const
void ReadData(size_t Z, const char *path=0)
G4ParticleDefinition * GetDefinition() const
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double Flor(G4double *, G4double)
#define G4MUTEX_INITIALIZER
G4double SetPsi(G4double, G4double)
G4double Finvlor(G4double *, G4double, G4double)
G4double Ftan(G4double *, G4double)
G4ThreeVector GetPerpendicularPolarization(const G4ThreeVector &direction0, const G4ThreeVector &polarization0) const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
size_t GetTableSize() const
G4double Fdlor(G4double *, G4double)
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
G4double Fintlor(G4double *, G4double)
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
std::vector< G4EmElementSelector * > * GetElementSelectors()
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double GetlogZ3() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double Fln(G4double, G4double, G4double)
void SystemOfRefChange(G4ThreeVector &direction0, G4ThreeVector &direction1, G4ThreeVector &polarization0)
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
G4double ScreenFunction2(G4double screenVariable)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double x[NPOINTSGL]
G4double ScreenFunction1(G4double screenVariable)
G4IonisParamElm * GetIonisation() const
G4double Glor(G4double *, G4double)
G4double SetPhi(G4double)
const G4ThreeVector & GetPolarization() const
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
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)
size_t GetNumberOfElements() const
G4ParticleChangeForGamma * fParticleChange
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4LPhysicsFreeVector * data[100]
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)
double epsilon(double density, double temperature)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const