improved rounding handling to get a more consistent expected chisq after fitting compared to pure chisq estimate of musrfit. Essentially get rid of rounding issues.

src/classes/PFitter.cpp

@@ -2171,15 +2171,22 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave)
 
   // calculate expected chisq
   std::vector<Double_t> param;
+  std::vector<Double_t> err;
   Double_t totalExpectedChisq = 0.0;
   std::vector<Double_t> expectedchisqPerRun;
   std::vector<UInt_t> ndfPerHisto;
 
-  for (UInt_t i=0; i<fParams.size(); i++)
+  for (UInt_t i=0; i<fParams.size(); i++) {
     param.push_back(fParams[i].fValue);
+    err.push_back(fParams[i].fPosError);
+  }
 
   // CalcExpectedChiSquare handles both, chisq and mlh
-  fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedchisqPerRun);
+  Bool_t ok;
+  PDoubleVector par_r = ParamRound(param, err, ok);
+  if (!ok)
+    par_r = param;
+  fFitterFcn->CalcExpectedChiSquare(par_r, totalExpectedChisq, expectedchisqPerRun);
 
   // calculate chisq per run
   std::vector<Double_t> chisqPerRun;
@@ -2905,6 +2912,46 @@ Double_t PFitter::MilliTime()
   return ((Double_t)now.tv_sec * 1.0e6 + (Double_t)now.tv_usec)/1.0e3;
 }
 
+//--------------------------------------------------------------------------
+// ParamRound (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Rounds the parameter vector value according to the given error estimate,
+ * so that the msr-file value and the fitter result are consistent with each
+ * other. This means that musrfit -c, and musrfit -e -t should give essentially
+ * the same values of expected chisq (up to small rounding values).
+ * 
+ * @param par parameter value vector
+ * @param err error value vector
+ * @param ok true if size of par and err are identically, otherwise false.
+ * 
+ * @return rounded parameter value vector, compatible with the msr-file output.
+ */
+PDoubleVector PFitter::ParamRound(const PDoubleVector &par, const PDoubleVector &err, Bool_t &ok)
+{
+  PDoubleVector par_r;
+  
+  par_r.resize(par.size());
+  ok = true;
+  
+  if (par.size() != err.size()) {
+    // error msg
+    ok = false;
+    return par_r;
+  }
+  
+  int exp;
+  double dval;
+  for (unsigned int i=0; i<par.size(); i++) {
+    exp = round(log10(fabs(err[i])))-2;
+    dval = round(par[i]*pow(10.0, -exp))/pow(10.0, -exp);
+    par_r[i] = dval;
+  }
+  
+  return par_r;
+}
+
+
 //-------------------------------------------------------------------------------------------------
 // end
 //-------------------------------------------------------------------------------------------------

src/include/PFitter.h

@@ -187,6 +187,7 @@ class PFitter
     Bool_t ExecuteSector(std::ofstream &fout);
 
     Double_t MilliTime();
+    PDoubleVector ParamRound(const PDoubleVector &par, const PDoubleVector &err, Bool_t &ok);
 };
 
 #endif // _PFITTER_H_