Handle ROOT TH1 with bin width !=1

Signed-off-by: Jan Mayer <jan.mayer@ikp.uni-koeln.de>

hdtv/histogram.py

@@ -27,11 +27,21 @@ import hdtv.rootext.mfile
 
 from hdtv.drawable import Drawable
 from hdtv.specreader import SpecReader, SpecReaderError
+from hdtv.cal import CalibrationFitter
 
 # Don't add created spectra to the ROOT directory
 ROOT.TH1.AddDirectory(ROOT.kFALSE)
 
 
+def HasPrimitiveBinning(hist):
+    if hist.GetNbinsX() != (hist.GetXaxis().GetXmax() - hist.GetXaxis().GetXmin()):
+        return False
+    for bin in range(0, hist.GetNbinsX()):
+        if hist.GetBinWidth(bin) != 1.:
+            return False
+    return True
+
+
 class Histogram(Drawable):
     """
     Histogram object
@@ -43,11 +53,24 @@ class Histogram(Drawable):
 
     def __init__(self, hist, color=hdtv.color.default, cal=None):
         Drawable.__init__(self, color, cal)
-        self._hist = hist
         self._norm = 1.0
         self._ID = None
         self.effCal = None
         self.typeStr = "spectrum"
+        self.cal = cal
+        if HasPrimitiveBinning(hist):
+            self._hist = hist
+        else:
+            hdtv.ui.info(hist.GetName() + " unconventional binning detected. Converting and creating calibration ...")
+            self._hist = ROOT.TH1D(hist.GetName(), hist.GetTitle(), hist.GetNbinsX(), 0, hist.GetNbinsX())
+            cf = CalibrationFitter()
+            # TODO: Slow
+            for bin in range(0, hist.GetNbinsX()):
+                cf.AddPair(bin, hist.GetXaxis().GetBinUpEdge(bin))
+                self._hist.SetBinContent(bin, hist.GetBinContent(bin))
+                # TODO: Copy Errors?
+            cf.FitCal(4)
+            self.cal = cf.calib
 
     def __str__(self):
         return self.name

hdtv/plugins/rootInterface.py

@@ -412,12 +412,15 @@ class RootFileInterface(object):
                     sid = self.spectra.Insert(spec)
                     spec.color = hdtv.color.ColorForID(sid.major)
                     loaded.append(sid)
-                    if args.load_cal:
-                        if spec.name in self.spectra.caldict:
-                            spec.cal = self.caldict[spec.name]
-                        else:
-                            print("Warning: no calibration found for %s" %
-                                  spec.name)
+                    if spec.cal:
+                        self.caldict[spec.name] = spec.cal
+                    else:
+                        if args.load_cal:
+                            if spec.name in self.spectra.caldict:
+                                spec.cal = self.caldict[spec.name]
+                            else:
+                                print("Warning: no calibration found for %s" %
+                                      spec.name)
                 else:
                     hdtv.ui.warn("%s is not a 1D histogram object" %
                                  obj.GetName())

test/fitxml/test_cal.py

@@ -147,6 +147,27 @@ def test_fit_in_calibrated_spectrum(region1, region2, peak, calfactor):
     markers_consistent(calfactor, region1, region2)
 
 
+# Check that bin sizes are handled in fitted peak volume (see also test_root_fit_volume_with_binning)
+@pytest.mark.parametrize("spectrum, calfactor", [("h", 1), ("h2", 0.5)])
+@pytest.mark.parametrize("region1, region2, peak, expected_volume", [(500., 1500., 1000., 1000000.)])
+def test_mfile_fit_volume_with_binning(region1, region2, peak, expected_volume, spectrum, calfactor):
+    spectra.Clear()
+
+    __main__.s.LoadSpectra(os.path.join("test", "share", "binning_root_" + spectrum + ".spc"))
+    spectra.ApplyCalibration("0", [0, calfactor])
+
+    spectra.SetMarker("region", region1)
+    spectra.SetMarker("region", region2)
+    spectra.SetMarker("peak", peak)
+    spectra.ExecuteFit()
+    fit = spectra.workFit
+    fit_volume = fit.ExtractParams()[0][0]["vol"].nominal_value
+    integral_volume= fit.ExtractIntegralParams()[0][0]["vol"].nominal_value
+
+    assert isclose(integral_volume, expected_volume, abs_tol=10)
+    assert isclose(fit_volume, expected_volume, abs_tol=200)
+
+
 @pytest.mark.parametrize("calfactor1", [
     2.1, 1.9])
 @pytest.mark.parametrize("calfactor2", [

test/plugins/test_rootInterface.py

@@ -22,31 +22,38 @@
 import os
 
 import pytest
+from test.helpers.utils import redirect_stdout, hdtvcmd, isclose, setup_io
 
-from test.helpers.utils import redirect_stdout, hdtvcmd
-
-import hdtv.cmdline
-import hdtv.options
 import hdtv.session
 
 import __main__
+
 # We don’t want to see the GUI. Can we prevent this?
 try:
     __main__.spectra = hdtv.session.Session()
 except RuntimeError:
     pass
+session = __main__.spectra
 
-#__main__.spectra.window.viewer.CloseWindow()
-
+import hdtv.plugins.specInterface
+import hdtv.plugins.fitInterface
 import hdtv.plugins.rootInterface
 
+import hdtv.cmdline
+import hdtv.options
+import hdtv.rfile_utils
+
+
 @pytest.fixture(autouse=True)
 def prepare(request):
-    #original_wd = os.path.abspath(os.path.join(__file__, os.pardir))
+    __main__.f.ResetFitterParameters()
+    # original_wd = os.path.abspath(os.path.join(__file__, os.pardir))
     original_wd = os.getcwd()
     os.chdir(original_wd)
     yield
     os.chdir(original_wd)
+    session.Clear()
+
 
 def test_cmd_root_pwd():
     f, ferr = hdtvcmd("root pwd")
@@ -68,3 +75,50 @@ def test_cmd_root_browse():
     hdtv.plugins.rootInterface.r.browser.SetName('Test')
     assert hdtv.plugins.rootInterface.r.browser.GetName() == 'Test'
     hdtv.plugins.rootInterface.r.browser = None
+
+
+# Check that bin sizes are handled in fitted peak volume (see also test_mfile_fit_volume_with_binning)
+@pytest.mark.parametrize("spectrum", ["h", "h2"])
+@pytest.mark.parametrize("region1, region2, peak, expected_volume", [(500., 1500., 1000., 1000000.)])
+def test_root_fit_volume_with_binning(region1, region2, peak, expected_volume, spectrum):
+    # Mock parser so RootGet can be used
+    # args = mock.Mock()
+    args = type('Test', (object,), {})
+    args.pattern = [os.path.join("test", "share", "binning.root", spectrum)]
+    args.replace = False
+    args.load_cal = False
+    args.invisible = False
+    hdtv.plugins.rootInterface.r.RootGet(args)
+
+    session.SetMarker("region", region1)
+    session.SetMarker("region", region2)
+    session.SetMarker("peak", peak)
+    session.ExecuteFit()
+    fit = session.workFit
+    fit_volume = fit.ExtractParams()[0][0]["vol"].nominal_value
+    integral_volume = fit.ExtractIntegralParams()[0][0]["vol"].nominal_value
+
+    assert isclose(integral_volume, expected_volume, abs_tol=10)
+    assert isclose(fit_volume, expected_volume, abs_tol=200)
+    hdtv.plugins.rootInterface.r.RootClose(None)
+
+
+def test_root_to_root_conversion_for_unconventional_binning():
+    # Mock parser so RootGet can be used
+    # args = mock.Mock()
+    args = type('Test', (object,), {})
+    args.replace = False
+    args.load_cal = False
+    args.invisible = False
+    args.pattern = [os.path.join("test", "share", "binning.root", "h")]
+    hdtv.plugins.rootInterface.r.RootGet(args)
+    args.pattern = [os.path.join("test", "share", "binning.root", "h2")]
+    hdtv.plugins.rootInterface.r.RootGet(args)
+
+    assert not get_spec(0).cal
+    assert isclose(get_spec(1).cal.GetCoeffs()[1], 0.5)
+
+
+def get_spec(specid):
+    return __main__.s.spectra.dict.get(
+        [x for x in list(__main__.s.spectra.dict) if x.major == specid][0])