feat: save histograms without re-binning (#165)

* save histograms with original binning as produced by coffea to root files
* re-bin histograms via cabinetry for statistical inference

Author: alexander-held

analyses/cms-open-data-ttbar/cabinetry_config.yml

@@ -9,11 +9,11 @@ Regions:
   - Name: "4j1b CR"
     RegionPath: "4j1b"
     Variable: "$H_T$ [GeV]"
-    Binning: [ 50.        ,  95.45454545, 140.90909091, 186.36363636, 231.81818182, 277.27272727, 322.72727273, 368.18181818,413.63636364, 459.09090909, 504.54545455, 550.        ]
+    Binning: [110, 150, 190, 230, 270, 310, 350, 390, 430, 470, 510, 550]
   - Name: "4j2b SR"
     RegionPath: "4j2b"
     Variable: "$m_{bjj}$ [GeV]"
-    Binning: [ 50.        ,  95.45454545, 140.90909091, 186.36363636, 231.81818182, 277.27272727, 322.72727273, 368.18181818,413.63636364, 459.09090909, 504.54545455, 550.        ]
+    Binning: [110, 150, 190, 230, 270, 310, 350, 390, 430, 470, 510, 550]
 
 Samples:
   - Name: "Pseudodata"

analyses/cms-open-data-ttbar/ttbar_analysis_pipeline.ipynb

@@ -503,12 +503,19 @@ def get_query(source: ObjectStream) -> ObjectStream:
 # %% [markdown]
 # ### Statistical inference
 #
+# We are going to perform a re-binning for the statistical inference.
+# This is planned to be conveniently provided via cabinetry (see [cabinetry#412](https://github.com/scikit-hep/cabinetry/issues/412), but in the meantime we can achieve this via [template building overrides](https://cabinetry.readthedocs.io/en/latest/advanced.html#overrides-for-template-building).
+# The implementation is provided in a function in `utils/`.
+#
 # A statistical model has been defined in `config.yml`, ready to be used with our output.
 # We will use `cabinetry` to combine all histograms into a `pyhf` workspace and fit the resulting statistical model to the pseudodata we built.
 
 # %% tags=[]
 config = cabinetry.configuration.load("cabinetry_config.yml")
-cabinetry.templates.collect(config)
+
+# rebinning: lower edge 110 GeV, merge bins 2->1
+rebinning_router = utils.get_cabinetry_rebinning_router(config, rebinning=slice(110j, None, hist.rebin(2)))
+cabinetry.templates.build(config, router=rebinning_router)
 cabinetry.templates.postprocess(config)  # optional post-processing (e.g. smoothing)
 ws = cabinetry.workspace.build(config)
 cabinetry.workspace.save(ws, "workspace.json")

analyses/cms-open-data-ttbar/ttbar_analysis_pipeline.py

@@ -1,6 +1,8 @@
 import asyncio
 import json
 
+import cabinetry
+from cabinetry.contrib import histogram_reader
 import hist
 import matplotlib as mpl
 import matplotlib.pyplot as plt
@@ -97,31 +99,31 @@ def save_histograms(all_histograms, fileset, filename):
 
     with uproot.recreate(filename) as f:
         for region in ["4j1b", "4j2b"]:
-            f[f"{region}_pseudodata"] = pseudo_data[120j::hist.rebin(2), region]
+            f[f"{region}_pseudodata"] = pseudo_data[:, region]
             for sample in nominal_samples:
                 sample_name = sample.split("__")[0]
-                f[f"{region}_{sample_name}"] = all_histograms[120j::hist.rebin(2), region, sample_name, "nominal"]
+                f[f"{region}_{sample_name}"] = all_histograms[:, region, sample_name, "nominal"]
 
                 # b-tagging variations
                 for i in range(4):
                     for direction in ["up", "down"]:
                         variation_name = f"btag_var_{i}_{direction}"
-                        f[f"{region}_{sample_name}_{variation_name}"] = all_histograms[120j::hist.rebin(2), region, sample_name, variation_name]
+                        f[f"{region}_{sample_name}_{variation_name}"] = all_histograms[:, region, sample_name, variation_name]
 
                 # jet energy scale variations
                 for variation_name in ["pt_scale_up", "pt_res_up"]:
-                    f[f"{region}_{sample_name}_{variation_name}"] = all_histograms[120j::hist.rebin(2), region, sample_name, variation_name]
+                    f[f"{region}_{sample_name}_{variation_name}"] = all_histograms[:, region, sample_name, variation_name]
 
             # ttbar modeling
-            f[f"{region}_ttbar_ME_var"] = all_histograms[120j::hist.rebin(2), region, "ttbar", "ME_var"]
-            f[f"{region}_ttbar_PS_var"] = all_histograms[120j::hist.rebin(2), region, "ttbar", "PS_var"]
+            f[f"{region}_ttbar_ME_var"] = all_histograms[:, region, "ttbar", "ME_var"]
+            f[f"{region}_ttbar_PS_var"] = all_histograms[:, region, "ttbar", "PS_var"]
 
-            f[f"{region}_ttbar_scaledown"] = all_histograms[120j :: hist.rebin(2), region, "ttbar", "scaledown"]
-            f[f"{region}_ttbar_scaleup"] = all_histograms[120j :: hist.rebin(2), region, "ttbar", "scaleup"]
+            f[f"{region}_ttbar_scaledown"] = all_histograms[:, region, "ttbar", "scaledown"]
+            f[f"{region}_ttbar_scaleup"] = all_histograms[:, region, "ttbar", "scaleup"]
 
             # W+jets scale
-            f[f"{region}_wjets_scale_var_down"] = all_histograms[120j :: hist.rebin(2), region, "wjets", "scale_var_down"]
-            f[f"{region}_wjets_scale_var_up"] = all_histograms[120j :: hist.rebin(2), region, "wjets", "scale_var_up"]
+            f[f"{region}_wjets_scale_var_down"] = all_histograms[:, region, "wjets", "scale_var_down"]
+            f[f"{region}_wjets_scale_var_up"] = all_histograms[:, region, "wjets", "scale_var_up"]
 
 
 class ServiceXDatasetGroup():
@@ -153,3 +155,25 @@ def get_data_rootfiles_uri(self, query, as_signed_url=True, title="Untitled"):
             files_per_process.update({process: all_files[parent_key[self.filelist[:,1]==process]]})
 
         return files_per_process
+
+
+def get_cabinetry_rebinning_router(config, rebinning):
+    # perform re-binning in cabinetry by providing a custom function reading histograms
+    # will eventually be replaced via https://github.com/scikit-hep/cabinetry/issues/412
+    rebinning_router = cabinetry.route.Router()
+
+    # this reimplements some of cabinetry.templates.collect
+    general_path = config["General"]["InputPath"]
+    variation_path = config["General"].get("VariationPath", None)
+
+    # define a custom template builder function that is executed for data samples
+    @rebinning_router.register_template_builder()
+    def build_data_hist(region, sample, systematic, template):
+        # get path to histogram
+        histo_path = cabinetry.templates.collector._histo_path(general_path, variation_path, region, sample, systematic, template)
+        h = hist.Hist(histogram_reader.with_uproot(histo_path))  # turn from boost-histogram into hist
+        # perform re-binning
+        h = h[rebinning]
+        return h
+
+    return rebinning_router