adding QCD estimation

uhh-cms · Jul 15, 2024 · 0470447 · 0470447
1 parent 27d3cbf
commit 0470447
Showing 1 changed file with 24 additions and 15 deletions.
diff --git a/hbt/config/hist_hooks.py b/hbt/config/hist_hooks.py
@@ -1,29 +1,38 @@
 import order as od
+import numpy as np
 import hist
 
 
-def create_regions(cat_id, regions, mc_hist, data_hist, cat_name, cat_parent):
+def create_regions(index, cat_id, regions, mc_hist, data_hist, cat_name, cat_parent):
     if "os" in cat_name:
         if "noniso" not in cat_name:
             regions[cat_parent]["A"] = {}
             regions[cat_parent]["A"]["mc"] = mc_hist[{"category": hist.loc(cat_id)}]
             regions[cat_parent]["A"]["data"] = data_hist[{"category": hist.loc(cat_id)}]
+            # regions[cat_parent]["A"]["mc"] = mc_hist.counts()[index]
+            # regions[cat_parent]["A"]["data"] = data_hist.counts()[index]
             print("... created region A for category " + cat_parent)
         else:
             regions[cat_parent]["B"] = {}
             regions[cat_parent]["B"]["mc"] = mc_hist[{"category": hist.loc(cat_id)}]
             regions[cat_parent]["B"]["data"] = data_hist[{"category": hist.loc(cat_id)}]
+            # regions[cat_parent]["B"]["mc"] = mc_hist.counts()[index]
+            # regions[cat_parent]["B"]["data"] = data_hist.counts()[index]
             print("... created region B for category " + cat_parent)
     elif "ss" in cat_name:
         if "noniso" not in cat_name:
             regions[cat_parent]["C"] = {}
             regions[cat_parent]["C"]["mc"] = mc_hist[{"category": hist.loc(cat_id)}]
             regions[cat_parent]["C"]["data"] = data_hist[{"category": hist.loc(cat_id)}]
+            # regions[cat_parent]["C"]["mc"] = mc_hist.counts()[index]
+            # regions[cat_parent]["C"]["data"] = data_hist.counts()[index]
             print("... created region C for category " + cat_parent)
         else:
             regions[cat_parent]["D"] = {}
             regions[cat_parent]["D"]["mc"] = mc_hist[{"category": hist.loc(cat_id)}]
             regions[cat_parent]["D"]["data"] = data_hist[{"category": hist.loc(cat_id)}]
+            # regions[cat_parent]["D"]["mc"] = mc_hist.counts()[index]
+            # regions[cat_parent]["D"]["data"] = data_hist.counts()[index]
             print("... created region D for category " + cat_parent)
 
 
@@ -52,28 +61,28 @@ def ana_ABCD(task, hists):
                 cat_id = mc_hist.axes[0].value(index)
                 cat_name = cfg.get_category(cat_id).name
                 cat_parent = cat_name.partition("_")[0]
-                print("... processing category " + cat_name)
+                print("... processing category " + str(cat_id) + " " + cat_name)
 
                 # define ABCD regions
                 if cat_parent not in regions:
                     regions[cat_parent] = {}
-                    create_regions(cat_id, regions, mc_hist, data_hist, cat_name, cat_parent)
+                    create_regions(index, cat_id, regions, mc_hist, data_hist, cat_name, cat_parent)
                 else:
-                    create_regions(cat_id, regions, mc_hist, data_hist, cat_name, cat_parent)
+                    create_regions(index, cat_id, regions, mc_hist, data_hist, cat_name, cat_parent)
 
-            from IPython import embed
-            embed()
-            quit()
+            # from IPython import embed
+            # embed()
+            # quit()
 
-            # calculate rate factor
-            #C_rate = regions[cat_parent]
+            # estimate QCD background
+            for cat_parent in regions.keys():
+                B_rate = regions[cat_parent]["B"]["data"] + ((-1) * regions[cat_parent]["B"]["mc"])
+                C_rate = regions[cat_parent]["C"]["data"] + ((-1) * regions[cat_parent]["C"]["mc"])
+                D_rate = regions[cat_parent]["D"]["data"] + ((-1) * regions[cat_parent]["D"]["mc"])
+                # per bin conversion factor
+                rate_factor = np.nan_to_num(C_rate.values() / D_rate.values())
 
-            #hists[cfg.processes.n.qcd] = B*(D/C)
-
-            # calculate rate factor
-            # C_rate =
-            # D_rate =
-            # rate_factor = C_rate / D_rate
+                hists[cfg.processes.n.qcd] = B_rate * rate_factor
 
             return hists