52 :
G4VEmModel(nam),fParticleChange(0),smallEnergy(2.*
MeV),isInitialised(false),
53 crossSectionHandler(0),meanFreePathTable(0)
56 highEnergyLimit = 100 *
GeV;
67 if(verboseLevel > 0) {
68 G4cout <<
"Livermore Gamma conversion is constructed " <<
G4endl
70 << lowEnergyLimit /
MeV <<
" MeV - "
71 << highEnergyLimit /
GeV <<
" GeV"
80 if (crossSectionHandler)
delete crossSectionHandler;
90 G4cout <<
"Calling G4LivermoreGammaConversionModelRC::Initialise()" <<
G4endl;
92 if (crossSectionHandler)
94 crossSectionHandler->
Clear();
95 delete crossSectionHandler;
101 crossSectionHandler->
Initialise(0,lowEnergyLimit,100.*
GeV,400);
102 G4String crossSectionFile =
"pair/pp-cs-";
103 crossSectionHandler->
LoadData(crossSectionFile);
107 if (verboseLevel > 2)
108 G4cout <<
"Loaded cross section files for Livermore Gamma Conversion model RC" <<
G4endl;
110 if (verboseLevel > 0) {
111 G4cout <<
"Livermore Gamma Conversion model is initialized " << G4endl
118 if(isInitialised)
return;
120 isInitialised =
true;
131 if (verboseLevel > 3) {
132 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModelRC"
135 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit)
return 0;
160 if (verboseLevel > 3)
161 G4cout <<
"Calling SampleSecondaries() of G4LivermoreGammaConversionModelRC" <<
G4endl;
175 if (photonEnergy < smallEnergy )
177 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
181 electronTotEnergy = (1. - epsilon) * photonEnergy;
182 positronTotEnergy = epsilon * photonEnergy;
186 positronTotEnergy = (1. - epsilon) * photonEnergy;
187 electronTotEnergy = epsilon * photonEnergy;
196 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries" <<
G4endl;
200 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
207 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
214 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
218 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
222 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
224 G4double epsilonRange = 0.5 - epsilonMin ;
230 G4double f10 = ScreenFunction1(screenMin) - fZ;
231 G4double f20 = ScreenFunction2(screenMin) - fZ;
234 G4double a=393.3750918,
b=115.3070201,
c=810.6428451, d=19.96497475, e=1016.874592, f=1.936685510,
235 gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
236 G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287, ff=-117602.840,
237 gg=-721455.467, hh=693957.8635, ii=156266.1085, jj=533209.9347;
239 G4double logepsMin = log(epsilonMin);
240 G4double NormaRC = a +
b*logepsMin +
c/logepsMin + d*pow(logepsMin,2.) + e/pow(logepsMin,2.) + f*pow(logepsMin,3.) +
241 gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
245 G4double HardPhotonThreshold = 0.08;
246 G4double r1, r2, r3, beta=0, gbeta, sigt = 582.068, sigh, rejet;
248 G4double cg = (11./2.)/(exp(-11.*HardPhotonThreshold/2.)-exp(-11./2.));
252 sigh = 1028.58*exp(-HardPhotonThreshold/0.09033) + 136.63;
253 if (r1 > 1.- sigh/sigt) {
258 beta = (-2./11.)*log(exp(-0.08*11./2.)-r3*11./(2.*cg));
259 gbeta = exp(-11.*beta/2.);
260 rejet = fbeta(beta)/(8000.*gbeta);
262 HardPhotonEnergy = beta * photonEnergy;
265 HardPhotonEnergy = 0.;
268 photonEnergy -= HardPhotonEnergy;
274 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.3333) ;
275 screen = screenFactor / (epsilon * (1. - epsilon));
276 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
281 screen = screenFactor / (epsilon * (1 - epsilon));
282 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
289 G4double deltaP_R1 = 1. + (a +
b*logepsilon +
c/logepsilon + d*pow(logepsilon,2.) + e/pow(logepsilon,2.) +
290 f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
291 j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
292 G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
293 / (1. + bb*logepsilon + dd*pow(logepsilon,2.) + ff*pow(logepsilon,3.) + hh*pow(logepsilon,4.)
294 + jj*pow(logepsilon,5.) ))/100.;
298 Rechazo = deltaP_R1/NormaRC;
302 Rechazo = deltaP_R2/NormaRC;
304 G4cout << Rechazo <<
" " << NormaRC <<
" " << epsilon <<
G4endl;
307 electronTotEnergy = (1. - epsilon) * photonEnergy;
308 positronTotEnergy = epsilon * photonEnergy;
337 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
338 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
350 electronDirection.
rotateUz(photonDirection);
362 positronDirection.
rotateUz(photonDirection);
370 fvect->push_back(particle1);
371 fvect->push_back(particle2);
373 if (HardPhotonEnergy > 0.)
377 G4double dxHardP= std::sin(thetaHardPhoton)*std::cos(phi);
378 G4double dyHardP= std::sin(thetaHardPhoton)*std::sin(phi);
379 G4double dzHardP =std::cos(thetaHardPhoton);
381 G4ThreeVector hardPhotonDirection (dxHardP, dyHardP, dzHardP);
382 hardPhotonDirection.
rotateUz(photonDirection);
386 fvect->push_back(particle3);
398 G4double G4LivermoreGammaConversionModelRC::ScreenFunction1(
G4double screenVariable)
404 if (screenVariable > 1.)
405 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
407 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
414 G4double G4LivermoreGammaConversionModelRC::ScreenFunction2(
G4double screenVariable)
420 if (screenVariable > 1.)
421 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
423 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
434 G4double Pi, gamma, eta, d, p1, p2, p3, p4, p5, p6, p7, ffbeta;
435 gamma = (1.-
x)*(1.-x)/
x;
437 d = Dilog(1./x)-Dilog(x);
439 p1 = -1.*(25528.*pow(gamma,2) + 116044.* gamma +151556.)/105.;
440 p2 = 256.* pow(gamma,3) + 1092.* pow(gamma,2) +1260.*gamma + 420.;
441 p3 = (676.*pow(gamma,3) + 9877.*pow(gamma,2) + 58415.*gamma + 62160.)/105.;
442 p4 = 64.*pow(gamma,3) + 305.*pow(gamma,2) + 475.*gamma + 269. - 276./gamma;
443 p5 = (676.*pow(gamma,3) + 38109.*pow(gamma,2) + 211637.*gamma + 266660. - 53632./gamma)/105.;
444 p6 = 32.*pow(gamma,2) + 416.*gamma + 1310. +1184./gamma;
445 p7 = 128.*pow(gamma,3) + 802.*pow(gamma,2) + 1028.*gamma - 470. - 1184./gamma;
446 ffbeta = (1.-
x) * (p1 + p2*Pi*Pi/6. + p3*log(gamma) +
447 p4*pow(log(x),2) + (p5 + p6*log(gamma))*eta*log(x) + p7*d*eta);
459 fdilog = pow(Pi,2)/6. + (1.-y)*(log(1-y)-1.)+pow((1.-y),2)*((1./2.)*log(1.-y)-1./4.)
460 +pow((1.-y),3)*((1./3.)*log(1.-y)-1./9.)+pow((1.-y),4)*((1./4.)*log(1.-y)-1./16.);
462 if (0.5 < y && y < 2.) {
463 fdilog = 1.-y+pow((1.-y),2)/4.+pow((1.-y),3)/9.+pow((1.-y),4)/16.+
464 pow((1.-y),5)/25.+pow((1.-y),6)/36.+pow((1.-y),7)/49.;
467 fdilog = -pow(log(y),2)/2. - pow(Pi,2)/6. + (log(y)+1.)/y +
468 (log(y)/2.+1./4.)/pow(y,2) + (log(y)/3.+1./9.)/pow(y,3);
G4double LowEnergyLimit() const
virtual ~G4LivermoreGammaConversionModelRC()
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
G4double GetfCoulomb() const
G4ParticleDefinition * GetDefinition() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void SetHighEnergyLimit(G4double)
G4double FindValue(G4int Z, G4double e) const
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
G4ParticleChangeForGamma * fParticleChange
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4LivermoreGammaConversionModelRC(const G4ParticleDefinition *p=0, const G4String &nam="LivermoreConversion")
G4double GetlogZ3() const
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()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4IonisParamElm * GetIonisation() const
void LoadData(const G4String &dataFile)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
const XML_Char int const XML_Char * value
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeTrackStatus(G4TrackStatus status)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()