Protections against arith. exceptions added.

diff --git a/PYTHIA8/pythia8175/src/BeamParticle.cxx b/PYTHIA8/pythia8175/src/BeamParticle.cxx
index 66e3f67c32025d9a6c175a692a3af6c93f164e5c..fe06947e45c9d8dad8c95cdd0db5362e9c19bd4d 100644 (file)
--- a/PYTHIA8/pythia8175/src/BeamParticle.cxx
+++ b/PYTHIA8/pythia8175/src/BeamParticle.cxx
@@ -400,6 +400,8 @@ double BeamParticle::xValFrac(int j, double Q2) {
 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: 
@@ -423,11 +425,18 @@ double BeamParticle::xCompFrac(double xs) {
         + 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;
+      }
   }
 }
 
@@ -438,6 +447,8 @@ double BeamParticle::xCompFrac(double xs) {
 // 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;
@@ -465,9 +476,17 @@ double BeamParticle::xCompDist(double xc, double 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;
+      }
   }
 }
 

diff --git a/PYTHIA8/pythia8175/src/FragmentationSystems.cxx b/PYTHIA8/pythia8175/src/FragmentationSystems.cxx
index e335445b2a5ad66e0fb3a184c180f2b8497f3fb9..da94b4de1155e7132c225a6bf21808cedeed9cb7 100644 (file)
--- a/PYTHIA8/pythia8175/src/FragmentationSystems.cxx
+++ b/PYTHIA8/pythia8175/src/FragmentationSystems.cxx
@@ -475,11 +475,22 @@ void StringRegion::setUp(Vec4 p1, Vec4 p2, bool isMassless) {
   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);
 

diff --git a/PYTHIA8/pythia8175/src/HelicityMatrixElements.cxx b/PYTHIA8/pythia8175/src/HelicityMatrixElements.cxx
index 20fdf3ca02606402e52bd13f9121db0919c9fc7c..8488efc2fe48ee25641f01ffa50c6695d4b3d953 100644 (file)
--- a/PYTHIA8/pythia8175/src/HelicityMatrixElements.cxx
+++ b/PYTHIA8/pythia8175/src/HelicityMatrixElements.cxx
@@ -2191,9 +2191,11 @@ double HMETau2FourPions::a1FormFactor(double s) {
 
 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
@@ -2207,10 +2209,12 @@ double HMETau2FourPions::rhoFormFactor1(double s) {
 // 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;