double BeamParticle::xCompFrac(double xs) {
// Select case by power of gluon (1-x_g) shape.
+ const double big = 1.e20;
+ const double tiny = 1.e-20;
switch (companionPower) {
case 0:
+ 6. * xs * log(xs) * (3. + 2. * xs * (3.+xs)) );
default:
- return ( -9. * xs * (xs*xs - 1.) * (5. + xs * (24. + xs)) + 12. * xs
- * log(xs) * (1. + 2. * xs) * (1. + 2. * xs * (5. + 2. * xs)) )
- / ( 8. * (1. + 2. * xs) * ((xs - 1.) * (1. + xs * (10. + xs))
- - 6. * xs * log(xs) * (1. + xs)) );
-
+ {
+ double p1 =
+ ( -9. * xs * (xs*xs - 1.) * (5. + xs * (24. + xs)) + 12. * xs
+ * log(xs) * (1. + 2. * xs) * (1. + 2. * xs * (5. + 2. * xs)) );
+ double q1 =
+ ( 8. * (1. + 2. * xs) * ((xs - 1.) * (1. + xs * (10. + xs))
+ - 6. * xs * log(xs) * (1. + xs)) );
+ if (abs(q1) > tiny)
+ return (p1/q1);
+ else
+ return p1 * big;
+ }
}
}
// The value corresponds to an unrescaled range between 0 and 1 - x_s.
double BeamParticle::xCompDist(double xc, double xs) {
+ const double big = 1.e20;
+ const double tiny = 1.e-20;
// Mother gluon momentum fraction. Check physical limit.
double xg = xc + xs;
+ 6. * xs * log(xs) * (3. + 2. * xs * (3. + xs)) );
default:
- return fac * pow4(1. - xg) / ( 2. * (1. + 2. * xs) * ((1. - xs)
- * (1. + xs * (10. + xs)) + 6. * xs * log(xs) * (1. + xs)) );
-
+ {
+ double p1 = fac * pow4(1. - xg);
+ double q1 = ( 2. * (1. + 2. * xs) * ((1. - xs)
+ * (1. + xs * (10. + xs))
+ + 6. * xs * log(xs) * (1. + xs)) );
+
+ if (abs(q1) > tiny)
+ return (p1/q1);
+ else
+ return p1 * big;
+ }
}
}
double pPosNeg = pPos * pNeg;
double kXPos = eX * pPos / pPosNeg;
double kXNeg = eX * pNeg / pPosNeg;
- double kXX = 1. / sqrt( 1. + 2. * kXPos * kXNeg * pPosNeg );
+
+ double arg = 1. + 2. * kXPos * kXNeg * pPosNeg;
+ double kXX = 1.e20;
+ if (arg > 0.)
+ kXX = 1. / sqrt( arg );
+
+
double kYPos = eY * pPos / pPosNeg;
double kYNeg = eY * pNeg / pPosNeg;
double kYX = kXX * (kXPos * kYNeg + kXNeg * kYPos) * pPosNeg;
- double kYY = 1. / sqrt(1. + 2. * kYPos * kYNeg * pPosNeg - pow2(kYX));
+
+ arg = 1. + 2. * kYPos * kYNeg * pPosNeg - pow2(kYX);
+ double kYY = 1.e20;
+ if (arg > 0.)
+ kYY = 1. / sqrt(arg);
+
eX = kXX * (eX - kXNeg * pPos - kXPos * pNeg);
eY = kYY * (eY - kYNeg * pPos - kYPos * pNeg - kYX * eX);
double HMETau2FourPions::rhoFormFactor1(double s) {
- double f = sqrtpos(1 - 4*picM*picM/s);
- if (s > 4*picM*picM)
+ double f = 0.;
+ if (s > 4*picM*picM) {
+ f = sqrtpos(1 - 4*picM*picM/s);
f = f * log((1 + f) / (1 - f)) * (s - 4*picM*picM) / M_PI;
+ }
else if (s < 0.0000001)
f = -8 * picM*picM / M_PI;
else
// Return the form factor for the rho(770) (equivalent to h(s) derivative).
double HMETau2FourPions::rhoFormFactor2(double s) {
-
- double f = sqrtpos(1 - 4*picM*picM/s);
- if (s > 4*picM*picM)
+ double f = 0.;
+
+ if (s > 4*picM*picM) {
+ f = sqrtpos(1 - 4*picM*picM/s);
f = f / (M_PI * s) * (s*f + (2*picM*picM + s)*log((1 + f) / (1 - f)));
+ }
else
f = 0;
return f;