#!/usr/bin/env python
import re
import numpy as np
import matplotlib.pyplot as plt
from filval.plotting import (decl_plot, render_plots, hist_plot, Plot, generate_dashboard, hists_to_table)
from filval.histogram import hist, hist_add
X_TICKS = ['CRZ', 'CRW']+[f'SR{i}' for i in range(1, 17)]
def read_yields(tag, postfit=False):
import ROOT
f = ROOT.TFile.Open(f'data/JetPtStudies/{tag}/mlfit.root')
hists = {'tttt': hist(f.Get('shapes_prefit/SS/tttt')),
'ttw': hist(f.Get('shapes_prefit/SS/ttw')),
'ttz': hist(f.Get('shapes_prefit/SS/ttz')),
'tth': hist(f.Get('shapes_prefit/SS/tth'))
}
try:
hists['fakes_mc'] = hist(f.Get('shapes_prefit/SS/fakes_mc'))
except:
print(tag, 'has no fakes_mc')
hists['fakes_mc'] = None
return hists
def set_xticks():
plt.xticks([x-0.5 for x in range(1, 19)], X_TICKS, rotation=60)
@decl_plot
def plot_yields(bJetPt):
"""\
[SR Definitions](https://github.com/cfangmeier/FTAnalysis/blob/new_baseline/analysis/misc/signal_regions.h#L127)
"""
_, ((ax_tttt, ax_ttw), (ax_ttz, ax_tth)) = plt.subplots(2, 2)
if bJetPt == 25:
cut_40 = read_yields('v1.0.4_JetPtCut40')
cut_35 = read_yields('v1.0.4_JetPtCut35')
cut_30 = read_yields('v1.0.4_JetPtCut30')
cut_25 = read_yields('v1.0.4_JetPtCut25')
cut_20 = read_yields('v1.0.4_JetPtCut20')
else: # bJetPt == 20
cut_40 = read_yields('v1.0.4_bJetPtCut20JetPtCut40')
cut_35 = read_yields('v1.0.4_bJetPtCut20JetPtCut35')
cut_30 = read_yields('v1.0.4_bJetPtCut20JetPtCut30')
cut_25 = read_yields('v1.0.4_bJetPtCut20JetPtCut25')
cut_20 = read_yields('v1.0.4_bJetPtCut20JetPtCut20')
def do_plot(axis, key):
plt.sca(axis)
hist_plot(cut_40[key], label="Jet $P_T>$40", include_errors=True)
hist_plot(cut_35[key], label="Jet $P_T>$35", include_errors=True)
hist_plot(cut_30[key], label="Jet $P_T>$30", include_errors=True)
hist_plot(cut_25[key], label="Jet $P_T>$25", include_errors=True)
hist_plot(cut_20[key], label="Jet $P_T>$20", include_errors=True)
plt.title(key)
set_xticks()
do_plot(ax_tttt, 'tttt')
do_plot(ax_ttw, 'ttw')
plt.legend()
do_plot(ax_ttz, 'ttz')
do_plot(ax_tth, 'tth')
tables = []
for cut, label in zip([cut_40, cut_35, cut_30, cut_25, cut_20],
["Jet $P_T>$40", "Jet $P_T>$35", "Jet $P_T>$30", "Jet $P_T>$25", "Jet $P_T>$20"]):
tables.append(f'
{label}
')
tables.append(hists_to_table([cut['tttt'], cut['ttw'], cut['tth'], cut['ttz']],
column_labels=X_TICKS, row_labels=['TTTT', 'TTW', 'TTH', 'TTZ']))
return ''.join(tables)
@decl_plot
def plot_yields_ratio():
_, ((ax_tttt, ax_ttw), (ax_ttz, ax_tth)) = plt.subplots(2, 2)
cut_40 = read_yields('v1.0.4_JetPtCut40')
cut_35 = read_yields('v1.0.4_JetPtCut35')
cut_30 = read_yields('v1.0.4_JetPtCut30')
cut_25 = read_yields('v1.0.4_JetPtCut25')
cut_20 = read_yields('v1.0.4_JetPtCut20')
def do_plot(axis, key):
plt.sca(axis)
c35 = (cut_35[key][0] / cut_40[key][0], *cut_35[key][1:])
c30 = (cut_30[key][0] / cut_40[key][0], *cut_30[key][1:])
c25 = (cut_25[key][0] / cut_40[key][0], *cut_25[key][1:])
c20 = (cut_20[key][0] / cut_40[key][0], *cut_25[key][1:])
hist_plot(c35, label="Jet $P_T>$35")
hist_plot(c30, label="Jet $P_T>$30")
hist_plot(c25, label="Jet $P_T>$25")
hist_plot(c20, label="Jet $P_T>$20")
plt.title(key)
axis.set_ylim((0, 10))
set_xticks()
do_plot(ax_tttt, 'tttt')
plt.legend()
do_plot(ax_ttw, 'ttw')
do_plot(ax_ttz, 'ttz')
do_plot(ax_tth, 'tth')
@decl_plot
def plot_signal_over_background(tags):
cuts = [read_yields(tag) for tag in tags]
def clean_label(label):
return label.replace('_', r'\_')
def do_plot_significance(hists, label):
background = hist_add(hists['ttw'], hists['ttz'], hists['tth'])
tttt = hists['tttt']
ttw = hists['ttw'][0]
ttz = hists['ttz'][0]
tth = hists['tth'][0]
sigma2_bg = (0.4*ttw)**2 + (0.4*ttz)**2 + (0.5*tth)**2
if hists['fakes_mc'] is not None:
fakes_mc = hists['fakes_mc']
background = hist_add(background, fakes_mc)
sigma2_bg += (0.2*fakes_mc[0])**2
ratio = tttt[0] / np.sqrt(tttt[0] + background[0] + sigma2_bg), *tttt[1:]
hist_plot(ratio, label=clean_label(label))
plt.sca(plt.subplot(221))
for cut, tag in zip(cuts, tags):
signal = cut['tttt']
hist_plot(signal, label=clean_label(tag))
set_xticks()
plt.ylabel('TTTT')
plt.ylim((0, 10))
plt.sca(plt.subplot(222))
for cut, tag in zip(cuts, tags):
background = hist_add(cut['ttw'], cut['ttz'], cut['tth'])
hist_plot(background, label=clean_label(tag))
set_xticks()
plt.ylabel('TTW+TTZ+TTH+MCFakes')
plt.ylim((0, 10))
plt.sca(plt.subplot(212))
for cut, tag in zip(cuts, tags):
do_plot_significance(cut, tag)
set_xticks()
plt.ylabel(r'SIG / $\sqrt{SIG+BG + \sigma_{BG}^2}$')
plt.legend()
tables = []
for cut, tag in zip(cuts, tags):
tables.append(f'{tag}
')
tables.append(hists_to_table([cut['tttt'], cut['ttw'], cut['tth'], cut['ttz']],
column_labels=X_TICKS, row_labels=['TTTT', 'TTW', 'TTH', 'TTZ']))
return ''.join(tables)
@decl_plot
def plot_sigs(include_fakes):
"""\
Significances above only include TTZ, TTW, TTH, and (optionally) MC-Fakes(TTBar) backgrounds.
"""
def get_sig(tag):
regex = re.compile('Sig: ([0-9.]+)')
ext = '' if include_fakes else '_nofakes'
with open(f'data/JetPtStudies/{tag}{ext}.txt') as f:
sig = float(regex.findall(f.read())[0])
return sig
sigs_b25 = {n: get_sig(f'v1.0.4_JetPtCut{n}') for n in [40, 35, 30, 25, 20]}
sigs_b20 = {n: get_sig(f'v1.0.4_bJetPtCut20JetPtCut{n}') for n in [40, 35, 30, 25, 20]}
plt.scatter(*zip(*sigs_b25.items()), label="bJetMinPt=25")
for key, val in sigs_b25.items():
plt.text(key, val-0.1, f'{key}: {val:0.3f}')
plt.scatter(*zip(*sigs_b20.items()), label="bJetMinPt=20")
for key, val in sigs_b20.items():
plt.text(key, val+0.1, f'{key}: {val:0.3f}')
plt.ylabel('Significance')
plt.xlabel('Jet $P_T$ Cut')
plt.xticks([20, 25, 30, 35, 40])
plt.xlim((18, 43))
plt.ylim((0, 3.0))
plt.legend()
def main():
sig_o_bg = plot_signal_over_background, (['v1.0.4_JetPtCut40', 'v1.0.4_JetPtCut25', 'v1.0.4_JetPtCut20'],)
sig_o_bg2 = plot_signal_over_background, (['v1.0.4_JetPtCut25', 'v1.0.4_bJetPtCut20JetPtCut25'],)
plots = [
Plot((plot_yields, (25,)), 'Prefit Yields (bJetPt>25)'),
Plot((plot_yields, (20,)), 'Prefit Yields (bJetPt>20)'),
Plot(plot_yields_ratio, 'Prefit Yields Relative to Baseline'),
Plot((plot_sigs, (True,)), 'Significances with Fakes'),
Plot((plot_sigs, (False,)), 'Significances w/o Fakes'),
Plot(sig_o_bg, 'Binned Yields with variable JetPt'),
Plot(sig_o_bg2, 'Binned Yields with variable bJetPt'),
]
render_plots(plots, to_disk=False)
generate_dashboard(plots, 'New Baseline Yields - Various JetPt Cuts',
output='new_baseline_yields_jet_pt_cuts.html',
source=__file__,
)
if __name__ == '__main__':
main()