protection against numerical instabilities

PwaDynamics/KPole.cc

@@ -50,7 +50,8 @@ void KPole::evalMatrix(const double mass, Spin OrbMom){
   double denom=_poleMass*_poleMass-mass*mass;
 
   if( std::abs(denom) < 1.e-10){
-     denom = 1.e-10;
+      if(denom<0.) denom = -1.e-10;
+      else denom = 1.e-10;
   }
 
   for (int i=0; i< int(_g_i.size()); ++i){

PwaDynamics/KPoleBarrier.cc

@@ -59,7 +59,8 @@ void KPoleBarrier::evalMatrix(const double mass, Spin OrbMom){
   
   double denom=_poleMass*_poleMass-mass*mass;
   if(std::abs(denom) < 1e-10){
-     denom = 1.E-10;
+    if(denom<0.) denom = -1.E-10;
+    else denom = 1.E-10;
   }
 
   for (int i=0; i< int(_g_i.size()); ++i){