Merge pull request #22 from donalrinho/master

Added hist plots of MC and reco vertices for paper

case-studies/flavour/Bc2TauNu/README.md

@@ -39,3 +39,4 @@ A few scripts are also used to generate plots and tables for the paper:
 - `exclusive_bkg_summary_table.py`: summarise the exclusive background statistics and efficiencies for the paper in a table.
 - `plot_max_hem_E.py`: plot charged vs neutral maximum hemisphere energies in signal and background, which are shown in paper.
 - `make_signal_yield_table.py`: make a table for the paper showing signal yields and uncertainties for different number of Z's.
+- `make_yield_BF_summary_tables.py`: make summary tables of yield and branching fraction precision as a function of number of Z's.

case-studies/flavour/Bc2TauNu/calc_BFs.py

@@ -8,6 +8,7 @@
 from userConfig import loc, train_vars, train_vars_vtx
 import plotting
 import utils as ut
+import matplotlib.ticker as plticker
 
 from matplotlib import rc
 rc('font',**{'family':'serif','serif':['Roman']})
@@ -147,9 +148,9 @@
     f.write("\\def\\BFBctoTauNuRelErr{" + "%.1f}\n" % x)
 
 #Make trend plots as a function of NZ for different variables
-params = {"N_Bc2TauNu": {"name": "$N(B_c^+ \\to \\tau^+ \\nu_\\tau)$ relative $\\sigma$"},
-          "BF_Bc2TauNu": {"name": "$\\mathcal{B}(B_c^+ \\to \\tau^+ \\nu_\\tau)$ relative $\\sigma$"},
-          "BF_ratio": {"name": "$R_c$ relative $\\sigma$"}
+params = {"N_Bc2TauNu": {"name": "$N(B_c^+ \\to \\tau^+ \\nu_\\tau)$ relative $\\sigma$","low": 0.02, "high": 0.12},
+          "BF_Bc2TauNu": {"name": "$\\mathcal{B}(B_c^+ \\to \\tau^+ \\nu_\\tau)$ relative $\\sigma$", "low": 0.08, "high": 0.15},
+          "BF_ratio": {"name": "$R_c$ relative $\\sigma$", "low": 0.02, "high": 0.12}
          }
 for v in params:
     fig, ax = plt.subplots(figsize=(9,8))
@@ -162,8 +163,16 @@
     for s in syst:
         plt.errorbar(x=number_of_zs,y=x[s],xerr=None,yerr=None,color=syst[s],fmt="o-",label="$\\sigma_{syst} = %s \\times \\sigma_{stat}$" % s)
     ax.tick_params(axis='both', which='major', labelsize=25)
+    #plt.grid(which="both",axis="y")
+    lmax = plticker.MultipleLocator(base=0.02)
+    lmin = plticker.MultipleLocator(base=0.01)
+    ax.yaxis.set_major_locator(lmax)
+    ax.yaxis.set_minor_locator(lmin)
+    # Add the grid
+    ax.grid(which='both', axis='y', linestyle='-')
     plt.ylabel(params[v]["name"],fontsize=30)
     plt.xlabel("$N_Z (\\times 10^{12})$",fontsize=30)
+    plt.ylim(params[v]["low"],params[v]["high"])
     if(v=="N_Bc2TauNu"):
         plt.legend(loc="upper right",fontsize=25)
     plt.tight_layout()

case-studies/flavour/Bc2TauNu/fit_MVA_dists.py

@@ -11,7 +11,7 @@
 from scipy import interpolate
 
 #Local code
-from userConfig import loc, mode, train_vars, train_vars_vtx, Ediff_cut
+from userConfig import loc, train_vars, train_vars_vtx, Ediff_cut
 import plotting
 import utils as ut
 
@@ -38,7 +38,7 @@
        "MVA2": 0.95
        }
 
-cut = f"EVT_ThrustDiff_E > {Ediff_cut} and EVT_MVA1 > {bdt['MVA1']} and EVT_MVA2 > {bdt['MVA2']}"
+cut = f"EVT_ThrustDiff_E {Ediff_cut} and EVT_MVA1 > {bdt['MVA1']} and EVT_MVA2 > {bdt['MVA2']}"
 
 #Load dataframes for each mode
 tree = {}

case-studies/flavour/Bc2TauNu/make_vertex_hist_plots.py

@@ -0,0 +1,63 @@
+import sys, os, argparse
+import json
+import numpy as np
+from uncertainties import *
+import matplotlib.pyplot as plt
+from ROOT import TFile
+
+#Local code
+from userConfig import loc, train_vars, train_vars_vtx
+import plotting
+import utils as ut
+import matplotlib.ticker as plticker
+
+from matplotlib import rc
+rc('font',**{'family':'serif','serif':['Roman']})
+rc('text', usetex=True)
+
+#Load ROOT file containing hists
+file = TFile.Open(f"{loc.EOS}/vertex_hists.root")
+
+hists = {"h_reco": {"hist": file.Get("h_reco"),"xname": "Number of vertices", "yname": "Density", "xrange": [0,12.1],"yrange": [0,0.4]},
+         "h_mc": {"hist": file.Get("h_mc"),"xname": "Number of vertices", "yname": "Density", "xrange": [0,12.1],"yrange": [0,0.4]},
+         "h_recoeff_SV_3trk": {"hist": file.Get("h_recoeff_SV_3trk"), "xname": "Number of MC tracks", "yname": "Fraction reco. as three-track vertex", "xrange": [0,7.1], "yrange": [0.,1.]}
+         }
+
+for h in hists:
+    hists[h]["bins"] = hists[h]["hist"].GetNbinsX()
+    hists[h]["contents"] = []
+    hists[h]["edges"] = []
+    for i in range(0,hists[h]["bins"]):
+        w = hists[h]["hist"].GetBinWidth(i)
+        hists[h]["contents"].append(hists[h]["hist"].GetBinContent(i))
+        hists[h]["edges"].append(hists[h]["hist"].GetBinCenter(i) - 0.5*w)
+    #Final upper bin edge
+    hists[h]["edges"].append(hists[h]["hist"].GetBinCenter(hists[h]["bins"]) + 0.5*w)
+
+#Plot reco MC comparison
+fig, ax = plt.subplots(figsize=(8,8))
+reco = plt.stairs(hists["h_reco"]["contents"], hists["h_reco"]["edges"], color="crimson", linewidth=2, label="Reco. vertices")
+mc = plt.stairs(hists["h_mc"]["contents"], hists["h_mc"]["edges"], color="dodgerblue", linewidth=2, label="MC vertices ($N_{\\mathrm{charged}} > 1$)")
+plt.xlim(hists["h_reco"]["xrange"])
+plt.ylim(hists["h_reco"]["yrange"])
+ax.tick_params(axis='both', which='major', labelsize=25)
+plt.xlabel(hists["h_reco"]["xname"],fontsize=30)
+plt.ylabel(hists["h_reco"]["yname"],fontsize=30)
+l = plticker.MultipleLocator(base=2.0) # this locator puts ticks at regular intervals
+ax.xaxis.set_major_locator(l)
+plt.legend(fontsize=25)
+plt.tight_layout()
+fig.savefig(f"{loc.PLOTS}/mc_vs_reco_vertices.pdf")
+
+#Plot the efficiency
+fig, ax = plt.subplots(figsize=(8,8))
+eff = plt.stairs(hists["h_recoeff_SV_3trk"]["contents"], hists["h_recoeff_SV_3trk"]["edges"], color="k", linewidth=2)
+plt.xlim(hists["h_recoeff_SV_3trk"]["xrange"])
+plt.ylim(hists["h_recoeff_SV_3trk"]["yrange"])
+ax.tick_params(axis='both', which='major', labelsize=25)
+plt.xlabel(hists["h_recoeff_SV_3trk"]["xname"],fontsize=30)
+plt.ylabel(hists["h_recoeff_SV_3trk"]["yname"],fontsize=30)
+l = plticker.MultipleLocator(base=1.0) # this locator puts ticks at regular intervals
+ax.xaxis.set_major_locator(l)
+plt.tight_layout()
+fig.savefig(f"{loc.PLOTS}/three_track_vertex_eff.pdf")

case-studies/flavour/Bc2TauNu/make_yield_BF_summary_tables.py

@@ -0,0 +1,49 @@
+import sys, os, argparse
+import json
+import numpy as np
+from uncertainties import *
+import matplotlib.pyplot as plt
+
+#Local code
+from userConfig import loc, train_vars, train_vars_vtx
+import plotting
+import utils as ut
+
+from matplotlib import rc
+rc('font',**{'family':'serif','serif':['Roman']})
+rc('text', usetex=True)
+
+number_of_zs = [0.5,1,2,3,4,5]
+syst = [0,0.25,0.5,1]
+params = {"N_Bc2TauNu": {"name": "$N(B_c^+ \\to \\tau^+ \\nu_\\tau)$","low": 0.02, "high": 0.12},
+          "BF_Bc2TauNu": {"name": "$\\mathcal{B}(B_c^+ \\to \\tau^+ \\nu_\\tau)$", "low": 0.08, "high": 0.15},
+          "BF_ratio": {"name": "$R_c$", "low": 0.02, "high": 0.12}
+         }
+
+with open(f'{loc.JSON}/BF_vals.json') as f:
+    vals = json.load(f)
+
+for p in params:
+    print("\\renewcommand{\\arraystretch}{1.4}{")
+    print("\\begin{table[h!]")
+    print("\\centering")
+    print("\\small")
+    print("\\begin{tabular}{ll}")
+    print("$N_Z (\\times 10^{12})$ & Relative $\\sigma$ ($\\sigma_{syst}^N = [0, 0.25, 0.5, 1] \\times \\sigma_{stat}^N$) \\\\ \\hline")
+    for nz in number_of_zs:
+        print(f"{nz} & ", end ="")
+        print("[", end="")
+        for s in syst:
+            x = round(vals[f"{p}_{nz}_{s}"][2],3)
+            if(s!=1):
+                print(x, end=", ")
+            else:
+                print(x, end="")
+        print("]", end="")
+        print(" \\\\")
+    print("\\hline")
+    print("\\end{tabular}")
+    print("\\caption{Estimated relative precision on %s as a function of $N_Z$, where four different levels of systematic uncertainty on the signal yield are shown.}" % params[p]['name'])
+    print("\\label{tab:%s_vs_NZ}" % p)
+    print("\\end{table}")
+    print("}")