update

check_yields.py

 import sys
 import numpy as np
 from collections import OrderedDict
+import matplotlib
+matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 import mplhep as hep
 plt.style.use(hep.style.ROOT)
@@ -22,10 +24,10 @@ if __name__ == '__main__':
     ])
 
     skip_vbf = [(0.5, 1, 1)] # for model_default
+    #skip_vbf = [(1, 2, 1)] # for model_no_vbf_C2V2
     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])
-
     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')
 
@@ -33,9 +35,11 @@ if __name__ == '__main__':
 
     content = read_datacard_structured(datacard)
     #all_vbf_names = [p["name"] for p in content["processes"] if "qqHH" in p["name"]]
+    all_vbf_labels = [sample.label for sample in vbf_formula.samples] # get sample label
     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"]]
 
+    zero_bin_names = [bin_name for bin_name in all_bin_names if content['observations'][bin_name]==0]
     shapes = OrderedDict()
     for vbf_name in all_vbf_names:
         shapes[vbf_name] = np.zeros(len(all_bin_names))
@@ -55,19 +59,35 @@ if __name__ == '__main__':
         plt.figure()
         hep.histplot(shapes[vbf_name], range(0,len(all_bin_names)+1))
         plt.legend(title=vbf_name)
-        #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()
 
+    print('zero_bins!')
+    print(zero_bin_names)
+    print("\n")
     for (C2V, kl) in [(0.5, 15), (2, 20), (1.5, -10), (2, 10), (0, 0)]:
+        print("formula for C2V=%s, kl=%s"%(C2V, kl))
+        print(vbf_formula.sigma[0].subs(zip(["CV", "C2V", "kl"]+list(map("xs{}".format, range(vbf_formula.n_samples))), [1, C2V, kl]+all_vbf_labels)))
+        shape = xsec_shape(C2V=C2V,kl=kl)
+        neg_bins = shape < 0
+        if np.sum(neg_bins): 
+            print('neg_bins!')
+            neg_bin_names = list(np.array(all_bin_names)[neg_bins])
+            print(neg_bin_names)
+            print("\n")
+
+            print('neg_bins & zero_bins!')
+            print(list(set(neg_bin_names) & set(zero_bin_names)))
+            print("\n")
+        else:
+            print("no neg bins!")
+            print("\n")
+
         plt.figure()
-        hep.histplot(xsec_shape(C2V=C2V,kl=kl), range(0,len(all_bin_names)+1))
+        hep.histplot(shape, 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()