Feature/qcd hook (#35)

* init new branch

* update submodules

* adding files

* fix linting

* updating inf_model

* adding limited dataset option

* adding hh producer

* deleting old qcd estimation

* test: add 2j cat and shape uncertanties

* adding nested dict for up/down shifts

* fixing uncertainties propagation

* under/overflow fix

* Linting.

* Remove test versions.

---------

Co-authored-by: Ana A <[email protected]>

hbt/config/categories.py

@@ -25,21 +25,8 @@ def add_categories(config: od.Config) -> None:
     add_category(config, name="noniso", id=13, selection="cat_noniso", label=r"$\tau_{h,2} non-isolated$", tags={"noniso"})  # noqa: E501
 
     # kinematic categories
-    add_category(
-        config,
-        name="incl",
-        id=100,
-        selection="cat_incl",
-        label="inclusive",
-    )
-
-    add_category(
-        config,
-        name="2j",
-        id=110,
-        selection="cat_2j",
-        label="2 jets",
-    )
+    add_category(config, name="incl", id=100, selection="cat_incl", label="inclusive")
+    add_category(config, name="2j", id=110, selection="cat_2j", label="2 jets")
 
     #
     # build groups

hbt/config/hist_hooks.py

@@ -141,21 +141,36 @@ def qcd_estimation(task, hists):
 
             # ABCD method
             # shape: (SHIFT, VAR)
-            os_iso_qcd = os_noniso_qcd * (int_ss_iso / int_ss_noniso)
+            os_iso_qcd = os_noniso_qcd * ((int_ss_iso / int_ss_noniso)[:, None])
 
             # combine uncertainties and store values in bare arrays
             os_iso_qcd_values = os_iso_qcd()
             os_iso_qcd_variances = os_iso_qcd(sn.UP, sn.ALL, unc=True)**2
 
+            # define uncertainties
+            unc_data = os_iso_qcd(sn.UP, ["os_noniso_data", "ss_iso_data", "ss_noniso_data"], unc=True)
+            unc_mc = os_iso_qcd(sn.UP, ["os_noniso_mc", "ss_iso_mc", "ss_noniso_mc"], unc=True)
+            unc_data_rel = abs(unc_data / os_iso_qcd_values)
+            unc_mc_rel = abs(unc_mc / os_iso_qcd_values)
+
+            # only keep the MC uncertainty if it is larger than the data uncertainty and larger than 15%
+            keep_variance_mask = (
+                np.isfinite(unc_mc_rel) &
+                (unc_mc_rel > unc_data_rel) &
+                (unc_mc_rel > 0.15)
+            )
+            os_iso_qcd_variances[keep_variance_mask] = unc_mc[keep_variance_mask]**2
+            os_iso_qcd_variances[~keep_variance_mask] = 0
+
             # retro-actively set values to zero for shifts that had negative integrals
             neg_int_mask = int_ss_iso_neg | int_ss_noniso_neg
             os_iso_qcd_values[neg_int_mask] = 1e-5
-            os_iso_qcd_variances[neg_int_mask] = 1e-5
+            os_iso_qcd_variances[neg_int_mask] = 0
 
             # residual zero filling
             zero_mask = os_iso_qcd_values <= 0
             os_iso_qcd_values[zero_mask] = 1e-5
-            os_iso_qcd_variances[zero_mask] = 1e-5
+            os_iso_qcd_variances[zero_mask] = 0
 
             # insert values into the qcd histogram
             cat_axis = qcd_hist.axes["category"]

hbt/inference/inf_model.py

@@ -0,0 +1,89 @@
+from columnflow.inference import inference_model, ParameterType, ParameterTransformation
+
+
+@inference_model
+def test_unc(self):
+    self.add_category(
+        "incl",
+        config_category="incl",
+        config_variable="hh_mass",
+        # fake data
+        data_from_processes=["TT", "dy"],
+        mc_stats=True,
+    )
+
+    self.add_category(
+        "2j",
+        config_category="2j",
+        config_variable="jet1_pt",
+        # fake data
+        data_from_processes=["TT", "dy"],
+        mc_stats=True,
+    )
+
+    #
+    # processes
+    #
+    self.add_process(
+        "dy",
+        config_process="dy",
+    )
+
+    self.add_process(
+        "TT",
+        config_process="tt_sl",
+    )
+
+    self.add_process(
+        "hh_ggf",
+        is_signal=True,
+        config_process="hh_ggf_hbb_htt_kl1_kt1",
+    )
+
+    #
+    # parameters
+    #
+
+    # groups
+    self.add_parameter_group("experiment")
+    self.add_parameter_group("theory")
+
+    # lumi
+    lumi = self.config_inst.x.luminosity
+    for unc_name in lumi.uncertainties:
+        self.add_parameter(
+            unc_name,
+            type=ParameterType.rate_gauss,
+            effect=lumi.get(names=unc_name, direction=("down", "up"), factor=True),
+        )
+        self.add_parameter_to_group(unc_name, "experiment")
+
+    # electron uncertainty
+    self.add_parameter(
+        "CMS_eff_e",  # this is the name of the uncertainty as it will show in the datacard. Let's use some variant of the official naming # noqa
+        process="*",
+        type=ParameterType.shape,
+        config_shift_source="e",  # this is the name of the shift (alias) in the config
+    )
+    self.add_parameter_to_group("CMS_eff_e", "experiment")
+
+    # a custom asymmetric uncertainty
+    self.add_parameter(
+        "QCDscale_ttbar",
+        process="TT",
+        type=ParameterType.shape,
+        transformations=[ParameterTransformation.effect_from_rate],
+        effect=(0.85, 1.1),
+    )
+    self.add_parameter_to_group("QCDscale_ttbar", "experiment")
+
+    """
+    # tune uncertainty
+    self.add_parameter(
+        "tune",
+        process="TT",
+        type=ParameterType.shape,
+        config_shift_source="tune",
+    )
+    self.add_parameter_to_group("tune", "experiment")
+    """

hbt/production/default.py

@@ -19,7 +19,6 @@
 from hbt.production.tau import tau_weights, trigger_weights
 from hbt.util import IF_DATASET_HAS_LHE_WEIGHTS
 
-
 ak = maybe_import("awkward")
 
 

law.cfg

@@ -47,7 +47,7 @@ production_modules: columnflow.production.{categories,normalization,processes},
 categorization_modules: hbt.categorization.default
 weight_production_modules: columnflow.weight.{empty,all_weights}, hbt.weight.default
 ml_modules: hbt.ml.test
-inference_modules: hbt.inference.{test,dl,base,derived,grav_model,test_unc}
+inference_modules: hbt.inference.{test,dl,base,derived,grav_model,inf_model}
 
 # whether MergeReducedEvents should keep its inputs from ReduceEvents by default
 # (otherwise they are removed after merging)