fixed numerical problem for calculating pipiS-wave dynamics at pole positions

PwaDynamics/FVector.cc

@@ -22,6 +22,7 @@
 //************************************************************************//
 
 #include "PwaDynamics/FVector.hh"
+#include "ErrLogger/ErrLogger.hh"
 //#include "qft++/relativistic-quantum-mechanics/Utils.hh"
 
 
@@ -88,6 +89,11 @@ complex<double> FVector::evalProjMatrix(const double mass, int index){
   for(int i=0;i<NumRows(); ++i){
     result+=denomMatrInv(index,i)*(*_Pvector)(i,0);
   }
+
+  if(result.real() != result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
   return result; 
 }
 

PwaDynamics/KMatrixSlowAdlerCorRel.cc

@@ -60,8 +60,21 @@ void KMatrixSlowAdlerCorRel::evalMatrix(const double mass){
     (*this)+= *(*it);
   }
 
-  double adlerFactor=(1.-_sAdler0)/(mass*mass-_sAdler0)*(mass*mass-_sAdler*0.1349766*0.1349766/2.);
-  double s0ScatFactor=(1.0-_s0Scat)/(mass*mass-_s0Scat);
+  double sAdler0Denom=mass*mass-_sAdler0;
+  if( fabs(sAdler0Denom)<1e-10){
+    if(sAdler0Denom<0) sAdler0Denom=-1e-10;
+    else sAdler0Denom=1e-10;
+  }
+
+ 
+  double adlerFactor=(1.-_sAdler0)/sAdler0Denom*(mass*mass-_sAdler*0.1349766*0.1349766/2.);
+
+  double s0ScatDenom=mass*mass-_s0Scat;
+  if( fabs(s0ScatDenom)<1e-10){
+    if(s0ScatDenom<0.) s0ScatDenom=-1e-10;
+    else s0ScatDenom=1e-10;
+  }
+  double s0ScatFactor=(1.0-_s0Scat)/s0ScatDenom;
 
   for (int i=0; i<NumRows(); ++i){
     for (int j=0; j<NumCols(); ++j){

PwaDynamics/PVectorSlowCorRel.cc

@@ -47,8 +47,14 @@ void PVectorSlowCorRel::evalMatrix(const double mass){
     thePVector += *(*it);
   }
 
+  double s0prodDenom=mass*mass-_s0prod;
+  if(fabs(s0prodDenom)<1e-10){
+    if (s0prodDenom<0.) s0prodDenom=-1e-10;
+    else s0prodDenom=1e-10;
+  }
+
   for (int i=0; i<thePVector.NumRows(); ++i){  
-    this->operator()(i,0)=thePVector(i,0)+ _fProdVec[i]*(1.0-_s0prod)/(mass*mass-_s0prod);
+    this->operator()(i,0)=thePVector(i,0)+ _fProdVec[i]*(1.0-_s0prod)/s0prodDenom;
   } 
 
 }

PwaUtils/HeliDecAmps.cc

@@ -129,7 +129,7 @@ complex<double> HeliDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs,
   int evtNo=theData->evtNo;
 
   if ( _cacheAmps && !_recalculate){
-    result=_cachedAmpMap[evtNo][lamX][lamFs];
+    result=_cachedAmpMap.at(evtNo).at(lamX).at(lamFs);
     result*=_absDyn->eval(theData, grandmaAmp);
     if(result.real()!=result.real()) DebugMsg << "result:\t" << result << endmsg;
     return result;
@@ -148,9 +148,9 @@ complex<double> HeliDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs,
     if(_enabledlamFsDaughter2 && lamFs!=lambda2) continue;
 
     double theMag=it->second;
-    double thePhi=_currentParamPhiLamLams[it->first];
+    double thePhi=_currentParamPhiLamLams.at(it->first);
     complex<double> expi(cos(thePhi), sin(thePhi));
-    complex<double> amp = it->first->parityFactor*theMag*expi*conj( theData->WignerDsString[_wignerDKey][it->first->J][lamX][lambda]);
+    complex<double> amp = it->first->parityFactor*theMag*expi*conj( theData->WignerDsString.at(_wignerDKey).at(it->first->J).at(lamX).at(lambda));
     result+=amp*daughterAmp(lambda1, lambda2, theData, lamFs);
   }
 
@@ -163,7 +163,10 @@ complex<double> HeliDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs,
 }
 
   result*=_absDyn->eval(theData, grandmaAmp);
-  if(result.real()!=result.real()) DebugMsg << "result:\t" << result << endmsg;
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
   return result;
 }
 

PwaUtils/LSDecAmps.cc

@@ -140,6 +140,13 @@ complex<double> LSDecAmps::XdecAmp(Spin& lamX, EvtData* theData, Spin& lamFs, Ab
   }
 
   result*=_absDyn->eval(theData, grandmaAmp);
+  if(result.real()!=result.real()){
+    Info << "dyn name: " << _absDyn->name() 
+	 << "\nname(): " << name()
+	 << endmsg;
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
   return result;
 }
 
@@ -169,6 +176,10 @@ complex<double> LSDecAmps::lsLoop(Spin& lamX, EvtData* theData, Spin& lam1Min, S
     }
   }
   result*=_preFactor*_isospinCG;
+  if(result.real()!=result.real()){
+    Alert << "result:\t" << result << endmsg;
+    exit(0);
+  }
   return result;
 }