70 : fPAIxscVector(nullptr),
71 fPAIdEdxVector(nullptr),
72 fPAIphotonVector(nullptr),
73 fPAIelectronVector(nullptr),
74 fChCosSqVector(nullptr),
75 fChWidthVector(nullptr)
96 for(j = 1; j < 5 ; j++)
137 energy1 = (*(*fMatSandiaMatrix)[i])[0];
138 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
140 if( energy2 - energy1 > 1.5*
fDelta*(energy1 + energy2) )
continue ;
143 for(j = i; j < fIntervalNumber-1; j++)
145 for( k = 0; k < 5; k++ )
174 energy1 = (*(*fMatSandiaMatrix)[i])[0];
175 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
189 for(j = 1; j < 5 ; j++)
237 a1 = (*(*fMatSandiaMatrix)[
k])[1];
238 a2 = (*(*fMatSandiaMatrix)[
k])[2];
239 a3 = (*(*fMatSandiaMatrix)[
k])[3];
240 a4 = (*(*fMatSandiaMatrix)[
k])[4];
242 c1 = (x2 -
x1)/x1/x2 ;
243 c2 = (x2 -
x1)*(x2 + x1)/x1/x1/x2/
x2 ;
244 c3 = (x2 -
x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/
x2 ;
247 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
258 G4double energy1, energy2, result = 0.;
264 energy1 = (*(*fMatSandiaMatrix)[i])[0];
271 energy1 = (*(*fMatSandiaMatrix)[i])[0];
277 energy1 = (*(*fMatSandiaMatrix)[i])[0];
278 energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
296 G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
298 a1 = (*(*fMatSandiaMatrix)[
k])[1];
299 a2 = (*(*fMatSandiaMatrix)[
k])[2];
300 a3 = (*(*fMatSandiaMatrix)[
k])[3];
301 a4 = (*(*fMatSandiaMatrix)[
k])[4];
303 energy2 = energy1*energy1;
304 energy3 = energy2*energy1;
305 energy4 = energy3*energy1;
307 result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
308 result *=
hbarc/energy1 ;
325 eIm2 = result*result;
328 eRe2 = result*result;
330 result = eIm2 + eRe2;
345 G4double x0, x02, x03, x04, x05,
x1,
x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
346 c1,
c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
353 x1 = (*(*fMatSandiaMatrix)[i])[0];
354 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
356 a1 = (*(*fMatSandiaMatrix)[i])[1];
357 a2 = (*(*fMatSandiaMatrix)[i])[2];
358 a3 = (*(*fMatSandiaMatrix)[i])[3];
359 a4 = (*(*fMatSandiaMatrix)[i])[4];
363 if(x0 >= x1) x0 = x1*(1+
fDelta);
368 if(x0 >= x2) x0 = x2*(1+
fDelta);
375 if( xx12 < 0 ) xx12 = -xx12;
379 xln3 = log((x2 + x0)/(x1 + x0)) ;
386 c1 = (x2 -
x1)/x1/x2 ;
387 c2 = (x2 -
x1)*(x2 +x1)/x1/x1/x2/
x2 ;
388 c3 = (x2 -
x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/
x2 ;
390 result -= (a1/x02 + a3/x04)*xln1 ;
391 result -= (a2/x02 + a4/x04)*c1 ;
392 result -= a3*c2/2/x02 ;
393 result -= a4*c3/3/x02 ;
395 cof1 = a1/x02 + a3/x04 ;
396 cof2 = a2/x03 + a4/x05 ;
398 result += 0.5*(cof1 +cof2)*xln2 ;
399 result += 0.5*(cof1 - cof2)*xln3 ;
418 G4double be2,cof,
x1,
x2,
x3,
x4,
x5,
x6,x7,x8,result ;
423 static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
424 be2 = betaGammaSq/(1 + betaGammaSq) ;
430 if( betaGammaSq < 0.01 ) x2 = log(be2) ;
433 x2 = -log( (1/betaGammaSq - epsilonRe)*
434 (1/betaGammaSq - epsilonRe) +
435 epsilonIm*epsilonIm )/2 ;
437 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
443 x3 = -epsilonRe + 1/betaGammaSq ;
444 x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
445 epsilonIm*epsilonIm) ;
447 x7 = atan2(epsilonIm,x3) ;
452 x4 = ((x1 +
x2)*epsilonIm + x6)/
hbarc ;
454 x8 = (1 + epsilonRe)*(1 + epsilonRe) +
457 result = (x4 + cof*integralTerm/omega/omega) ;
458 if(result < 1.0
e-8) result = 1.0e-8 ;
459 result *= fine_structure_const/be2/
pi ;
462 result *= (1-exp(-be4/betaBohr4)) ;
497 static const G4double cofBetaBohr = 4.0 ;
499 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
501 be2 = betaGammaSq/(1 + betaGammaSq) ;
504 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ;
507 logarithm = -log( (1/betaGammaSq - epsilonRe)*
508 (1/betaGammaSq - epsilonRe) +
509 epsilonIm*epsilonIm )*0.5 ;
510 logarithm += log(1+1.0/betaGammaSq) ;
513 if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
519 x3 = -epsilonRe + 1.0/betaGammaSq ;
520 x5 = -1.0 - epsilonRe +
521 be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
522 epsilonIm*epsilonIm) ;
523 if( x3 == 0.0 ) argument = 0.5*
pi;
524 else argument = atan2(epsilonIm,x3) ;
527 dNdxC = ( logarithm*epsilonIm + argument )/
hbarc ;
529 if(dNdxC < 1.0
e-8) dNdxC = 1.0e-8 ;
531 dNdxC *= fine_structure_const/be2/
pi ;
533 dNdxC *= (1-exp(-be4/betaBohr4)) ;
537 modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
558 G4double cof, resonance, modul2, dNdxP ;
562 static const G4double cofBetaBohr = 4.0 ;
564 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
566 be2 = betaGammaSq/(1 + betaGammaSq) ;
570 resonance *= epsilonIm/
hbarc ;
573 dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
575 if( dNdxP < 1.0
e-8 ) dNdxP = 1.0e-8 ;
577 dNdxP *= fine_structure_const/be2/
pi ;
578 dNdxP *= (1-exp(-be4/betaBohr4)) ;
582 modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
599 G4double energy1, energy2, result = 0.;
619 for( k =
fPAIbin - 2; k >= 0; k-- )
647 for( i = i2; i >= i1; i-- )
651 if( i==i2 ) result += integral.
Legendre10(
this,
655 else if( i == i1 ) result += integral.
Legendre10(
this,
680 G4double energy1, energy2, result = 0.;
700 for( k =
fPAIbin - 2; k >= 0; k-- )
728 for( i = i2; i >= i1; i-- )
732 if( i==i2 ) result += integral.
Legendre10(
this,
736 else if( i == i1 ) result += integral.
Legendre10(
this,
760 G4double energy1, energy2, beta2, module2, cos2,
width, result = 0.;
788 for( k =
fPAIbin - 2; k >= 0; k-- )
824 for( i = i2; i >= i1; i-- )
828 if( i==i2 ) result += integral.
Legendre10(
this,
832 else if( i == i1 ) result += integral.
Legendre10(
this,
856 G4double energy1, energy2, result = 0.;
876 for( k =
fPAIbin - 2; k >= 0; k-- )
904 for( i = i2; i >= i1; i-- )
908 if( i==i2 ) result += integral.
Legendre10(
this,
912 else if( i == i1 ) result += integral.
Legendre10(
this,
937 omega2 = omega*omega ;
938 omega3 = omega2*omega ;
939 omega4 = omega2*omega2 ;
947 G4cout<<
"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<
G4endl;
951 a1 = (*(*fMatSandiaMatrix)[i])[1];
952 a2 = (*(*fMatSandiaMatrix)[i])[2];
953 a3 = (*(*fMatSandiaMatrix)[i])[3];
954 a4 = (*(*fMatSandiaMatrix)[i])[4];
956 lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);