fix: handle bad ctof_tdcadc fits (#223)

JeffersonLab · Aug 15, 2024 · 4f0ac7d · 4f0ac7d
1 parent f9a2487
commit 4f0ac7d
Showing 1 changed file with 13 additions and 10 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
@@ -37,22 +37,25 @@ class CTOFFitter {
 
   static def tdcadcdifffit(H1F h1) {
 
+    // set min/max allowed values
+    def valOrMinAllowed = { val -> [val, h1.getAxis().min()].max() }
+    def valOrMaxAllowed = { val -> [val, h1.getAxis().max()].min() }
+
     def fit_peak = { peakbin, prefix ->
       def f1 = new F1D("${prefix}:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", -5.0, 5.0)
       double hAmp  = h1.getBinContent(peakbin);
       double hMean = h1.getAxis().getBinCenter(peakbin)
       double hRMS  = h1.getRMS();
       double factor1 = 3.0
       double factor2 = 1.57
-      double rangeMin = (hMean - factor1*0.5);
-      double rangeMax = (hMean + factor2*0.5);
+      double rangeMin = valOrMinAllowed(hMean - factor1*0.5);
+      double rangeMax = valOrMaxAllowed(hMean + factor2*0.5);
       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, 1);
-      // f1.setParLimits(2, 0.5*hRMS, 1.5*hRMS);
+      f1.setParLimits(1, rangeMin, rangeMax);
+      f1.setParLimits(2, 0, [5*hRMS, h1.getAxis().max()-h1.getAxis().min()].min())
 
       def makefit = {func->
         hMean = func.getParameter(1)
@@ -70,18 +73,18 @@ class CTOFFitter {
     }
 
     // find the highest peak, and fit it
-    def bins = (1..h1.getAxis().getNBins())
+    def bins = (1..h1.getAxis().getNBins()-1)
     def peakbin1 = bins.max{h1.getBinContent(it)}
     def func1 = fit_peak(peakbin1,'fit1')
 
     // find the 2nd highest peak by excluding the region around the 1st highest
     // peak, and searching for the new highest max, then fit it
-    def peak1Start = func1.getParameter(1) - 2*func1.getParameter(2)
-    def peak1End   = func1.getParameter(1) + 2*func1.getParameter(2)
-    def peakbin2 = bins
+    def peak1Start = valOrMinAllowed(func1.getParameter(1) - 2*func1.getParameter(2))
+    def peak1End   = valOrMaxAllowed(func1.getParameter(1) + 2*func1.getParameter(2))
+    def dataWithoutPeak1 = bins
       .collect{ [ it, h1.getBinContent(it) ] }
       .findAll{ h1.getAxis().getBinCenter(it[0]) < peak1Start || h1.getAxis().getBinCenter(it[0]) > peak1End }
-      .max{ it[1] }[0]
+    def peakbin2 = dataWithoutPeak1.size() > 0 ? dataWithoutPeak1.max{it[1]}[0] : peakbin1 // if peak1 fit was bad, just repeat the fit for peak 2 and leave it bad
     def func2 = fit_peak(peakbin2,'fit2')
 
     // decide which fit result is upstream and downstream