diff --git a/python/makePostFitPlot_unblind_cpws.py b/python/makePostFitPlot_unblind_cpws.py index d8bdfb1..23204b6 100644 --- a/python/makePostFitPlot_unblind_cpws.py +++ b/python/makePostFitPlot_unblind_cpws.py @@ -1,13 +1,11 @@ from __future__ import print_function +import argparse import numpy as np import math import ROOT as r -import shlex -import sys import os -from ROOT import TLorentzVector,TGraphAsymmErrors,TMath -from ROOT import Double -from colors import * +from ROOT import TGraphAsymmErrors, Math +from colors import mycolors, default_colors leftMargin = 0.17 rightMargin = 0.04 @@ -21,6 +19,24 @@ r.gStyle.SetStatH(0.15) r.gROOT.SetBatch(1) +def get_tgraph(data, denom=None): + alpha = 1-0.6827 + TGraph_data = TGraphAsymmErrors(data) + for i in range(TGraph_data.GetN()): + idenom = 1 + if denom: + idenom = denom.GetBinContent(i+1) + N = TGraph_data.GetY()[i] + L = 0 + if N>0: + L = Math.gamma_quantile(alpha/2,N,1.) + U = Math.gamma_quantile_c(alpha/2,N+1,1) + TGraph_data.SetPointEYlow(i, (N-L)/idenom) + TGraph_data.SetPointEYhigh(i, (U-N)/idenom) + TGraph_data.SetPoint(i, TGraph_data.GetX()[i], N/idenom) + TGraph_data.SetPointEXlow(i, 0) + TGraph_data.SetPointEXhigh(i, 0) + return TGraph_data def makeplot_single( h1_sig=None, @@ -34,7 +50,8 @@ def makeplot_single( sig_scale_=1.0, dir_name_="plots", output_name_=None, - extraoptions=None + extraoptions=None, + cms_label="Internal" ): if h1_sig == None or h1_bkg == None: @@ -87,27 +104,18 @@ def makeplot_single( stack = r.THStack("stack", "stack") nS = np.zeros(h1_bkg[0].GetNbinsX()) eS = np.zeros(h1_bkg[0].GetNbinsX()) - #hist_all is used to make the data/mc ratio. remove signal for the moment due to signal is scaled right now + # hist_all is used to make the data/mc ratio. remove signal for the moment due to signal is scaled right now hist_all = h1_bkg[0].Clone("hist_all") hist_all.Scale(0.0) - #hist_s = h1_sig[0].Clone("hist_s") hist_b = h1_bkg[0].Clone("hist_b") for idx in range(len(h1_bkg)): stack.Add(h1_bkg[idx]) for ib in range(h1_bkg[0].GetNbinsX()): nS[ib] += h1_bkg[idx].GetBinContent(ib+1) eS[ib] = math.sqrt(eS[ib]*eS[ib] + h1_bkg[idx].GetBinError(ib+1)*h1_bkg[idx].GetBinError(ib+1)) - hist_all.Add(h1_bkg[idx]) + hist_all.Add(h1_bkg[idx]) if idx > 0: hist_b.Add(h1_bkg[idx]) - - #for idx in range(len(h1_sig)): - # print("ggH signal yield: ", hist_s.Integral()) - # if idx > 0: - # hist_temp = h1_sig[idx].Clone(h1_sig[idx].GetName()+"_temp") - #hist_all.Add(hist_temp) - #hist_s.Add(h1_sig[idx]) - #print("all signal yield: ", hist_s.Integral()) stack.SetTitle("") @@ -121,7 +129,6 @@ def makeplot_single( eS[ib] = math.sqrt(eS[ib]*eS[ib] + h1_sig[idx].GetBinError(ib+1)*h1_sig[idx].GetBinError(ib+1)) if stack.GetMaximum() > maxY: maxY = stack.GetMaximum() - #if "SR" in h.GetTitle(): stack.Draw("hist") else: stack.Draw("hist") @@ -130,11 +137,8 @@ def makeplot_single( for idx in range(len(h1_sig)): if h1_sig[idx].GetMaximum() > maxY: maxY = h1_sig[idx].GetMaximum() - #if "SR" in h1_bkg[0].GetTitle(): - #h1_sig[idx].Draw("samehist") - #hist_s.Draw("samehist") - ##draw stack total unc on top of total histogram + # draw stack total unc on top of total histogram box = r.TBox(0,0,1,1,) box.SetFillStyle(3002) box.SetLineWidth(0) @@ -144,40 +148,18 @@ def makeplot_single( if h1_data: if h1_data.GetMaximum() > maxY: - maxY = h1_data.GetMaximum()+np.sqrt(h1_data.GetMaximum()) - #if not "SR" in h1_data.GetTitle() or "fail" in h1_data.GetTitle(): - if True: - #print("debug h1_data.GetName()",h1_data.GetName(), h1_data.GetTitle()) - TGraph_data = TGraphAsymmErrors(h1_data) - for i in range(TGraph_data.GetN()): - #data point - var_x, var_y = Double(0.), Double(0.) - TGraph_data.GetPoint(i,var_x,var_y) - if np.fabs(var_y) < 1e-5: - TGraph_data.SetPoint(i,var_x,-1.0) - TGraph_data.SetPointEYlow(i,-1) - TGraph_data.SetPointEYhigh(i,-1) - #print("zero bins in the data TGraph: bin",i+1) - else: - TGraph_data.SetPoint(i,var_x,var_y) - err_low = var_y - (0.5*TMath.ChisquareQuantile(0.1586555,2.*var_y)) - TGraph_data.SetPointEYlow(i, var_y - (0.5*TMath.ChisquareQuantile(0.1586555,2.*var_y))) - TGraph_data.SetPointEYhigh(i, (0.5*TMath.ChisquareQuantile(1.-0.1586555,2.*(var_y+1))) - var_y) - - TGraph_data.SetMarkerColor(1) - TGraph_data.SetMarkerSize(1) - TGraph_data.SetMarkerStyle(20) - TGraph_data.Draw("same P") + maxY = h1_data.GetMaximum()+np.sqrt(h1_data.GetMaximum()) + TGraph_data = get_tgraph(h1_data) + TGraph_data.SetMarkerColor(1) + TGraph_data.SetMarkerSize(1) + TGraph_data.SetMarkerStyle(20) + TGraph_data.Draw("same ZP") stack.GetYaxis().SetTitle("Events") stack.GetYaxis().SetTitleOffset(1.05) stack.GetYaxis().SetTitleSize(0.08) stack.GetYaxis().SetLabelSize(0.06) - #stack.GetYaxis().CenterTitle() stack.GetXaxis().SetLabelSize(0.) - #stack.GetXaxis().SetLabelOffset(0.013) - #if "xaxis_range" in extraoptions: - # stack.GetXaxis().SetRangeUser(float(extraoptions["xaxis_range"][0]),float(extraoptions["xaxis_range"][1])) leg = r.TLegend(0.2, 0.60, 0.9, 0.88) leg.SetNColumns(3) @@ -187,8 +169,6 @@ def makeplot_single( leg.SetTextSize(0.035) for idx in range(len(h1_bkg)): leg.AddEntry(h1_bkg[idx], bkg_legends_[idx], "F") - #if "SR" in hist_s.GetTitle(): - # leg.AddEntry(hist_s, 'HH #times {:1.2}'.format(sig_scale_), "L") leg.AddEntry(box, "Total unc", "F") if h1_data: @@ -203,40 +183,21 @@ def makeplot_single( ratio_High = 2.0 if h1_data: - ratio = TGraphAsymmErrors(h1_data) - for i in range(ratio.GetN()): - - #bkg prediction - imc = Double(hist_all.GetBinContent(i+1)) - #data point - var_x, var_y = Double(0.), Double(0.) - #if not ("SR" in h1_data.GetTitle() and (i>5 and i<9)): - ratio.GetPoint(i,var_x,var_y) - if var_y == 0.: - ratio.SetPoint(i,var_x,-1.0) - ratio.SetPointEYlow(i,-1) - ratio.SetPointEYhigh(i,-1) - continue - ratio.SetPoint(i,var_x,var_y/imc) - err_low = (var_y - (0.5*TMath.ChisquareQuantile(0.1586555,2.*var_y)))/imc - err_high = ((0.5*TMath.ChisquareQuantile(1.-0.1586555,2.*(var_y+1))) - var_y)/imc - ratio.SetPointEYlow(i, err_low) - ratio.SetPointEYhigh(i, err_high) - + ratio = get_tgraph(h1_data, denom=hist_all) ratio.SetMarkerColor(1) ratio.SetMarkerSize(1) ratio.SetMarkerStyle(20) ratio.GetXaxis().SetTitle("j_{2} regressed mass [GeV]") - #myC.Update() if "ratio_range" in extraoptions: ratio_Low = extraoptions["ratio_range"][0] ratio_High = extraoptions["ratio_range"][1] - ratio.GetYaxis().SetTitle("data/mc") + ratio.GetYaxis().SetTitle("Data/pred.") + ratio.GetYaxis().SetTitleSize(0.01) ratio.GetYaxis().SetRangeUser(ratio_Low, ratio_High) ratio.GetXaxis().SetRangeUser(50, 220) ratio.SetTitle("") - ratio.Draw("same AP") + ratio.Draw("same APZ") pad2.Update() print(ratio.GetTitle(),ratio.GetName(),"debug") @@ -244,7 +205,7 @@ def makeplot_single( ratio = h1_sig[0].Clone("ratio") ratio_High = 0.0 for ibin in range(1,ratio.GetNbinsX()+1): - s = 0#hist_s.GetBinContent(ibin) + s = 0 b = hist_b.GetBinContent(ibin) L = 0.0 if b > 0.0: @@ -271,10 +232,8 @@ def makeplot_single( ratio.GetYaxis().SetLabelSize(0.13) ratio.GetYaxis().SetTickLength(0.01) ratio.GetYaxis().SetNdivisions(505) - #if "xaxis_range" in extraoptions: - # ratio.GetXaxis().SetRangeUser(float(extraoptions["xaxis_range"][0]),float(extraoptions["xaxis_range"][1])) - #draw stack total unc on the ratio plot to present the background uncertainty + # draw stack total unc on the ratio plot to present the background uncertainty box_ratio = r.TBox(0,0,1,1,) box_ratio.SetFillStyle(3002) box_ratio.SetLineWidth(0) @@ -290,7 +249,7 @@ def makeplot_single( ratio.GetXaxis().SetTitle(x_title) ratio.GetYaxis().CenterTitle() - ##########draw CMS preliminary + # draw CMS preliminary pad1.cd() tex1 = r.TLatex(leftMargin, 0.91, "CMS") tex1.SetNDC() @@ -298,7 +257,7 @@ def makeplot_single( tex1.SetTextSize(0.070) tex1.SetLineWidth(2) tex1.Draw() - tex2 = r.TLatex(leftMargin+0.12,0.912,"Internal") + tex2 = r.TLatex(leftMargin+0.12,0.912,cms_label) tex2.SetNDC() tex2.SetTextFont(52) tex2.SetTextSize(0.055) @@ -320,10 +279,9 @@ def makeplot_single( else: outFile = outFile + "/" + hist_name_ - #print("maxY = "+str(maxY)) stack.SetMaximum(maxY*1.7) - #print everything into txt file + # print everything into txt file text_file = open(outFile+"_linY.txt", "w") text_file.write("bin | x ") for idx in range(len(h1_bkg)): @@ -340,7 +298,7 @@ def makeplot_single( text_file.write("-------") text_file.write("\n") - #print yield table for AN + # print yield table for AN text_file.write("print yield table for AN\n") bkg_all = 0 bkg_all_errsq = 0 @@ -374,7 +332,8 @@ def makeplot_single( myC.SaveAs(outFile+"_logY.png") myC.SaveAs(outFile+"_logY.pdf") myC.SaveAs(outFile+"_logY.C") - #save histogram and ratio to root file + + # save histogram and ratio to root file outFile_root = r.TFile(outFile+".root", "recreate") outFile_root.cd() for idx in range(len(h1_bkg)): @@ -384,14 +343,13 @@ def makeplot_single( if h1_data: h1_data.Write() ratio.Write() - #outFile_root.Write() outFile_root.Close() -def main(vbdt, HH_limit): +def main(vbdt, HH_limit, cms_label): dirNameSig = "shapes_fit_s" - debug = False; + debug = False for dirName in ["shapes_prefit", "shapes_fit_s", "shapes_fit_b"]: pnames_sig = [] @@ -400,16 +358,18 @@ def main(vbdt, HH_limit): pnames_bkg = ["ttH_hbb", "VH_hbb", "bbbb_boosted_ggf_others", "ttbar", "bbbb_boosted_ggf_qcd_datadriven","qqHH_CV_1_C2V_1_kl_1_hbbhbb","ggHH_kl_1_kt_1_hbbhbb"] bkg_legends = ["t#bar{t}H", "VH", "V+jets,VV", "t#bar{t}+jets", "QCD+ggH+VBFH", "VBFHH #mu = {:1.2}".format(HH_limit) ,"ggHH #mu = {:1.2}".format(HH_limit)] sig_legends = [] - pname_data = "data" + + # color scheme for an sig_colors = [2, 839, 800, 1, 632] bkg_colors = [2003, 2011, 2001, 2005, 2007, 7, 2] - - #color scheme for paper draft - #sig_colors = [2, 4, 800, 1, 632] - #bkg_colors = [616, 2011, 2001, 601, 797, 2007, 839] + + if cms_label in ["", "Supplementary"]: + # color scheme for paper draft + sig_colors = [2, 4, 800, 1, 632] + bkg_colors = [616, 2011, 2001, 601, 797, 2, 839] region = "SBplusfail" - #the first three are the sideband region, the next three is the full AK8 jet 2 mass region, the last one is the common fail region + # the first three are the sideband region, the next three is the full AK8 jet 2 mass region, the last one is the common fail region bdtbins = ["SRBin1", "SRBin2", "SRBin3","fitfail"] for bdtbin in bdtbins: print("looking at file: ",vbdt+"/fitDiagnostics"+region+".root") @@ -471,14 +431,9 @@ def main(vbdt, HH_limit): if dirthis_data.GetListOfKeys().Contains("data"): g = fin.Get(dirName+"/"+bdtbin_data+"/data") - x,y = r.Double(0), r.Double(0) for idx in range(nBins): - g.GetPoint(idx, x, y) + y = g.GetY()[idx] h1_data.SetBinContent(idx+1, y*10) - #if "SR" in h1_data.GetTitle(): - # h1_data.SetBinContent(7,0) - # h1_data.SetBinContent(8,0) - # h1_data.SetBinContent(9,0) h1_data.GetXaxis().SetTitle("j_{2} regressed mass (GeV)") makeplot_single( @@ -492,12 +447,15 @@ def main(vbdt, HH_limit): sig_scale_=HH_limit, output_name_=dirName+"_"+bdtbin+"_BDTv8p2", dir_name_= "plots/postfitplots/output_all_"+vbdt+"/", - extraoptions={"stack_signal": False} + extraoptions={"stack_signal": False}, + cms_label=cms_label ) if __name__ == "__main__": - vbdt = "v8p2yield_AN_sr_sys_1123v1_kl" - #"v8p2yield_AN_sr_sys_0830_fix2017trigSF0908_SDv1" - #"v8p2yield_AN_sr_sys_0830_fix2017trigSF0908" - mu_HH = 3.0; - main(vbdt,mu_HH) + parser = argparse.ArgumentParser() + parser.add_argument( '--muHH', help='signal strength best-fit' , type=float, default=1.0) + parser.add_argument( '--vbdt', help='dir for postfit histogram', default='v8p2yield_AN_sr_sys_1123v1_kl') + parser.add_argument( '--label', choices = ["Internal", "Preliminary", "Supplementary", ""], help='CMS label', default='Internal') + args = parser.parse_args() + + main(args.vbdt, args.muHH, args.label)