Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

update plot style #9

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

update plot style #9

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 64 additions & 106 deletions python/makePostFitPlot_unblind_cpws.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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,
Expand All @@ -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:
Expand Down Expand Up @@ -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("")

Expand All @@ -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")
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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:
Expand All @@ -203,48 +183,29 @@ 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")
else:
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:
Expand All @@ -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)
Expand All @@ -290,15 +249,15 @@ 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()
tex1.SetTextFont(61)
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)
Expand All @@ -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)):
Expand All @@ -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
Expand Down Expand Up @@ -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)):
Expand All @@ -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 = []
Expand All @@ -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")
Expand Down Expand Up @@ -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(
Expand All @@ -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)