|
| 1 | +import os |
| 2 | +import sys |
| 3 | +import time |
| 4 | +import argparse |
| 5 | +import json |
| 6 | + |
| 7 | +import ROOT as r |
| 8 | +import os |
| 9 | +import numpy as np |
| 10 | +import matplotlib |
| 11 | +matplotlib.use("Agg") |
| 12 | +import matplotlib.pyplot as plt |
| 13 | +import matplotlib.transforms as transforms |
| 14 | +import mplhep as hep |
| 15 | + |
| 16 | +plt.style.use(hep.style.ROOT) |
| 17 | + |
| 18 | +BASE = os.getcwd() |
| 19 | + |
| 20 | +def get_vals(fname): |
| 21 | + rfile = r.TFile.Open(fname) |
| 22 | + rtree = rfile.Get("limit") |
| 23 | + vals = [] |
| 24 | + for i in range(rtree.GetEntries()): |
| 25 | + rtree.GetEntry(i) |
| 26 | + mu = rtree.limit |
| 27 | + vals.append(mu) |
| 28 | + return vals |
| 29 | + |
| 30 | +def fval(lambda1, lambda2, p1, p2, nbins): |
| 31 | + |
| 32 | + numerator = -2.0 * np.log(np.array(lambda1) / np.array(lambda2)) / (p2 - p1) |
| 33 | + denominator = -2.0 * np.log(np.array(lambda2)) / (nbins - p2) |
| 34 | + |
| 35 | + return numerator / denominator |
| 36 | + |
| 37 | +def fplot(dname, ref, alt, year=2017, savename='fplotX', nbins=130): |
| 38 | + ref_pt, ref_rho = ref |
| 39 | + alt_pt, alt_rho = alt |
| 40 | + p1 = (ref_pt + 1) * (ref_rho + 1) |
| 41 | + p2 = (alt_pt + 1) * (alt_rho + 1) |
| 42 | + |
| 43 | + alt = get_vals('{dname}/bkgtest_{ref_pt}-{ref_rho}_{alt_pt}-{alt_rho}/gofsalt.root' |
| 44 | + .format(dname=dname, ref_pt=ref_pt, ref_rho=ref_rho, alt_pt=alt_pt, alt_rho=alt_rho)) |
| 45 | + base = get_vals('{dname}/bkgtest_{ref_pt}-{ref_rho}_{alt_pt}-{alt_rho}/gofsbase.root' |
| 46 | + .format(dname=dname, ref_pt=ref_pt, ref_rho=ref_rho, alt_pt=alt_pt, alt_rho=alt_rho)) |
| 47 | + if len(alt) != len(base): |
| 48 | + raise ValueError("Number of toys for base and ref does not match.") |
| 49 | + fvals = fval(base, alt, 2, 3, nbins) |
| 50 | + f_data = fval(get_vals('{dname}/bkgtest_{ref_pt}-{ref_rho}_{alt_pt}-{alt_rho}/refbase.root' |
| 51 | + .format(dname=dname, ref_pt=ref_pt, ref_rho=ref_rho, alt_pt=alt_pt, alt_rho=alt_rho)), |
| 52 | + get_vals('{dname}/bkgtest_{ref_pt}-{ref_rho}_{alt_pt}-{alt_rho}/refalt.root' |
| 53 | + .format(dname=dname, ref_pt=ref_pt, ref_rho=ref_rho, alt_pt=alt_pt, alt_rho=alt_rho)), |
| 54 | + 2, 3, nbins)[0] |
| 55 | + |
| 56 | + from scipy.stats import f |
| 57 | + x_lim = max(np.percentile(fvals, 90), f_data*1.2) |
| 58 | + x = np.linspace(0, x_lim, 200) |
| 59 | + bins = np.linspace(0, x_lim, 30) |
| 60 | + width = bins[1] - bins[0] |
| 61 | + |
| 62 | + fig, ax = plt.subplots() |
| 63 | + trans = transforms.blended_transform_factory(ax.transData, ax.transAxes) |
| 64 | + ax.plot(x, len(base) * width * f.pdf(x, p2 - p1, nbins - p2), color='red', label='F-dist, ndf({},{})'.format(p2-p1, nbins-p2)) |
| 65 | + ax.hist(fvals, bins, facecolor='none', edgecolor='black', histtype='stepfilled', lw=2, |
| 66 | + label="Toys, N = {}".format(len(fvals))) |
| 67 | + ax.hist(fvals[fvals > f_data], bins, facecolor='steelblue', edgecolor='gray', histtype='stepfilled', alpha=0.3, |
| 68 | + label='p-value = {}'.format(round(float(len(fvals[fvals > f_data]))/len(fvals), 3))); |
| 69 | + ax.annotate("", xy=(f_data, 0), xycoords=trans, |
| 70 | + xytext=(f_data, 0.25), textcoords=trans, |
| 71 | + arrowprops=dict(lw='4', color='b', arrowstyle="->,head_length=1.5,head_width=0.5"), |
| 72 | + ) |
| 73 | + ax.plot([], [], color='blue', lw=2, label="Observed = {:.3f}".format(f_data)) |
| 74 | + |
| 75 | + title = "TF({},{}) x TF({},{})".format(ref_pt, ref_rho, alt_pt, alt_rho) |
| 76 | + ax.legend(title=title) |
| 77 | + hep.cms.label(data=True, year=year, ax=ax) |
| 78 | + ax.set_xlim(0, x_lim) |
| 79 | + |
| 80 | + fig.savefig('{}.pdf'.format(savename), dpi=300, transparent=True, bbox_inches='tight') |
| 81 | + fig.savefig('{}.png'.format(savename), dpi=300, transparent=True, bbox_inches='tight') |
| 82 | + |
| 83 | +if __name__ == "__main__": |
| 84 | + parser = argparse.ArgumentParser(description='Make plots for ftest/gof files.') |
| 85 | + parser.add_argument('dname', nargs='?', default=os.getcwd()) |
| 86 | + parser.add_argument("--degs", |
| 87 | + type=str, |
| 88 | + default=None, |
| 89 | + help="Polynomial degrees in the shape 'pt,rho' e.g. '2,2'") |
| 90 | + parser.add_argument('-d', '--savedir', type=str, default=None, help="Plots out") |
| 91 | + parser.add_argument('--year', type=str, default=None, help="Year to display on plots.") |
| 92 | + args = parser.parse_args() |
| 93 | + |
| 94 | + configs = json.load(open(os.path.join(BASE, args.dname, "config.json"))) |
| 95 | + if args.year is None: |
| 96 | + args.year = str(configs['year']) |
| 97 | + if args.degs is None: |
| 98 | + args.degs = str(configs['degs']) |
| 99 | + ref_pt, ref_rho = tuple([int(s) for s in args.degs.split(',')]) |
| 100 | + nbins = int(configs['NBINS']) + int(configs['NBINSMU']) |
| 101 | + |
| 102 | + sname = "{}_{}_{}".format(os.path.join(args.savedir, args.dname), |
| 103 | + str(ref_pt) + str(ref_rho), |
| 104 | + str(ref_pt + 1) + str(ref_rho), |
| 105 | + ) |
| 106 | + fplot(os.path.join(BASE, args.dname), |
| 107 | + (ref_pt, ref_rho), (ref_pt + 1, ref_rho), args.year, sname) |
| 108 | + sname = "{}_{}_{}".format(os.path.join(args.savedir, args.dname), |
| 109 | + str(ref_pt) + str(ref_rho), |
| 110 | + str(ref_pt) + str(ref_rho + 1), |
| 111 | + ) |
| 112 | + fplot(os.path.join(BASE, args.dname), |
| 113 | + (ref_pt, ref_rho), (ref_pt, ref_rho + 1), args.year, sname) |
| 114 | + sname = "{}_{}_{}".format(os.path.join(args.savedir, args.dname), |
| 115 | + str(ref_pt) + str(ref_rho), |
| 116 | + str(ref_pt + 1) + str(ref_rho + 1), |
| 117 | + ) |
| 118 | + fplot(os.path.join(BASE, args.dname), |
| 119 | + (ref_pt, ref_rho), (ref_pt + 1, ref_rho + 1), args.year, sname) |
0 commit comments