From: morsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Thu, 9 May 2013 08:51:24 +0000 (+0000)
Subject: Protections against arith. exceptions added.
X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=461c66a409c782ddc45c4e52aff6435b28f596f9;hp=8bc1dd564fae313f29b0aa48a37adab89b35d5e7

Protections against arith. exceptions added.
---

diff --git a/PYTHIA8/pythia8175/src/BeamParticle.cxx b/PYTHIA8/pythia8175/src/BeamParticle.cxx
index 66e3f67c320..fe06947e45c 100644
--- 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 e335445b2a5..da94b4de115 100644
--- 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 20fdf3ca026..8488efc2fe4 100644
--- 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;