fix: undershooting Gaussian fits of ctof_time (#206)

JeffersonLab · May 29, 2024 · 100cfbe · 100cfbe
1 parent 22d89c2
commit 100cfbe
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 34 deletions.
diff --git a/detectors/src/main/java/org/jlab/clas/timeline/fitter/CTOFFitter.groovy b/detectors/src/main/java/org/jlab/clas/timeline/fitter/CTOFFitter.groovy
@@ -2,29 +2,27 @@ package org.jlab.clas.timeline.fitter
 import org.jlab.groot.fitter.DataFitter
 import org.jlab.groot.data.H1F
 import org.jlab.groot.math.F1D
+import org.jlab.clas.timeline.util.HistoUtil
 
 
 class CTOFFitter {
   static F1D timefit(H1F h1) {
     def f1 = new F1D("fit:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", -5.0, 5.0)
-    double hAmp  = h1.getBinContent(h1.getMaximumBin());
-    double hMean = h1.getAxis().getBinCenter(h1.getMaximumBin())
-    double hRMS  = Math.min(h1.getRMS(), 0.2)
-    double rangeMin = (hMean - 2.5*hRMS);
-    double rangeMax = (hMean + 2.5*hRMS);
+    double hAmp   = h1.getBinContent(h1.getMaximumBin());
+    double hMu    = h1.getAxis().getBinCenter(h1.getMaximumBin())
+    double hSigma = Math.min(HistoUtil.getHistoIQR(h1) / 2.0, 0.2)
+    def rangeFactor = 1.5
+    double rangeMin = (hMu - rangeFactor*hSigma);
+    double rangeMax = (hMu + rangeFactor*hSigma);
     f1.setRange(rangeMin, rangeMax);
     f1.setParameter(0, hAmp);
-    // f1.setParLimits(0, hAmp*0.8, hAmp*1.2);
-    f1.setParameter(1, hMean);
-    // f1.setParLimits(1, 0, 0.5);
-    f1.setParameter(2, hRMS);
-    // f1.setParLimits(2, 0.5*hRMS, 1.5*hRMS);
-
+    f1.setParameter(1, hMu);
+    f1.setParameter(2, hSigma);
 
     def makefit = {func->
-      hMean = func.getParameter(1)
-      hRMS = func.getParameter(2).abs()
-      func.setRange(hMean-2*hRMS,hMean+2*hRMS)
+      def fMu = func.getParameter(1)
+      def fSigma = func.getParameter(2).abs()
+      func.setRange(fMu-rangeFactor*fSigma, fMu+rangeFactor*fSigma)
       DataFitter.fit(func,h1,"Q")
       return [func.getChiSquare(), (0..<func.getNPars()).collect{func.getParameter(it)}]
     }

diff --git a/detectors/src/main/java/org/jlab/clas/timeline/fitter/HTCCFitter.groovy b/detectors/src/main/java/org/jlab/clas/timeline/fitter/HTCCFitter.groovy
@@ -9,6 +9,7 @@ package org.jlab.clas.timeline.fitter
 import org.jlab.groot.fitter.DataFitter
 import org.jlab.groot.data.H1F
 import org.jlab.groot.math.F1D
+import org.jlab.clas.timeline.util.HistoUtil
 
 
 class HTCCFitter{
@@ -17,26 +18,7 @@ class HTCCFitter{
 
     // compute the IQR of `h1`, which is a better estimate of the initial sigma parameter than
     // the RMS, since it is not (typically) biased by outliers
-    def l1 = []
-    h1.getAxis().getNBins().times { bin ->
-      def counts = h1.getBinContent(bin).toInteger() // FIXME: assumes the histogram is unweighted
-      def value  = h1.getAxis().getBinCenter(bin)
-      counts.times { l1 += value }
-    }
-    def listMedian = { d ->
-      if(d.size()==0) {
-        System.err.println("WARNING: attempt to calculate median of an empty list")
-        return 0
-      }
-      d.sort()
-      def m = d.size().intdiv(2)
-      d.size() % 2 ? d[m] : (d[m-1]+d[m]) / 2
-    }
-    def mq = listMedian(l1)
-    def lq = listMedian(l1.findAll{it<mq})
-    def uq = listMedian(l1.findAll{it>mq})
-    def iqr = uq - lq
-
+    def iqr = HistoUtil.getHistoIQR(h1)
 
     def f1 = new F1D("fit:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", maxV-2, maxV+3);
 

diff --git a/detectors/src/main/java/org/jlab/clas/timeline/util/HistoUtil.groovy b/detectors/src/main/java/org/jlab/clas/timeline/util/HistoUtil.groovy
@@ -57,4 +57,28 @@ class HistoUtil {
     nBinsIn.times{ histOut.fill(histIn.getXaxis().getBinCenter(it), histIn.getBinContent(it)) }
     return histOut
   }
+
+  /// @returns the Interquartile Range (IQR) of a histogram
+  /// @param the histogram; FIXME: only unweighted histograms are supported at the moment
+  static def getHistoIQR(H1F hist) {
+    def hist_list = []
+    hist.getAxis().getNBins().times { bin ->
+      def counts = hist.getBinContent(bin).toInteger() // FIXME: assumes the histogram is unweighted
+      def value  = hist.getAxis().getBinCenter(bin)
+      counts.times { hist_list += value }
+    }
+    def listMedian = { d ->
+      if(d.size()==0) {
+        System.err.println("WARNING in HistoUtil.getHistoIQR: attempt to calculate median of an empty list")
+        return 0
+      }
+      d.sort()
+      def m = d.size().intdiv(2)
+      d.size() % 2 ? d[m] : (d[m-1]+d[m]) / 2
+    }
+    def mq = listMedian(hist_list)
+    def lq = listMedian(hist_list.findAll{it<mq})
+    def uq = listMedian(hist_list.findAll{it>mq})
+    return uq - lq
+  }
 }