54 fParticleChange =
nullptr;
66 G4cout <<
"G4LivermoreGammaConversionModelRC is constructed " <<
G4endl;
92 G4cout <<
"Calling Initialise() of G4LivermoreGammaConversionModelRC."
95 << LowEnergyLimit() /
MeV <<
" MeV - "
96 << HighEnergyLimit() /
GeV <<
" GeV"
105 InitialiseElementSelectors(particle, cuts);
109 char* path = getenv(
"G4LEDATA");
116 for(
G4int i=0; i<numOfCouples; ++i)
123 for (
G4int j=0; j<nelm; ++j)
128 if(!
data[Z]) { ReadData(Z, path); }
132 if(isInitialised) {
return; }
133 fParticleChange = GetParticleChangeForGamma();
134 isInitialised =
true;
152 return lowEnergyLimit;
159 if (verboseLevel > 1)
161 G4cout <<
"Calling ReadData() of G4LivermoreGammaConversionModelRC"
165 if(
data[Z]) {
return; }
167 const char* datadir = path;
171 datadir = getenv(
"G4LEDATA");
174 G4Exception(
"G4LivermoreGammaConversionModelRC::ReadData()",
176 "Environment variable G4LEDATA not defined");
187 std::ostringstream ost;
188 ost << datadir <<
"/livermore/pair/pp-cs-" << Z <<
".dat";
189 std::ifstream
fin(ost.str().c_str());
194 ed <<
"G4LivermoreGammaConversionModelRC data file <" << ost.str().c_str()
195 <<
"> is not opened!" <<
G4endl;
196 G4Exception(
"G4LivermoreGammaConversionModelRC::ReadData()",
198 ed,
"G4LEDATA version should be G4EMLOW6.27 or later.");
205 if(verboseLevel > 3) {
G4cout <<
"File " << ost.str()
206 <<
" is opened by G4LivermoreGammaConversionModelRC" <<
G4endl;}
212 data[
Z] ->SetSpline(
true);
224 if (verboseLevel > 1)
226 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC"
230 if (GammaEnergy < lowEnergyLimit) {
return 0.0; }
236 if(intZ < 1 || intZ >
maxZ) {
return xs; }
244 InitialiseForElement(0, intZ);
246 if(!pv) {
return xs; }
249 xs = pv->
Value(GammaEnergy);
254 G4cout <<
"****** DEBUG: tcs value for Z=" << Z <<
" at energy (MeV)="
256 G4cout <<
" cs (Geant4 internal unit)=" << xs <<
G4endl;
257 G4cout <<
" -> first cs value in EADL data file (iu) =" << (*pv)[0] <<
G4endl;
258 G4cout <<
" -> last cs value in EADL data file (iu) =" << (*pv)[
n] <<
G4endl;
259 G4cout <<
"*********************************************************" <<
G4endl;
269 std::vector<G4DynamicParticle*>* fvect,
285 if (verboseLevel > 1) {
286 G4cout <<
"Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC"
302 if (photonEnergy < smallEnergy )
304 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
307 electronTotEnergy = (1. -
epsilon) * photonEnergy;
308 positronTotEnergy = epsilon * photonEnergy;
312 positronTotEnergy = (1. -
epsilon) * photonEnergy;
313 electronTotEnergy = epsilon * photonEnergy;
321 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
325 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
332 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
339 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
347 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
349 G4double epsilonRange = 0.5 - epsilonMin ;
355 G4double f10 = ScreenFunction1(screenMin) - fZ;
356 G4double f20 = ScreenFunction2(screenMin) - fZ;
362 G4double a=393.3750918, b=115.3070201, c=810.6428451,
d=19.96497475,
e=1016.874592,
f=1.936685510,
363 gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
364 G4double aa=-18.6371131,
bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287,
ff=-117602.840,
365 gg=-721455.467,
hh=693957.8635, ii=156266.1085, jj=533209.9347;
367 G4double logepsMin = log(epsilonMin);
368 G4double NormaRC = a + b*logepsMin + c/logepsMin +
d*pow(logepsMin,2.) +
e/pow(logepsMin,2.) +
f*pow(logepsMin,3.) +
369 gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
372 G4double HardPhotonThreshold = 0.08;
373 G4double r1, r2, r3,
beta=0, gbeta, sigt = 582.068, sigh, rejet;
378 sigh = 1028.58*
G4Exp(-HardPhotonThreshold/0.09033) + 136.63;
381 if (r1 > 1.- sigh/sigt) {
386 beta = (-2./11.)*log(
G4Exp(-0.08*11./2.)-r3*11./(2.*cg));
387 gbeta =
G4Exp(-11.*beta/2.);
388 rejet = fbeta(beta)/(8000.*gbeta);
390 HardPhotonEnergy = beta * photonEnergy;
393 HardPhotonEnergy = 0.;
396 photonEnergy -= HardPhotonEnergy;
404 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.333333) ;
405 screen = screenFactor / (epsilon * (1. -
epsilon));
406 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
411 screen = screenFactor / (epsilon * (1 -
epsilon));
412 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
421 G4double deltaP_R1 = 1. + (a + b*logepsilon + c/logepsilon +
d*pow(logepsilon,2.) +
e/pow(logepsilon,2.) +
422 f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
423 j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
424 G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
425 / (1. +
bb*logepsilon + dd*pow(logepsilon,2.) +
ff*pow(logepsilon,3.) +
hh*pow(logepsilon,4.)
426 + jj*pow(logepsilon,5.) ))/100.;
430 Rechazo = deltaP_R1/NormaRC;
434 Rechazo = deltaP_R2/NormaRC;
440 electronTotEnergy = (1. -
epsilon) * photonEnergy;
441 positronTotEnergy = epsilon * photonEnergy;
471 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
472 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
482 electronDirection.
rotateUz(photonDirection);
492 positronDirection.
rotateUz(photonDirection);
499 fvect->push_back(particle1);
500 fvect->push_back(particle2);
502 if (HardPhotonEnergy > 0.)
506 G4double dxHardP= std::sin(thetaHardPhoton)*std::cos(phi);
507 G4double dyHardP= std::sin(thetaHardPhoton)*std::sin(phi);
508 G4double dzHardP =std::cos(thetaHardPhoton);
510 G4ThreeVector hardPhotonDirection (dxHardP, dyHardP, dzHardP);
511 hardPhotonDirection.
rotateUz(photonDirection);
515 fvect->push_back(particle3);
519 fParticleChange->SetProposedKineticEnergy(0.);
532 if (screenVariable > 1.)
533 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
535 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
548 if (screenVariable > 1.)
549 value = 42.24 - 8.368 *
G4Log(screenVariable + 0.952);
551 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
561 G4double Pi, gamma, eta,
d, p1, p2, p3, p4, p5, p6, p7, ffbeta;
562 gamma = (1.-
x)*(1.-x)/
x;
564 d = Dilog(1./x)-Dilog(x);
566 p1 = -1.*(25528.*pow(gamma,2) + 116044.* gamma +151556.)/105.;
567 p2 = 256.* pow(gamma,3) + 1092.* pow(gamma,2) +1260.*gamma + 420.;
568 p3 = (676.*pow(gamma,3) + 9877.*pow(gamma,2) + 58415.*gamma + 62160.)/105.;
569 p4 = 64.*pow(gamma,3) + 305.*pow(gamma,2) + 475.*gamma + 269. - 276./gamma;
570 p5 = (676.*pow(gamma,3) + 38109.*pow(gamma,2) + 211637.*gamma + 266660. - 53632./gamma)/105.;
571 p6 = 32.*pow(gamma,2) + 416.*gamma + 1310. +1184./gamma;
572 p7 = 128.*pow(gamma,3) + 802.*pow(gamma,2) + 1028.*gamma - 470. - 1184./gamma;
573 ffbeta = (1.-
x) * (p1 + p2*Pi*Pi/6. + p3*log(gamma) +
574 p4*pow(log(x),2) + (p5 + p6*log(gamma))*eta*log(x) + p7*d*eta);
585 fdilog = pow(Pi,2)/6. + (1.-
y)*(log(1-y)-1.)+pow((1.-y),2)*((1./2.)*log(1.-y)-1./4.)
586 +pow((1.-y),3)*((1./3.)*log(1.-y)-1./9.)+pow((1.-y),4)*((1./4.)*log(1.-y)-1./16.);
588 if (0.5 < y && y < 2.) {
589 fdilog = 1.-y+pow((1.-y),2)/4.+pow((1.-y),3)/9.+pow((1.-y),4)/16.+
590 pow((1.-y),5)/25.+pow((1.-y),6)/36.+pow((1.-y),7)/49.;
593 fdilog = -pow(log(y),2)/2. - pow(Pi,2)/6. + (log(y)+1.)/y +
594 (log(y)/2.+1./4.)/pow(y,2) + (log(y)/3.+1./9.)/pow(y,3);
608 G4AutoLock l(&LivermoreGammaConversionModelRCMutex);
611 if(!
data[Z]) { ReadData(Z); }
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
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
static constexpr double MeV
G4double fbeta(G4double x)
const G4ThreeVector & GetMomentumDirection() const
const G4ElementVector * GetElementVector() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4double Value(G4double theEnergy, size_t &lastidx) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
void ReadData(size_t Z, const char *path=0)
G4double G4Log(G4double x)
static G4LPhysicsFreeVector * data[100]
const XML_Char const XML_Char * data
G4LivermoreGammaConversionModelRC(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreGammaConversionRC_1")
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
#define G4MUTEX_INITIALIZER
G4ParticleDefinition * GetDefinition() const
const XML_Char int const XML_Char * value
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double)
static constexpr double electron_mass_c2
virtual ~G4LivermoreGammaConversionModelRC()
static constexpr double twopi
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double ScreenFunction2(G4double screenVariable)
static G4Positron * Positron()
static G4Electron * Electron()
size_t GetTableSize() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
double epsilon(double density, double temperature)
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::vector< G4Element * > G4ElementVector
G4IonisParamElm * GetIonisation() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4double ScreenFunction1(G4double screenVariable)
const G4Material * GetMaterial() const
G4double Dilog(G4double x)
G4double GetfCoulomb() const
G4double GetlogZ3() const
static constexpr double GeV
size_t GetNumberOfElements() const
size_t GetVectorLength() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments