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;
173 if (photonEnergy < smallEnergy )
175 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
179 electronTotEnergy = (1. - epsilon) * photonEnergy;
180 positronTotEnergy = epsilon * photonEnergy;
184 positronTotEnergy = (1. - epsilon) * photonEnergy;
185 electronTotEnergy = epsilon * photonEnergy;
194 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries" <<
G4endl;
198 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - element = 0"
205 G4cout <<
"G4LivermoreGammaConversionModelRC::SampleSecondaries - ionisation = 0"
212 if (photonEnergy > 50. *
MeV) fZ += 8. * (element->
GetfCoulomb());
216 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
217 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
220 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
221 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
222 G4double epsilonRange = 0.5 - epsilonMin ;
228 G4double f10 = ScreenFunction1(screenMin) - fZ;
229 G4double f20 = ScreenFunction2(screenMin) - fZ;
230 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
231 G4double normF2 = std::max(1.5 * f20,0.);
232 G4double a=393.3750918,
b=115.3070201,
c=810.6428451,
d=19.96497475,
e=1016.874592,
f=1.936685510,
233 gLocal=751.2140962, h=0.099751048, i=299.9466339, j=0.002057250, k=49.81034926;
234 G4double aa=-18.6371131, bb=-1729.95248, cc=9450.971186, dd=106336.0145, ee=55143.09287,
ff=-117602.840,
235 gg=-721455.467,
hh=693957.8635, ii=156266.1085, jj=533209.9347;
237 G4double logepsMin = log(epsilonMin);
238 G4double NormaRC = a +
b*logepsMin +
c/logepsMin +
d*pow(logepsMin,2.) +
e/pow(logepsMin,2.) +
f*pow(logepsMin,3.) +
239 gLocal/pow(logepsMin,3.) + h*pow(logepsMin,4.) + i/pow(logepsMin,4.) + j*pow(logepsMin,5.) +
246 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.3333) ;
247 screen = screenFactor / (epsilon * (1. - epsilon));
248 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
253 screen = screenFactor / (epsilon * (1 - epsilon));
254 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
261 G4double deltaP_R1 = 1. + (a +
b*logepsilon +
c/logepsilon +
d*pow(logepsilon,2.) +
e/pow(logepsilon,2.) +
262 f*pow(logepsilon,3.) + gLocal/pow(logepsilon,3.) + h*pow(logepsilon,4.) + i/pow(logepsilon,4.) +
263 j*pow(logepsilon,5.) + k/pow(logepsilon,5.))/100.;
264 G4double deltaP_R2 = 1.+((aa + cc*logepsilon + ee*pow(logepsilon,2.) + gg*pow(logepsilon,3.) + ii*pow(logepsilon,4.))
265 / (1. + bb*logepsilon + dd*pow(logepsilon,2.) +
ff*pow(logepsilon,3.) +
hh*pow(logepsilon,4.)
266 + jj*pow(logepsilon,5.) ))/100.;
270 Rechazo = deltaP_R1/NormaRC;
274 Rechazo = deltaP_R2/NormaRC;
276 G4cout << Rechazo <<
" " << NormaRC <<
" " << epsilon <<
G4endl;
279 electronTotEnergy = (1. - epsilon) * photonEnergy;
280 positronTotEnergy = epsilon * photonEnergy;
309 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
310 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
322 electronDirection.
rotateUz(photonDirection);
334 positronDirection.
rotateUz(photonDirection);
338 positronDirection, positronKineEnergy);
341 fvect->push_back(particle1);
342 fvect->push_back(particle2);
352 G4double G4LivermoreGammaConversionModelRC::ScreenFunction1(
G4double screenVariable)
358 if (screenVariable > 1.)
359 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
361 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
368 G4double G4LivermoreGammaConversionModelRC::ScreenFunction2(
G4double screenVariable)
374 if (screenVariable > 1.)
375 value = 42.24 - 8.368 * std::log(screenVariable + 0.952);
377 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);