update to scan different values

check_yields.py

@@ -4,6 +4,7 @@ from collections import OrderedDict
 import matplotlib.pyplot as plt
 import mplhep as hep
 plt.style.use(hep.style.ROOT)
+import sympy
 
 from dhi.datacard_tools import read_datacard_structured
 from dhi.models.hh_model import VBFSample, VBFFormula
@@ -22,14 +23,17 @@ if __name__ == '__main__':
 
     skip_vbf = [(0.5, 1, 1)] # for model_default
     vbf_keys = [key for key in vbf_samples if skip_vbf is None or key not in skip_vbf]
+
     vbf_formula = VBFFormula([vbf_samples[key] for key in vbf_keys])
 
-    print('formula (plug in templates for xsi)', vbf_formula.sigma)
+    symbol_names = ["C2V", "CV", "kl"] + list(map("xs{}".format, range(vbf_formula.n_samples)))
+    xsec_func = sympy.lambdify(sympy.symbols(symbol_names), vbf_formula.sigma, 'numpy')
 
     datacard = '/afs/cern.ch/user/w/woodson/work/HH/datacards_run2/bbbb_boosted_vbf/v2/datacard.txt'
 
     content = read_datacard_structured(datacard)
-    all_vbf_names = [p["name"] for p in content["processes"] if "qqHH" in p["name"]]
+    #all_vbf_names = [p["name"] for p in content["processes"] if "qqHH" in p["name"]]
+    all_vbf_names = [sample.label+"_hbbhbb" for sample in vbf_formula.samples] # get sample names in same order as formula
     all_bin_names = [p["name"] for p in content["bins"]]
 
     shapes = OrderedDict()
@@ -42,6 +46,11 @@ if __name__ == '__main__':
             if vbf_name in rates:
                 shapes[vbf_name][bin_index] = rates[vbf_name]
 
+    def xsec_shape(C2V=1, CV=1, kl=1):
+        xsec = xsec_func(C2V, CV, kl,
+                         *(shapes[vbf_name] for vbf_name in all_vbf_names))[0][0]
+        return xsec
+
     for vbf_name in all_vbf_names:
         plt.figure()
         hep.histplot(shapes[vbf_name], range(0,len(all_bin_names)+1))
@@ -49,4 +58,16 @@ if __name__ == '__main__':
         #plt.xticks([])
         #plt.xticks(range(0,len(all_bin_names)), all_bin_names, rotation=45)
         plt.savefig('VBF/{}.pdf'.format(vbf_name))
+        plt.savefig('VBF/{}.png'.format(vbf_name))
+        plt.close()
+
+    for (C2V, kl) in [(0.5, 15), (2, 20), (1.5, -10), (2, 10), (0, 0)]:
+        plt.figure()
+        hep.histplot(xsec_shape(C2V=C2V,kl=kl), range(0,len(all_bin_names)+1))
+        vbf_name = 'qqHH_CV_{}_C2V_{}_kl_{}_hbbhbb'.format(1, C2V, kl).replace('.','p')
+        plt.legend(title=vbf_name+'\n(interp.)')
+        #plt.xticks([])
+        #plt.xticks(range(0,len(all_bin_names)), all_bin_names, rotation=45)
+        plt.savefig('VBF/{}.pdf'.format(vbf_name))
+        plt.savefig('VBF/{}.png'.format(vbf_name))
         plt.close()