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 index 64cce3b8..fa46461f 100644 --- 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.. - 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{itmq}) - 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 index bf77c728..ce5c5534 100644 --- 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{itmq}) + return uq - lq + } }