|
| 1 | +import matplotlib |
| 2 | +matplotlib.use("Agg") |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +from matplotlib import transforms |
| 5 | +import numpy as np |
| 6 | +import ROOT as r |
| 7 | +import mplhep as hep |
| 8 | +import argparse |
| 9 | + |
| 10 | +plt.style.use(hep.style.ROOT) |
| 11 | + |
| 12 | +def get_vals(fname): |
| 13 | + rfile = r.TFile.Open(fname) |
| 14 | + rtree = rfile.Get("limit") |
| 15 | + vals = [] |
| 16 | + for i in range(rtree.GetEntries()): |
| 17 | + rtree.GetEntry(i) |
| 18 | + mu = rtree.limit |
| 19 | + vals.append(mu) |
| 20 | + return vals |
| 21 | + |
| 22 | +def gofplot(datafile, mcfile, year=2017, savename='fplotX', nbins=130, algo='saturated'): |
| 23 | + gofs = np.array(get_vals(mcfile)) |
| 24 | + gof_data = get_vals(datafile)[0] |
| 25 | + |
| 26 | + print("XXXXXXX") |
| 27 | + print(gof_data) |
| 28 | + print(np.array([np.around(np.mean(gofs) + x * np.std(gofs), 3) for x in [-3,-2,-1,0,1,2,3]])) |
| 29 | + print("XXXXXXX") |
| 30 | + |
| 31 | + from scipy.stats import chi2 |
| 32 | + x_lim = np.max(gofs) * 1.2 |
| 33 | + x_low = np.min(gofs) * 0.9 |
| 34 | + x = np.linspace(x_low, x_lim, 200) |
| 35 | + bins = np.linspace(0, x_lim, 50) |
| 36 | + width = bins[1] - bins[0] |
| 37 | + |
| 38 | + fig, ax = plt.subplots() |
| 39 | + trans = transforms.blended_transform_factory(ax.transData, ax.transAxes) |
| 40 | + if algo == 'saturated': |
| 41 | + ax.plot(x, len(gofs) * width * chi2.pdf(x, np.mean(gofs)), color='red', label='$\chi^2 fit$, ndf = {:.2f}'.format(np.mean(gofs))) |
| 42 | + h, _, _ = ax.hist(gofs, bins, facecolor='none', edgecolor='black', histtype='stepfilled', lw=2, |
| 43 | + label="Toys, N = {}".format(len(gofs))) |
| 44 | + ax.hist(gofs[gofs > gof_data], bins, facecolor='steelblue', edgecolor='gray', histtype='stepfilled', alpha=0.3, |
| 45 | + label='p-value = {}'.format(round(float(len(gofs[gofs > gof_data]))/len(gofs), 3))); |
| 46 | + print("P-value", round(float(len(gofs[gofs > gof_data]))/len(gofs), 3)) |
| 47 | + ax.annotate("", xy=(gof_data, 0), xycoords=trans, |
| 48 | + xytext=(gof_data, 0.25), textcoords=trans, |
| 49 | + arrowprops=dict(lw='4', color='b', arrowstyle="->,head_length=1.5,head_width=0.5"), |
| 50 | + ) |
| 51 | + ax.plot([], [], color='blue', lw=2, label="Observed = {:.2f}".format(gof_data)) |
| 52 | + |
| 53 | + ax.legend() |
| 54 | + hep.cms.label(llabel='Private Work', data=True, year=year, ax=ax) |
| 55 | + ax.set_xlim(np.mean(gofs)-np.std(gofs)*4, np.mean(gofs)+np.std(gofs)*5) |
| 56 | + ax.set_ylim(0, max(h) * 1.4) |
| 57 | + if algo == 'saturated': |
| 58 | + xlab = r"$-2log(\lambda)$" |
| 59 | + else: |
| 60 | + xlab = "KS" |
| 61 | + ax.set_xlabel(xlab , x=1, ha='right') |
| 62 | + ax.set_ylabel("Pseudoexperiments", y=1, ha='right') |
| 63 | + fig.savefig('{}.pdf'.format(savename), dpi=300, transparent=True, bbox_inches='tight') |
| 64 | + fig.savefig('{}.png'.format(savename), dpi=300, transparent=True, bbox_inches='tight') |
| 65 | + |
| 66 | + |
| 67 | +if __name__ == "__main__": |
| 68 | + parser = argparse.ArgumentParser(description='Make plots for ftest/gof files.') |
| 69 | + parser.add_argument('datagof') |
| 70 | + parser.add_argument('mcgofs') |
| 71 | + parser.add_argument("--year", choices=['2016', '2017', '2018', ""], default="") |
| 72 | + parser.add_argument("--algo", choices=['saturated', 'KS'], default="saturated") |
| 73 | + args = parser.parse_args() |
| 74 | + |
| 75 | +sname = "plots/gof_" + args.year + "_" + args.algo |
| 76 | +gofplot(args.datagof, args.mcgofs, year=args.year, savename=sname, algo=args.algo) |
| 77 | + |
| 78 | + |
0 commit comments