From 22d89c28b4a7bae4ed33d99bbef22754c17e2431 Mon Sep 17 00:00:00 2001
From: Christopher Dilks <c-dilks@users.noreply.github.com>
Date: Wed, 1 May 2024 19:25:20 -0400
Subject: [PATCH] fix: use IQR for initial `sigma` for HTCC vtimediff fit
 stabilization (#205)

---
 .../clas/timeline/fitter/HTCCFitter.groovy    | 27 ++++++++++++++++++-
 1 file changed, 26 insertions(+), 1 deletion(-)

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
index 38296f1f..e94f3aa3 100644
--- a/detectors/src/main/java/org/jlab/clas/timeline/fitter/HTCCFitter.groovy
+++ b/detectors/src/main/java/org/jlab/clas/timeline/fitter/HTCCFitter.groovy
@@ -14,6 +14,30 @@ import org.jlab.groot.math.F1D
 class HTCCFitter{
   static F1D timeIndPMT(H1F h1) {
     double maxV = h1.getDataX(h1.getMaximumBin());
+
+    // 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 f1 = new F1D("fit:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", maxV-2, maxV+3);
 
     f1.setLineColor(2);
@@ -22,7 +46,8 @@ class HTCCFitter{
 
     f1.setParameter(0, h1.getMax());
     f1.setParameter(1, maxV);
-    f1.setParameter(2, 0.7);
+    f1.setParameter(2, Math.max(iqr / 2.0, 0.1));
+    f1.setParLimits(2, Math.max(iqr / 5.0, 0.05), Math.min(iqr * 5.0, h1.getAxis().getBinCenter(h1.getAxis().getNBins()-1)));
 
     DataFitter.fit(f1, h1, "");