From 9e8d53463581fb917bc506d6de73d561fc01b59b Mon Sep 17 00:00:00 2001
From: Andrew Zhen Ning Yong <andrew.yong@ed.ac.uk>
Date: Thu, 20 May 2021 11:04:21 +0100
Subject: [PATCH] New methods in XYSampleData class to skip fits that are
 non-converging/large chi2PerDof. This is tracked by boolean goodFit_.

---
 lib/Statistics/XYSampleData.cpp | 55 ++++++++++++++++++++++++---------
 lib/Statistics/XYSampleData.hpp |  3 ++
 2 files changed, 43 insertions(+), 15 deletions(-)

diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp
index 6e2d2d6..19342c9 100644
--- a/lib/Statistics/XYSampleData.cpp
+++ b/lib/Statistics/XYSampleData.cpp
@@ -234,6 +234,21 @@ DVec XYSampleData::getYError(const Index j)
     return data_.getYError(j);
 }
 
+bool XYSampleData::checkFit()
+{
+    return goodFit_;
+}
+
+void XYSampleData::checkChi2PerDof(double Chi2PerDof)
+{
+    if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) 
+    {
+        goodFit_ = false;
+        cerr << "chi2PerDof = " << Chi2PerDof << ". Aborting fit now." << endl;
+    }
+
+}
+
 // get total fit variance matrix and its pseudo-inverse ////////////////////////
 const DMat & XYSampleData::getFitVarMat(void)
 {
@@ -292,24 +307,34 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
     result.resize(nSample_);
     result.chi2_.resize(nSample_);
     result.model_.resize(v.size());
+    double chi2PerDof;
+    goodFit_ = true;
     FOR_STAT_ARRAY(result, s)
     {
-        setDataToSample(s);
-        if (s == central)
+        if(goodFit_)
         {
-            sampleResult = data_.fit(minimizer, initCopy, v);
-            initCopy     = sampleResult.segment(0, initCopy.size());
-        }
-        else
-        {
-            sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
-        }
-        result[s]       = sampleResult;
-        result.chi2_[s] = sampleResult.getChi2();
-        for (unsigned int j = 0; j < v.size(); ++j)
-        {
-            result.model_[j].resize(nSample_);
-            result.model_[j][s] = sampleResult.getModel(j);
+            setDataToSample(s);
+            if (s == central)
+            {
+                sampleResult = data_.fit(minimizer, initCopy, v);
+                initCopy     = sampleResult.segment(0, initCopy.size());
+                chi2PerDof   = sampleResult.getChi2PerDof();
+                checkChi2PerDof(chi2PerDof);
+            }
+            else
+            {
+                sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
+                chi2PerDof   = sampleResult.getChi2PerDof();
+                checkChi2PerDof(chi2PerDof);
+            }
+            result[s]       = sampleResult;
+            result.chi2_[s] = sampleResult.getChi2();
+            for (unsigned int j = 0; j < v.size(); ++j)
+            {
+                result.model_[j].resize(nSample_);
+                result.model_[j][s] = sampleResult.getModel(j);
+
+            }
         }
     }
     result.nPar_    = sampleResult.getNPar();
diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp
index 6cba855..a66c9dd 100644
--- a/lib/Statistics/XYSampleData.hpp
+++ b/lib/Statistics/XYSampleData.hpp
@@ -91,6 +91,8 @@ public:
     const DMat &       getXYVar(const Index i, const Index j);
     DVec               getXError(const Index i);
     DVec               getYError(const Index j);
+    bool               checkFit(); // check fit candidate based on chi2PerDof
+    void               checkChi2PerDof(double Chi2PerDof);
     // get total fit variance matrix and its pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
@@ -133,6 +135,7 @@ private:
     Index                                 nSample_, dataSample_{central};
     bool                                  initData_{true}, computeVarMat_{true};
     bool                                  initXMap_{true};
+    bool                                  goodFit_{true}; // used to break minimisation if central sample chi2PerDof is bad
 };
 
 /******************************************************************************