From ccdc55f64215fe5853447f71b9c7a0359ef3163d Mon Sep 17 00:00:00 2001 From: "[Minxi" <[yhzz598@gmail.com]> Date: Thu, 30 Jan 2020 20:13:33 -0500 Subject: [PATCH 1/4] leptonFlavor_plot --- helpers.py | 16 +- plotAllMCFlavor.py | 955 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 965 insertions(+), 6 deletions(-) create mode 100644 plotAllMCFlavor.py diff --git a/helpers.py b/helpers.py index 73e94ab..c600adf 100755 --- a/helpers.py +++ b/helpers.py @@ -39,13 +39,13 @@ def negWeightFractions(path,muon=True): def binning(channel='muon'): - if channel == 'muon': + #if channel == 'muon': - return [60., 65.43046396, 71.3524269, 77.81037328, 84.85281374, 92.53264952, 100.90756983, 110.04048518, 120., 129.95474058, 140.73528833, 152.41014904, 165.0535115, 178.74571891, 193.57377942, 209.63191906, 227.02218049, 245.85507143, 266.2502669, 288.3373697, - 312.25673399, 338.16035716, 366.21284574, 396.59246138, 429.49225362, 465.12128666, 503.70596789, 545.49148654, 590.74337185, 639.74918031, 692.82032303, 750.29404456, 812.53556599, 879.94040575, 952.93689296, - 1031.98888927, 1117.59873655, 1210.310449, 1310.71317017, 1419.4449167, - 1537.19663264, 1664.71658012, 1802.81509423, 1952.36973236, 2114.3308507, 2289.72764334, 2479.6746824, 2685.37900061, 2908.14776151, 3149.39656595, 3410.65844758, 3693.59361467, 4000. ] - if channel == 'electron': + #return [60., 65.43046396, 71.3524269, 77.81037328, 84.85281374, 92.53264952, 100.90756983, 110.04048518, 120., 129.95474058, 140.73528833, 152.41014904, 165.0535115, 178.74571891, 193.57377942, 209.63191906, 227.02218049, 245.85507143, 266.2502669, 288.3373697, + #312.25673399, 338.16035716, 366.21284574, 396.59246138, 429.49225362, 465.12128666, 503.70596789, 545.49148654, 590.74337185, 639.74918031, 692.82032303, 750.29404456, 812.53556599, 879.94040575, 952.93689296, + #1031.98888927, 1117.59873655, 1210.310449, 1310.71317017, 1419.4449167, + #1537.19663264, 1664.71658012, 1802.81509423, 1952.36973236, 2114.3308507, 2289.72764334, 2479.6746824, 2685.37900061, 2908.14776151, 3149.39656595, 3410.65844758, 3693.59361467, 4000. ] + #if channel == 'electron': # ~ return ([j for j in range(50, 120, 5)] + # ~ [j for j in range(120, 150, 5)] + # ~ [j for j in range(150, 200, 10)] + @@ -373,6 +373,10 @@ def __init__(self,processes,lumi,plot,zScaleFac): self.theHistogram = temphist.Clone() else: self.theHistogram.Add(temphist.Clone()) + def Add(self,addstack): + for h in addstack.theStack.GetHists(): + self.theStack.Add(h.Clone()) + self.theHistogram.Add(h.Clone()) class TheStackRun2: from ROOT import THStack diff --git a/plotAllMCFlavor.py b/plotAllMCFlavor.py new file mode 100644 index 0000000..2e5ec1c --- /dev/null +++ b/plotAllMCFlavor.py @@ -0,0 +1,955 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT, TGaxis +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy +import numpy as np +import root_numpy + +# Muon sys uncertainty (%) +# as a function of mass +def znorm(hist): + lowLimit=hist.FindBin(60) + upLimit=hist.FindBin(120) + znorm=0 + for i in range(lowLimit+1,upLimit+1): + znorm+=hist.GetBinContent(i) + return znorm + +def getMuErr(mass, chann, norm=False): + lumi = 0.0 + znorm = 0.0 + pileup = 0.0 # we don't use pileup 0.046 + dybkg = 0.0 # 0.07 + #pdf = 0.01*(0.433+0.003291*mass-2.159e-6*mass**2+9.044e-10*mass**3-1.807e-13*mass**4+1.51e-17*mass**5) + pdf = 0.0 + + # muons only next + if chann: mutrig = 0.003 + else: mutrig = 0.007 + resolution = 0.01 + muid = 0.05 + + if norm: + lumi = 0.0 + znorm = 0.0 + dybkg = 0 + return math.sqrt(lumi**2+znorm**2+pileup**2+dybkg**2+pdf**2+mutrig**2+resolution**2+muid**2) + + +# chann = True if BB +# chann = False if BE +def getElErr(mass, chann, norm=False): + + lumi = 0.0 + znorm = 0.0 + pileup = 0.046 # we don't use pileup 0.046 + dybkg = 0.0 # 0.07 + + # poly values are in % + #pdf = 0.01*(0.433 + 0.003291*mass - 2.159e-6*mass**2 + 9.044e-10*mass**3 - 1.807e-13*mass**4 + 1.51e-17*mass**5) + pdf = 0.0 + + # the following two are electrons only + if chann: energyscale = 0.02 + else: energyscale = 0.01 + + if chann: + if mass < 90: idscale = 0.01 + elif mass < 1000: idscale = 0.00002198 * mass + 0.008 + else: idscale = 0.03 + else: + if mass < 90: idscale = 0.01 + elif mass < 300: idscale = 0.00014286 * mass - 0.00285 + else: idscale = 0.04 + + if chann: scalefac = 0.03 + else: scalefac = 0.05 + + if norm: + lumi = 0.0 + znorm = 0.0 + dybkg = 0 + return math.sqrt(lumi**2+znorm**2+ pileup**2 + dybkg**2 + pdf**2 + energyscale**2 + idscale**2 + scalefac**2) + +def getErrors(default, others, keys): + dfarr=root_numpy.hist2array(default) + errs={} + for other, key in zip(others,keys): + if type(other)==list: + err1=root_numpy.hist2array(other[0])-dfarr + err1=abs(err1) + err2=root_numpy.hist2array(other[1])-dfarr + err2=abs(err2) + err=np.maximum(err1,err2) + errs[key]=err + else: + err=root_numpy.hist2array(other)-dfarr + err=abs(err) + errs[key]=err + return errs + +# multiply hist by 1/(Acceptance x Efficiency) +def inverseAE(hist, plotObj, year): + # muon and electron + # BB and BE + if year == 2016: + if plotObj.muon: + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 600: + ae = 2.129-0.1268*math.exp(-(mass-119.2)/22.35)-2.386*mass**(-0.03619) + else: + ae = 2.891-2.291e+04/(mass+8294.)-0.0001247*mass + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 450: + ae = 13.56-6.672*math.exp((mass+4.879e+06)/7.233e+06)-826*mass**(-1.567) + else: + ae = 0.2529+0.06511*mass**0.8755*math.exp(-(mass+4601.)/1147) + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + else: # is electron + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = 0.6386-497.7/(mass+348.7) + 69570.0/(mass**2+115400.0) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = -0.03377+735.1/(mass+1318)-88890.0/(mass**2+75720)+14240000.0/(mass**3+23420000) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + + elif year == 2017: + if plotObj.muon: + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 600: + ae = 2.13-0.1313*math.exp(-(mass-110.9)/20.31)-2.387*mass**(-0.03616) + else: + ae = 4.931-55500.0/(mass+11570.0)-0.0002108*mass + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 450: + ae = 13.39-6.696*math.exp((mass+4855000.0)/7431000.0)-108.8*mass**(-1.138) + else: + ae = 0.3148+0.04447*mass**1.42*math.exp(-(mass+5108.0)/713.5) + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + else: # is electron + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = 0.585-404.6/(mass+279.5) + 56180.0/(mass**2+91430.0) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = 0.02066+429.7/(mass+729)-90960.0/(mass**2+71900)+13780000.0/(mass**3+22050000) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif year == 2018: + if plotObj.muon: + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 600: + ae = 2.14-0.1286*math.exp(-(mass-110.6)/22.44)-2.366*mass**(-0.03382) + else: + ae = 5.18-58450.0/(mass+11570.0)-0.0002255*mass + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + if mass < 450: + ae = 13.4-6.693*math.exp((mass+4852000.0)/7437000.0)-81.43*mass**(-1.068) + else: + ae = 0.3154+0.04561*mass**1.362*math.exp(-(mass+4927.0)/727.5) + #print mass, ae + if mass < 120: ae = float("inf") + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + else: # is electron + if "BB" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = 0.576-417.7/(mass+381.8) + 46070.0/(mass**2+107200) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + elif "BE" in plotObj.fileName: + for i in range(1, hist.GetSize()-1): + mass = hist.GetBinCenter(i) + ae = 0.01443+475.7/(mass+639.1)-105600.0/(mass**2+82810)+12890000.0/(mass**3+23170000) + #print mass, ae + hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + +def Stacks(processes,lumi,plot,zScale): + stacks=[] + for i in range(3): + stacks.append(TheStack(processes[i],lumi[i],plot,zScale[i])) + return stacks +def Addhist(histlist): + tempHist=histlist[0] + for i in range(1,3): + tempHist.Add(histlist[i]) + return tempHist +def Addstack(Stacklist): + tempStack=Stacklist[0] + for i in range(1,3): + tempStack.Add(Stacklist[i]) + return tempStack +def plotDataMC(args,plot_mu,plot_el, category): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if args.ratio: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.5,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.5) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw() + ratioPad.Draw() + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + + # Data load processes + colors = createMyColors() + if args.use2016: + data_mu = Process(Data2016, normalized=True) + data_el = Process(Data2016, normalized=True) + elif args.use2018: + data_mu = Process(Data2018, normalized=True) + data_el = Process(Data2018, normalized=True) + elif args.useall: + data_all=[Process(Data2016, normalized=True),Process(Data, normalized=True),Process(Data2018, normalized=True)] + else: + data_mu = Process(Data, normalized=True) + data_el = Process(Data, normalized=True) + + eventCounts_mu = totalNumberOfGeneratedEvents(path,plot_mu["default"].muon) + eventCounts_el = totalNumberOfGeneratedEvents(path,plot_el["default"].muon) + negWeights_mu = negWeightFractions(path,plot_mu["default"].muon) + negWeights_el = negWeightFractions(path,plot_el["default"].muon) + + # Background load processes + backgrounds = copy(args.backgrounds) + if plot_mu["default"].useJets: + if "Wjets" in backgrounds: + backgrounds.remove("Wjets") + backgrounds.insert(0,"Jets") + processes_mu = [] + processes_el = [] + processes_mu2016=[] + processes_mu2017=[] + processes_mu2018=[] + processes_el2016=[] + processes_el2017=[] + processes_el2018=[] + for background in backgrounds: + if args.use2016: + if background == "Jets": + processes_mu.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el,normalized=True)) + else: + processes_mu.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) + processes_el.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el)) + elif args.use2018: + if background == "Jets": + processes_mu.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el,normalized=True)) + else: + processes_mu.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu)) + processes_el.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el)) + elif args.useall: + if background == "Jets": + processes_mu2016.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el2016.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el,normalized=True)) + processes_mu2017.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el2017.append(Process(getattr(Backgrounds,background),eventCounts_el,negWeights_el,normalized=True)) + processes_mu2018.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el2018.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el,normalized=True)) + processes_mu=[processes_mu2016,processes_mu2017,processes_mu2018] + processes_el=[processes_mu2016,processes_mu2017,processes_mu2018] + else: + processes_mu2016.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) + processes_el2016.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el)) + processes_mu2017.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu)) + processes_el2017.append(Process(getattr(Backgrounds,background),eventCounts_el,negWeights_el)) + processes_mu2018.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu)) + processes_el2018.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el)) + processes_mu=[processes_mu2016,processes_mu2017,processes_mu2018] + processes_el=[processes_mu2016,processes_mu2017,processes_mu2018] + else: + if background == "Jets": + processes_mu.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu,normalized=True)) + processes_el.append(Process(getattr(Backgrounds,background),eventCounts_el,negWeights_el,normalized=True)) + else: + processes_mu.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu)) + processes_el.append(Process(getattr(Backgrounds,background),eventCounts_el,negWeights_el)) + + + legend = TLegend(0.55, 0.75, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + + + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + + # Modify legend information + legendHistData_mu = ROOT.TH1F() + legendHistData_el = ROOT.TH1F() + dy_mu = ROOT.TH1F() + dy_el = ROOT.TH1F() + if args.data: + legendHistData_mu.SetMarkerColor(ROOT.kViolet) + legendHistData_el.SetMarkerColor(ROOT.kOrange) + dy_mu.SetLineColor(ROOT.kBlue-3) + dy_el.SetLineColor(ROOT.kRed-3) + legend.AddEntry(legendHistData_mu,"Data #rightarrow #mu^{+}#mu^{-}","pe") + legend.AddEntry(legendHistData_el,"Data #rightarrow e^{+}e^{-}", "pe") + legend.AddEntry(dy_mu, "MC Inclusive #rightarrow #mu^{+}#mu^{-}", "l") + legend.AddEntry(dy_el, "MC Inclusive #rightarrow e^{+}e^{-}", "l") + + if args.useall: + for i in range(3): + for process in reversed(processes_mu[i]): + if not plot_mu["default"].muon and "#mu^{+}#mu^{-}" in process.label: + process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+}e^{-}") + process.theColor = ROOT.kBlue + process.theLineColor = ROOT.kBlue + temphist = ROOT.TH1F() + temphist.SetFillColor(process.theColor) + + for process in reversed(processes_el[i]): + if not plot_el["default"].muon and "#mu^{+}#mu^{-}" in process.label: + process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+}e^{-}") + process.theColor = ROOT.kRed + process.theLineColor = ROOT.kRed + temphist = ROOT.TH1F() + temphist.SetFillColor(process.theColor) + else: + for process in reversed(processes_mu): + if not plot_mu["default"].muon and "#mu^{+}#mu^{-}" in process.label: + process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+}e^{-}") + process.theColor = ROOT.kBlue + process.theLineColor = ROOT.kBlue + temphist = ROOT.TH1F() + temphist.SetFillColor(process.theColor) + + for process in reversed(processes_el): + if not plot_el["default"].muon and "#mu^{+}#mu^{-}" in process.label: + process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+}e^{-}") + process.theColor = ROOT.kRed + process.theLineColor = ROOT.kRed + temphist = ROOT.TH1F() + temphist.SetFillColor(process.theColor) + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + logScale = plot_mu["default"].log + + if logScale == True: + plotPad.SetLogy() + + if args.use2016: + lumi_el = 35.9*1000 + lumi_mu = 36.3*1000 + elif args.use2018: + lumi_el = 59.97*1000 + lumi_mu = 61.608*1000 + elif args.useall: + lumi_el = [35.9*1000,41.529*1000,59.97*1000] + lumi_mu = [36.3*1000,42.135*1000,61.608*1000] + else: + lumi_el = 41.529*1000 + lumi_mu = 42.135*1000 + if args.use2016: + zScaleFac_mu = zScale2016["muons"] + if "BE" in plot_el["default"].fileName: + zScaleFac_el = zScale2016["electrons"][0] + elif "BB" in plot_el["default"].fileName: + zScaleFac_el = zScale2016["electrons"][1] + elif "BE" in plot_el["default"].fileName: + zScaleFac_el = zScale2016["electrons"][2] + else: + zScaleFac_el = zScale2016["electrons"][0] + elif args.use2018: + zScaleFac_mu = zScale2018["muons"] + if "BBBE" in plot_el["default"].fileName: + zScaleFac_el = zScale2018["electrons"][0] + elif "BB" in plot_el["default"].fileName: + zScaleFac_el = zScale2018["electrons"][1] + elif "BE" in plot_el["default"].fileName: + zScaleFac_el = zScale2018["electrons"][2] + else: + zScaleFac_el = zScale2018["electrons"][0] + + elif args.useall: + zScaleFac_mu = [zScale2016["muons"],zScale["muons"],zScale2018["muons"]] + if "BBBE" in plot_el["default"].fileName: + zScaleFac_el = [zScale2016["electrons"][0],zScale["electrons"][0],zScale2018["electrons"][0]] + elif "BB" in plot_el["default"].fileName: + zScaleFac_el = [zScale2016["electrons"][1],zScale["electrons"][1],zScale2018["electrons"][1]] + elif "BE" in plot_el["default"].fileName: + zScaleFac_el = [zScale2016["electrons"][2],zScale["electrons"][2],zScale2018["electrons"][2]] + else: + zScaleFac_el = [zScale2016["electrons"][0],zScale["electrons"][0],zScale2018["electrons"][0]] + else: + zScaleFac_mu = zScale["muons"] + if "BBBE" in plot_el["default"].fileName: + zScaleFac_el = zScale["electrons"][0] + elif "BB" in plot_el["default"].fileName: + zScaleFac_el = zScale["electrons"][1] + elif "BE" in plot_el["default"].fileName: + zScaleFac_el = zScale["electrons"][2] + else: + zScaleFac_el = zScale["electrons"][0] + + + + # Data and background loading + if args.useall: + datamu=[] + datael=[] + for i in range(3): + datamu.append(data_all[i].loadHistogram(plot_mu["default"],lumi_mu[i],zScaleFac_mu[i])) + datael.append(data_all[i].loadHistogram(plot_el["default"],lumi_el[i],zScaleFac_el[i])) + stackmu = Stacks(processes_mu,lumi_mu,plot_mu["default"],zScaleFac_mu) + mu_scaleup=Stacks(processes_mu,lumi_mu,plot_mu["scale_up"],zScaleFac_mu) + mu_scaledown=Stacks(processes_mu,lumi_mu,plot_mu["scale_down"],zScaleFac_mu) + mu_ID=Stacks(processes_mu,lumi_mu,plot_mu["ID"],zScaleFac_mu) + mu_reso=Stacks(processes_mu,lumi_mu,plot_mu["reso"],zScaleFac_mu) + stackel = Stacks(processes_el,lumi_el,plot_el["default"],zScaleFac_el) + el_scaleup=Stacks(processes_el,lumi_el,plot_el["scale_up"],zScaleFac_el) + el_scaledown=Stacks(processes_el,lumi_el,plot_el["scale_down"],zScaleFac_el) + el_PUup=Stacks(processes_el,lumi_el,plot_el["PU_up"],zScaleFac_el) + el_PUdown=Stacks(processes_el,lumi_el,plot_el["PU_down"],zScaleFac_el) + else: + datamu = data_mu.loadHistogram(plot_mu["default"],lumi_mu,zScaleFac_mu) + datael = data_el.loadHistogram(plot_el["default"],lumi_el,zScaleFac_el) + stackmu = TheStack(processes_mu,lumi_mu,plot_mu["default"],zScaleFac_mu) + mu_scaleup=TheStack(processes_mu,lumi_mu,plot_mu["scale_up"],zScaleFac_mu) + mu_scaledown=TheStack(processes_mu,lumi_mu,plot_mu["scale_down"],zScaleFac_mu) + mu_ID=TheStack(processes_mu,lumi_mu,plot_mu["ID"],zScaleFac_mu) + mu_reso=TheStack(processes_mu,lumi_mu,plot_mu["reso"],zScaleFac_mu) + + stackel = TheStack(processes_el,lumi_el,plot_el["default"],zScaleFac_el) + el_scaleup=TheStack(processes_el,lumi_el,plot_el["scale_up"],zScaleFac_el) + el_scaledown=TheStack(processes_el,lumi_el,plot_el["scale_down"],zScaleFac_el) + el_PUup=TheStack(processes_el,lumi_el,plot_el["PU_up"],zScaleFac_el) + el_PUdown=TheStack(processes_el,lumi_el,plot_el["PU_down"],zScaleFac_el) + + if args.ae: + if args.useall: + i=0 + for year in range(2016,2019): + for h in stackmu[i].theStack.GetHists(): inverseAE(h, plot_mu["default"], year) + for h in stackel[i].theStack.GetHists(): inverseAE(h, plot_el["default"], year) + inverseAE(stackmu[i].theHistogram, plot_mu["default"], year) + inverseAE(stackel[i].theHistogram, plot_el["default"], year) + inverseAE(mu_scaleup[i].theHistogram, plot_mu["scale_up"], year) + inverseAE(mu_scaledown[i].theHistogram, plot_mu["scale_down"], year) + inverseAE(mu_ID[i].theHistogram, plot_mu["ID"], year) + inverseAE(mu_reso[i].theHistogram, plot_mu["reso"], year) + inverseAE(el_scaleup[i].theHistogram, plot_el["scale_up"], year) + inverseAE(el_scaledown[i].theHistogram, plot_el["scale_down"], year) + inverseAE(el_PUup[i].theHistogram, plot_el["PU_up"], year) + inverseAE(el_PUdown[i].theHistogram, plot_el["PU_down"], year) + inverseAE(datamu[i], plot_mu["default"], year) + inverseAE(datael[i], plot_el["default"], year) + i+=1 + else: + if args.use2016: year = 2016 + elif args.use2018: year = 2018 + else: year =2017 + for h in stackmu.theStack.GetHists(): inverseAE(h, plot_mu["default"], year) + for h in stackel.theStack.GetHists(): inverseAE(h, plot_el["default"], year) + inverseAE(stackmu.theHistogram, plot_mu["default"], year) + inverseAE(stackel.theHistogram, plot_el["default"], year) + inverseAE(mu_scaleup.theHistogram, plot_mu["scale_up"], year) + inverseAE(mu_scaledown.theHistogram, plot_mu["scale_down"], year) + inverseAE(mu_ID.theHistogram, plot_mu["ID"], year) + inverseAE(mu_reso.theHistogram, plot_mu["reso"], year) + inverseAE(el_scaleup.theHistogram, plot_el["scale_up"], year) + inverseAE(el_scaledown.theHistogram, plot_el["scale_down"], year) + inverseAE(el_PUup.theHistogram, plot_el["PU_up"], year) + inverseAE(el_PUdown.theHistogram, plot_el["PU_down"], year) + inverseAE(datamu, plot_mu["default"], year) + inverseAE(datael, plot_el["default"], year) + + if args.znorm: + if args.useall: + for i in range(3): + znum=znorm(stackmu[i].theHistogram) + stackmu[i].theHistogram.Scale(1./znum) + znum=znorm(stackel[i].theHistogram) + stackel[i].theHistogram.Scale(1./znum) + znum=znorm(mu_scaleup[i].theHistogram) + mu_scaleup[i].theHistogram.Scale(1./znum) + znum=znorm(mu_scaledown[i].theHistogram) + mu_scaledown[i].theHistogram.Scale(1./znum) + znum=znorm(mu_ID[i].theHistogram) + mu_ID[i].theHistogram.Scale(1./znum) + znum=znorm(mu_reso[i].theHistogram) + mu_reso[i].theHistogram.Scale(1./znum) + znum=znorm(el_scaleup[i].theHistogram) + el_scaleup[i].theHistogram.Scale(1./znum) + znum=znorm(el_scaledown[i].theHistogram) + el_scaledown[i].theHistogram.Scale(1./znum) + znum=znorm(el_PUup[i].theHistogram) + el_PUup[i].theHistogram.Scale(1./znum) + znum=znorm(el_PUdown[i].theHistogram) + el_PUdown[i].theHistogram.Scale(1./znum) + znum=znorm(datamu[i]) + datamu[i].Scale(1./znum) + znum=znorm(datael[i]) + datael[i].Scale(1./znum) + else: + znum=znorm(stackmu.theHistogram) + stackmu.theHistogram.Scale(1./znum) + znum=znorm(stackel.theHistogram) + stackel.theHistogram.Scale(1./znum) + znum=znorm(mu_scaleup.theHistogram) + mu_scaleup.theHistogram.Scale(1./znum) + znum=znorm(mu_scaledown.theHistogram) + mu_scaledown.theHistogram.Scale(1./znum) + znum=znorm(mu_ID.theHistogram) + mu_ID.theHistogram.Scale(1./znum) + znum=znorm(mu_reso.theHistogram) + mu_reso.theHistogram.Scale(1./znum) + znum=znorm(el_scaleup.theHistogram) + el_scaleup.theHistogram.Scale(1./znum) + znum=znorm(el_scaledown.theHistogram) + el_scaledown.theHistogram.Scale(1./znum) + znum=znorm(el_PUup.theHistogram) + el_PUup.theHistogram.Scale(1./znum) + znum=znorm(el_PUdown.theHistogram) + el_PUdown.theHistogram.Scale(1./znum) + znum=znorm(datamu) + datamu.Scale(1./znum) + znum=znorm(datael) + datael.Scale(1./znum) + if args.useall: + i=0 + Errs_mu=[] + Errs_el=[] + for year in range(2016,2019): + lis_mu=[[mu_scaleup[i].theHistogram,mu_scaledown[i].theHistogram],mu_ID[i].theHistogram,mu_reso[i].theHistogram] + key_mu=["massScale","ID","resolution"] + lis_el=[[el_scaleup[i].theHistogram,el_scaledown[i].theHistogram],[el_PUup[i].theHistogram,el_PUdown[i].theHistogram]] + key_el=["massScale","pileUp"] + Errs_mu.append(getErrors(stackmu[i].theHistogram,lis_mu,key_mu)) + Errs_el.append(getErrors(stackel[i].theHistogram,lis_el,key_el)) + i+=1 + if category=="bb": + errmu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2+Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2+Errs_mu[0]["ID"]**2+(Errs_mu[1]["ID"]+Errs_mu[2]["ID"])**2 + else: + errmu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2+Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2+Errs_mu[0]["ID"]**2+Errs_mu[1]["ID"]**2+Errs_mu[2]["ID"]**2 + errel=Errs_el[0]["massScale"]**2+Errs_el[0]["pileUp"]**2+Errs_el[1]["massScale"]**2+Errs_el[1]["pileUp"]**2+Errs_el[2]["massScale"]**2+Errs_el[2]["pileUp"]**2 + scale_el=Errs_el[0]["massScale"]**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 + scale_mu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 + stackmu=Addstack(stackmu) + stackel=Addstack(stackel) + mu_scaleup=Addstack(mu_scaleup) + mu_scaledown=Addstack(mu_scaledown) + mu_ID=Addstack(mu_ID) + mu_reso=Addstack(mu_reso) + el_scaleup=Addstack(el_scaleup) + el_scaledown=Addstack(el_scaledown) + el_PUup=Addstack(el_PUup) + el_PUdown=Addstack(el_PUdown) + datamu=Addhist(datamu) + datael=Addhist(datael) + else: + lis_mu=[[mu_scaleup.theHistogram,mu_scaledown.theHistogram],mu_ID.theHistogram,mu_reso.theHistogram] + lis_el=[[el_scaleup.theHistogram,el_scaledown.theHistogram],[el_PUup.theHistogram,el_PUdown.theHistogram]] + key_mu=["massScale","ID","resolution"] + Errs_mu=getErrors(stackmu.theHistogram,lis_mu,key_mu) + errmu=Errs_mu["massScale"]**2+Errs_mu["resolution"]**2+Errs_mu["ID"]**2 + key_el=["massScale","pileUp"] + Errs_el=getErrors(stackel.theHistogram,lis_el,key_el) + errel=Errs_el["massScale"]**2+Errs_el["pileUp"]**2 + scale_el=Errs_el["massScale"]**2 + scale_mu=Errs_mu["massScale"]**2 + if args.data: + yMax = datamu.GetBinContent(datamu.GetMaximumBin()) + if "Mass" in plot_mu["default"].fileName: + yMin = 0.00001 + else: + yMin = 0.01 + xMax = datamu.GetXaxis().GetXmax() + xMin = datael.GetXaxis().GetXmin() + else: + yMax = stackmu.theHistogram.GetBinContent(datamu.GetMaximumBin()) + yMin = 0.01 + xMax = stackmu.theHistogram.GetXaxis().GetXmax() + xMin = stackmu.theHistogram.GetXaxis().GetXmin() + if plot_mu["default"].yMax == None: + if logScale: + yMax = yMax*10000 + else: + yMax = yMax*1.5 + else: yMax = plot_mu["default"].yMax + + if "Mass" in plot_mu["default"].fileName: + yMax = 20000000 + + if not plot_mu["default"].yMin == None: + yMin = plot_mu.yMin + if not plot_mu["default"].xMin == None: + xMin = plot_mu["default"].xMin + if not plot_mu["default"].xMax == None: + xMax = plot_mu["default"].xMax + + + xMin = 200 + xMax = 2000 + yMin = 1e-3 + yMax = 1e4 + + if args.ae: + yMin = 0.00001 / 40 + yMax = 200000000.0 / 40 + xMin = 200 + xMax = 2000 + yMin *= 10000 + yMax /= 10 + if args.ae: + # ~ vh = plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; %s ; %s" %("m(l^{+}l^{-}) [GeV]","3 years data")) + vh = plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; %s ; %s" %("m(l^{+}l^{-}) [GeV]","Events / GeV * 1/(acc. x eff)")) + else: + # ~ vh = plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; %s ; %s" %("m(l^{+}l^{-}) [GeV]","Lumi #times d#sigma(pp#rightarrow ll)")) + vh = plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; %s ; %s" %("m(l^{+}l^{-}) [GeV]","Events / GeV")) + vh.GetXaxis().SetMoreLogLabels() + + drawStack_mu = stackmu + drawStack_el = stackel + + + # Draw background from stack + drawStack_mu.theHistogram.SetLineColor(ROOT.kBlue-3) + drawStack_el.theHistogram.SetLineColor(ROOT.kRed-3) + drawStack_mu.theHistogram.Draw("same hist") + for i in range(1, drawStack_mu.theHistogram.GetSize()-1): + print(drawStack_mu.theHistogram.GetBinContent(i)) + print(drawStack_el.theHistogram.GetBinContent(i)) + drawStack_el.theHistogram.Draw("same hist") + + + # Draw data + datamu.SetMinimum(0.0001) + if args.data: + datamu.SetMarkerColor(ROOT.kViolet) + datael.SetMarkerColor(ROOT.kOrange) + for i in range(1, datamu.GetSize()-1): + print(datamu.GetBinContent(i)) + print(datael.GetBinContent(i)) + datamu.Draw("samep") + datael.Draw("samep") + + # Draw legend + if "Eta" in plot_mu["default"].fileName or "CosTheta" in plot_mu["default"].fileName: + legendEta.Draw() + else: + legend.Draw() + + plotPad.SetLogx(plot_mu["default"].logX) + if args.useall: + latex.DrawLatex(0.95, 0.96, "three years data") + else: + latex.DrawLatex(0.95, 0.96, "%.1f fb^{-1} (13 TeV, #mu#mu), %.1f fb^{-1} (13 TeV, ee)"%(lumi_mu*0.001, lumi_el*0.001)) + yLabelPos = 0.85 + cmsExtra = "Private Work" + if not args.data: + cmsExtra = "#splitline{Private Work}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + + if args.ratio: + ratioPad.cd() + ratioPad.SetLogx(plot_mu["default"].logX) + + hhmu = drawStack_mu.theHistogram + hhel = drawStack_el.theHistogram + h_mu_scaleup=mu_scaleup.theHistogram + h_mu_scaledown=mu_scaledown.theHistogram + h_mu_ID=mu_ID.theHistogram + h_mu_reso=mu_reso.theHistogram + h_el_scaleup=el_scaleup.theHistogram + h_el_scaledown=el_scaledown.theHistogram + h_el_PUup=el_PUup.theHistogram + h_el_PUdown=el_PUdown.theHistogram + ratioGraphs = ROOT.TGraphAsymmErrors(hhmu.GetSize()-2) + chann = True if "BB" in plot_mu["default"].fileName else False + #print(errel) + #print(errmu) + for i in range(1, hhmu.GetSize()-1): + xval = hhmu.GetBinCenter(i) + xerr = hhmu.GetBinWidth(i)/2 + if hhel.GetBinContent(i) == 0: continue + if hhmu.GetBinContent(i) == 0: continue + print(hhmu.GetBinContent(i)) + print(hhel.GetBinContent(i)) + yval = hhmu.GetBinContent(i)*1.0/hhel.GetBinContent(i) + yerr = yval*math.sqrt((errel[i-1]**0.5/hhel.GetBinContent(i))**2+(errmu[i-1]**0.5/hhmu.GetBinContent(i))**2+(hhmu.GetBinError(i)/hhmu.GetBinContent(i))**2+(hhel.GetBinError(i)/hhel.GetBinContent(i))**2) + ratioGraphs.SetPoint(i, xval, yval) + ratioGraphs.SetPointError(i, xerr, xerr, yerr, yerr) + print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) + ratioData = ROOT.TGraphAsymmErrors(datamu.GetSize()-2) + for i in range(1, datamu.GetSize()-1): + xval = datamu.GetBinCenter(i) + xerr = datamu.GetBinWidth(i)/2 + if datael.GetBinContent(i) == 0: continue + if datamu.GetBinContent(i) == 0: continue + #print(datael.GetBinContent(i)) + #print(datael.GetBinError(i)) + yval = datamu.GetBinContent(i)*1.0/datael.GetBinContent(i) + yerr = yval*math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))**2+scale_mu[i-1]/(datamu.GetBinContent(i)**2)+scale_el[i-1]/(datael.GetBinContent(i)**2)) + ratioData.SetPoint(i, xval, yval) + ratioData.SetPointError(i, xerr, xerr, yerr, yerr) + print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) + nBinsX = 20 + nBinsY = 10 + if args.ae: + hAxis = ROOT.TH2F("hAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 2.5) + else: + hAxis = ROOT.TH2F("hAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 5) + hAxis.Draw("AXIS") + + hAxis.GetYaxis().SetNdivisions(408) + hAxis.SetTitleOffset(0.4, "Y") + hAxis.SetTitleSize(0.09, "Y") + hAxis.SetTitleSize(0.06, "X") + hAxis.SetYTitle("R_{#mu#mu/ee}") + hAxis.SetXTitle("m(l^{+}l^{-}) [GeV]") + hAxis.GetXaxis().SetLabelSize(0.048) + hAxis.GetYaxis().SetLabelSize(0.048) + #hAxis.GetXaxis().SetTicks("+") + #hAxis.SetTitleSize(0.15, "Y") + hAxis.GetXaxis().SetMoreLogLabels() + oneLine = ROOT.TLine(xMin, 1.0, xMax, 1.0) + oneLine.SetLineStyle(2) + oneLine.Draw() + + ratioGraphs.SetFillColor(ROOT.kBlue-3) + ratioGraphs.SetMarkerColor(ROOT.kBlue-3) + ratioGraphs.GetXaxis().SetLabelSize(0.0) + ratioGraphs.SetFillStyle(3002) + ratioGraphs.Draw("same p") + ratioData.SetMarkerColor(ROOT.kViolet) + ratioData.Draw("same p") + + rlegend = TLegend(0.2, 0.65, 0.5, 0.925) + rlegend.SetFillStyle(0) + rlegend.SetBorderSize(1) + rlegend.SetTextFont(42) + rlegend.AddEntry(ratioGraphs, "e/mu in MC inclusive", "pe") + rlegend.AddEntry(ratioData, "e/mu in data", "pe") + rlegend.Draw("same") + + ratioPad.Update() + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if args.ratio: + ratioPad.RedrawAxis() + if not os.path.exists("lepFlavor"): + os.makedirs("lepFlavor") + + if args.use2016: year = "2016" + elif args.useall: year = "2016_to_2018" + elif args.use2018: year = "2018" + else: year = "2017" + + if args.ae: year += "_inverseAE" + if args.znorm: year += "_znorm" + + hCanvas.Print("lepFlavor/%s_%s_datamc.pdf"%(plot_mu["default"].fileName, year)) + if args.datadiviMC: + ers=0 + dmCanvas = TCanvas("hCanvas", "Distribution", 800,800) + ratioDataMC = ROOT.TGraphAsymmErrors(datamu.GetSize()-2) + for i in range(1, datamu.GetSize()-1): + xval = datamu.GetBinCenter(i) + xerr = datamu.GetBinWidth(i)/2 + if datael.GetBinContent(i) == 0: continue + if datamu.GetBinContent(i) == 0: continue + dataval = datamu.GetBinContent(i)*1.0/datael.GetBinContent(i) + dataerr = math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))) + if hhel.GetBinContent(i) == 0: continue + if hhmu.GetBinContent(i) == 0: continue + MCval = hhmu.GetBinContent(i)*1.0/hhel.GetBinContent(i) + MCerr = math.sqrt((errel[i-1]**0.5/hhel.GetBinContent(i))**2+(errmu[i-1]**0.5/hhmu.GetBinContent(i))**2+(hhmu.GetBinError(i)/hhmu.GetBinContent(i))**2+(hhel.GetBinError(i)/hhel.GetBinContent(i))**2) + yval = dataval/MCval + yerr=yval*math.sqrt((MCerr)**2+(dataerr)**2) + ratioDataMC.SetPoint(i, xval, yval) + ratioDataMC.SetPointError(i, xerr, xerr, yerr, yerr) + print ("M = %f, r+-e = %f +- %f"%(xval, yval, abs(yval-1)/yerr)) + + + #print("max err") + #print(ers) + dmPad = ROOT.TPad("dmPad","dmPad",0,0,1,1) + dmPad.SetLogx(True) + #setTDRStyle() + #dmPad.UseCurrentStyle() + dmPad.Draw() + dmPad.cd() + dmAxis = ROOT.TH2F("dmAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 2.5) + dmAxis.Draw("AXIS") + dmAxis.GetYaxis().SetNdivisions(408) + dmAxis.SetTitleOffset(1.1, "Y") + dmAxis.SetTitleOffset(1.2, "X") + dmAxis.SetTitleSize(0.06, "Y") + dmAxis.SetTitleSize(0.04, "X") + dmAxis.SetYTitle("Data/MC") + dmAxis.SetXTitle("m(l^{+}l^{-}) [GeV]") + dmAxis.GetXaxis().SetLabelSize(0.048) + dmAxis.GetYaxis().SetLabelSize(0.048) + #hAxis.GetXaxis().SetTicks("+") + #hAxis.SetTitleSize(0.15, "Y") + dmAxis.GetXaxis().SetMoreLogLabels() + dmLine = ROOT.TLine(xMin, 1.0, xMax, 1.0) + dmLine.SetLineStyle(2) + dmLine.Draw() + + ratioDataMC.SetFillColor(ROOT.kBlue-3) + ratioDataMC.SetMarkerColor(ROOT.kBlue-3) + ratioDataMC.GetXaxis().SetLabelSize(0.0) + ratioDataMC.SetFillStyle(3002) + ratioDataMC.Draw("SAME p") + dmPad.Update() + dmPad.RedrawAxis() + dmCanvas.Print("lepFlavor/%s_%s_datadivimc.pdf"%(plot_mu["default"].fileName, year)) + +if __name__ == "__main__": + + + parser = argparse.ArgumentParser(description='Process some integers.') + + parser.add_argument("-d", "--data", action="store_true", dest="data", default=False, + help="plot data points.") + parser.add_argument("-m", "--mc", action="store_true", dest="mc", default=False, + help="plot mc backgrounds.") + parser.add_argument("-p", "--plot", dest="plot", nargs=1, default="", + help="plot to plot.") + parser.add_argument("-n", "--norm", action="store_true", dest="norm", default=False, + help="normalize to data.") + parser.add_argument("-2016", "--2016", action="store_true", dest="use2016", default=False, + help="use 2016 data and MC.") + parser.add_argument("-2018", "--2018", action="store_true", dest="use2018", default=False, + help="use 2018 data with 2018 MC.") + parser.add_argument("-r", "--ratio", action="store_true", dest="ratio", default=False, + help="plot ratio plot") + parser.add_argument("-l", "--log", action="store_true", dest="log", default=False, + help="plot with log scale for y axis") + parser.add_argument("-b", "--backgrounds", dest="backgrounds", action="append", default=[], + help="backgrounds to plot.") + parser.add_argument("--ae", action="store_true", dest="ae", default=False,help="times inverse Acceptance x Efficiency") + parser.add_argument("--znorm", action="store_true", dest="znorm", default=False, help="normalize to z peak") + parser.add_argument("--all", action="store_true", dest="useall", default=False, help="add the data from 2016 to 2018") + parser.add_argument("-c", action="store_true", dest="datadiviMC", default=False, help="get data/MC") + args = parser.parse_args() + if len(args.backgrounds) == 0: + args.backgrounds = ["Wjets","Other","DrellYan"] + + + muplots_bb={"default":"massPlotBB","scale_up":"massPlotBBScaleUpNoLog","scale_down":"massPlotBBScaleDownNoLog","reso":"massPlotBBSmearNoLog","ID":"massPlotBBMuonIDNoLog"} + muplots_be={"default":"massPlotBE","scale_up":"massPlotBEScaleUpNoLog","scale_down":"massPlotBEScaleDownNoLog","reso":"massPlotBESmearNoLog","ID":"massPlotBEMuonIDNoLog"} + eleplots_bb={"default":"massPlotEleBB","scale_up":"massPlotEleBBScaleUpNoLog","scale_down":"massPlotEleBBScaleDownNoLog","PU_up":"massPlotEleBBPUScaleUpNoLog","PU_down":"massPlotEleBBPUScaleDownNoLog"} + eleplots_be={"default":"massPlotEleBE","scale_up":"massPlotEleBEScaleUpNoLog","scale_down":"massPlotEleBEScaleDownNoLog","PU_up":"massPlotEleBEPUScaleUpNoLog","PU_down":"massPlotEleBEPUScaleDownNoLog"} + plot_mu_bb={} + plot_el_bb={} + plot_mu_be={} + plot_el_be={} + for key in muplots_bb.keys(): + plot_mu_bb[key] = getPlot(muplots_bb[key]) + plot_mu_bb[key].logX=True + for key in eleplots_bb.keys(): + plot_el_bb[key] = getPlot(eleplots_bb[key]) + plot_el_bb[key].logX=True + for key in muplots_be.keys(): + plot_mu_be[key] = getPlot(muplots_be[key]) + plot_mu_be[key].logX=True + for key in eleplots_be.keys(): + plot_el_be[key] = getPlot(eleplots_be[key]) + plot_el_be[key].logX=True + + plotDataMC(args,plot_mu_bb,plot_el_bb,"bb") + plotDataMC(args,plot_mu_be,plot_el_be,"be") From d7adec0f866740c19115182e6e7498bfc37161d7 Mon Sep 17 00:00:00 2001 From: "[Minxi" <[yhzz598@gmail.com]> Date: Fri, 27 Mar 2020 19:02:30 -0400 Subject: [PATCH 2/4] update --- plotAllMCFlavor.py | 574 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 473 insertions(+), 101 deletions(-) diff --git a/plotAllMCFlavor.py b/plotAllMCFlavor.py index 2e5ec1c..cc4358f 100644 --- a/plotAllMCFlavor.py +++ b/plotAllMCFlavor.py @@ -13,7 +13,7 @@ from copy import copy import numpy as np import root_numpy - +from copy import deepcopy # Muon sys uncertainty (%) # as a function of mass def znorm(hist): @@ -23,7 +23,14 @@ def znorm(hist): for i in range(lowLimit+1,upLimit+1): znorm+=hist.GetBinContent(i) return znorm - + +def scale(hist1,hist2): + lowLimit=hist1.FindBin(200) + upLimit=hist1.FindBin(400) + sFac1=hist1.Integral(lowLimit,upLimit) + sFac2=hist2.Integral(lowLimit,upLimit) + print(sFac2) + return sFac1/sFac2 def getMuErr(mass, chann, norm=False): lumi = 0.0 znorm = 0.0 @@ -98,9 +105,34 @@ def getErrors(default, others, keys): return errs # multiply hist by 1/(Acceptance x Efficiency) +def Rebin(hist): + bng=[50, 120,150,200,300,400,500,690,900,1250,1610, 3490, 6070] + #bng=([j for j in range(50, 120, 5)] + + # [j for j in range(120, 150, 5)] + + # [j for j in range(150, 200, 10)] + + # [j for j in range(200, 600, 20)]+ + # [j for j in range(600, 900, 30) ]+ + # [j for j in range(900, 1250, 50)] + + # [j for j in range(1250, 1610, 60) ] + + # [j for j in range(1610, 1890, 70) ] + + # [j for j in range(1890, 3970, 80) ] + + # [j for j in range(3970, 6070, 100) ] + + # [6070]) + #bng=[50,200,1500,6070] + hist.Sumw2() + for i in range(0,hist.GetNbinsX()): + hist.SetBinContent(i,hist.GetBinContent(i)*hist.GetBinWidth(i)) + hist.SetBinError(i,hist.GetBinError(i)*hist.GetBinWidth(i)) + hist=hist.Rebin(len(bng) - 1, 'hist_' + uuid.uuid4().hex, array('d', bng)) + for i in range(0,hist.GetNbinsX()): + hist.SetBinContent(i,hist.GetBinContent(i)/hist.GetBinWidth(i)) + hist.SetBinError(i,hist.GetBinError(i)/hist.GetBinWidth(i)) + return deepcopy(hist) + def inverseAE(hist, plotObj, year): # muon and electron # BB and BE + if year == 2016: if plotObj.muon: if "BB" in plotObj.fileName: @@ -113,6 +145,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2016 BB mu ae = %f" %(ae)) elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) @@ -123,6 +158,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2016 BE mu ae = %f" %(ae)) else: # is electron if "BB" in plotObj.fileName: for i in range(1, hist.GetSize()-1): @@ -130,12 +168,19 @@ def inverseAE(hist, plotObj, year): ae = 0.6386-497.7/(mass+348.7) + 69570.0/(mass**2+115400.0) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2016 BB e ae = %f" %(ae)) + elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) ae = -0.03377+735.1/(mass+1318)-88890.0/(mass**2+75720)+14240000.0/(mass**3+23420000) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2016 BE e ae = %f" %(ae)) elif year == 2017: if plotObj.muon: @@ -149,6 +194,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2017 BB mu ae = %f" %(ae)) elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) @@ -159,6 +207,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2017 BE mu ae = %f" %(ae)) else: # is electron if "BB" in plotObj.fileName: for i in range(1, hist.GetSize()-1): @@ -166,12 +217,18 @@ def inverseAE(hist, plotObj, year): ae = 0.585-404.6/(mass+279.5) + 56180.0/(mass**2+91430.0) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2017 BB e ae = %f" %(ae)) elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) ae = 0.02066+429.7/(mass+729)-90960.0/(mass**2+71900)+13780000.0/(mass**3+22050000) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2017 BE e ae = %f" %(ae)) elif year == 2018: if plotObj.muon: if "BB" in plotObj.fileName: @@ -184,6 +241,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2018 BB mu ae = %f" %(ae)) elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) @@ -194,6 +254,9 @@ def inverseAE(hist, plotObj, year): #print mass, ae if mass < 120: ae = float("inf") hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2018 BE mu ae = %f" %(ae)) else: # is electron if "BB" in plotObj.fileName: for i in range(1, hist.GetSize()-1): @@ -201,13 +264,18 @@ def inverseAE(hist, plotObj, year): ae = 0.576-417.7/(mass+381.8) + 46070.0/(mass**2+107200) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2018 BB e ae = %f" %(ae)) elif "BE" in plotObj.fileName: for i in range(1, hist.GetSize()-1): mass = hist.GetBinCenter(i) ae = 0.01443+475.7/(mass+639.1)-105600.0/(mass**2+82810)+12890000.0/(mass**3+23170000) #print mass, ae hist.SetBinContent(i, hist.GetBinContent(i)*1.0/ae) - + hist.SetBinError(i, hist.GetBinError(i)*1.0/ae) + #print(mass) + #print("2018 BE e ae = %f" %(ae)) def Stacks(processes,lumi,plot,zScale): stacks=[] for i in range(3): @@ -267,7 +335,7 @@ def plotDataMC(args,plot_mu,plot_el, category): if plot_mu["default"].useJets: if "Wjets" in backgrounds: backgrounds.remove("Wjets") - backgrounds.insert(0,"Jets") + backgrounds.insert(0,"Jets") processes_mu = [] processes_el = [] processes_mu2016=[] @@ -281,6 +349,9 @@ def plotDataMC(args,plot_mu,plot_el, category): if background == "Jets": processes_mu.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu,normalized=True)) processes_el.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el,normalized=True)) + elif background == "Other": + processes_el.append(Process(getattr(Backgrounds2016,"OtherEle"),eventCounts_el,negWeights_el)) + processes_mu.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) else: processes_mu.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) processes_el.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el)) @@ -301,6 +372,16 @@ def plotDataMC(args,plot_mu,plot_el, category): processes_el2018.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el,normalized=True)) processes_mu=[processes_mu2016,processes_mu2017,processes_mu2018] processes_el=[processes_mu2016,processes_mu2017,processes_mu2018] + elif background == "Other": + processes_el2016.append(Process(getattr(Backgrounds2016,"OtherEle"),eventCounts_el,negWeights_el)) + processes_mu2016.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) + processes_mu2017.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu)) + processes_el2017.append(Process(getattr(Backgrounds,background),eventCounts_el,negWeights_el)) + processes_mu2018.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu)) + processes_el2018.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el)) + processes_mu=[processes_mu2016,processes_mu2017,processes_mu2018] + processes_el=[processes_el2016,processes_el2017,processes_el2018] + else: processes_mu2016.append(Process(getattr(Backgrounds2016,background),eventCounts_mu,negWeights_mu)) processes_el2016.append(Process(getattr(Backgrounds2016,background),eventCounts_el,negWeights_el)) @@ -309,7 +390,7 @@ def plotDataMC(args,plot_mu,plot_el, category): processes_mu2018.append(Process(getattr(Backgrounds2018,background),eventCounts_mu,negWeights_mu)) processes_el2018.append(Process(getattr(Backgrounds2018,background),eventCounts_el,negWeights_el)) processes_mu=[processes_mu2016,processes_mu2017,processes_mu2018] - processes_el=[processes_mu2016,processes_mu2017,processes_mu2018] + processes_el=[processes_el2016,processes_el2017,processes_el2018] else: if background == "Jets": processes_mu.append(Process(getattr(Backgrounds,background),eventCounts_mu,negWeights_mu,normalized=True)) @@ -324,7 +405,7 @@ def plotDataMC(args,plot_mu,plot_el, category): legend.SetBorderSize(0) legend.SetTextFont(42) - + latex = ROOT.TLatex() latex.SetTextFont(42) latex.SetTextAlign(31) @@ -437,6 +518,7 @@ def plotDataMC(args,plot_mu,plot_el, category): else: lumi_el = 41.529*1000 lumi_mu = 42.135*1000 + if args.use2016: zScaleFac_mu = zScale2016["muons"] if "BE" in plot_el["default"].fileName: @@ -513,96 +595,250 @@ def plotDataMC(args,plot_mu,plot_el, category): el_PUup=TheStack(processes_el,lumi_el,plot_el["PU_up"],zScaleFac_el) el_PUdown=TheStack(processes_el,lumi_el,plot_el["PU_down"],zScaleFac_el) - if args.ae: + if args.alter: + if args.useall: + for i in range(3): + sFac=scale(datael[i],datamu[i]) + datamu[i].Scale(sFac) + sFac=scale(stackel[i].theHistogram,stackmu[i].theHistogram) + stackmu[i].theHistogram.Scale(sFac) + mu_scaleup[i].theHistogram.Scale(sFac) + mu_scaledown[i].theHistogram.Scale(sFac) + mu_ID[i].theHistogram.Scale(sFac) + mu_reso[i].theHistogram.Scale(sFac) + func=ROOT.TF1("func","[0]+[1]*log(x)+[2]*x^2",200,3500) + + i=0 + Errs_mu_fit=[] + Errs_el_fit=[] + for year in range(2016,2019): + lis_mu_fit=[[mu_scaleup[i].theHistogram,mu_scaledown[i].theHistogram],mu_ID[i].theHistogram,mu_reso[i].theHistogram] + key_mu_fit=["massScale","ID","resolution"] + lis_el_fit=[[el_scaleup[i].theHistogram,el_scaledown[i].theHistogram],[el_PUup[i].theHistogram,el_PUdown[i].theHistogram]] + key_el_fit=["massScale","pileUp"] + Errs_mu_fit.append(getErrors(stackmu[i].theHistogram,lis_mu_fit,key_mu_fit)) + Errs_el_fit.append(getErrors(stackel[i].theHistogram,lis_el_fit,key_el_fit)) + i+=1 + if category=="bb": + errmu_fit=Errs_mu_fit[0]["massScale"]**2+(Errs_mu_fit[1]["massScale"]+Errs_mu_fit[2]["massScale"])**2+Errs_mu_fit[0]["resolution"]**2+Errs_mu_fit[1]["resolution"]**2+Errs_mu_fit[2]["resolution"]**2+Errs_mu_fit[0]["ID"]**2+(Errs_mu_fit[1]["ID"]+Errs_mu_fit[2]["ID"])**2 + else: + errmu_fit=Errs_mu_fit[0]["massScale"]**2+(Errs_mu_fit[1]["massScale"]+Errs_mu_fit[2]["massScale"])**2+Errs_mu_fit[0]["resolution"]**2+Errs_mu_fit[1]["resolution"]**2+Errs_mu_fit[2]["resolution"]**2+Errs_mu_fit[0]["ID"]**2+Errs_mu_fit[1]["ID"]**2+Errs_mu_fit[2]["ID"]**2 + errel_fit=Errs_el_fit[0]["massScale"]**2+Errs_el_fit[0]["pileUp"]**2+Errs_el_fit[1]["massScale"]**2+Errs_el_fit[1]["pileUp"]**2+Errs_el_fit[2]["massScale"]**2+Errs_el_fit[2]["pileUp"]**2 + scale_el_fit=Errs_el_fit[0]["massScale"]**2+Errs_el_fit[1]["massScale"]**2+Errs_el_fit[2]["massScale"]**2 + scale_mu_fit=Errs_mu_fit[0]["massScale"]**2+(Errs_mu_fit[1]["massScale"]+Errs_mu_fit[2]["massScale"])**2 + stackmu_fit=Addstack(stackmu) + stackel_fit=Addstack(stackel) + mu_scaleup_fit=Addstack(mu_scaleup) + mu_scaledown_fit=Addstack(mu_scaledown) + mu_ID_fit=Addstack(mu_ID) + mu_reso_fit=Addstack(mu_reso) + el_scaleup_fit=Addstack(el_scaleup) + el_scaledown_fit=Addstack(el_scaledown) + el_PUup_fit=Addstack(el_PUup) + el_PUdown_fit=Addstack(el_PUdown) + datamu_fit=Addhist(datamu) + datael_fit=Addhist(datael) + if args.ratio: + hhmu_fit = stackmu_fit.theHistogram + hhel_fit = stackel_fit.theHistogram + h_mu_scaleup_fit=mu_scaleup_fit.theHistogram + h_mu_scaledown_fit=mu_scaledown_fit.theHistogram + h_mu_ID_fit=mu_ID_fit.theHistogram + h_mu_reso_fit=mu_reso_fit.theHistogram + h_el_scaleup_fit=el_scaleup_fit.theHistogram + h_el_scaledown_fit=el_scaledown_fit.theHistogram + h_el_PUup_fit=el_PUup_fit.theHistogram + h_el_PUdown_fit=el_PUdown_fit.theHistogram + ratioGraphs_fit = ROOT.TGraphAsymmErrors(hhmu_fit.GetSize()-2) + for i in range(1, hhmu_fit.GetSize()-1): + xval_fit = hhmu_fit.GetBinCenter(i) + xerr_fit = hhmu_fit.GetBinWidth(i)/2 + if hhel_fit.GetBinContent(i) == 0: continue + if hhmu_fit.GetBinContent(i) == 0: continue + #print(hhmu.GetBinContent(i)) + #print(hhel.GetBinContent(i)) + yval_fit = hhmu_fit.GetBinContent(i)*1.0/hhel_fit.GetBinContent(i) + hhmu_fit.SetBinError(i,math.sqrt(errmu_fit[i-1]+hhmu_fit.GetBinError(i)**2)) + hhel_fit.SetBinError(i,math.sqrt(errel_fit[i-1]+hhel_fit.GetBinError(i)**2)) + yerr_fit = yval_fit*math.sqrt((hhmu_fit.GetBinError(i)/hhmu_fit.GetBinContent(i))**2+(hhel_fit.GetBinError(i)/hhel_fit.GetBinContent(i))**2) + ratioGraphs_fit.SetPoint(i, xval_fit, yval_fit) + ratioGraphs_fit.SetPointError(i, xerr_fit, xerr_fit, yerr_fit, yerr_fit) + #func.SetParLimits(2,0.8,1.2) + #func.SetParLimits(1,1.5e-1,1.9e-1) + ratioGraphs_fit.Fit(func,"","",200,3500) + + #ratioGraphs_fit.Fit(func,"","",200,3500) + #ratioGraphs_fit.Fit(func,"","",200,3500) + #ratioGraphs_fit.Fit(func,"","",200,3500) + #ratioGraphs_fit.Fit(func,"","",200,3500) + print(func.GetNDF()) + print(func.GetChisquare()) + print(func.GetChisquare()/func.GetNDF()) + + if args.ae: + if args.useall: i=0 for year in range(2016,2019): - for h in stackmu[i].theStack.GetHists(): inverseAE(h, plot_mu["default"], year) - for h in stackel[i].theStack.GetHists(): inverseAE(h, plot_el["default"], year) + print(year) + for h in stackmu[i].theStack.GetHists(): + inverseAE(h, plot_mu["default"], year) + h=Rebin(h) + for h in stackel[i].theStack.GetHists(): + inverseAE(h, plot_el["default"], year) + h=Rebin(h) inverseAE(stackmu[i].theHistogram, plot_mu["default"], year) + stackmu[i].theHistogram=Rebin(stackmu[i].theHistogram) inverseAE(stackel[i].theHistogram, plot_el["default"], year) + stackel[i].theHistogram=Rebin(stackel[i].theHistogram) inverseAE(mu_scaleup[i].theHistogram, plot_mu["scale_up"], year) + mu_scaleup[i].theHistogram=Rebin(mu_scaleup[i].theHistogram) inverseAE(mu_scaledown[i].theHistogram, plot_mu["scale_down"], year) + mu_scaledown[i].theHistogram=Rebin(mu_scaledown[i].theHistogram) inverseAE(mu_ID[i].theHistogram, plot_mu["ID"], year) + mu_ID[i].theHistogram=Rebin(mu_ID[i].theHistogram) inverseAE(mu_reso[i].theHistogram, plot_mu["reso"], year) + mu_reso[i].theHistogram=Rebin(mu_reso[i].theHistogram) inverseAE(el_scaleup[i].theHistogram, plot_el["scale_up"], year) + el_scaleup[i].theHistogram=Rebin(el_scaleup[i].theHistogram) inverseAE(el_scaledown[i].theHistogram, plot_el["scale_down"], year) + el_scaledown[i].theHistogram=Rebin(el_scaledown[i].theHistogram) inverseAE(el_PUup[i].theHistogram, plot_el["PU_up"], year) + el_PUup[i].theHistogram=Rebin(el_PUup[i].theHistogram) inverseAE(el_PUdown[i].theHistogram, plot_el["PU_down"], year) + el_PUdown[i].theHistogram=Rebin(el_PUdown[i].theHistogram) inverseAE(datamu[i], plot_mu["default"], year) + datamu[i]=Rebin(datamu[i]) inverseAE(datael[i], plot_el["default"], year) + datael[i]=Rebin(datael[i]) i+=1 else: if args.use2016: year = 2016 elif args.use2018: year = 2018 else: year =2017 - for h in stackmu.theStack.GetHists(): inverseAE(h, plot_mu["default"], year) - for h in stackel.theStack.GetHists(): inverseAE(h, plot_el["default"], year) + for h in stackmu.theStack.GetHists(): + inverseAE(h, plot_mu["default"], year) + h=Rebin(h) + for h in stackel.theStack.GetHists(): + inverseAE(h, plot_el["default"], year) + h=Rebin(h) inverseAE(stackmu.theHistogram, plot_mu["default"], year) + stackmu.theHistogram=Rebin(stackmu.theHistogram) inverseAE(stackel.theHistogram, plot_el["default"], year) + stackel.theHistogram=Rebin(stackel.theHistogram) inverseAE(mu_scaleup.theHistogram, plot_mu["scale_up"], year) + mu_scaleup.theHistogram=Rebin(mu_scaleup.theHistogram) inverseAE(mu_scaledown.theHistogram, plot_mu["scale_down"], year) + mu_scaledown.theHistogram=Rebin(mu_scaledown.theHistogram) inverseAE(mu_ID.theHistogram, plot_mu["ID"], year) + mu_ID.theHistogram=Rebin(mu_ID.theHistogram) inverseAE(mu_reso.theHistogram, plot_mu["reso"], year) + mu_reso.theHistogram=Rebin(mu_reso.theHistogram) inverseAE(el_scaleup.theHistogram, plot_el["scale_up"], year) + el_scaleup.theHistogram=Rebin(el_scaleup.theHistogram) inverseAE(el_scaledown.theHistogram, plot_el["scale_down"], year) + el_scaledown.theHistogram=Rebin(el_scaledown.theHistogram) inverseAE(el_PUup.theHistogram, plot_el["PU_up"], year) + el_PUup.theHistogram=Rebin(el_PUup.theHistogram) inverseAE(el_PUdown.theHistogram, plot_el["PU_down"], year) + el_PUdown.theHistogram=Rebin(el_PUdown.theHistogram) inverseAE(datamu, plot_mu["default"], year) + datamu=Rebin(datamu) inverseAE(datael, plot_el["default"], year) - - if args.znorm: - if args.useall: - for i in range(3): - znum=znorm(stackmu[i].theHistogram) - stackmu[i].theHistogram.Scale(1./znum) - znum=znorm(stackel[i].theHistogram) - stackel[i].theHistogram.Scale(1./znum) - znum=znorm(mu_scaleup[i].theHistogram) - mu_scaleup[i].theHistogram.Scale(1./znum) - znum=znorm(mu_scaledown[i].theHistogram) - mu_scaledown[i].theHistogram.Scale(1./znum) - znum=znorm(mu_ID[i].theHistogram) - mu_ID[i].theHistogram.Scale(1./znum) - znum=znorm(mu_reso[i].theHistogram) - mu_reso[i].theHistogram.Scale(1./znum) - znum=znorm(el_scaleup[i].theHistogram) - el_scaleup[i].theHistogram.Scale(1./znum) - znum=znorm(el_scaledown[i].theHistogram) - el_scaledown[i].theHistogram.Scale(1./znum) - znum=znorm(el_PUup[i].theHistogram) - el_PUup[i].theHistogram.Scale(1./znum) - znum=znorm(el_PUdown[i].theHistogram) - el_PUdown[i].theHistogram.Scale(1./znum) - znum=znorm(datamu[i]) - datamu[i].Scale(1./znum) - znum=znorm(datael[i]) - datael[i].Scale(1./znum) - else: - znum=znorm(stackmu.theHistogram) - stackmu.theHistogram.Scale(1./znum) - znum=znorm(stackel.theHistogram) - stackel.theHistogram.Scale(1./znum) - znum=znorm(mu_scaleup.theHistogram) - mu_scaleup.theHistogram.Scale(1./znum) - znum=znorm(mu_scaledown.theHistogram) - mu_scaledown.theHistogram.Scale(1./znum) - znum=znorm(mu_ID.theHistogram) - mu_ID.theHistogram.Scale(1./znum) - znum=znorm(mu_reso.theHistogram) - mu_reso.theHistogram.Scale(1./znum) - znum=znorm(el_scaleup.theHistogram) - el_scaleup.theHistogram.Scale(1./znum) - znum=znorm(el_scaledown.theHistogram) - el_scaledown.theHistogram.Scale(1./znum) - znum=znorm(el_PUup.theHistogram) - el_PUup.theHistogram.Scale(1./znum) - znum=znorm(el_PUdown.theHistogram) - el_PUdown.theHistogram.Scale(1./znum) - znum=znorm(datamu) - datamu.Scale(1./znum) - znum=znorm(datael) - datael.Scale(1./znum) + datael=Rebin(datael) + #else: + #if args.useall: + #i=0 + #for year in range(2016,2019): + #for h in stackmu[i].theStack.GetHists(): + # h=Rebin(h) + #for h in stackel[i].theStack.GetHists(): + # h=Rebin(h) + #stackmu[i].theHistogram=Rebin(stackmu[i].theHistogram) + #stackel[i].theHistogram=Rebin(stackel[i].theHistogram) + #mu_scaleup[i].theHistogram=Rebin(mu_scaleup[i].theHistogram) + #mu_scaledown[i].theHistogram=Rebin(mu_scaledown[i].theHistogram) + #mu_ID[i].theHistogram=Rebin(mu_ID[i].theHistogram) + #mu_reso[i].theHistogram=Rebin(mu_reso[i].theHistogram) + #el_scaleup[i].theHistogram=Rebin(el_scaleup[i].theHistogram) + #el_scaledown[i].theHistogram=Rebin(el_scaledown[i].theHistogram) + #el_PUup[i].theHistogram=Rebin(el_PUup[i].theHistogram) + #el_PUdown[i].theHistogram=Rebin(el_PUdown[i].theHistogram) + #datamu[i]=Rebin(datamu[i]) + #datael[i]=Rebin(datael[i]) + #i+=1 + #else: + #if args.use2016: year = 2016 + #elif args.use2018: year = 2018 + #else: year =2017 + #for h in stackmu.theStack.GetHists(): + # h=Rebin(h) + #for h in stackel.theStack.GetHists(): + # h=Rebin(h) + #stackmu.theHistogram=Rebin(stackmu.theHistogram) + #stackel.theHistogram=Rebin(stackel.theHistogram) + #mu_scaleup.theHistogram=Rebin(mu_scaleup.theHistogram) + #mu_scaledown.theHistogram=Rebin(mu_scaledown.theHistogram) + #mu_ID.theHistogram=Rebin(mu_ID.theHistogram) + #mu_reso.theHistogram=Rebin(mu_reso.theHistogram) + #el_scaleup.theHistogram=Rebin(el_scaleup.theHistogram) + #el_scaledown.theHistogram=Rebin(el_scaledown.theHistogram) + #el_PUup.theHistogram=Rebin(el_PUup.theHistogram) + #el_PUdown.theHistogram=Rebin(el_PUdown.theHistogram) + #datamu=Rebin(datamu) + #datael=Rebin(datael) + #if args.znorm: + # if args.useall: + # for i in range(3): + # znum=znorm(stackmu[i].theHistogram) + # stackmu[i].theHistogram.Scale(1./znum) + # znum=znorm(stackel[i].theHistogram) + # stackel[i].theHistogram.Scale(1./znum) + # znum=znorm(mu_scaleup[i].theHistogram) + # mu_scaleup[i].theHistogram.Scale(1./znum) + # znum=znorm(mu_scaledown[i].theHistogram) + # mu_scaledown[i].theHistogram.Scale(1./znum) + # znum=znorm(mu_ID[i].theHistogram) + # mu_ID[i].theHistogram.Scale(1./znum) + # znum=znorm(mu_reso[i].theHistogram) + # mu_reso[i].theHistogram.Scale(1./znum) + # znum=znorm(el_scaleup[i].theHistogram) + # el_scaleup[i].theHistogram.Scale(1./znum) + # znum=znorm(el_scaledown[i].theHistogram) + # el_scaledown[i].theHistogram.Scale(1./znum) + # znum=znorm(el_PUup[i].theHistogram) + # el_PUup[i].theHistogram.Scale(1./znum) + # znum=znorm(el_PUdown[i].theHistogram) + # el_PUdown[i].theHistogram.Scale(1./znum) + # znum=znorm(datamu[i]) + # datamu[i].Scale(1./znum) + # znum=znorm(datael[i]) + # datael[i].Scale(1./znum) + # else: + # znum=znorm(stackmu.theHistogram) + # stackmu.theHistogram.Scale(1./znum) + # znum=znorm(stackel.theHistogram) + # stackel.theHistogram.Scale(1./znum) + # znum=znorm(mu_scaleup.theHistogram) + # mu_scaleup.theHistogram.Scale(1./znum) + # znum=znorm(mu_scaledown.theHistogram) + # mu_scaledown.theHistogram.Scale(1./znum) + # znum=znorm(mu_ID.theHistogram) + # mu_ID.theHistogram.Scale(1./znum) + # znum=znorm(mu_reso.theHistogram) + # mu_reso.theHistogram.Scale(1./znum) + # znum=znorm(el_scaleup.theHistogram) + # el_scaleup.theHistogram.Scale(1./znum) + # znum=znorm(el_scaledown.theHistogram) + # el_scaledown.theHistogram.Scale(1./znum) + # znum=znorm(el_PUup.theHistogram) + # el_PUup.theHistogram.Scale(1./znum) + # znum=znorm(el_PUdown.theHistogram) + # el_PUdown.theHistogram.Scale(1./znum) + # znum=znorm(datamu) + # datamu.Scale(1./znum) + # znum=znorm(datael) + # datael.Scale(1./znum) if args.useall: i=0 Errs_mu=[] @@ -617,9 +853,21 @@ def plotDataMC(args,plot_mu,plot_el, category): i+=1 if category=="bb": errmu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2+Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2+Errs_mu[0]["ID"]**2+(Errs_mu[1]["ID"]+Errs_mu[2]["ID"])**2 + ErrMu={} + ErrMu["massScale"]=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 + ErrMu["resolution"]=Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2 + ErrMu["ID"]=Errs_mu[0]["ID"]**2+(Errs_mu[1]["ID"]+Errs_mu[2]["ID"])**2 + else: errmu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2+Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2+Errs_mu[0]["ID"]**2+Errs_mu[1]["ID"]**2+Errs_mu[2]["ID"]**2 + ErrMu={} + ErrMu["massScale"]=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 + ErrMu["resolution"]=Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2 + ErrMu["ID"]=Errs_mu[0]["ID"]**2+Errs_mu[1]["ID"]**2+Errs_mu[2]["ID"]**2 errel=Errs_el[0]["massScale"]**2+Errs_el[0]["pileUp"]**2+Errs_el[1]["massScale"]**2+Errs_el[1]["pileUp"]**2+Errs_el[2]["massScale"]**2+Errs_el[2]["pileUp"]**2 + ErrEl={} + ErrEl["massScale"]=Errs_el[0]["massScale"]**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 + ErrEl["pileUp"]=Errs_el[0]["pileUp"]**2+Errs_el[1]["pileUp"]**2+Errs_el[2]["pileUp"]**2 scale_el=Errs_el[0]["massScale"]**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 scale_mu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 stackmu=Addstack(stackmu) @@ -677,7 +925,7 @@ def plotDataMC(args,plot_mu,plot_el, category): xMin = 200 - xMax = 2000 + xMax = 3490 yMin = 1e-3 yMax = 1e4 @@ -685,7 +933,7 @@ def plotDataMC(args,plot_mu,plot_el, category): yMin = 0.00001 / 40 yMax = 200000000.0 / 40 xMin = 200 - xMax = 2000 + xMax = 3490 yMin *= 10000 yMax /= 10 if args.ae: @@ -703,24 +951,55 @@ def plotDataMC(args,plot_mu,plot_el, category): # Draw background from stack drawStack_mu.theHistogram.SetLineColor(ROOT.kBlue-3) drawStack_el.theHistogram.SetLineColor(ROOT.kRed-3) - drawStack_mu.theHistogram.Draw("same hist") - for i in range(1, drawStack_mu.theHistogram.GetSize()-1): - print(drawStack_mu.theHistogram.GetBinContent(i)) - print(drawStack_el.theHistogram.GetBinContent(i)) - drawStack_el.theHistogram.Draw("same hist") + drawStack_mu.theHistogram.SetMarkerStyle(6) + drawStack_el.theHistogram.SetMarkerStyle(6) + drawStack_mu.theHistogram.SetFillStyle(0) + drawStack_el.theHistogram.SetFillStyle(0) + #for i in range(1, drawStack_mu.theHistogram.GetSize()-1): + # val=drawStack_mu.theHistogram.GetBinContent(i) + # drawStack_mu.theHistogram.SetBinError(i,math.sqrt(errmu[i-1]+drawStack_mu.theHistogram.GetBinError(i)**2)) + # drawStack_mu.theHistogram.SetBinContent(i,val) + # val=drawStack_el.theHistogram.GetBinContent(i) + # drawStack_el.theHistogram.SetBinError(i,math.sqrt(errel[i-1]+drawStack_el.theHistogram.GetBinError(i)**2)) + # drawStack_el.theHistogram.SetBinContent(i,val) + # val=datamu.GetBinContent(i) + # datamu.SetBinError(i,math.sqrt(scale_mu[i-1]+datamu.GetBinError(i)**2)) + # datamu.SetBinContent(i,val) + # val=datael.GetBinContent(i) + # datael.SetBinError(i,math.sqrt(scale_el[i-1]+datael.GetBinError(i)**2)) + # datael.SetBinContent(i,val) + #drawStack_mu.theHistogram.SetLineWidth(2) + #drawStack_el.theHistogram.SetLineWidth(2) + drawStack_mu.theHistogram.Draw("same histe ") + #for i in range(1, drawStack_mu.theHistogram.GetSize()-1): + #print(drawStack_mu.theHistogram.GetBinContent(i)) + #print(drawStack_el.theHistogram.GetBinContent(i)) + drawStack_el.theHistogram.Draw("same histe ") + drawStack_mu.theHistogram.SetLineColor(ROOT.kBlue) + drawStack_el.theHistogram.SetLineColor(ROOT.kRed) + drawStack_mu.theHistogram.SetFillStyle(0) + drawStack_el.theHistogram.SetFillStyle(0) # Draw data datamu.SetMinimum(0.0001) + Box = ROOT.TBox(0.39,0.79,0.5,0.9) + Box.SetLineColor(ROOT.kRed) + Box.SetFillStyle(0) + #Box.SetFillAlphaColor(0,0) + Box.SetLineWidth(2) if args.data: datamu.SetMarkerColor(ROOT.kViolet) datael.SetMarkerColor(ROOT.kOrange) - for i in range(1, datamu.GetSize()-1): - print(datamu.GetBinContent(i)) - print(datael.GetBinContent(i)) - datamu.Draw("samep") - datael.Draw("samep") - + #for i in range(1, datamu.GetSize()-1): + #print(datamu.GetBinContent(i)) + #print(datael.GetBinContent(i)) + datamu.Draw("ESAMEp") + datael.Draw("ESAMEp") + Box.Draw("same") + #a1=datael.Integral(60,120) + #a2=drawStack_el.theHistogram.Integral(60,120) + #print(a1/a2) # Draw legend if "Eta" in plot_mu["default"].fileName or "CosTheta" in plot_mu["default"].fileName: legendEta.Draw() @@ -729,17 +1008,38 @@ def plotDataMC(args,plot_mu,plot_el, category): plotPad.SetLogx(plot_mu["default"].logX) if args.useall: - latex.DrawLatex(0.95, 0.96, "three years data") + latex.DrawLatex(0.95, 0.96, "%.1f fb^{-1} (13 TeV, #mu#mu), %.1f fb^{-1} (13 TeV, ee)"%((lumi_mu[0]+lumi_mu[1]+lumi_mu[2])*0.001, (lumi_el[0]+lumi_el[1]+lumi_el[2])*0.001)) else: latex.DrawLatex(0.95, 0.96, "%.1f fb^{-1} (13 TeV, #mu#mu), %.1f fb^{-1} (13 TeV, ee)"%(lumi_mu*0.001, lumi_el*0.001)) yLabelPos = 0.85 - cmsExtra = "Private Work" - if not args.data: - cmsExtra = "#splitline{Private Work}{Simulation}" - yLabelPos = 0.82 + cmsExtra = "Preliminary" + #if not args.data: + #cmsExtra = "#splitline{Private Work}{Simulation}" + #yLabelPos = 0.82 + if category=="bb": + if args.useall: + dataLabel = "Full Run 2 BB" + elif args.use2016: + dataLabel = "2016 BB" + elif args.use2018: + dataLabel = "2018 BB" + else: + dataLabel = "2017 BB" + else: + if args.useall: + dataLabel = "Full Run 2 BE" + elif args.use2016: + dataLabel = "2016 BE" + elif args.use2018: + dataLabel = "2018 BE" + else: + dataLabel = "2017 BE" + pave=ROOT.TPaveText(0.39,0.79,0.5,0.9) + pave.AddBox(0.39,0.79,0.5,0.9) + pave.Draw() latexCMS.DrawLatex(0.19,0.89,"CMS") latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) - + latexCMSExtra.DrawLatex(0.39,0.79,"%s"%(dataLabel)) if args.ratio: ratioPad.cd() ratioPad.SetLogx(plot_mu["default"].logX) @@ -758,18 +1058,39 @@ def plotDataMC(args,plot_mu,plot_el, category): chann = True if "BB" in plot_mu["default"].fileName else False #print(errel) #print(errmu) + #print(hhmu.GetSize()) for i in range(1, hhmu.GetSize()-1): xval = hhmu.GetBinCenter(i) xerr = hhmu.GetBinWidth(i)/2 if hhel.GetBinContent(i) == 0: continue if hhmu.GetBinContent(i) == 0: continue - print(hhmu.GetBinContent(i)) - print(hhel.GetBinContent(i)) + muScale=math.sqrt(ErrMu["massScale"][i-1])/hhmu.GetBinContent(i) + muReso=math.sqrt(ErrMu["resolution"][i-1])/hhmu.GetBinContent(i) + muID=math.sqrt(ErrMu["ID"][i-1])/hhmu.GetBinContent(i) + muStat=hhmu.GetBinError(i)/hhmu.GetBinContent(i) + elScale=math.sqrt(ErrEl["massScale"][i-1])/hhel.GetBinContent(i) + elPile=math.sqrt(ErrEl["pileUp"][i-1])/hhel.GetBinContent(i) + elStat=hhel.GetBinError(i)/hhel.GetBinContent(i) + print "mass scale %f, resolution %f, ID %f, stat %f, mass %f, dimuon channel" %(muScale, muReso, muID, muStat, xval) + print "mass scale %f, pile up %f, stat %f, mass %f, dielectron channel" %(elScale, elPile, elStat, xval) + #print(hhmu.GetBinContent(i)) + #print(hhel.GetBinContent(i)) yval = hhmu.GetBinContent(i)*1.0/hhel.GetBinContent(i) - yerr = yval*math.sqrt((errel[i-1]**0.5/hhel.GetBinContent(i))**2+(errmu[i-1]**0.5/hhmu.GetBinContent(i))**2+(hhmu.GetBinError(i)/hhmu.GetBinContent(i))**2+(hhel.GetBinError(i)/hhel.GetBinContent(i))**2) + hhmu.SetBinError(i,math.sqrt(errmu[i-1]+hhmu.GetBinError(i)**2)) + hhel.SetBinError(i,math.sqrt(errel[i-1]+hhel.GetBinError(i)**2)) + + yerr = yval*math.sqrt((hhmu.GetBinError(i)/hhmu.GetBinContent(i))**2+(hhel.GetBinError(i)/hhel.GetBinContent(i))**2) ratioGraphs.SetPoint(i, xval, yval) - ratioGraphs.SetPointError(i, xerr, xerr, yerr, yerr) - print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) + ratioGraphs.SetPointError(i, xerr, xerr, yerr, yerr) + if args.alter: + #h_fit=deepcopy(hhmu) + #hhmu.Sumw2() + #hhel.Sumw2() + #h_fit.Divide(hhmu,hhel) + #func1=ROOT.TF1("func1","[0]*log(x)+[1]",200,3500) + #ratioGraphs.Fit(func1,"","",200,3500) + #print(func1.GetChisquare()/func1.GetNDF()) + print("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) ratioData = ROOT.TGraphAsymmErrors(datamu.GetSize()-2) for i in range(1, datamu.GetSize()-1): xval = datamu.GetBinCenter(i) @@ -780,15 +1101,17 @@ def plotDataMC(args,plot_mu,plot_el, category): #print(datael.GetBinError(i)) yval = datamu.GetBinContent(i)*1.0/datael.GetBinContent(i) yerr = yval*math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))**2+scale_mu[i-1]/(datamu.GetBinContent(i)**2)+scale_el[i-1]/(datael.GetBinContent(i)**2)) + print(yval) + print(yerr) ratioData.SetPoint(i, xval, yval) ratioData.SetPointError(i, xerr, xerr, yerr, yerr) - print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) + print("M = %f, r+-e = %f +- %f"%(xval, yval, yerr/yval)) nBinsX = 20 nBinsY = 10 if args.ae: hAxis = ROOT.TH2F("hAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 2.5) else: - hAxis = ROOT.TH2F("hAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 5) + hAxis = ROOT.TH2F("hAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 2.5) hAxis.Draw("AXIS") hAxis.GetYaxis().SetNdivisions(408) @@ -811,15 +1134,19 @@ def plotDataMC(args,plot_mu,plot_el, category): ratioGraphs.GetXaxis().SetLabelSize(0.0) ratioGraphs.SetFillStyle(3002) ratioGraphs.Draw("same p") + #func.Draw("same") ratioData.SetMarkerColor(ROOT.kViolet) ratioData.Draw("same p") - + if args.alter: + Expression="f(x)=p_{0}+p_{1}lnx+p_{2}x^{2}" + latexCMSExtra.DrawLatex(0.59,yLabelPos,"%s"%(Expression)) + func.Draw("same") rlegend = TLegend(0.2, 0.65, 0.5, 0.925) rlegend.SetFillStyle(0) rlegend.SetBorderSize(1) rlegend.SetTextFont(42) - rlegend.AddEntry(ratioGraphs, "e/mu in MC inclusive", "pe") - rlegend.AddEntry(ratioData, "e/mu in data", "pe") + rlegend.AddEntry(ratioGraphs, "mumu/ee in MC inclusive", "pe") + rlegend.AddEntry(ratioData, "mumu/ee in data", "pe") rlegend.Draw("same") ratioPad.Update() @@ -838,7 +1165,7 @@ def plotDataMC(args,plot_mu,plot_el, category): else: year = "2017" if args.ae: year += "_inverseAE" - if args.znorm: year += "_znorm" + #if args.znorm: year += "_znorm" hCanvas.Print("lepFlavor/%s_%s_datamc.pdf"%(plot_mu["default"].fileName, year)) if args.datadiviMC: @@ -851,16 +1178,49 @@ def plotDataMC(args,plot_mu,plot_el, category): if datael.GetBinContent(i) == 0: continue if datamu.GetBinContent(i) == 0: continue dataval = datamu.GetBinContent(i)*1.0/datael.GetBinContent(i) - dataerr = math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))) + dataerr = math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))**2+scale_mu[i-1]/(datamu.GetBinContent(i)**2)+scale_el[i-1]/(datael.GetBinContent(i)**2)) if hhel.GetBinContent(i) == 0: continue if hhmu.GetBinContent(i) == 0: continue MCval = hhmu.GetBinContent(i)*1.0/hhel.GetBinContent(i) MCerr = math.sqrt((errel[i-1]**0.5/hhel.GetBinContent(i))**2+(errmu[i-1]**0.5/hhmu.GetBinContent(i))**2+(hhmu.GetBinError(i)/hhmu.GetBinContent(i))**2+(hhel.GetBinError(i)/hhel.GetBinContent(i))**2) yval = dataval/MCval yerr=yval*math.sqrt((MCerr)**2+(dataerr)**2) - ratioDataMC.SetPoint(i, xval, yval) ratioDataMC.SetPointError(i, xerr, xerr, yerr, yerr) - print ("M = %f, r+-e = %f +- %f"%(xval, yval, abs(yval-1)/yerr)) + print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr)) + if args.alter: + ers=0 + dmCanvas = TCanvas("hCanvas", "Distribution", 800,800) + ratioDataMC = ROOT.TGraphAsymmErrors(datamu.GetSize()-2) + for i in range(1, datamu.GetSize()-1): + xval=datamu.GetBinCenter(i) + fval=func.Eval(xval) + print "At the bin center with the mass %f, scaling function value is %f, dimuon event yield is %f, dielectron event yield is %f"%(datamu.GetBinCenter(i), fval, datamu.GetBinContent(i)*datamu.GetBinWidth(i), datael.GetBinContent(i)*datael.GetBinWidth(i)) + muval=datamu.GetBinContent(i)/func.Eval(xval) + ParErr0=func.GetParError(0) + ParErr1=func.GetParError(1) + ParErr2=func.GetParError(2) + fErr=math.sqrt(ParErr0**2+(ParErr1*math.log(xval))**2+(ParErr2*xval**2)**2) + #fErr=0 + #err=math.sqrt(datamu.GetBinError(i)**2+scale_mu[i-1]+muval**2*fErr**2)/fval + err=muval*fErr/fval + datamu.SetBinContent(i,muval) + datamu.SetBinError(i,err) + #elerr=math.sqrt(datael.GetBinError(i)**2+scale_el[i-1]) + #datael.SetBinError(i,elerr) + datamu.Sumw2() + datamu=Rebin(datamu) + datael.Sumw2() + datael=Rebin(datael) + for i in range(1, datamu.GetSize()-1): + xval=datamu.GetBinCenter(i) + xerr=datamu.GetBinWidth(i) + if datael.GetBinContent(i) == 0: continue + if datamu.GetBinContent(i) == 0: continue + yval=datamu.GetBinContent(i)/datael.GetBinContent(i) + yerr=yval*math.sqrt((datamu.GetBinError(i)/datamu.GetBinContent(i))**2+(datael.GetBinError(i)/datael.GetBinContent(i))**2) + ratioDataMC.SetPoint(i, xval, yval) + ratioDataMC.SetPointError(i, xerr/2, xerr/2, yerr, yerr) + print ("M = %f, r+-e = %f +- %f"%(xval, yval, yerr)) #print("max err") @@ -869,6 +1229,7 @@ def plotDataMC(args,plot_mu,plot_el, category): dmPad.SetLogx(True) #setTDRStyle() #dmPad.UseCurrentStyle() + Box.Draw("same") dmPad.Draw() dmPad.cd() dmAxis = ROOT.TH2F("dmAxis", "", nBinsX, xMin, xMax, nBinsY, 0.5, 2.5) @@ -878,7 +1239,10 @@ def plotDataMC(args,plot_mu,plot_el, category): dmAxis.SetTitleOffset(1.2, "X") dmAxis.SetTitleSize(0.06, "Y") dmAxis.SetTitleSize(0.04, "X") - dmAxis.SetYTitle("Data/MC") + if args.datadiviMC: + dmAxis.SetYTitle("R^{data}_{#mu#mu/ee}/R^{mc}_{#mu#mu/ee}") + else: + dmAxis.SetYTitle("R^{data}_{#mu#mu/ee}") dmAxis.SetXTitle("m(l^{+}l^{-}) [GeV]") dmAxis.GetXaxis().SetLabelSize(0.048) dmAxis.GetYaxis().SetLabelSize(0.048) @@ -896,6 +1260,14 @@ def plotDataMC(args,plot_mu,plot_el, category): ratioDataMC.Draw("SAME p") dmPad.Update() dmPad.RedrawAxis() + latexCMS.DrawLatex(0.19,0.89,"CMS") + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + latexCMSExtra.DrawLatex(0.59,yLabelPos,"%s"%(dataLabel)) + + if args.useall: + latex.DrawLatex(0.95, 0.96, "%.1f fb^{-1} (13 TeV, #mu#mu), %.1f fb^{-1} (13 TeV, ee)"%((lumi_mu[0]+lumi_mu[1]+lumi_mu[2])*0.001, (lumi_el[0]+lumi_el[1]+lumi_el[2])*0.001)) + else: + latex.DrawLatex(0.95, 0.96, "%.1f fb^{-1} (13 TeV, #mu#mu), %.1f fb^{-1} (13 TeV, ee)"%(lumi_mu*0.001, lumi_el*0.001)) dmCanvas.Print("lepFlavor/%s_%s_datadivimc.pdf"%(plot_mu["default"].fileName, year)) if __name__ == "__main__": @@ -922,13 +1294,13 @@ def plotDataMC(args,plot_mu,plot_el, category): parser.add_argument("-b", "--backgrounds", dest="backgrounds", action="append", default=[], help="backgrounds to plot.") parser.add_argument("--ae", action="store_true", dest="ae", default=False,help="times inverse Acceptance x Efficiency") - parser.add_argument("--znorm", action="store_true", dest="znorm", default=False, help="normalize to z peak") + parser.add_argument("--alter", action="store_true", dest="alter", default=False, help="Another approach") parser.add_argument("--all", action="store_true", dest="useall", default=False, help="add the data from 2016 to 2018") parser.add_argument("-c", action="store_true", dest="datadiviMC", default=False, help="get data/MC") args = parser.parse_args() if len(args.backgrounds) == 0: - args.backgrounds = ["Wjets","Other","DrellYan"] - + #args.backgrounds = ["Wjets","Other","DrellYan"] + args.backgrounds = ["DrellYan","Other"] muplots_bb={"default":"massPlotBB","scale_up":"massPlotBBScaleUpNoLog","scale_down":"massPlotBBScaleDownNoLog","reso":"massPlotBBSmearNoLog","ID":"massPlotBBMuonIDNoLog"} muplots_be={"default":"massPlotBE","scale_up":"massPlotBEScaleUpNoLog","scale_down":"massPlotBEScaleDownNoLog","reso":"massPlotBESmearNoLog","ID":"massPlotBEMuonIDNoLog"} From 8b7e2b02324c0a474e5cdfc4b974a965689d97d2 Mon Sep 17 00:00:00 2001 From: "[Minxi" <[yhzz598@gmail.com]> Date: Tue, 29 Sep 2020 17:04:05 -0400 Subject: [PATCH 3/4] update unfolding --- .plotAllMCFlavor.py.swo | Bin 0 -> 24576 bytes RooUnfold | 1 + dataCollection.py | 1017 +++++++++++++++++++++++++++++++++++++++ defs.py | 4 +- helpers.py | 61 +-- histHandle.py | 116 +++++ plotAllMCFlavor.py | 20 +- plotTreeStudy.py | 567 ++++++++++++++++++++++ run.sh | 11 + unfoldDataPlot.py | 251 ++++++++++ unfoldDataPlot_V2.py | 190 ++++++++ unfoldPlot.py | 303 ++++++++++++ unfoldPlot_V2.py | 196 ++++++++ unfoldTestPlot.py | 214 ++++++++ unfoldingMC.py | 189 ++++++++ unfoldingSample.py | 260 ++++++++++ unfoldingSample_V2.py | 266 ++++++++++ 17 files changed, 3626 insertions(+), 40 deletions(-) create mode 100644 .plotAllMCFlavor.py.swo create mode 160000 RooUnfold create mode 100644 dataCollection.py create mode 100644 histHandle.py create mode 100644 plotTreeStudy.py create mode 100644 run.sh create mode 100755 unfoldDataPlot.py create mode 100755 unfoldDataPlot_V2.py create mode 100755 unfoldPlot.py create mode 100755 unfoldPlot_V2.py create mode 100755 unfoldTestPlot.py create mode 100644 unfoldingMC.py create mode 100644 unfoldingSample.py create mode 100644 unfoldingSample_V2.py diff --git a/.plotAllMCFlavor.py.swo b/.plotAllMCFlavor.py.swo new file mode 100644 index 0000000000000000000000000000000000000000..da2a7c1617cf877fa9c8ef26e52897f9f10cf687 GIT binary patch literal 24576 zcmeI4dyFJUea9PYlWgKgu#J)+1jP&|&+N|hPS5Vn?p|&bZg%&~qP@G#?cN!&KA(GM zYIggsr+eJpbGLJcZ5}d4M>KTUrNQEyce_aKH;s2w;9TomSD*Qzi{$>^aY%2UpD!ix)6zYFLg}+tB_fz3F zD*vBR&u3EMH!1(OsOLYV!hfLj->br>Zjdiw`v0xM*HrkUsc_}geo=*gB^9o`vqxru z%mSGOG7DrD$SjaqAhSSbfy@G#1u_d{7Pu}8m@Q5FCB!VzdKTyZ;rjo3w`tn9!BgN# zV1f-W4PFgy1Z%fy+8aR)FmMZa@fJTS#B$IK2_Y5kv6jHs2WzOZPNt9o_OfQe?OrS%r7 zi?e#MjH;?sjOy&XFq$+=vL2KMqq+!Dqgq>tYeCjy*)SIMYH_++)o134NM>dst^jDC zsTqqyO{`XVsjEO#Kx1qLAULd-PS15vC-CWW(_{5It5%E})l!}{8q9QBtg#|oBq`j{ zC)VO^RPlN4`TX!ofvU6=>MaR%Xb*YT>DhtRwfmxBvF$!HrVEFgP2LSmr^(qKTWPng zrp29Rzlcu2<3Z1J823D%jsH>AkjU0K-alQ(oARie+kOE~VPib&p6do@oL;Bf*C~Tl z)AP)J!3}n~S1h;Ma3ZGcKjPHJYInMN$K10zy$&T}l<=n-CRGj7PWbGT`j&Az&H$Cg z2;&=}Mq|I+Z&ecUQGG$b%T>ARbsI%R9nu#oidtf@6Mc>st*#XEOYrgOgGJJQ{nCTl zgs{e#BvI8S>anO%o$v5q#q(UxFDUPlOvM+;kkY6cMK|X^-~|P{*Rduit5X*p*Xz&= z>!M5Ob@9rkSgXH%Zb!a$TW$Hm9q#S){ibPi>0z}J4)NMHn+PI%2?gk+D|i=4(7R49 zN;q`DPLw=KjdF`YBUB?L=ERrK+%+BNAR7#tp*Sj%%7~0JE4-mu4Hz@{n%RH^x|G@K zLXAyIkuqhhQe8k|^xdgwi|JT~%Hrr$ksV|vj&7ji8cMcOH5M1bWT!;33s+4PI+X!p zD3O>eML6DBQXG#q|HmBdDf4ZXq%$*6g7d=gb#zs(HZBZEPOY;9v7Ys z%#FgHbDvKx`r@X-?`=wcZBw0I;92wmPI=sS?H(?{I^^o~tk!jXVGJ3f$kI}1>akww z6*g~-50td@Bn!vwY9Is+p2iZ_ku8Tgd@s0a^P#y{+;}h?4>p~tp@xG=QmWZZjc7E3 zO=D!+5KZPfnuh9MkM{O&8EBtspcI`QK5Frts*jx8%8?*3cf&`yW_7G!$|&`^>W#(% zp`?ypQ~64DM=Ym>8FT?tVSF)R+`+IKf++7ogP045_G{vsV(p&Oqz1)Y+1sH7AEn8J z-}f1UxEGk1M)z6aKoQNx>U3Q%VDZoy`iN1SL8H^v+Me5CO}C3dO^VZdDCqj(8+{8> z7J2hSS&U!nn3k4U!)%`0@!X!%!qAm)&=e@`u~1?g)(&okLsY#;v4rFB@^XnSn}I3+ z%5ZYy4e`6kjhrE^iMi6n4I0}@Og^P_^u8eTz02(`CRkzCC$*gu$B%D{HH{fGcMH@7 z6s4%ogU#g=8$sVjry#A8@Cr>LF1Y(L2B``~Hh_9SS`viX!p?*WH*C`I;tjw%*4|t8MM*t2-EK`I6h|U=q1zIsAB+8iB9+tuQ|n zDSS7*ovw-5QdqB@_q7>BLnec3g#1*eZ64BE?RTu5iZ1 z;`K{8wey-gn4|r}b=2UB&2i>Z%86~uGrev)NGvu_c8|Gh?oK4uaxXX|1(y*Jimsm= zp>F+TN+G01X;X5sKI}BdMcE!(ltz@tp`tWWDa#U%NL&iWNtbe?iqg2oqJ(m>J`^RY z<}Pyo_Byh!Jfb|@EY)C zto=U%P60Z%|2M4lKMURgPJmwn>tG2Kz%PUU!al$=;N!pphrwOoR`7kCw?7U35PSrz zf*P0xuLG|I&*7Z?+h7)yz-N&8kAg?Q1KW!Rx>koUQ*E zcsuBV`@tf3HMoL*DDXGncf)@$pz47hID*4o?e1=;r>76o=^PQq5lx)ATaFljq$szq zMRX5GQs`w9n(wh%!7k3Br*X0{F|mi;`oZv}#6k6lV``Pffn#cUc2CE6ogPMriP+6; za)&$b;U0SqpGSD0jS+h*>VZ#3XWKa3z!8Ahs)qux(RR0Ic6VvZOYDtp3g=(_o-tGuImU}h zbwY5IMhr&O|J2oUSPueLmnnF4d1PRe1(Jyz<3*95P;f+kur~_JwmkEK*q?;ihBTq( z2Sm_T2Mih36ZSl5PZdMsUO|n@$r`4q6g3J_rzG2K72rnmr&*jXyv3oSZ3cXA>FCDF zUf`MfGP1TtJ~+%W7EAp#bDP`iu3yR>@N)+W_#c!Axgu1B5z}#SU1YJ8gGS8p;sffz zw@?pCn!_&gJVG(Pcg48lXNmy!}chSD~!6A`QXG$YZ%P3zmbSRQI z%oNs?dT|FqN)G=`>OIo4ZFZvPuqx!BeJD~Qx2A2QlB5;TjU?5CW^3e7%i|=3Nt8rF ztu(^0LQ&oZ>Oq~<4i+vYKUfZtNFoWDb9gk-ns!{TU$4`&pXTEEs8yl*(ORg7%qNSM z_q(e-Ne4L%ZV!vO%}7M1k@FJ;}(BPbEu- z*`f^(G~mU#67sZ&!<5_s-DtD+_I;Pi#)GV2%&<*Jf09{wwq*yv!+i}+6Hl zu1?`3y~Ko*+tNwjj;Jp(^G_=SVGeNTVoGHO*Egcq>gi&bHz~>OND{B1Bwiy)%JxwG zQpa;m9T62(oWuxYQ%@8&S?RqJ$i7GQQ8~D;I(A{rms@iE&0T(VVQLzTIyWZMqb4V= z!WB$o+7yFW-mKQg_RRV!?v-I6vkSFKNuDZBxsu*F8n%+6?ExESq6~vt2 z+;MT8j`XR>Bc)DRuFY}(CTcY#jCx4Ku^e$RiN2Bw-62v{L=rSx8($kqPKB6PE z(jk7b=Q_b#<(aGLPe*FtqLqlHx~J{v4h9pKI()inAkP2KV$JzotU+o0AO8FHFJQfY z68s43{Kvq3;7)KEYyA%b3zWeRvDW`0_yG7#@P6`#LwjZz~Ih9wu6%Gp!^iJ zaIzf~{Og%qd{mI_ponGZsHK0jh)pbd10: +# k=11 +# else: +# k=i +# if j>10: +# l=11 +# else: +# l=j +# val1=hist.GetBinContent(i,j) +# val2=tmp.GetBinContent(k,l) +# tmp.SetBinContent(k,l,val1+val2) +# +# hist=tmp.Clone() +#histContainer.processes(reBin2D,["response"]) +#histContainer.processes(reBin1D,["recMass","data","bkg"]) +#histContainer.saveHists("dataCollectionDY2016.root") + + + + + + + + + + + diff --git a/defs.py b/defs.py index 9df41fe..f428703 100755 --- a/defs.py +++ b/defs.py @@ -4,7 +4,8 @@ import sys import copy -path = "../files/" +path = "/depot/cms/users/schul105/files/" +#path="/depot/cms/users/minxi/" @@ -22,6 +23,7 @@ # ~ "muons":1.0 } + crossSections = { "dyInclusive50":5765.4, "dy50to120":2112.90, diff --git a/helpers.py b/helpers.py index c600adf..b454bd3 100755 --- a/helpers.py +++ b/helpers.py @@ -39,35 +39,35 @@ def negWeightFractions(path,muon=True): def binning(channel='muon'): - #if channel == 'muon': - - #return [60., 65.43046396, 71.3524269, 77.81037328, 84.85281374, 92.53264952, 100.90756983, 110.04048518, 120., 129.95474058, 140.73528833, 152.41014904, 165.0535115, 178.74571891, 193.57377942, 209.63191906, 227.02218049, 245.85507143, 266.2502669, 288.3373697, + if channel == 'muon' or 'electron': + #return [60,120,400,600,900,1300,1800,6000] + return [50, 120,150,200,300,400,500,690,900,1250,1610, 2000, 3000,4000, 6070]# 152.41014904, 165.0535115, 178.74571891, 193.57377942, 209.63191906, 227.02218049, 245.85507143, 266.2502669, 288.3373697, #312.25673399, 338.16035716, 366.21284574, 396.59246138, 429.49225362, 465.12128666, 503.70596789, 545.49148654, 590.74337185, 639.74918031, 692.82032303, 750.29404456, 812.53556599, 879.94040575, 952.93689296, #1031.98888927, 1117.59873655, 1210.310449, 1310.71317017, 1419.4449167, #1537.19663264, 1664.71658012, 1802.81509423, 1952.36973236, 2114.3308507, 2289.72764334, 2479.6746824, 2685.37900061, 2908.14776151, 3149.39656595, 3410.65844758, 3693.59361467, 4000. ] #if channel == 'electron': - # ~ return ([j for j in range(50, 120, 5)] + - # ~ [j for j in range(120, 150, 5)] + - # ~ [j for j in range(150, 200, 10)] + - # ~ [j for j in range(200, 600, 20)]+ - # ~ [j for j in range(600, 900, 30) ]+ - # ~ [j for j in range(900, 1250, 50)] + - # ~ [j for j in range(1250, 1600, 60) ] + - # ~ [j for j in range(1600, 1900, 70) ] + - # ~ [j for j in range(1900, 4000, 80) ] + - # ~ [j for j in range(4000, 5000, 100) ] + - # ~ [5000]) - return ([j for j in range(50, 120, 5)] + - [j for j in range(120, 150, 5)] + - [j for j in range(150, 200, 10)] + - [j for j in range(200, 600, 20)]+ - [j for j in range(600, 900, 30) ]+ - [j for j in range(900, 1250, 50)] + - [j for j in range(1250, 1610, 60) ] + - [j for j in range(1610, 1890, 70) ] + - [j for j in range(1890, 3970, 80) ] + - [j for j in range(3970, 6070, 100) ] + - [6070]) + # return ([j for j in range(50, 120, 5)] + + # [j for j in range(120, 150, 5)] + + # [j for j in range(150, 200, 10)] + + # [j for j in range(200, 600, 20)]+ + # [j for j in range(600, 900, 30) ]+ + # [j for j in range(900, 1250, 50)] + + # [j for j in range(1250, 1600, 60) ] + + # [j for j in range(1600, 1900, 70) ] + + # [j for j in range(1900, 4000, 80) ] + + # [j for j in range(4000, 5000, 100) ] + + # [5000]) + #return ([j for j in range(50, 120, 5)] + + # [j for j in range(120, 150, 5)] + + # [j for j in range(150, 200, 10)] + + # [j for j in range(200, 600, 20)]+ + # [j for j in range(600, 900, 30) ]+ + # [j for j in range(900, 1250, 50)] + + # [j for j in range(1250, 1610, 60) ] + + # [j for j in range(1610, 1890, 70) ] + + # [j for j in range(1890, 3970, 80) ] + + # [j for j in range(3970, 6070, 100) ] + + # [6070]) # Calculate logarithmic bins # ~ width = (math.log(m_max) - math.log(m_min)) / nbins # ~ logbins = [] # ~ # Exceed m_max to start with Z' binning, but reach 5 TeV @@ -75,7 +75,7 @@ def binning(channel='muon'): # ~ logbins.append(int(math.exp(math.log(m_min) + width * i))) # ~ return logbins - + # return [j for j in range(50,6071,5)] def loadHistoFromFile(fileName,histName,rebin,muon=True,logBins=False): """ returns histogram from file @@ -128,7 +128,10 @@ def loadHistoFromFile(fileName,histName,rebin,muon=True,logBins=False): bng = binning("electron") else: bng = binning("muon") - + #for i in range(0,result.GetNbinsX()): + # print(result.GetBinLowEdge(i)) + #print(histName) + #print(len(bng) - 1) result = result.Rebin(len(bng) - 1, 'hist_' + uuid.uuid4().hex, array('d', bng)) for i in range(0,result.GetNbinsX()): @@ -374,9 +377,7 @@ def __init__(self,processes,lumi,plot,zScaleFac): else: self.theHistogram.Add(temphist.Clone()) def Add(self,addstack): - for h in addstack.theStack.GetHists(): - self.theStack.Add(h.Clone()) - self.theHistogram.Add(h.Clone()) + self.theHistogram.Add(addstack.theHistogram.Clone()) class TheStackRun2: from ROOT import THStack diff --git a/histHandle.py b/histHandle.py new file mode 100644 index 0000000..5a65c75 --- /dev/null +++ b/histHandle.py @@ -0,0 +1,116 @@ +import ROOT +from copy import deepcopy + + +class histHandle(object): + + def __init__(self,name,layers): + self.name=name + self.layers=layers + self.nodes={} + self.__tree={} + ROOT.TH1.AddDirectory(0) + ROOT.TH2.AddDirectory(0) + for layer in layers: + self.nodes[layer]=set([]) + + def __addLeaf(self,histName,nodes,h,tree="none"): + if tree=="none": + tree=self.__tree + if len(nodes)>1: + if nodes[0] not in tree.keys(): + #print nodes[0] + tree[nodes[0]]={} + self.__addLeaf(histName,nodes[1:],h,tree[nodes[0]]) + else: + if nodes[0] in tree.keys(): + #print nodes[0] + if "hist" in tree[nodes[0]].keys(): + tree[nodes[0]]["hist"].Add(h.Clone()) + else: + tree[nodes[0]]["hist"]=h.Clone(histName) + else: + tree[nodes[0]]={"hist":h.Clone(histName)} + + def __loopOver(self,node,func,selection=[]): + flag=True + for key in node.keys(): + #print key + if key =="hist": + func(node[key]) + elif key in selection: + self.__loopOver(node[key],func,selection) + flag=False + if flag: + for key in node.keys(): + #print key + if key !="hist": + self.__loopOver(node[key],func,selection) + + def loadHist(self,filename,histName,leafPath): + f=ROOT.TFile.Open(filename) + hist=f.Get(histName) + layerNum=len(leafPath) + for key in leafPath.keys(): + self.nodes[key].add(leafPath[key]) + subLayers=self.layers[0:layerNum] + nodes=[] + for layer in subLayers: + nodes.append(leafPath[layer]) + s="_" + histName=s.join(nodes) + #print histName + self.__addLeaf(histName,nodes,histm) + f.Close() + + def addHist(self,hist,leafPath): + layerNum=len(leafPath) + for key in leafPath.keys(): + self.nodes[key].add(leafPath[key]) + subLayers=self.layers[0:layerNum] + nodes=[] + for layer in subLayers: + nodes.append(leafPath[layer]) + s="_" + histName=s.join(nodes) + #print histName + self.__addLeaf(histName,nodes,hist) + + def loadHists(self,filename,histList): + for dual in histList: + self.loadHist(filename,dual[0],dual[1]) + + def getHist(self,leafPath): + tree=self.__tree.copy() + nodes=[] + num=len(leafPath) + layers=self.layers[0:num] + for layer in layers: + nodes.append(leafPath[layer]) + #print nodes + for node in nodes: + #print tree + tree=tree[node].copy() + h=tree["hist"].Clone() + return h + def process(self,func,leafPath): + tree=self.__tree.copy() + num=len(leafPath) + layers=self.layers[0:num] + nodes=[] + for layer in layers: + nodes.append(leafPath[layer]) + for node in nodes: + tree=tree[node].copy() + func(tree["hist"]) + def processes(self,func,selection=""): + tree=self.__tree.copy() + if selection=="": + self.__loopOver(tree,func) + else: + self.__loopOver(tree,func,selection) + def saveHists(self,filename,selection=""): + f=ROOT.TFile(filename,"RECREATE") + def saveHist(hist): + hist.Write() + self.processes(saveHist,selection) diff --git a/plotAllMCFlavor.py b/plotAllMCFlavor.py index cc4358f..a425833 100644 --- a/plotAllMCFlavor.py +++ b/plotAllMCFlavor.py @@ -30,6 +30,7 @@ def scale(hist1,hist2): sFac1=hist1.Integral(lowLimit,upLimit) sFac2=hist2.Integral(lowLimit,upLimit) print(sFac2) + print sFac1/sFac2 return sFac1/sFac2 def getMuErr(mass, chann, norm=False): lumi = 0.0 @@ -328,6 +329,7 @@ def plotDataMC(args,plot_mu,plot_el, category): eventCounts_mu = totalNumberOfGeneratedEvents(path,plot_mu["default"].muon) eventCounts_el = totalNumberOfGeneratedEvents(path,plot_el["default"].muon) negWeights_mu = negWeightFractions(path,plot_mu["default"].muon) + print negWeights_mu negWeights_el = negWeightFractions(path,plot_el["default"].muon) # Background load processes @@ -746,7 +748,7 @@ def plotDataMC(args,plot_mu,plot_el, category): datamu=Rebin(datamu) inverseAE(datael, plot_el["default"], year) datael=Rebin(datael) - #else: + #els: #if args.useall: #i=0 #for year in range(2016,2019): @@ -864,10 +866,10 @@ def plotDataMC(args,plot_mu,plot_el, category): ErrMu["massScale"]=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 ErrMu["resolution"]=Errs_mu[0]["resolution"]**2+Errs_mu[1]["resolution"]**2+Errs_mu[2]["resolution"]**2 ErrMu["ID"]=Errs_mu[0]["ID"]**2+Errs_mu[1]["ID"]**2+Errs_mu[2]["ID"]**2 - errel=Errs_el[0]["massScale"]**2+Errs_el[0]["pileUp"]**2+Errs_el[1]["massScale"]**2+Errs_el[1]["pileUp"]**2+Errs_el[2]["massScale"]**2+Errs_el[2]["pileUp"]**2 + errel=Errs_el[0]["massScale"]**2+(Errs_el[0]["pileUp"]+Errs_el[1]["pileUp"]+Errs_el[2]["pileUp"])**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 ErrEl={} ErrEl["massScale"]=Errs_el[0]["massScale"]**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 - ErrEl["pileUp"]=Errs_el[0]["pileUp"]**2+Errs_el[1]["pileUp"]**2+Errs_el[2]["pileUp"]**2 + ErrEl["pileUp"]=(Errs_el[0]["pileUp"]+Errs_el[1]["pileUp"]+Errs_el[2]["pileUp"])**2 scale_el=Errs_el[0]["massScale"]**2+Errs_el[1]["massScale"]**2+Errs_el[2]["massScale"]**2 scale_mu=Errs_mu[0]["massScale"]**2+(Errs_mu[1]["massScale"]+Errs_mu[2]["massScale"])**2 stackmu=Addstack(stackmu) @@ -1066,11 +1068,11 @@ def plotDataMC(args,plot_mu,plot_el, category): if hhmu.GetBinContent(i) == 0: continue muScale=math.sqrt(ErrMu["massScale"][i-1])/hhmu.GetBinContent(i) muReso=math.sqrt(ErrMu["resolution"][i-1])/hhmu.GetBinContent(i) - muID=math.sqrt(ErrMu["ID"][i-1])/hhmu.GetBinContent(i) + muID=math.sqrt(ErrMu["ID"][i-1])/(hhmu.GetBinContent(i)) muStat=hhmu.GetBinError(i)/hhmu.GetBinContent(i) elScale=math.sqrt(ErrEl["massScale"][i-1])/hhel.GetBinContent(i) elPile=math.sqrt(ErrEl["pileUp"][i-1])/hhel.GetBinContent(i) - elStat=hhel.GetBinError(i)/hhel.GetBinContent(i) + elStat=hhel.GetBinError(i)/(hhel.GetBinContent(i)) print "mass scale %f, resolution %f, ID %f, stat %f, mass %f, dimuon channel" %(muScale, muReso, muID, muStat, xval) print "mass scale %f, pile up %f, stat %f, mass %f, dielectron channel" %(elScale, elPile, elStat, xval) #print(hhmu.GetBinContent(i)) @@ -1205,8 +1207,8 @@ def plotDataMC(args,plot_mu,plot_el, category): err=muval*fErr/fval datamu.SetBinContent(i,muval) datamu.SetBinError(i,err) - #elerr=math.sqrt(datael.GetBinError(i)**2+scale_el[i-1]) - #datael.SetBinError(i,elerr) + elerr=math.sqrt(datael.GetBinError(i)**2+scale_el[i-1]) + datael.SetBinError(i,elerr) datamu.Sumw2() datamu=Rebin(datamu) datael.Sumw2() @@ -1299,8 +1301,8 @@ def plotDataMC(args,plot_mu,plot_el, category): parser.add_argument("-c", action="store_true", dest="datadiviMC", default=False, help="get data/MC") args = parser.parse_args() if len(args.backgrounds) == 0: - #args.backgrounds = ["Wjets","Other","DrellYan"] - args.backgrounds = ["DrellYan","Other"] + args.backgrounds = ["Wjets","Other"] + #args.backgrounds = ["DrellYan"] muplots_bb={"default":"massPlotBB","scale_up":"massPlotBBScaleUpNoLog","scale_down":"massPlotBBScaleDownNoLog","reso":"massPlotBBSmearNoLog","ID":"massPlotBBMuonIDNoLog"} muplots_be={"default":"massPlotBE","scale_up":"massPlotBEScaleUpNoLog","scale_down":"massPlotBEScaleDownNoLog","reso":"massPlotBESmearNoLog","ID":"massPlotBEMuonIDNoLog"} diff --git a/plotTreeStudy.py b/plotTreeStudy.py new file mode 100644 index 0000000..004c950 --- /dev/null +++ b/plotTreeStudy.py @@ -0,0 +1,567 @@ +import ROOT +import os +import numpy +import subprocess +import root_numpy +import math +import copy + + +def loopOver(path): + pathList=[] + for filename in os.listdir(path): + pathNew=path+"/"+filename + if os.path.isfile(pathNew) and pathNew.endswith(".root"): + pathList.append(pathNew) + elif os.path.isdir(pathNew): + return loopOver(pathNew) + return pathList + +def getBkg(path,name): + pathList=[] + for filename in os.listdir(path): + pathNew=path+"/"+filename + if os.path.isdir(pathNew) and filename.startswith(name): + pathList+=loopOver(pathNew) + return pathList + + +ROOT.gStyle.SetOptStat(0) +path="/mnt/hadoop/store/user/minxi" +pathList=getBkg(path,"ZToMuMu") +pathList+=getBkg(path,"DY") +pathList2016=[] +pathList2017=[] +for path in pathList: + if(path.find("2016")==-1): + pathList2016.append(path) + else: + pathList2017.append(path) + +#chain=ROOT.TChain("SimpleNtupler") +mass2016=y2016=pt2016=CS2016=CS_label2016=numpy.ones(0) +mass2017=y2017=pt2017=CS2017=CS_label2017=numpy.ones(0) +for path in pathList2016: + #path=path.replace("/mnt/hadoop/","root://xrootd.rcac.purdue.edu//") + splitPath=path.split("/") + pathMap="/home/yang1452/lepratio/CIAnalysis/crabFile/"+splitPath[-3]+"/"+splitPath[-2] + if not os.path.exists(pathMap): + os.makedirs(pathMap) + cmd="cp "+path+" "+pathMap + os.system(cmd) + pathMap=pathMap+"/"+splitPath[-1] + f2016=ROOT.TFile.Open(pathMap) + tree2016=f2016.Get("SimpleNtupler/t") + mass2016=numpy.concatenate([mass2016,root_numpy.tree2array(tree2016,"dil_mass")]) + y2016=numpy.concatenate([y2016,root_numpy.tree2array(tree2016,"dil_rap")]) + pt2016=numpy.concatenate([pt2016,root_numpy.tree2array(tree2016,"dil_pt")]) + CS2016=numpy.concatenate([CS2016,root_numpy.tree2array(tree2016,"cos_cs")]) + CS_label2016=numpy.concatenate([CS_label2016,root_numpy.tree2array(tree2016,"CS_label")]) + + + +for path in pathList2017: + #path=path.replace("/mnt/hadoop/","root://xrootd.rcac.purdue.edu//") + splitPath=path.split("/") + pathMap="/home/yang1452/lepratio/CIAnalysis/crabFile/"+splitPath[-3]+"/"+splitPath[-2] + if not os.path.exists(pathMap): + os.makedirs(pathMap) + cmd="cp "+path+" "+pathMap + os.system(cmd) + pathMap=pathMap+"/"+splitPath[-1] + f2017=ROOT.TFile.Open(pathMap) + tree2017=f2017.Get("SimpleNtupler/t") + mass2017=numpy.concatenate([mass2017,root_numpy.tree2array(tree2017,"dil_mass")]) + y2017=numpy.concatenate([y2017,root_numpy.tree2array(tree2017,"dil_rap")]) + pt2017=numpy.concatenate([pt2017,root_numpy.tree2array(tree2017,"dil_pt")]) + CS2017=numpy.concatenate([CS2017,root_numpy.tree2array(tree2017,"cos_cs")]) + CS_label2017=numpy.concatenate([CS_label2017,root_numpy.tree2array(tree2017,"CS_label")]) + + + + + +#h_CS800to1000=ROOT.TH1D("h_CS","CS angle",100,-1.,1.) +#h_CS800to1000.SetLineColor(2) +#h_CS800to1000.GetXaxis().SetTitle("cos #theta") +#h_CS5500to6000=ROOT.TH1D("h_CS","CS angle",100,-1.,1.) +#h_CS5500to6000.SetLineColor(3) +#h_CS5500to6000T=ROOT.TH1D("h_CS","CS angle",100,-1.,1.) +#h_CS5500to6000T.SetLineColor(4) +bng=[i for i in range(200,400,20)]+[i for i in range(400,800,40)]+[i for i in range(800,1400,60)]+[i for i in range(1400,2300,90)]+[i for i in range(2300,3500,120)]+[i for i in range(3500,4500,100)]+[i for i in range(4500,6000,150)] +bng=numpy.asarray(bng,dtype=numpy.float64) +CS_2016true=copy.deepcopy(CS2016) +CS_2016true[CS_label2016==0]=-CS_2016true[CS_label2016==0] +CS_2017true=copy.deepcopy(CS2017) +CS_2017true[CS_label2017==0]=-CS_2017true[CS_label2017==0] +h_prob_yVSmassn2016=ROOT.TH2D("h_prob_yVSmassn2016","y_vs_mass",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmassd2016=ROOT.TH2D("h_prob_yVSmassd2016","y_vs_mass",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmass2016=ROOT.TH2D("h_prob_yVSmass2016","Probability to identify quark with a correct direction",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmassErr2016=ROOT.TH2D("h_prob_yVSmassErr2016","Relative Error",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmassn2017=ROOT.TH2D("h_prob_yVSmassn2017","y_vs_mass",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmassd2017=ROOT.TH2D("h_prob_yVSmassd2017","y_vs_mass",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmass2017=ROOT.TH2D("h_prob_yVSmass2017","Probability to identify quark with a correct direction",len(bng)-1,bng,30,0,2.4) +h_prob_yVSmassErr2017=ROOT.TH2D("h_prob_yVSmassErr2017","Relative Error",len(bng)-1,bng,30,0,2.4) + +h_prob_yVSmassn2016.FillN(len(y2016[CS_label2016==1]),mass2016[CS_label2016==1],numpy.absolute(y2016[CS_label2016==1]),numpy.ones(len(y2016[CS_label2016==1]))) +h_prob_yVSmassn2016.Sumw2() +h_prob_yVSmassd2016.FillN(len(y2016),mass2016,numpy.absolute(y2016),numpy.ones(len(y2016))) +h_prob_yVSmassd2016.Sumw2() +h_prob_yVSmass2016.Divide(h_prob_yVSmassn2016,h_prob_yVSmassd2016) +h_prob_yVSmassn2017.FillN(len(y2017[CS_label2017==1]),mass2017[CS_label2017==1],numpy.absolute(y2017[CS_label2017==1]),numpy.ones(len(y2017[CS_label2017==1]))) +h_prob_yVSmassn2017.Sumw2() +h_prob_yVSmassd2017.FillN(len(y2017),mass2017,numpy.absolute(y2017),numpy.ones(len(y2017))) +h_prob_yVSmassd2017.Sumw2() +h_prob_yVSmass2017.Divide(h_prob_yVSmassn2017,h_prob_yVSmassd2017) + +for i in range(0,h_prob_yVSmass2016.GetNbinsX()): + for j in range(0,h_prob_yVSmass2016.GetNbinsY()): + if h_prob_yVSmassn2016.GetBinContent(i+1,j+1)!=0: + err=h_prob_yVSmass2016.GetBinError(i+1,j+1)/h_prob_yVSmass2016.GetBinContent(i+1,j+1) + else: + err=0 + h_prob_yVSmassErr2016.SetBinContent(i+1,j+1,err) + +for i in range(0,h_prob_yVSmass2017.GetNbinsX()): + for j in range(0,h_prob_yVSmass2017.GetNbinsY()): + if h_prob_yVSmassn2017.GetBinContent(i+1,j+1)!=0: + err=h_prob_yVSmass2017.GetBinError(i+1,j+1)/h_prob_yVSmass2017.GetBinContent(i+1,j+1) + else: + err=0 + h_prob_yVSmassErr2017.SetBinContent(i+1,j+1,err) + +c1=ROOT.TCanvas("c1","c1",800,800) +c1.Divide(2,2) +c1.cd(1) +h_prob_yVSmass2016.GetXaxis().SetTitle("mll[GeV]") +h_prob_yVSmass2016.GetYaxis().SetTitle("y") +h_prob_yVSmass2016.Draw("COLZ") +c1.cd(2) +h_prob_yVSmassErr2016.GetXaxis().SetTitle("mll[GeV]") +h_prob_yVSmassErr2016.GetYaxis().SetTitle("y") +h_prob_yVSmassErr2016.Draw("COLZ") +c1.cd(3) +h_prob_yVSmass2017.GetXaxis().SetTitle("mll[GeV]") +h_prob_yVSmass2017.GetYaxis().SetTitle("y") +h_prob_yVSmass2017.Draw("COLZ") +c1.cd(4) +h_prob_yVSmassErr2017.GetXaxis().SetTitle("mll[GeV]") +h_prob_yVSmassErr2017.GetYaxis().SetTitle("y") +h_prob_yVSmassErr2017.Draw("COLZ") +c1.Print("Study/Result1.pdf") + + + + + + + + + +h_CS1_16=ROOT.TH1D("h_CS1_16","CS angle",100,-1.,1.) +h_CS2_16=ROOT.TH1D("h_CS2_16","CS angle",100,-1.,1.) +h_CS3_16=ROOT.TH1D("h_CS3_16","CS angle",100,-1.,1.) +h_CS4_16=ROOT.TH1D("h_CS4_16","CS angle",100,-1.,1.) +h_CS5_16=ROOT.TH1D("h_CS5_16","CS angle",100,-1.,1.) +h_CS1_16t=ROOT.TH1D("h_CS1_16t","CS angle",100,-1.,1.) +h_CS2_16t=ROOT.TH1D("h_CS2_16t","CS angle",100,-1.,1.) +h_CS3_16t=ROOT.TH1D("h_CS3_16t","CS angle",100,-1.,1.) +h_CS4_16t=ROOT.TH1D("h_CS4_16t","CS angle",100,-1.,1.) +h_CS5_16t=ROOT.TH1D("h_CS5_16t","CS angle",100,-1.,1.) +h_CS1_17=ROOT.TH1D("h_CS1_17","CS angle",100,-1.,1.) +h_CS2_17=ROOT.TH1D("h_CS2_17","CS angle",100,-1.,1.) +h_CS3_17=ROOT.TH1D("h_CS3_17","CS angle",100,-1.,1.) +h_CS4_17=ROOT.TH1D("h_CS4_17","CS angle",100,-1.,1.) +h_CS5_17=ROOT.TH1D("h_CS5_17","CS angle",100,-1.,1.) +h_CS1_17t=ROOT.TH1D("h_CS1_17t","CS angle",100,-1.,1.) +h_CS2_17t=ROOT.TH1D("h_CS2_17t","CS angle",100,-1.,1.) +h_CS3_17t=ROOT.TH1D("h_CS3_17t","CS angle",100,-1.,1.) +h_CS4_17t=ROOT.TH1D("h_CS4_17t","CS angle",100,-1.,1.) +h_CS5_17t=ROOT.TH1D("h_CS5_17t","CS angle",100,-1.,1.) +#h_CS800to1000.SetLineColor(2) +#h_CS800to1000.GetXaxis().SetTitle("cos #theta") +#h_CS5500to6000=ROOT.TH1D("h_CS","CS angle",100,-1.,1.) +#h_CS5500to6000.SetLineColor(3) +#h_CS5500to6000T=ROOT.TH1D("h_CS","CS angle",100,-1.,1.) +#h_CS5500to6000T.SetLineColor(4) +h_16=[h_CS1_16,h_CS2_16,h_CS3_16,h_CS4_16,h_CS5_16] +h_16t=[h_CS1_16t,h_CS2_16t,h_CS3_16t,h_CS4_16t,h_CS5_16t] +h_17=[h_CS1_17,h_CS2_17,h_CS3_17,h_CS4_17,h_CS5_17] +h_17t=[h_CS1_17t,h_CS2_17t,h_CS3_17t,h_CS4_17t,h_CS5_17t] +massLimit=[1000.,2000.,3000.,4000.,5000.,6000.] +h_CS1_16t.FillN(len(CS2016[(mass2016massLimit[0]) & (y2016<1.)]),CS2016[(mass2016massLimit[0]) & (y2016<1.)],numpy.ones(len(CS2016[(mass2016massLimit[0]) & (y2016<1.)]))) +c4=ROOT.TCanvas("c4","c4",800,800) +h_CS1_16t.Draw() +c4.Print("Study/result4.pdf") +i=0 +for h in h_16t: + h_16[i].FillN(len(CS2016[(mass2016massLimit[i]) & (y2016<1.)]),CS2016[(mass2016massLimit[i]) & (y2016<1.)],numpy.ones(len(CS2016[(mass2016massLimit[i]) & (y2016<1.)]))) + #print len(CS2016[(mass2016massLimit[i]) & (y2016<1.)]) + h_16[i].Scale(1./len(CS2016[(mass2016massLimit[i]) & (y2016<1.)])) + h.FillN(len(CS_2016true[(mass2016massLimit[i]) & (y2016<1.)]),CS_2016true[(mass2016massLimit[i]) & (y2016<1.)],numpy.ones(len(CS_2016true[(mass2016massLimit[i]) & (y2016<1.)]))) + h.Scale(1./len(CS_2016true[(mass2016massLimit[i]) & (y2016<1.)])) + h_17[i].FillN(len(CS2017[(mass2017massLimit[i]) & (y2017<1.)]),CS2017[(mass2017massLimit[i]) & (y2017<1.)],numpy.ones(len(CS2017[(mass2017massLimit[i]) & (y2017<1.)]))) + h_17[i].Scale(1./len(CS2017[(mass2017massLimit[i]) & (y2017<1.)])) + h_17t[i].FillN(len(CS_2017true[(mass2017massLimit[i]) & (y2017<1.)]),CS_2017true[(mass2017massLimit[i]) & (y2017<1.)],numpy.ones(len(CS_2017true[(mass2017massLimit[i]) & (y2017<1.)]))) + h_17t[i].Scale(1./len(CS_2017true[(mass2017massLimit[i]) & (y2017<1.)])) + i+=1 + + +c2=ROOT.TCanvas("c2","c2",800,800) +c2.Divide(2,3) +c2.cd(1) +lyt1=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt1.SetTextFont(28) +leg1=["1000 to 2000 GeV","2000 to 3000 GeV"," 3000 to 4000 GeV"," 4000 to 5000 GeV"," 5000 to 6000 GeV"] +#h_plott1=ROOT.TH1D("h_plott1","2016 truth of the CS angle",100,-1.,1.) +#h_plott1.GetXaxis().SetTitle("cos #theta") +#h_plott1.Draw() +h_16t[0].SetTitle("2016 truth of the CS angle") +for i in range(len(h_16)): + h_16t[i].SetLineColor(i+2) + lyt1.AddEntry(h_16t[i],leg1[i]) + if i>0: + h_16t[i].Draw("samehist") + else: + h_16t[i].Draw("hist") +lyt1.Draw() + +c2.cd(2) +lyt2=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt2.SetTextFont(28) +#h_plott2=ROOT.TH1D("h_plott2","2017 truth of the CS angle",100,-1.,1.) +#h_plott2.GetXaxis().SetTitle("cos #theta") +#h_plott2.Draw() +h_17t[0].SetTitle("2017 truth of the CS angle") +for i in range(len(h_17)): + h_17t[i].SetLineColor(i+2) + lyt2.AddEntry(h_17t[i],leg1[i]) + if i>0: + h_17t[i].Draw("samehist") + else: + h_17t[i].Draw("hist") +lyt2.Draw() + +c2.cd(3) +lyt3=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt3.SetTextFont(28) +#h_plott3=ROOT.TH1D("h_plott3","Comparsion of the CS truth between 2016 and 2017 at 2 TeV to 3 TeV",100,-1.,1.) +#h_plott3.GetXaxis().SetTitle("cos #theta") +h_16t[1].SetTitle("Comparsion of the CS truth between 2016 and 2017 at 2 TeV to 3 TeV") +#h_plott3.Draw() +h_16t[1].SetLineColor(2) +h_17t[1].SetLineColor(3) +lyt3.AddEntry(h_16t[1],"2016") +lyt3.AddEntry(h_17t[1],"2017") +h_16t[1].Draw("hist") +h_17t[1].Draw("samehist") +lyt3.Draw() + + +c2.cd(4) +lyt4=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt4.SetTextFont(28) +#h_plott4=ROOT.TH1D("h_plott4","Comparsion of the CS truth between 2016 and 2017 at 3 TeV to 4 TeV",100,-1.,1.) +#h_plott4.GetXaxis().SetTitle("cos #theta") +#h_plott4.Draw() +h_16t[2].SetTitle("Comparsion of the CS truth between 2016 and 2017 at 3 TeV to 4 TeV") +h_16t[2].SetLineColor(2) +h_17t[2].SetLineColor(3) +lyt4.AddEntry(h_16t[2],"2016") +lyt4.AddEntry(h_17t[2],"2017") +h_16t[2].Draw("hist") +h_17t[2].Draw("samehist") +lyt4.Draw() + +c2.cd(5) +lyt5=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt5.SetTextFont(28) +#h_plott5=ROOT.TH1D("h_plott5","Comparsion of the CS truth between 2016 and 2017 at 4 TeV to 5 TeV",100,-1.,1.) +#h_plott5.GetXaxis().SetTitle("cos #theta") +h_16t[3].SetTitle("Comparsion of the CS truth between 2016 and 2017 at 4 TeV to 5 TeV") +#h_plott5.Draw() +h_16t[3].SetLineColor(2) +h_17t[3].SetLineColor(3) +lyt5.AddEntry(h_16t[3],"2016") +lyt5.AddEntry(h_17t[3],"2017") +h_16t[3].Draw("hist") +h_17t[3].Draw("samehist") +lyt5.Draw() + +c2.cd(6) +lyt6=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +lyt6.SetTextFont(28) +#h_plott6=ROOT.TH1D("h_plott6","Comparsion of the CS truth between 2016 and 2017 at 5 TeV to 6 TeV",100,-1.,1.) +#h_plott6.GetXaxis().SetTitle("cos #theta") +h_16t[4].SetTitle("Comparsion of the CS truth between 2016 and 2017 at 5 TeV to 6 TeV") +#h_plott6.Draw() +h_16t[4].SetLineColor(2) +h_17t[4].SetLineColor(3) +#lyt6.AddEntry(h_plott6,"") +lyt6.AddEntry(h_16t[4],"2016") +lyt6.AddEntry(h_17t[4],"2017") +h_16t[4].Draw("hist") +h_17t[4].Draw("samehist") +lyt6.Draw() +c2.Update() +c2.Print("Study/Result2.pdf") + + + + +c3=ROOT.TCanvas("c3","c3",800,800) +c3.Divide(2,3) +c3.cd(1) +ly1=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly1.SetTextFont(28) +#h_plot1=ROOT.TH1D("h_plot1","2016 CS angle",100,-1.,1.) +#h_plot1.GetXaxis().SetTitle("cos #theta") +#h_plot1.Draw() +h_16[0].SetTitle("2016 CS angle") +for i in range(len(h_16)): + h_16[i].SetLineColor(i+2) + ly1.AddEntry(h_16[i],leg1[i]) + if i>0: + h_16[i].Draw("samehist") + else: + h_16[i].Draw("hist") +ly1.Draw() + +c3.cd(2) +ly2=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly2.SetTextFont(28) +#h_plot2=ROOT.TH1D("h_plot2","2017 CS angle",100,-1.,1.) +#h_plot2.GetXaxis().SetTitle("cos #theta") +#h_plot2.Draw() +h_17[0].SetTitle("2017 CS angle") +for i in range(len(h_17)): + h_17[i].SetLineColor(i+2) + ly2.AddEntry(h_17[i],leg1[i]) + if i>0: + h_17[i].Draw("samehist") + else: + h_17[i].Draw("hist") +ly2.Draw() + +c3.cd(3) +ly3=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly3.SetTextFont(28) +#h_plot3=ROOT.TH1D("h_plot3","Comparsion of the CS between 2016 and 2017 at 2 TeV to 3 TeV",100,-1.,1.) +#h_plot3.GetXaxis().SetTitle("cos #theta") +#h_plot3.Draw() +h_16[1].SetTitle("Comparsion of the CS between 2016 and 2017 at 2 TeV to 3 TeV") +h_16[1].SetLineColor(2) +h_17[1].SetLineColor(3) +ly3.AddEntry(h_16[1],"2016") +ly3.AddEntry(h_17[1],"2017") +h_16[1].Draw("hist") +h_17[1].Draw("samehist") +ly3.Draw() + + +c3.cd(4) +ly4=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly4.SetTextFont(28) +#h_plot4=ROOT.TH1D("h_plot4","Comparsion of the CS between 2016 and 2017 at 3 TeV to 4 TeV",100,-1.,1.) +#h_plot4.GetXaxis().SetTitle("cos #theta") +#h_plot4.Draw() +h_16[2].SetTitle("Comparsion of the CS between 2016 and 2017 at 3 TeV to 4 TeV") +h_16[2].SetLineColor(2) +h_17[2].SetLineColor(3) +ly4.AddEntry(h_16[2],"2016") +ly4.AddEntry(h_17[2],"2017") +h_16[2].Draw("hist") +h_17[2].Draw("samehist") +ly4.Draw() + +c3.cd(5) +ly5=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly5.SetTextFont(28) +#h_plot5=ROOT.TH1D("h_plot5","Comparsion of the CS between 2016 and 2017 at 4 TeV to 5 TeV",100,-1.,1.) +#h_plot5.GetXaxis().SetTitle("cos #theta") +#h_plot5.Draw() +h_16[3].SetTitle("Comparsion of the CS between 2016 and 2017 at 4 TeV to 5 TeV") +h_16[3].SetLineColor(2) +h_17[3].SetLineColor(3) +ly5.AddEntry(h_16[3],"2016") +ly5.AddEntry(h_17[3],"2017") +h_16[3].Draw("hist") +h_17[3].Draw("samehist") +ly5.Draw() + + +c3.cd(6) +ly6=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +ly6.SetTextFont(28) +#h_plot6=ROOT.TH1D("h_plot6","Comparsion of the CS between 2016 and 2017 at 5 TeV to 6 TeV",100,-1.,1.) +#h_plot6.GetXaxis().SetTitle("cos #theta") +#h_plot6.Draw() +h_16[4].SetTitle("Comparsion of the CS between 2016 and 2017 at 5 TeV to 6 TeV") +h_16[4].SetLineColor(2) +h_17[4].SetLineColor(3) +ly6.AddEntry(h_16[4],"2016") +ly6.AddEntry(h_17[4],"2017") +h_16[4].Draw("hist") +h_17[4].Draw("samehist") +ly6.Draw() + +c3.Print("Study/Result3.pdf") + +#masses +#bng1=[i for i in range(200,1000,5)]+[i for i in range(1000,2000,10)]+[i for i in range(2000,3500,20)] +#bng1=numpy.asarray(bng1,dtype=numpy.float64) +#h_Corr1=ROOT.TH1D("h_Corr1", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_Corr1.SetLineColor(2) +#h_Corr2=ROOT.TH1D("h_Corr2", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_Corr2.SetLineColor(3) +#h_Corr3=ROOT.TH1D("h_Corr3", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_Corr3.SetLineColor(4) +#h_Corr4=ROOT.TH1D("h_Corr4", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_Corr4.SetLineColor(5) +#h_Corr4.GetXaxis().SetTitle("mll[GeV]") +#h_Corr4.GetYaxis().SetRangeUser(-0.15,0.15) +#h_CorrT=ROOT.TH1D("h_CorrT", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_CorrT.SetLineColor(2) +#h_Corr=ROOT.TH1D("h_Corr", "Pearson Correlation Coefficient",len(bng1)-1,bng1) +#h_Corr.SetLineColor(3) +#h_Corr.GetXaxis().SetTitle("mll[GeV]") +#h_Corr.GetYaxis().SetRangeUser(-0.02,0.06) +#ly2=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) +#ly3=ROOT.TLegend(0.1, 0.75, 0.4, 0.925) + + +#for i in range(1,len(bng1)): +# if bng1[i]<1000: +# lowerLimit=0.9*bng1[i] +# upperLimit=1.1*bng1[i] +# elif bng1[i]<2000: +# lowerLimit=0.85*bng1[i] +# upperLimit=1.15*bng1[i] +# else: +# lowerLimit=0.8*bng1[i] +# upperLimit=1.2*bng1[i] +# obsmass1=mass[(masslowerLimit) &(numpy.absolute(y)<0.25)] +# obsCS1=CS[(masslowerLimit)&(numpy.absolute(y)<0.25)] +# corr1=numpy.corrcoef(obsmass1,obsCS1) +# h_Corr1.SetBinContent(i,corr1[0,1]) +# obsmass2=mass[(masslowerLimit) &(numpy.absolute(y)>0.25)&(numpy.absolute(y)<0.5)] +# obsCS2=CS[(masslowerLimit)&(numpy.absolute(y)>0.25)&(numpy.absolute(y)<0.5)] +# corr2=numpy.corrcoef(obsmass2,obsCS2) +# h_Corr2.SetBinContent(i,corr2[0,1]) +# obsmass3=mass[(masslowerLimit) &(numpy.absolute(y)>0.5)&(numpy.absolute(y)<0.75)] +# obsCS3=CS[(masslowerLimit)&(numpy.absolute(y)>0.5)&(numpy.absolute(y)<0.75)] +# corr3=numpy.corrcoef(obsmass3,obsCS3) +# h_Corr3.SetBinContent(i,corr3[0,1]) +# obsmass4=mass[(masslowerLimit) &(numpy.absolute(y)>0.75)&(numpy.absolute(y)<1.)] +# obsCS4=CS[(masslowerLimit)&(numpy.absolute(y)>0.75)&(numpy.absolute(y)<1.)] +# corr4=numpy.corrcoef(obsmass4,obsCS4) +# h_Corr4.SetBinContent(i,corr4[0,1]) +# obsmass=mass[(masslowerLimit)] +# obsCST=CS_true[(masslowerLimit)] +# obsCS=CS[(masslowerLimit)] +# corrT=numpy.corrcoef(obsmass,obsCST) +# h_CorrT.SetBinContent(i,corrT[0,1]) +# corr=numpy.corrcoef(obsmass,obsCS) +# h_Corr.SetBinContent(i,corr[0,1]) +#plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) +#plotPad.Draw() +#plotPad.cd() +#c2=ROOT.TCanvas("c2","c2",800,800) +#c2.Divide(2,1) +#c2.cd(1) +#ly2.AddEntry(h_Corr1,"0<|y|<0.25") +#ly2.AddEntry(h_Corr2,"0.25<|y|<0.5") +#ly2.AddEntry(h_Corr3,"0.5<|y|<0.75") +#ly2.AddEntry(h_Corr4,"0.75<|y|<1") +#h_Corr4.Draw() +#h_Corr3.Draw("same") +#h_Corr2.Draw("same") +#h_Corr1.Draw("same") +#ly2.SetTextFont(32) +#ly2.Draw() +#ROOT.gPad.SetLogx() +#c2.cd(2) +#ly3.AddEntry(h_Corr,"CS") +#ly3.AddEntry(h_CorrT,"True CS") +#h_Corr.Draw() +#h_CorrT.Draw("same") +#ly3.SetTextFont(32) +#ly3.Draw() +#plotPad.SetLogx() +#ROOT.gPad.SetLogx() +#plotPad.Update() +#c2.Print("Study/Result2.pdf") +#plotPad.Update() + +#h_CS1000to1200=ROOT.TH1D("h_CS1","CS angle",100,-1.,1.) +#h_CS1000to1200.SetLineColor(2) +#h_CS1000to1200.GetXaxis().SetTitle("cos #theta") +#h_CS3000to3400=ROOT.TH1D("h_CS2","CS angle",100,-1.,1.) +#h_CS3000to3400.SetLineColor(3) + +#h_CS1000to1200T=ROOT.TH1D("h_CST1","CS angle truth",100,-1.,1.) +#h_CS1000to1200T.SetLineColor(2) +#h_CS3000to3400T=ROOT.TH1D("h_CST2","CS angle truth",100,-1.,1.) +#h_CS3000to3400T.SetLineColor(3) +#h_CS3000to3400T.GetXaxis().SetTitle("cos #theta") + +#h_CS1000to1200r1=ROOT.TH1D("h_CS1r1","CS angle |y|<0.5",100,-1.,1.) +#h_CS1000to1200r1.SetLineColor(2) +#h_CS3000to3400r1=ROOT.TH1D("h_CS2r1","CS angle |y|<0.5",100,-1.,1.) +#h_CS3000to3400r1.SetLineColor(3) +#h_CS3000to3400r1.GetXaxis().SetTitle("cos #theta") + +#h_CS1000to1200r2=ROOT.TH1D("h_CS1r2","CS angle",100,-1.,1.) +#h_CS1000to1200r2.SetLineColor(2) +#h_CS3000to3400r2=ROOT.TH1D("h_CS2r2","CS angle 0.5<|y|<1",100,-1.,1.) +#h_CS3000to3400r2.SetLineColor(3) +#h_CS3000to3400r2.GetXaxis().SetTitle("cos #theta") + +#h_CS1000to1200Tr=ROOT.TH1D("h_CS1tr","CS angle truth |y|<1",100,-1.,1.) +#h_CS1000to1200Tr.SetLineColor(2) +#h_CS3000to3400Tr=ROOT.TH1D("h_CS2tr","CS angle truth |y|<1",100,-1.,1.) +#h_CS3000to3400Tr.SetLineColor(3) +#h_CS3000to3400Tr.GetXaxis().SetTitle("cos #theta") + +#h_CS1000to1200.FillN(len(CS[(mass<1200) & (mass>1000)]),CS[(mass<1200) & (mass>1000)],numpy.ones(len(CS[(mass<1200) & (mass>1000)]))) +#h_CS1000to1200.Scale(1./len(CS[(mass<1200) & (mass>1000)])) +#h_CS3000to3400.FillN(len(CS[(mass<3400) & (mass>3000)]),CS[(mass<3400) & (mass>3000)],numpy.ones(len(CS[(mass<3400) & (mass>3000)]))) +#h_CS3000to3400.Scale(1./len(CS[(mass<3400) & (mass>3000)])) +#h_CS1000to1200T.FillN(len(CS_true[(mass<1200) & (mass>1000)]),CS_true[(mass<1200) & (mass>1000)],numpy.ones(len(CS_true[(mass<1200) & (mass>1000)]))) +#h_CS1000to1200T.Scale(1./len(CS_true[(mass<1200) & (mass>1000)])) +#h_CS3000to3400T.FillN(len(CS_true[(mass<3400) & (mass>3000)]),CS_true[(mass<3400) & (mass>3000)],numpy.ones(len(CS_true[(mass<3400) & (mass>3000)]))) +#h_CS3000to3400T.Scale(1./len(CS_true[(mass<3400) & (mass>3000)])) +#h_CS1000to1200r1.FillN(len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)<0.5)]),CS[(mass<1200) & (mass>1000)&(numpy.absolute(y)<0.5)],numpy.ones(len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)<0.5)]))) +#h_CS1000to1200r1.Scale(1./len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)<0.5)])) +#h_CS3000to3400r1.FillN(len(CS[(mass<3400) & (mass>3000) & (numpy.absolute(y)<0.5)]),CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)<0.5)],numpy.ones(len(CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)<0.5)]))) +#h_CS3000to3400r1.Scale(1./len(CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)<0.5)])) +#h_CS1000to1200r2.FillN(len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)>0.5)&(numpy.absolute(y)<1)]),CS[(mass<1200) & (mass>1000)&(numpy.absolute(y)>0.5)&(numpy.absolute(y)<1)],numpy.ones(len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)>0.5)&(numpy.absolute(y)<1)]))) +#h_CS1000to1200r2.Scale(1./len(CS[(mass<1200) & (mass>1000) & (numpy.absolute(y)>0.5)&(numpy.absolute(y)<1)])) +#h_CS3000to3400r2.FillN(len(CS[(mass<3400) & (mass>3000) & (numpy.absolute(y)>0.5)& (numpy.absolute(y)<1)]),CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)>0.5) & (numpy.absolute(y)<1)],numpy.ones(len(CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)>0.5) & (numpy.absolute(y)<1)]))) +#h_CS3000to3400r2.Scale(1./len(CS[(mass<3400) & (mass>3000)& (numpy.absolute(y)>0.5) & (numpy.absolute(y)<1)])) +#h_CS1000to1200Tr.FillN(len(CS_true[(mass<1200) & (mass>1000) &(numpy.absolute(y)<1)]),CS_true[(mass<1200) & (mass>1000)&(numpy.absolute(y)<1)],numpy.ones(len(CS_true[(mass<1200) & (mass>1000) &(numpy.absolute(y)<1)]))) +#h_CS1000to1200Tr.Scale(1./len(CS_true[(mass<1200) & (mass>1000) &(numpy.absolute(y)<1)])) +#h_CS3000to3400Tr.FillN(len(CS_true[(mass<3400) & (mass>3000) & (numpy.absolute(y)<1)]),CS_true[(mass<3400) & (mass>3000)& (numpy.absolute(y)<1)],numpy.ones(len(CS_true[(mass<3400) & (mass>3000)& (numpy.absolute(y)<1)]))) +#h_CS3000to3400Tr.Scale(1./len(CS_true[(mass<3400) & (mass>3000)& (numpy.absolute(y)<1)])) + +#c3=ROOT.TCanvas("c3","c3",800,800) +#c3.Divide(2,3) +#c3.cd(1) +#h_CS3000to3400.Draw() +#h_CS1000to1200.Draw("same") +#c3.cd(2) +#h_CS3000to3400T.Draw() +#h_CS1000to1200T.Draw("same") +#c3.cd(3) +#h_CS3000to3400r1.Draw() +#h_CS1000to1200r1.Draw("same") +#c3.cd(4) +#h_CS3000to3400r2.Draw() +#h_CS1000to1200r2.Draw("same") +#c3.cd(5) +#h_CS3000to3400Tr.Draw() +#h_CS1000to1200Tr.Draw("same") + +#c3.Print("Study/Result3.pdf") + + diff --git a/run.sh b/run.sh new file mode 100644 index 0000000..a45bed6 --- /dev/null +++ b/run.sh @@ -0,0 +1,11 @@ +python plotAllMCFlavor.py -d -r -c --all --ae +python plotAllMCFlavor.py -d -r -c -2018 --ae +python plotAllMCFlavor.py -d -r -c -2016 --ae +python plotAllMCFlavor.py -d -r -c --ae +python plotAllMCFlavor.py -d -r -c --all +python plotAllMCFlavor.py -d -r -c -2018 +python plotAllMCFlavor.py -d -r -c -2016 +python plotAllMCFlavor.py -d -r -c + + + diff --git a/unfoldDataPlot.py b/unfoldDataPlot.py new file mode 100755 index 0000000..e1dda48 --- /dev/null +++ b/unfoldDataPlot.py @@ -0,0 +1,251 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, SignalsADD, Signals2018ADD, Signals2018, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy + + +def plotData(datahist,mchist,usedata, label1, label2,name,filename): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if usedata==True: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.3,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.3) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw("hist") + ratioPad.Draw("hist") + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + colors = createMyColors() + + legend = TLegend(0.55, 0.6, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legendEta = TLegend(0.45, 0.75, 0.925, 0.925) + legendEta.SetFillStyle(0) + legendEta.SetBorderSize(0) + legendEta.SetTextFont(42) + legendEta.SetNColumns(2) + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + legendHistData = ROOT.TH1F() + category = ROOT.TLatex() + category.SetNDC(True) + category.SetTextSize(0.04) + if usedata==True: + legend.AddEntry(legendHistData,label1,"pe") + legendEta.AddEntry(legendHistData,label1,"pe") + + #process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+^{*}}e^{-}") + temphist = ROOT.TH1F() + temphist.SetFillColor(3) + legendHists.append(temphist.Clone) + #legend.AddEntry(temphist,label2,"l") + #legendEta.AddEntry(temphist,process.label,"f") + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + plotPad.SetLogy() + if usedata==True: + yMax = datahist.GetBinContent(datahist.GetMaximumBin())*1000 + yMin = 0.00000001 + xMax = datahist.GetXaxis().GetXmax() + xMin = datahist.GetXaxis().GetXmin() + else: + yMax = mchist.GetBinContent(datahist.GetMaximumBin()) + yMin = 0.00000001 + xMax = mchist.GetXaxis().GetXmax() + xMin = mchist.GetXaxis().GetXmin() + yMax = yMax*10000 + if name.find("dimuon")!=-1: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{#mu#mu}[GeV] ;Events/GeV") + else: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{ee}[GeV] ;Events/GeV") + # Draw signal information + + # Draw background from stack + mchist.SetFillColor(0) + mchist.SetLineColor(2) + legend.AddEntry(mchist,label2,"l") + mchist.Draw("samehist") + + # Draw data + datahist.SetMinimum(0.0001) + if usedata==True: + datahist.SetMarkerStyle(8) + datahist.Draw("samepehist") + + # Draw legend + legend.Draw() + plotPad.SetLogx() + latex.DrawLatex(0.95,0.96,"13 TeV") + yLabelPos = 0.85 + cmsExtra = "Preliminary" + if not usedata==True: + cmsExtra = "#splitline{Preliminary}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + category.DrawLatex(0.3,0.7,name) + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + #~ print datahist.Integral() + if usedata==True: + try: + ratioPad.cd() + ratioPad.SetLogx() + except AttributeError: + print ("Plot fails. Look up in errs/failedPlots.txt") + outFile =open("errs/failedPlots.txt","a") + outFile.write('%s\n'%plot.filename%("_"+run.label+"_"+dilepton)) + outFile.close() + #plot.cuts=baseCut + return 1 + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title="After/Before",yMin=0,yMax=2.0,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if usedata==True: + + ratioPad.RedrawAxis() + if not os.path.exists("plots"): + os.makedirs("plots") + hCanvas.Print("%s.pdf"%filename) + + + + +f=ROOT.TFile.Open("RooUnfold/UnfoldedSample.root") +recBB_17_mu=f.Get("datamubb2017") +unfoldBB_17_mu=f.Get("BB2017mu_invert") +unfoldBB_17_muEM=f.Get("BB2017mu_reg") +recBE_17_mu=f.Get("datamube2017") +unfoldBE_17_mu=f.Get("BE2017mu_invert") +unfoldBE_17_muEM=f.Get("BE2017mu_reg") +recBB_16_mu=f.Get("datamubb2016") +unfoldBB_16_mu=f.Get("BB2016mu_invert") +unfoldBB_16_muEM=f.Get("BB2016mu_reg") +recBE_16_mu=f.Get("datamube2016") +unfoldBE_16_mu=f.Get("BE2016mu_invert") +unfoldBE_16_muEM=f.Get("BE2016mu_reg") +recBB_18_mu=f.Get("datamubb2018") +for i in range(recBB_18_mu.GetNbinsX()+2): + val=recBB_18_mu.GetBinContent(i) + print val +unfoldBB_18_mu=f.Get("BB2018mu_invert") +unfoldBB_18_muEM=f.Get("BB2018mu_reg") +unfoldBB_18_muiter400=f.Get("BB2018muiter_400") +recBE_18_mu=f.Get("datamube2018") +unfoldBE_18_mu=f.Get("BE2018mu_invert") +unfoldBE_18_muEM=f.Get("BE2018mu_reg") + +recBB_17_e=f.Get("dataebb2017") +unfoldBB_17_e=f.Get("BB2017el_invert") +unfoldBB_17_eEM=f.Get("BB2017el_reg") +recBE_17_e=f.Get("dataebe2017") +unfoldBE_17_e=f.Get("BE2017el_invert") +unfoldBE_17_eEM=f.Get("BE2017el_reg") +recBB_16_e=f.Get("dataebb2016") +unfoldBB_16_e=f.Get("BB2016el_invert") +unfoldBB_16_eEM=f.Get("BB2016el_reg") +recBE_16_e=f.Get("dataebe2016") +unfoldBE_16_e=f.Get("BE2016el_invert") +unfoldBE_16_eEM=f.Get("BE2016el_reg") +recBB_18_e=f.Get("dataebb2018") +unfoldBB_18_e=f.Get("BB2018el_invert") +unfoldBB_18_eEM=f.Get("BB2018el_reg") +recBE_18_e=f.Get("dataebe2018") +unfoldBE_18_e=f.Get("BE2018el_invert") +unfoldBE_18_eEM=f.Get("BE2018el_reg") + + + +plotDataMC(unfoldBB_17_mu,recBB_17_mu,True,"Unfolded","Reconstructed","2017 BB dimuon","data_BB_17_mu" ) +plotDataMC(unfoldBE_17_mu,recBE_17_mu,True,"Unfolded","Reconstructed","2017 BE dimuon","data_BE_17_mu" ) +plotDataMC(unfoldBB_17_muEM,recBB_17_mu,True,"Unfolded","Reconstruted","2017 BB dimuon","data_BB_17_muEM" ) +plotDataMC(unfoldBE_17_muEM,recBE_17_mu,True,"Unfolded","Reconstruted","2017 BE dimuon","data_BE_17_muEM" ) + +plotDataMC(unfoldBB_16_mu,recBB_16_mu,True,"Unfolded","Reconstructed","2016 BB dimuon","data_BB_16_mu" ) +plotDataMC(unfoldBE_16_mu,recBE_16_mu,True,"Unfolded","Reconstructed","2016 BE dimuon","data_BE_16_mu" ) +plotDataMC(unfoldBB_16_muEM,recBB_16_mu,True,"Unfolded","Reconstructed","2016 BB dimuon","data_BB_16_muEM" ) +plotDataMC(unfoldBE_16_muEM,recBE_16_mu,True,"Unfolded","Reconstructed","2016 BE Dimuon","data_BE_16_muEM" ) + +plotDataMC(unfoldBB_18_e,recBB_18_e,True,"Unfolded","Reconstructed","2018 BB dielectron","data_BB_18_e" ) +plotDataMC(unfoldBE_18_e,recBE_18_e,True,"Unfolded","Reconstructed","2018 BE dielectron","data_BE_18_e" ) +plotDataMC(unfoldBB_18_eEM,recBB_18_e,True,"Unfolded","Reconstructed","2018 BB dielectron","data_BB_18_eEM" ) +plotDataMC(unfoldBE_18_eEM,recBE_18_e,True,"Unfolded","Reconstructed","2018 BE dielectron","data_BE_18_eEM" ) + +plotDataMC(unfoldBB_17_e,recBB_17_e,True,"Unfolded","Reconstructed","2017 BB dielectron","data_BB_17_e" ) +plotDataMC(unfoldBE_17_e,recBE_17_e,True,"Unfolded","Reconstructed","2017 BE dielectron","data_BE_17_e" ) +plotDataMC(unfoldBB_17_eEM,recBB_17_e,True,"Unfolded","Reconstructed","2017 BB dielectron","data_BB_17_eEM" ) +plotDataMC(unfoldBE_17_eEM,recBE_17_e,True,"Unfolded","Reconstructed","2017 BE dielectron","data_BE_17_eEM" ) + +plotDataMC(unfoldBB_16_e,recBB_16_e,True,"Unfolded","Reconstructed","2016 BB dielectron","data_BB_16_e" ) +plotDataMC(unfoldBE_16_e,recBE_16_e,True,"Unfolded","Reconstructed","2016 BE dielectron","data_BE_16_e" ) +plotDataMC(unfoldBB_16_eEM,recBB_16_e,True,"Unfolded","Reconstructed","2016 BB dielectron","data_BB_16_eEM" ) +plotDataMC(unfoldBE_16_eEM,recBE_16_e,True,"Unfolded","Reconstructed","2016 BE dielectron","data_BE_16_eEM" ) + +plotDataMC(unfoldBB_18_mu,recBB_18_mu,True,"Unfolded","Reconsructed","2018 BB dimuon","data_BB_18_mu" ) +plotDataMC(unfoldBE_18_mu,recBE_18_mu,True,"Unfolded","Reconstructed","2018 BE dimuon","data_BE_18_mu" ) +plotDataMC(unfoldBB_18_muEM,recBB_18_mu,True,"Unfolded","Reconstructed","2018 BB dimuon","data_BB_18_muEM" ) +plotDataMC(unfoldBE_18_muEM,recBE_18_mu,True,"Unfolded","Reconstructed","2018 BE Dimuon","data_BE_18_muEM" ) +plotDataMC(unfoldBB_18_muiter400,recBB_18_mu,True,"Unfolded","Reconstructed","2018 BB dimuon","data_BB_18_muiter400" ) + diff --git a/unfoldDataPlot_V2.py b/unfoldDataPlot_V2.py new file mode 100755 index 0000000..9bd2608 --- /dev/null +++ b/unfoldDataPlot_V2.py @@ -0,0 +1,190 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, SignalsADD, Signals2018ADD, Signals2018, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy + + +def plotDataMC(datahist,mchist,usedata, label1, label2,name,filename): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if usedata==True: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.3,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.3) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw("hist") + ratioPad.Draw("hist") + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + colors = createMyColors() + + legend = TLegend(0.55, 0.6, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legendEta = TLegend(0.45, 0.75, 0.925, 0.925) + legendEta.SetFillStyle(0) + legendEta.SetBorderSize(0) + legendEta.SetTextFont(42) + legendEta.SetNColumns(2) + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + legendHistData = ROOT.TH1F() + category = ROOT.TLatex() + category.SetNDC(True) + category.SetTextSize(0.04) + if usedata==True: + legend.AddEntry(legendHistData,label1,"pe") + legendEta.AddEntry(legendHistData,label1,"pe") + + #process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+^{*}}e^{-}") + temphist = ROOT.TH1F() + temphist.SetFillColor(3) + legendHists.append(temphist.Clone) + #legend.AddEntry(temphist,label2,"l") + #legendEta.AddEntry(temphist,process.label,"f") + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + plotPad.SetLogy() + if usedata==True: + yMax = datahist.GetBinContent(datahist.GetMaximumBin())*1000 + yMin = 0.00000001 + xMax = datahist.GetXaxis().GetXmax() + xMin = datahist.GetXaxis().GetXmin() + else: + yMax = mchist.GetBinContent(datahist.GetMaximumBin()) + yMin = 0.00000001 + xMax = mchist.GetXaxis().GetXmax() + xMin = mchist.GetXaxis().GetXmin() + yMax = yMax*10000 + if name.find("dimuon")!=-1: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{#mu#mu}[GeV] ;Events/GeV") + else: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{ee}[GeV] ;Events/GeV") + # Draw signal information + + # Draw background from stack + mchist.SetFillColor(0) + mchist.SetLineColor(2) + legend.AddEntry(mchist,label2,"l") + mchist.Draw("samehist") + + # Draw data + datahist.SetMinimum(0.0001) + if usedata==True: + datahist.SetMarkerStyle(8) + datahist.Draw("samepehist") + + # Draw legend + legend.Draw() + plotPad.SetLogx() + latex.DrawLatex(0.95,0.96,"13 TeV") + yLabelPos = 0.85 + cmsExtra = "Preliminary" + if not usedata==True: + cmsExtra = "#splitline{Preliminary}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + category.DrawLatex(0.3,0.7,name) + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + #~ print datahist.Integral() + if usedata==True: + try: + ratioPad.cd() + ratioPad.SetLogx() + except AttributeError: + print ("Plot fails. Look up in errs/failedPlots.txt") + outFile =open("errs/failedPlots.txt","a") + outFile.write('%s\n'%plot.filename%("_"+run.label+"_"+dilepton)) + outFile.close() + #plot.cuts=baseCut + return 1 + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title="R_{mumu}/R_{ee}",yMin=0,yMax=2.0,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if usedata==True: + + ratioPad.RedrawAxis() + if not os.path.exists("plots"): + os.makedirs("plots") + hCanvas.Print("plots/%s.pdf"%filename) + + + + +fkFac=ROOT.TFile.Open("RooUnfold/kFacTest.root") +unfold_mu=fkFac.Get("response_mu_invert") +unfold_mu_kFac=fkFac.Get("response_mu_kFac_invert") +unfold_el=fkFac.Get("response_el_invert") +unfold_el_kFac=fkFac.Get("response_el_kFac_invert") + +plotDataMC(unfold_mu_kFac,unfold_mu,True,"kFac/o","kFac/w","dimuon","kFac_mu" ) +plotDataMC(unfold_el_kFac,unfold_el,True,"kFac/o","kFac/w","dielectron","kFac_el" ) + +fUnscaled=ROOT.TFile.Open("RooUnfold/unfoldedMass.root") +unfoldedMass_mu=fUnscaled.Get("response_mu_invert") +unfoldedMass_el=fUnscaled.Get("response_el_invert") + +plotDataMC(unfoldedMass_mu,unfoldedMass_el,True,"muon","electron","uncorrected ratio","Uncorrected" ) + diff --git a/unfoldPlot.py b/unfoldPlot.py new file mode 100755 index 0000000..b060bb8 --- /dev/null +++ b/unfoldPlot.py @@ -0,0 +1,303 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, SignalsADD, Signals2018ADD, Signals2018, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy + + +def plotDataMC(datahist,mchist,usedata, label1, label2,name,filename): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if usedata==True: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.3,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.3) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw("hist") + ratioPad.Draw("hist") + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + colors = createMyColors() + legend = TLegend(0.55, 0.6, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legendEta = TLegend(0.45, 0.75, 0.925, 0.925) + legendEta.SetFillStyle(0) + legendEta.SetBorderSize(0) + legendEta.SetTextFont(42) + legendEta.SetNColumns(2) + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + legendHistData = ROOT.TH1F() + category = ROOT.TLatex() + category.SetNDC(True) + category.SetTextSize(0.04) + if usedata==True: + legend.AddEntry(legendHistData,label1,"pe") + legendEta.AddEntry(legendHistData,label1,"pe") + + #process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+^{*}}e^{-}") + temphist = ROOT.TH1F() + if label2=="Generated": + temphist.SetFillColor(3) + mchist.SetFillColor(3) + else: + temphist.SetFillColor(2) + mchist.SetFillColor(2) + legendHists.append(temphist.Clone) + legend.AddEntry(temphist,label2,"f") + #legendEta.AddEntry(temphist,process.label,"f") + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + plotPad.SetLogy() + if usedata==True: + yMax = datahist.GetBinContent(datahist.GetMaximumBin())*1000 + yMin = 0.00000001 + xMax = datahist.GetXaxis().GetXmax() + xMin = datahist.GetXaxis().GetXmin() + else: + yMax = mchist.GetBinContent(datahist.GetMaximumBin()) + yMin = 0.00000001 + xMax = mchist.GetXaxis().GetXmax() + xMin = mchist.GetXaxis().GetXmin() + yMax = yMax*10000 + if name.find("dimuon")!=-1: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{#mu#mu}[GeV] ;fb/GeV") + else: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{ee}[GeV] ;fb/GeV") + + + # Draw signal information + + # Draw background from stack + mchist.Draw("samehist") + + # Draw data + datahist.SetMinimum(0.0001) + if usedata==True: + datahist.SetMarkerStyle(8) + datahist.Draw("samepehist") + + # Draw legend + legend.Draw() + plotPad.SetLogx() + latex.DrawLatex(0.95,0.96,"13 TeV") + yLabelPos = 0.85 + cmsExtra = "Preliminary" + if not usedata==True: + cmsExtra = "#splitline{Preliminary}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + category.DrawLatex(0.3,0.7,name) + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + #~ print datahist.Integral() + if usedata==True: + try: + ratioPad.cd() + ratioPad.SetLogx() + except AttributeError: + print ("Plot fails. Look up in errs/failedPlots.txt") + outFile =open("errs/failedPlots.txt","a") + outFile.write('%s\n'%plot.filename%("_"+run.label+"_"+dilepton)) + outFile.close() + #plot.cuts=baseCut + return 1 + if label2=="Generated": + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title=" Unfolded/Gen",yMin=0.7,yMax=1.3,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + else: + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title="Unfolded/Reco",yMin=0.7,yMax=1.3,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if usedata==True: + + ratioPad.RedrawAxis() + if not os.path.exists("plots"): + os.makedirs("plots") + hCanvas.Print("%s.pdf"%filename) + + + + +f=ROOT.TFile.Open("RooUnfold/ClosureTest.root") +genBB_17_mu=f.Get("genBB2017_mu") +recBB_17_mu=f.Get("recBB2017_mu") +unfoldBB_17_mu=f.Get("BB2017mu_invert") +unfoldBB_17_muEM=f.Get("BB2017mu_reg") +genBE_17_mu=f.Get("genBE2017_mu") +recBE_17_mu=f.Get("recBE2017_mu") +unfoldBE_17_mu=f.Get("BE2017mu_invert") +unfoldBE_17_muEM=f.Get("BE2017mu_reg") +genBB_16_mu=f.Get("genBB2016_mu") +recBB_16_mu=f.Get("recBB2016_mu") +unfoldBB_16_mu=f.Get("BB2016mu_invert") +unfoldBB_16_muEM=f.Get("BB2016mu_reg") +genBE_16_mu=f.Get("genBE2016_mu") +recBE_16_mu=f.Get("recBE2016_mu") +unfoldBE_16_mu=f.Get("BE2016mu_invert") +unfoldBE_16_muEM=f.Get("BE2016mu_reg") +genBB_18_mu=f.Get("genBB2018_mu") +recBB_18_mu=f.Get("recBB2018_mu") +unfoldBB_18_mu=f.Get("BB2018mu_invert") +unfoldBB_18_muEM=f.Get("BB2018mu_reg") +genBE_18_mu=f.Get("genBE2018_mu") +recBE_18_mu=f.Get("recBE2018_mu") +unfoldBE_18_mu=f.Get("BE2018mu_invert") +unfoldBE_18_muEM=f.Get("BE2018mu_reg") +genBB_17_e=f.Get("genBB2017_e") +recBB_17_e=f.Get("recBB2017_e") +unfoldBB_17_e=f.Get("BB2017el_invert") +unfoldBB_17_eEM=f.Get("BB2017el_reg") +genBE_17_e=f.Get("genBE2017_e") +recBE_17_e=f.Get("recBE2017_e") +unfoldBE_17_e=f.Get("BE2017el_invert") +unfoldBE_17_eEM=f.Get("BE2017el_reg") +genBB_16_e=f.Get("genBB2016_e") +recBB_16_e=f.Get("recBB2016_e") +unfoldBB_16_e=f.Get("BB2016el_invert") +unfoldBB_16_eEM=f.Get("BB2016el_reg") +genBE_16_e=f.Get("genBE2016_e") +recBE_16_e=f.Get("recBE2016_e") +unfoldBE_16_e=f.Get("BE2016el_invert") +unfoldBE_16_eEM=f.Get("BE2016el_reg") +genBB_18_e=f.Get("genBB2018_e") +recBB_18_e=f.Get("recBB2018_e") +unfoldBB_18_e=f.Get("BB2018el_invert") +unfoldBB_18_eEM=f.Get("BB2018el_reg") +genBE_18_e=f.Get("genBE2018_e") +recBE_18_e=f.Get("recBE2018_e") +unfoldBE_18_e=f.Get("BE2018el_invert") +unfoldBE_18_eEM=f.Get("BE2018el_reg") + + + +plotDataMC(unfoldBB_17_mu,genBB_17_mu,True,"Unfolded","Generated","2017 BB dimuon","gen_closure_BB_17_mu" ) +plotDataMC(unfoldBB_17_mu,recBB_17_mu,True,"Unfolded","Reconstructed","2017 BB dimuon","rec_closure_BB_17_mu" ) +plotDataMC(unfoldBE_17_mu,genBE_17_mu,True,"Unfolded","Generated","2017 BE dimuon","gen_closure_BE_17_mu" ) +plotDataMC(unfoldBE_17_mu,recBE_17_mu,True,"Unfolded","Reconstructed","2017 BE dimuon","rec_closure_BE_17_mu" ) + + +plotDataMC(unfoldBB_17_e,genBB_17_e,True,"Unfolded","Generated","2017 BB dielectron","gen_closure_BB_17_el" ) +plotDataMC(unfoldBB_17_e,recBB_17_e,True,"Unfolded","Reconstructed","2017 BB dielectron","rec_closure_BB_17_el" ) +plotDataMC(unfoldBE_17_e,genBE_17_e,True,"Unfolded","Generated","2017 BE dielectron","gen_closure_BE_17_el" ) +plotDataMC(unfoldBE_17_e,recBE_17_e,True,"Unfolded","Reconstructed","2017 BE dielectron","rec_closure_BE_17_el" ) + +plotDataMC(unfoldBB_16_mu,genBB_16_mu,True,"Unfolded","Generated","2016 BB dimuon","gen_closure_BB_16_mu" ) +plotDataMC(unfoldBB_16_mu,recBB_16_mu,True,"Unfolded","Reconstructed","2016 BB dimuon","rec_closure_BB_16_mu" ) +plotDataMC(unfoldBE_16_mu,genBE_16_mu,True,"Unfolded","Generated","2016 BE dimuon","gen_closure_BE_16_mu" ) +plotDataMC(unfoldBE_16_mu,recBE_16_mu,True,"Unfolded","Reconstructed","2016 BE dimuon","rec_closure_BE_16_mu" ) + + +plotDataMC(unfoldBB_16_e,genBB_16_e,True,"Unfolded","Generated","2016 BB dielectron","gen_closure_BB_16_el" ) +plotDataMC(unfoldBB_16_e,recBB_16_e,True,"Unfolded","Reconstructed","2016 BB dielectron","rec_closure_BB_16_el" ) +plotDataMC(unfoldBE_16_e,genBE_16_e,True,"Unfolded","Generated","2016 BE dielectron","gen_closure_BE_16_el" ) +plotDataMC(unfoldBE_16_e,recBE_16_e,True,"Unfolded","Recoconstructed","2016 BE dielectron","rec_closure_BE_16_el" ) + +plotDataMC(unfoldBB_18_mu,genBB_18_mu,True,"Unfolded","Generated","2018 BB dimuon","gen_closure_BB_18_mu" ) +plotDataMC(unfoldBB_18_mu,recBB_18_mu,True,"Unfolded","Reconstructed","2018 BB dimuon","rec_closure_BB_18_mu" ) +plotDataMC(unfoldBE_18_mu,genBE_18_mu,True,"Unfolded","Generated","2018 BE dimuon","gen_closure_BE_18_mu" ) +plotDataMC(unfoldBE_18_mu,recBE_18_mu,True,"Unfolded","Reconstructed","2018 BE dimuon","rec_closure_BE_18_mu" ) + + +plotDataMC(unfoldBB_18_e,genBB_18_e,True,"Unfolded","Generated","2018 BB dielectron","gen_closure_BB_18_el" ) +plotDataMC(unfoldBB_18_e,recBB_18_e,True,"Unfolded","Reconstructed","2018 BB dielectron","rec_closure_BB_18_el" ) +plotDataMC(unfoldBE_18_e,genBE_18_e,True,"Unfolded","Generated","2018 BE dielectron","gen_closure_BE_18_el" ) +plotDataMC(unfoldBE_18_e,recBE_18_e,True,"Unfolded","Reconstructed","2018 BE dielectron","rec_closure_BE_18_el" ) + + + +plotDataMC(unfoldBB_17_muEM,genBB_17_mu,True,"Unfolded","Generated","2017 BB dimuon","gen_closure_BB_17_muEM" ) +plotDataMC(unfoldBB_17_muEM,recBB_17_mu,True,"Unfolded","Reconstructed","2017 BB dimuon","rec_closure_BB_17_muEM" ) +plotDataMC(unfoldBE_17_muEM,genBE_17_mu,True,"Unfolded","Generated","2017 BE dimuon","gen_closure_BE_17_muEM" ) +plotDataMC(unfoldBE_17_muEM,recBE_17_mu,True,"Unfolded","Reconstructed","2017 BE dimuon","rec_closure_BE_17_muEM" ) + + +plotDataMC(unfoldBB_17_eEM,genBB_17_e,True,"Unfolded","Generated","2017 BB dielectron","gen_closure_BB_17_elEM" ) +plotDataMC(unfoldBB_17_eEM,recBB_17_e,True,"Unfolded","Reconstructed","2017 BB dielectron","rec_closure_BB_17_elEM" ) +plotDataMC(unfoldBE_17_eEM,genBE_17_e,True,"Unfolded","Generated","2017 BE dielectron","gen_closure_BE_17_elEM" ) +plotDataMC(unfoldBE_17_eEM,recBE_17_e,True,"Unfolded","Reconstructed","2017 BE dielectron","rec_closure_BE_17_elEM" ) + +plotDataMC(unfoldBB_16_muEM,genBB_16_mu,True,"Unfolded","Generated","2016 BB dimuon","gen_closure_BB_16_muEM" ) +plotDataMC(unfoldBB_16_muEM,recBB_16_mu,True,"Unfolded","Reconstructed","2016 BB dimuon","rec_closure_BB_16_muEM" ) +plotDataMC(unfoldBE_16_muEM,genBE_16_mu,True,"Unfolded","Generated","2016 BE dimuon","gen_closure_BE_16_muEM" ) +plotDataMC(unfoldBE_16_muEM,recBE_16_mu,True,"Unfolded","Reconstructed","2016 BE dimuon","rec_closure_BE_16_muEM" ) + +plotDataMC(unfoldBB_16_eEM,genBB_16_e,True,"Unfolded","Generated","2016 BB dielectron","gen_closure_BB_16_elEM" ) +plotDataMC(unfoldBB_16_eEM,recBB_16_e,True,"Unfolded","Reconstructed","2016 BB dielectron","rec_closure_BB_16_elEM" ) +plotDataMC(unfoldBE_16_eEM,genBE_16_e,True,"Unfolded","Generated","2016 BE dielectron","gen_closure_BE_16_elEM" ) +plotDataMC(unfoldBE_16_eEM,recBE_16_e,True,"Unfolded","Reconstructed","2016 BE dielectron","rec_closure_BE_16_elEM" ) + +plotDataMC(unfoldBB_18_muEM,genBB_18_mu,True,"Unfolded","Generated","2018 BB dimuon","gen_closure_BB_18_muEM" ) +plotDataMC(unfoldBB_18_muEM,recBB_18_mu,True,"Unfolded","Reconstructed","2018 BB dimuon","rec_closure_BB_18_muEM" ) +plotDataMC(unfoldBE_18_muEM,genBE_18_mu,True,"Unfolded","Generated","2018 BE dimuon","gen_closure_BE_18_muEM" ) +plotDataMC(unfoldBE_18_muEM,recBE_18_mu,True,"Unfolded","Reconstructed","2018 BE dimuon","rec_closure_BE_18_muEM" ) + + +plotDataMC(unfoldBB_18_eEM,genBB_18_e,True,"Unfolded","Generated","2018 BB dielectron","gen_closure_BB_18_elEM" ) +plotDataMC(unfoldBB_18_eEM,recBB_18_e,True,"Unfolded","Reconstructed","2018 BB dielectron","rec_closure_BB_18_elEM" ) +plotDataMC(unfoldBE_18_eEM,genBE_18_e,True,"Unfolded","Generated","2018 BE dielectron","gen_closure_BE_18_elEM" ) +plotDataMC(unfoldBE_18_eEM,recBE_18_e,True,"Unfolded","Reconstructed","2018 BE dielectron","rec_closure_BE_18_elEM" ) + + diff --git a/unfoldPlot_V2.py b/unfoldPlot_V2.py new file mode 100755 index 0000000..9e3348a --- /dev/null +++ b/unfoldPlot_V2.py @@ -0,0 +1,196 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, SignalsADD, Signals2018ADD, Signals2018, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy + + +def plotDataMC(datahist,mchist,usedata, label1, label2,name,filename): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if usedata==True: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.3,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.3) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw("hist") + ratioPad.Draw("hist") + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + colors = createMyColors() + legend = TLegend(0.55, 0.6, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legendEta = TLegend(0.45, 0.75, 0.925, 0.925) + legendEta.SetFillStyle(0) + legendEta.SetBorderSize(0) + legendEta.SetTextFont(42) + legendEta.SetNColumns(2) + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + legendHistData = ROOT.TH1F() + category = ROOT.TLatex() + category.SetNDC(True) + category.SetTextSize(0.04) + if usedata==True: + legend.AddEntry(legendHistData,label1,"pe") + legendEta.AddEntry(legendHistData,label1,"pe") + + #process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+^{*}}e^{-}") + temphist = ROOT.TH1F() + if label2=="Generated": + temphist.SetFillColor(3) + mchist.SetFillColor(3) + else: + temphist.SetFillColor(2) + mchist.SetFillColor(2) + legendHists.append(temphist.Clone) + legend.AddEntry(temphist,label2,"f") + #legendEta.AddEntry(temphist,process.label,"f") + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + plotPad.SetLogy() + if usedata==True: + yMax = datahist.GetBinContent(datahist.GetMaximumBin())*1000 + yMin = 0.00000001 + xMax = datahist.GetXaxis().GetXmax() + xMin = datahist.GetXaxis().GetXmin() + else: + yMax = mchist.GetBinContent(datahist.GetMaximumBin()) + yMin = 0.00000001 + xMax = mchist.GetXaxis().GetXmax() + xMin = mchist.GetXaxis().GetXmin() + yMax = yMax*10000 + if name.find("dimuon")!=-1: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{#mu#mu}[GeV] ;fb/GeV") + else: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{ee}[GeV] ;fb/GeV") + + + # Draw signal information + + # Draw background from stack + mchist.Draw("samehist") + + # Draw data + datahist.SetMinimum(0.0001) + if usedata==True: + datahist.SetMarkerStyle(8) + datahist.Draw("samepehist") + + # Draw legend + legend.Draw() + plotPad.SetLogx() + latex.DrawLatex(0.95,0.96,"13 TeV") + yLabelPos = 0.85 + cmsExtra = "Preliminary" + if not usedata==True: + cmsExtra = "#splitline{Preliminary}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + category.DrawLatex(0.3,0.7,name) + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + #~ print datahist.Integral() + if usedata==True: + try: + ratioPad.cd() + ratioPad.SetLogx() + except AttributeError: + print ("Plot fails. Look up in errs/failedPlots.txt") + outFile =open("errs/failedPlots.txt","a") + outFile.write('%s\n'%plot.filename%("_"+run.label+"_"+dilepton)) + outFile.close() + #plot.cuts=baseCut + return 1 + if label2=="Generated": + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title=" nokFac/kFac",yMin=0.7,yMax=1.3,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + else: + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title="nokFac/kFac",yMin=0.7,yMax=1.3,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if usedata==True: + + ratioPad.RedrawAxis() + if not os.path.exists("plots"): + os.makedirs("plots") + hCanvas.Print("plots/%s.pdf"%filename) + + + + +f=ROOT.TFile.Open("RooUnfold/kFacTest.root") +h_mu=f.Get("muon_combine_combine_recoMass_combine") +h_mu_kFac=f.Get("muon_combine_combine_recoMass_combine_kFac") +h_el=f.Get("electron_combine_combine_recoMass_combine") +h_el_kFac=f.Get("electron_combine_combine_recoMass_combine_kFac") + + +plotDataMC(h_mu_kFac,h_mu,True,"kFac/o","kFac/w","reco dimuon","rec_MukFac") +plotDataMC(h_el_kFac,h_el,True,"kFac/o","kFac/w","reco dielectron","rec_elkFac") + + + + diff --git a/unfoldTestPlot.py b/unfoldTestPlot.py new file mode 100755 index 0000000..7ca6cf9 --- /dev/null +++ b/unfoldTestPlot.py @@ -0,0 +1,214 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True + +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, SignalsADD, Signals2018ADD, Signals2018, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy + + +def plotDataMC(datahist,mchist,usedata, label1, label2,name,filename): + + + hCanvas = TCanvas("hCanvas", "Distribution", 800,800) + if usedata==True: + plotPad = ROOT.TPad("plotPad","plotPad",0,0.3,1,1) + ratioPad = ROOT.TPad("ratioPad","ratioPad",0,0.,1,0.3) + setTDRStyle() + plotPad.UseCurrentStyle() + ratioPad.UseCurrentStyle() + plotPad.Draw("hist") + ratioPad.Draw("hist") + plotPad.cd() + else: + plotPad = ROOT.TPad("plotPad","plotPad",0,0,1,1) + setTDRStyle() + plotPad.UseCurrentStyle() + plotPad.Draw() + plotPad.cd() + colors = createMyColors() + + legend = TLegend(0.55, 0.6, 0.925, 0.925) + legend.SetFillStyle(0) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legendEta = TLegend(0.45, 0.75, 0.925, 0.925) + legendEta.SetFillStyle(0) + legendEta.SetBorderSize(0) + legendEta.SetTextFont(42) + legendEta.SetNColumns(2) + latex = ROOT.TLatex() + latex.SetTextFont(42) + latex.SetTextAlign(31) + latex.SetTextSize(0.04) + latex.SetNDC(True) + latexCMS = ROOT.TLatex() + latexCMS.SetTextFont(61) + latexCMS.SetTextSize(0.06) + latexCMS.SetNDC(True) + latexCMSExtra = ROOT.TLatex() + latexCMSExtra.SetTextFont(52) + latexCMSExtra.SetTextSize(0.045) + latexCMSExtra.SetNDC(True) + legendHists = [] + legendHistData = ROOT.TH1F() + category = ROOT.TLatex() + category.SetNDC(True) + category.SetTextSize(0.04) + if usedata==True: + legend.AddEntry(legendHistData,label1,"pe") + legendEta.AddEntry(legendHistData,label1,"pe") + + #process.label = process.label.replace("#mu^{+}#mu^{-}","e^{+^{*}}e^{-}") + temphist = ROOT.TH1F() + temphist.SetFillColor(3) + legendHists.append(temphist.Clone) + #legend.AddEntry(temphist,label2,"l") + #legendEta.AddEntry(temphist,process.label,"f") + + # Modify plot pad information + nEvents=-1 + + ROOT.gStyle.SetOptStat(0) + + intlumi = ROOT.TLatex() + intlumi.SetTextAlign(12) + intlumi.SetTextSize(0.045) + intlumi.SetNDC(True) + intlumi2 = ROOT.TLatex() + intlumi2.SetTextAlign(12) + intlumi2.SetTextSize(0.07) + intlumi2.SetNDC(True) + scalelabel = ROOT.TLatex() + scalelabel.SetTextAlign(12) + scalelabel.SetTextSize(0.03) + scalelabel.SetNDC(True) + metDiffLabel = ROOT.TLatex() + metDiffLabel.SetTextAlign(12) + metDiffLabel.SetTextSize(0.03) + metDiffLabel.SetNDC(True) + chi2Label = ROOT.TLatex() + chi2Label.SetTextAlign(12) + chi2Label.SetTextSize(0.03) + chi2Label.SetNDC(True) + hCanvas.SetLogy() + + + # Luminosity information + plotPad.cd() + plotPad.SetLogy(0) + plotPad.SetLogy() + if usedata==True: + yMax = datahist.GetBinContent(datahist.GetMaximumBin())*1000 + yMin = 0.00000001 + xMax = datahist.GetXaxis().GetXmax() + xMin = datahist.GetXaxis().GetXmin() + else: + yMax = mchist.GetBinContent(datahist.GetMaximumBin()) + yMin = 0.00000001 + xMax = mchist.GetXaxis().GetXmax() + xMin = mchist.GetXaxis().GetXmin() + yMax = yMax*10000 + if name.find("dimuon")!=-1: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{#mu#mu}[GeV] ;Events/GeV") + else: + plotPad.DrawFrame(xMin,yMin,xMax,yMax,"; m_{ee}[GeV] ;Events/GeV") + # Draw signal information + + # Draw background from stack + mchist.SetFillColor(0) + mchist.SetLineColor(2) + legend.AddEntry(mchist,label2,"l") + mchist.Draw("samehist") + + # Draw data + datahist.SetMinimum(0.0001) + if usedata==True: + datahist.SetMarkerStyle(8) + datahist.Draw("samepehist") + + # Draw legend + legend.Draw() + plotPad.SetLogx() + latex.DrawLatex(0.95,0.96,"13 TeV") + yLabelPos = 0.85 + cmsExtra = "Preliminary" + if not usedata==True: + cmsExtra = "#splitline{Preliminary}{Simulation}" + yLabelPos = 0.82 + latexCMS.DrawLatex(0.19,0.89,"CMS") + category.DrawLatex(0.3,0.7,name) + latexCMSExtra.DrawLatex(0.19,yLabelPos,"%s"%(cmsExtra)) + #~ print datahist.Integral() + if usedata==True: + try: + ratioPad.cd() + ratioPad.SetLogx() + except AttributeError: + print ("Plot fails. Look up in errs/failedPlots.txt") + outFile =open("errs/failedPlots.txt","a") + outFile.write('%s\n'%plot.filename%("_"+run.label+"_"+dilepton)) + outFile.close() + #plot.cuts=baseCut + return 1 + ratioGraphs = ratios.RatioGraph(datahist,mchist, xMin=xMin, xMax=xMax,title="After/Before",yMin=0,yMax=2.0,ndivisions=10,color=ROOT.kBlack,adaptiveBinning=10000000000000,labelSize=0.125,pull=False) + ratioGraphs.draw(ROOT.gPad,True,False,True,chi2Pos=0.8) + + + ROOT.gPad.RedrawAxis() + plotPad.RedrawAxis() + if usedata==True: + + ratioPad.RedrawAxis() + if not os.path.exists("plots"): + os.makedirs("plots") + hCanvas.Print("plots/%s.jpg"%filename) + + + + +f=ROOT.TFile.Open("RooUnfold/invVsiter.root") + +recBB_18_mu=f.Get("datamubb2018") +unfoldBB_18_mu=f.Get("BB2018mu_invert") +unfoldBB_18_muiter1=f.Get("BB2018muiter_1") +unfoldBB_18_muiter3=f.Get("BB2018muiter_3") +unfoldBB_18_muiter5=f.Get("BB2018muiter_5") +unfoldBB_18_muiter10=f.Get("BB2018muiter_10") +unfoldBB_18_muiter100=f.Get("BB2018muiter_100") +unfoldBB_18_muiter1000=f.Get("BB2018muiter_1000") +unfoldBB_18_muiter10000=f.Get("BB2018muiter_10000") + +plotDataMC(unfoldBB_18_mu,recBB_18_mu,True,"Unfolded","Reconsructed","2018 BB dimuon","data_BB_18_mu_invt" ) +plotDataMC(unfoldBB_18_muiter1,recBB_18_mu,True,"Unfolded_iter1","Reconstructed","2018 BB dimuon","data_BE_18_muiter1" ) +plotDataMC(unfoldBB_18_muiter3,recBB_18_mu,True,"Unfolded_iter3","Reconstructed","2018 BB dimuon","data_BE_18_muiter3" ) +plotDataMC(unfoldBB_18_muiter5,recBB_18_mu,True,"Unfolded_iter5","Reconstructed","2018 BB dimuon","data_BE_18_muiter5" ) +plotDataMC(unfoldBB_18_muiter10,recBB_18_mu,True,"Unfolded_iter10","Reconstructed","2018 BB dimuon","data_BE_18_muiter10" ) +plotDataMC(unfoldBB_18_muiter100,recBB_18_mu,True,"Unfolded_iter100","Reconstructed","2018 BB dimuon","data_BE_18_muiter100" ) +plotDataMC(unfoldBB_18_muiter1000,recBB_18_mu,True,"Unfolded_iter1000","Reconstructed","2018 BB dimuon","data_BE_18_muiter1000" ) +plotDataMC(unfoldBB_18_muiter10000,recBB_18_mu,True,"Unfolded_iter10000","Reconstructed","2018 BB dimuon","data_BE_18_muiter10000" ) + +f2=ROOT.TFile.Open("RooUnfold/invVsiter_obkg.root") +bkgBB_18_mu=f2.Get("bkgmubb2018") +recBB_18_mu1=f2.Get("datamubb2018") +unfolditer=f2.Get("BB2018muiter_10000") +unfoldInv=f2.Get("BB2018mu_invert") +recBB_18_mu1.Add(bkgBB_18_mu,-1) +for i in range(recBB_18_mu.GetNbinsX()+2): + val=recBB_18_mu1.GetBinContent(i) + err=recBB_18_mu1.GetBinError(i) + wid=recBB_18_mu1.GetBinWidth(i) + val=val/wid + err=err/wid + recBB_18_mu1.SetBinContent(i,val) + recBB_18_mu1.SetBinError(i,err) +plotDataMC(unfoldInv,recBB_18_mu1,True,"Unfolded","Reconsructed","2018 BB dimuon","BB_18_mu_invt" ) +plotDataMC(unfolditer,recBB_18_mu1,True,"Unfoldediter10000","Reconstructed","2018 BB dimuon","BB_18_muiter" ) + + diff --git a/unfoldingMC.py b/unfoldingMC.py new file mode 100644 index 0000000..a31a696 --- /dev/null +++ b/unfoldingMC.py @@ -0,0 +1,189 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True +import numpy +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT, TGaxis +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy +import numpy as np +import root_numpy +from copy import deepcopy +from ROOT import TUnfold +def Addhist(histlist): + tempHist=histlist[0] + for i in range(1,3): + tempHist.Add(histlist[i]) + return tempHist +def Addstack(Stacklist): + tempStack=Stacklist[0] + for i in range(1,3): + tempStack.Add(Stacklist[i]) + return tempStack +def Stacks(processes,lumi,plot,zScale): + stacks=[] + for i in range(3): + stacks.append(TheStack(processes[i],lumi[i],plot,zScale[i])) + return stacks +def average(hist): + for i in range(hist.GetNbinsX()+2): + val=hist.GetBinContent(i) + err=hist.GetBinError(i) + wid=hist.GetBinWidth(i) + val=val/wid + err=err/wid + if i==hist.GetNbinsX()+1: + print wid + hist.SetBinContent(i,val) + hist.SetBinError(i,err) +def reAverage(hist): + for i in range(hist.GetNbinsX()+2): + val=hist.GetBinContent(i) + err=hist.GetBinError(i) + wid=hist.GetBinWidth(i) + val=val*wid + err=err*wid + hist.SetBinContent(i,val) + hist.SetBinError(i,err) + +plot_mu_bb = getPlot("massPlotBB") +plot_mu_be = getPlot("massPlotBE") +lumi_mu = [36.3*1000,42.135*1000,61.608*1000] +zScaleFac_mu = [zScale2016["muons"],zScale["muons"],zScale2018["muons"]] +bng=[50, 120,150,200,300,400,500,690,900,1250,1610, 2000, 3000,4000, 6070] +bng1=[50, 120,150,200,300,400,500,690,900,1250,1610, 1970, 3010,3970, 6070] +bng=numpy.asarray(bng,dtype=numpy.float64) +bng1=numpy.asarray(bng1,dtype=numpy.float64) +plot_e_bb = getPlot("massPlotEleBB") +plot_e_be = getPlot("massPlotEleBE") +lumi_e = [35.9*1000,41.529*1000,59.97*1000] +zScaleFac_e = [zScale2016["electrons"],zScale["electrons"],zScale2018["electrons"]] +eventCounts_e = totalNumberOfGeneratedEvents(path,False) +eventCounts_mu = totalNumberOfGeneratedEvents(path,True) +negWeights_mu = negWeightFractions(path,True) +negWeights_e = negWeightFractions(path,False) +processes_mu2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +processes_mu2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +processes_mu2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +stackmu_bb2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_bb,zScaleFac_mu[0]) +stackmu_bb2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_bb,zScaleFac_mu[1]) +stackmu_bb2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_bb,zScaleFac_mu[2]) +stackmu_be2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_be,zScaleFac_mu[0]) +stackmu_be2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_be,zScaleFac_mu[1]) +stackmu_be2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_be,zScaleFac_mu[2]) +stacke_bb2016 = TheStack(processes_e2016,lumi_e[0],plot_e_bb,zScaleFac_e[0][1]) +stacke_bb2017 = TheStack(processes_e2017,lumi_e[1],plot_e_bb,zScaleFac_e[1][1]) +stacke_bb2018 = TheStack(processes_e2018,lumi_e[2],plot_e_bb,zScaleFac_e[2][1]) +stacke_be2016 = TheStack(processes_e2016,lumi_e[0],plot_e_be,zScaleFac_e[0][2]) +stacke_be2017 = TheStack(processes_e2017,lumi_e[1],plot_e_be,zScaleFac_e[1][2]) +stackmu_be2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_be,zScaleFac_mu[2]) +stacke_bb2016 = TheStack(processes_e2016,lumi_e[0],plot_e_bb,zScaleFac_e[0][1]) +stacke_bb2017 = TheStack(processes_e2017,lumi_e[1],plot_e_bb,zScaleFac_e[1][1]) +stacke_bb2018 = TheStack(processes_e2018,lumi_e[2],plot_e_bb,zScaleFac_e[2][1]) +stacke_be2016 = TheStack(processes_e2016,lumi_e[0],plot_e_be,zScaleFac_e[0][2]) +stacke_be2017 = TheStack(processes_e2017,lumi_e[1],plot_e_be,zScaleFac_e[1][2]) +stacke_be2018 = TheStack(processes_e2018,lumi_e[2],plot_e_be,zScaleFac_e[2][2]) +bkgmubb=[stackmu_bb2016.theHistogram,stackmu_bb2017.theHistogram,stackmu_bb2018.theHistogram] +bkgmube=[stackmu_be2016.theHistogram,stackmu_be2017.theHistogram,stackmu_be2018.theHistogram] +bkgebb=[stacke_bb2016.theHistogram,stacke_bb2017.theHistogram,stacke_bb2018.theHistogram] +bkgebe=[stacke_be2016.theHistogram,stacke_be2017.theHistogram,stacke_be2018.theHistogram] +bkgebb2016=ROOT.TH1D("bkgebb2016","bkgebb2016",len(bng)-1,bng) +bkgebb2017=ROOT.TH1D("bkgebb2017","bkgebb2017",len(bng)-1,bng) +bkgebb2018=ROOT.TH1D("bkgebb2018","bkgebb2018",len(bng)-1,bng) +bkgebe2016=ROOT.TH1D("bkgebe2016","bkgebe2016",len(bng)-1,bng) +bkgebe2017=ROOT.TH1D("bkgebe2017","bkgebe2017",len(bng)-1,bng) +bkgebe2018=ROOT.TH1D("bkgebe2018","bkgebe2018",len(bng)-1,bng) +reAverage(bkgebb[0]) +bkgebb0=bkgebb[0].Rebin(len(bng1)-1,"",bng1) +reAverage(bkgebb[1]) +bkgebb1=bkgebb[1].Rebin(len(bng1)-1,"",bng1) +reAverage(bkgebb[2]) +bkgebb2=bkgebb[2].Rebin(len(bng1)-1,"",bng1) +reAverage(bkgebe[0]) +bkgebe0=bkgebe[0].Rebin(len(bng1)-1,"",bng1) +reAverage(bkgebe[1]) +bkgebe1=bkgebe[1].Rebin(len(bng1)-1,"",bng1) +reAverage(bkgebe[2]) +bkgebe2=bkgebe[2].Rebin(len(bng1)-1,"",bng1) +for i in range(len(bng)+1): + val1=0 + val2=0 + #val=bkgmubb[0].GetBinContent(i) + #err=bkgmubb[0].GetBinError(i) + #bkgmubb2016.SetBinContent(i,val) + #bkgmubb2016.SetBinError(i,err) + #val=bkgmubb[1].GetBinContent(i) + #err=bkgmubb[1].GetBinError(i) + #bkgmubb2017.SetBinContent(i,val) + #bkgmubb2017.SetBinError(i,err) + #val=bkgmubb[2].GetBinContent(i) + #err=bkgmubb[2].GetBinError(i) + #bkgmubb2018.SetBinContent(i,val) + #bkgmubb2018.SetBinError(i,err) + #val=bkgmube[0].GetBinContent(i) + #err=bkgmube[0].GetBinError(i) + #bkgmube2016.SetBinContent(i,val) + #bkgmube2016.SetBinError(i,err) + #val=bkgmube[1].GetBinContent(i) + #err=bkgmube[1].GetBinError(i) + #bkgmube2017.SetBinContent(i,val) + #bkgmube2017.SetBinError(i,err) + #val=bkgmube[2].GetBinContent(i) + #err=bkgmube[2].GetBinError(i) + #bkgmube2018.SetBinContent(i,val) + #bkgmube2018.SetBinError(i,err) + val=bkgebb0.GetBinContent(i) + err=bkgebb0.GetBinError(i) + val1+=val + bkgebb2016.SetBinContent(i,val) + bkgebb2016.SetBinError(i,err) + val=bkgebb1.GetBinContent(i) + err=bkgebb1.GetBinError(i) + val1+=val + bkgebb2017.SetBinContent(i,val) + bkgebb2017.SetBinError(i,err) + val=bkgebb2.GetBinContent(i) + err=bkgebb2.GetBinError(i) + val1+=val + bkgebb2018.SetBinContent(i,val) + bkgebb2018.SetBinError(i,err) + val=bkgebe0.GetBinContent(i) + err=bkgebe0.GetBinError(i) + val2+=val + bkgebe2016.SetBinContent(i,val) + bkgebe2016.SetBinError(i,err) + val=bkgebe1.GetBinContent(i) + err=bkgebe1.GetBinError(i) + val2+=val + bkgebe2017.SetBinContent(i,val) + bkgebe2017.SetBinError(i,err) + val=bkgebe2.GetBinContent(i) + err=bkgebe2.GetBinError(i) + val2+=val + bkgebe2018.SetBinContent(i,val) + bkgebe2018.SetBinError(i,err) + print val1 + print val2 +average(bkgebb2016) +average(bkgebb2017) +average(bkgebb2018) +average(bkgebe2016) +average(bkgebe2017) +average(bkgebe2018) +f=ROOT.TFile("unfoldingMC_V3.root","RECREATE") +bkgebb2016.Write() +bkgebb2017.Write() +bkgebb2018.Write() +bkgebe2016.Write() +bkgebe2017.Write() +bkgebe2018.Write() +f.Close() + diff --git a/unfoldingSample.py b/unfoldingSample.py new file mode 100644 index 0000000..16889d4 --- /dev/null +++ b/unfoldingSample.py @@ -0,0 +1,260 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True +import numpy +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT, TGaxis +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy +import numpy as np +import root_numpy +from copy import deepcopy +from ROOT import TUnfold +def Addhist(histlist): + tempHist=histlist[0] + for i in range(1,3): + tempHist.Add(histlist[i]) + return tempHist +def Addstack(Stacklist): + tempStack=Stacklist[0] + for i in range(1,3): + tempStack.Add(Stacklist[i]) + return tempStack +def Stacks(processes,lumi,plot,zScale): + stacks=[] + for i in range(3): + stacks.append(TheStack(processes[i],lumi[i],plot,zScale[i])) + return stacks + +datamubb=[] +datamube=[] +data_mu=[Process(Data2016, normalized=True),Process(Data, normalized=True),Process(Data2018, normalized=True)] +plot_mu_bb = getPlot("massPlotBB") +plot_mu_be = getPlot("massPlotBE") +lumi_mu = [36.3*1000,42.135*1000,61.608*1000] +zScaleFac_mu = [zScale2016["muons"],zScale["muons"],zScale2018["muons"]] +for i in range(3): + datamubb.append(data_mu[i].loadHistogram(plot_mu_bb,lumi_mu[i],zScaleFac_mu[i])) + datamube.append(data_mu[i].loadHistogram(plot_mu_be,lumi_mu[i],zScaleFac_mu[i])) + +bng=[50, 120,150,200,300,400,500,690,900,1250,1610, 2000, 4000, 6070] +#bng=[60,120,400,600,900,1300,1800,6000] +bng=numpy.asarray(bng,dtype=numpy.float64) +datamubb2016=ROOT.TH1D("datamubb2016","datamubb2016",len(bng)-1,bng) +datamubb2017=ROOT.TH1D("datamubb2017","datamubb2017",len(bng)-1,bng) +datamubb2018=ROOT.TH1D("datamubb2018","datamubb2018",len(bng)-1,bng) +datamube2016=ROOT.TH1D("datamube2016","datamube2016",len(bng)-1,bng) +datamube2017=ROOT.TH1D("datamube2017","datamube2017",len(bng)-1,bng) +datamube2018=ROOT.TH1D("datamube2018","datamube2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val=datamubb[0].GetBinContent(i) + err=datamubb[0].GetBinError(i) + datamubb2016.SetBinContent(i,val) + datamubb2016.SetBinError(i,err) + val=datamubb[1].GetBinContent(i) + err=datamubb[1].GetBinError(i) + datamubb2017.SetBinContent(i,val) + datamubb2017.SetBinError(i,err) + val=datamubb[2].GetBinContent(i) + err=datamubb[2].GetBinError(i) + datamubb2018.SetBinContent(i,val) + datamubb2018.SetBinError(i,err) + val=datamube[0].GetBinContent(i) + err=datamube[0].GetBinError(i) + datamube2016.SetBinContent(i,val) + datamube2016.SetBinError(i,err) + val=datamube[1].GetBinContent(i) + err=datamube[1].GetBinError(i) + datamube2017.SetBinContent(i,val) + datamube2017.SetBinError(i,err) + val=datamube[2].GetBinContent(i) + err=datamube[2].GetBinError(i) + datamube2018.SetBinContent(i,val) + datamube2018.SetBinError(i,err) + +dataebb=[] +dataebe=[] +data_e=[Process(Data2016, normalized=True),Process(Data, normalized=True),Process(Data2018, normalized=True)] +plot_e_bb = getPlot("massPlotEleBB") +plot_e_be = getPlot("massPlotEleBE") +lumi_e = [35.9*1000,41.529*1000,59.97*1000] +zScaleFac_e = [zScale2016["electrons"],zScale["electrons"],zScale2018["electrons"]] +for i in range(3): + dataebb.append(data_e[i].loadHistogram(plot_e_bb,lumi_e[i],zScaleFac_e[i][0])) + dataebe.append(data_e[i].loadHistogram(plot_e_be,lumi_e[i],zScaleFac_e[i][0])) + +dataebb2016=ROOT.TH1D("dataebb2016","dataebb2016",len(bng)-1,bng) +dataebb2017=ROOT.TH1D("dataebb2017","dataebb2017",len(bng)-1,bng) +dataebb2018=ROOT.TH1D("dataebb2018","dataebb2018",len(bng)-1,bng) +dataebe2016=ROOT.TH1D("dataebe2016","dataebe2016",len(bng)-1,bng) +dataebe2017=ROOT.TH1D("dataebe2017","dataebe2017",len(bng)-1,bng) +dataebe2018=ROOT.TH1D("dataebe2018","dataebe2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val=dataebb[0].GetBinContent(i) + err=dataebb[0].GetBinError(i) + dataebb2016.SetBinContent(i,val) + dataebb2016.SetBinError(i,err) + val=dataebb[1].GetBinContent(i) + err=dataebb[1].GetBinError(i) + dataebb2017.SetBinContent(i,val) + dataebb2017.SetBinError(i,err) + val=dataebb[2].GetBinContent(i) + err=dataebb[2].GetBinError(i) + dataebb2018.SetBinContent(i,val) + dataebb2018.SetBinError(i,err) + val=dataebe[0].GetBinContent(i) + err=dataebe[0].GetBinError(i) + dataebe2016.SetBinContent(i,val) + dataebe2016.SetBinError(i,err) + val=dataebe[1].GetBinContent(i) + err=dataebe[1].GetBinError(i) + dataebe2017.SetBinContent(i,val) + dataebe2017.SetBinError(i,err) + val=dataebe[2].GetBinContent(i) + err=dataebe[2].GetBinError(i) + dataebe2018.SetBinContent(i,val) + dataebe2018.SetBinError(i,err) + +backgrounds=["Wjets","Other"] + +f=ROOT.TFile("unfoldingData.root","RECREATE") +datamubb2016.Write() +datamubb2017.Write() +datamubb2018.Write() +datamube2016.Write() +datamube2017.Write() +datamube2018.Write() +dataebb2016.Write() +dataebb2017.Write() +dataebb2018.Write() +dataebe2016.Write() +dataebe2017.Write() +dataebe2018.Write() +f.Close() + +print plot_mu_be.useJets +print plot_e_be.useJets +print plot_e_bb.fileName +print plot_e_be.fileName +eventCounts_e = totalNumberOfGeneratedEvents(path,False) +eventCounts_mu = totalNumberOfGeneratedEvents(path,True) +print eventCounts_e +negWeights_mu = negWeightFractions(path,True) +negWeights_e = negWeightFractions(path,False) +processes_mu2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_mu,negWeights_mu,normalized=True),Process(getattr(Backgrounds2016,"Other"),eventCounts_mu,negWeights_mu)] +processes_e2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_e,negWeights_e,normalized=True),Process(getattr(Backgrounds2016,"OtherEle"),eventCounts_e,negWeights_e)] +processes_mu2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_mu,negWeights_mu,normalized=True),Process(getattr(Backgrounds,"Other"),eventCounts_mu,negWeights_mu)] +processes_e2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_e,negWeights_e,normalized=True),Process(getattr(Backgrounds,"Other"),eventCounts_e,negWeights_e)] +processes_mu2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_mu,negWeights_mu,normalized=True),Process(getattr(Backgrounds2018,"Other"),eventCounts_mu,negWeights_mu)] +processes_e2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_e,negWeights_e,normalized=True),Process(getattr(Backgrounds2018,"Other"),eventCounts_e,negWeights_e)] +stackmu_bb2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_bb,zScaleFac_mu[0]) +stackmu_bb2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_bb,zScaleFac_mu[1]) +stackmu_bb2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_bb,zScaleFac_mu[2]) +stackmu_be2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_be,zScaleFac_mu[0]) +stackmu_be2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_be,zScaleFac_mu[1]) +stackmu_be2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_be,zScaleFac_mu[2]) +stacke_bb2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_bb,zScaleFac_e[0][1]) +stacke_bb2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_bb,zScaleFac_e[1][1]) +stacke_bb2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_bb,zScaleFac_e[2][1]) +stacke_be2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_be,zScaleFac_e[0][2]) +stacke_be2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_be,zScaleFac_e[1][2]) +stacke_be2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_be,zScaleFac_e[2][2]) +bkgmubb=[stackmu_bb2016.theHistogram,stackmu_bb2017.theHistogram,stackmu_bb2018.theHistogram] +bkgmube=[stackmu_be2016.theHistogram,stackmu_be2017.theHistogram,stackmu_be2018.theHistogram] +bkgebb=[stacke_bb2016.theHistogram,stacke_bb2017.theHistogram,stacke_bb2018.theHistogram] +bkgebe=[stacke_be2016.theHistogram,stacke_be2017.theHistogram,stacke_be2018.theHistogram] +bkgebb2016=ROOT.TH1D("bkgebb2016","bkgebb2016",len(bng)-1,bng) +bkgebb2017=ROOT.TH1D("bkgebb2017","bkgebb2017",len(bng)-1,bng) +bkgebb2018=ROOT.TH1D("bkgebb2018","bkgebb2018",len(bng)-1,bng) +bkgebe2016=ROOT.TH1D("bkgebe2016","bkgebe2016",len(bng)-1,bng) +bkgebe2017=ROOT.TH1D("bkgebe2017","bkgebe2017",len(bng)-1,bng) +bkgebe2018=ROOT.TH1D("bkgebe2018","bkgebe2018",len(bng)-1,bng) +bkgmubb2016=ROOT.TH1D("bkgmubb2016","bkgmubb2016",len(bng)-1,bng) +bkgmubb2017=ROOT.TH1D("bkgmubb2017","bkgmubb2017",len(bng)-1,bng) +bkgmubb2018=ROOT.TH1D("bkgmubb2018","bkgmubb2018",len(bng)-1,bng) +bkgmube2016=ROOT.TH1D("bkgmube2016","bkgmube2016",len(bng)-1,bng) +bkgmube2017=ROOT.TH1D("bkgmube2017","bkgmube2017",len(bng)-1,bng) +bkgmube2018=ROOT.TH1D("bkgmube2018","bkgmube2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val1=0 + val2=0 + val=bkgmubb[0].GetBinContent(i) + err=bkgmubb[0].GetBinError(i) + bkgmubb2016.SetBinContent(i,val) + bkgmubb2016.SetBinError(i,err) + val=bkgmubb[1].GetBinContent(i) + err=bkgmubb[1].GetBinError(i) + bkgmubb2017.SetBinContent(i,val) + bkgmubb2017.SetBinError(i,err) + val=bkgmubb[2].GetBinContent(i) + err=bkgmubb[2].GetBinError(i) + bkgmubb2018.SetBinContent(i,val) + bkgmubb2018.SetBinError(i,err) + val=bkgmube[0].GetBinContent(i) + err=bkgmube[0].GetBinError(i) + bkgmube2016.SetBinContent(i,val) + bkgmube2016.SetBinError(i,err) + val=bkgmube[1].GetBinContent(i) + err=bkgmube[1].GetBinError(i) + bkgmube2017.SetBinContent(i,val) + bkgmube2017.SetBinError(i,err) + val=bkgmube[2].GetBinContent(i) + err=bkgmube[2].GetBinError(i) + bkgmube2018.SetBinContent(i,val) + bkgmube2018.SetBinError(i,err) + val=bkgebb[0].GetBinContent(i) + val1+=val + err=bkgebb[0].GetBinError(i) + bkgebb2016.SetBinContent(i,val) + bkgebb2016.SetBinError(i,err) + val=bkgebb[1].GetBinContent(i) + val1+=val + err=bkgebb[1].GetBinError(i) + bkgebb2017.SetBinContent(i,val) + bkgebb2017.SetBinError(i,err) + val=bkgebb[2].GetBinContent(i) + err=bkgebb[2].GetBinError(i) + val1+=val + bkgebb2018.SetBinContent(i,val) + bkgebb2018.SetBinError(i,err) + val=bkgebe[0].GetBinContent(i) + err=bkgebe[0].GetBinError(i) + val2+=val + bkgebe2016.SetBinContent(i,val) + bkgebe2016.SetBinError(i,err) + val=bkgebe[1].GetBinContent(i) + err=bkgebe[1].GetBinError(i) + val2+=val + bkgebe2017.SetBinContent(i,val) + bkgebe2017.SetBinError(i,err) + val=bkgebe[2].GetBinContent(i) + err=bkgebe[2].GetBinError(i) + val2+=val + bkgebe2018.SetBinContent(i,val) + bkgebe2018.SetBinError(i,err) + print val1 + print val2 +f=ROOT.TFile("unfoldingMC.root","RECREATE") +bkgmubb2016.Write() +bkgmubb2017.Write() +bkgmubb2018.Write() +bkgmube2016.Write() +bkgmube2017.Write() +bkgmube2018.Write() +bkgebb2016.Write() +bkgebb2017.Write() +bkgebb2018.Write() +bkgebe2016.Write() +bkgebe2017.Write() +bkgebe2018.Write() +f.Close() + + + + + diff --git a/unfoldingSample_V2.py b/unfoldingSample_V2.py new file mode 100644 index 0000000..07d5346 --- /dev/null +++ b/unfoldingSample_V2.py @@ -0,0 +1,266 @@ +import argparse +import ROOT +ROOT.PyConfig.IgnoreCommandLineOptions = True +import numpy +from ROOT import TCanvas, TPad, TH1F, TH1I, THStack, TLegend, TMath, gROOT, TGaxis +import ratios +from setTDRStyle import setTDRStyle +gROOT.SetBatch(True) +from helpers import * +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, Data, Data2016, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import math +import os +from copy import copy +import numpy as np +import root_numpy +from copy import deepcopy +from ROOT import TUnfold +def Addhist(histlist): + tempHist=histlist[0] + for i in range(1,3): + tempHist.Add(histlist[i]) + return tempHist +def Addstack(Stacklist): + tempStack=Stacklist[0] + for i in range(1,3): + tempStack.Add(Stacklist[i]) + return tempStack +def Stacks(processes,lumi,plot,zScale): + stacks=[] + for i in range(3): + stacks.append(TheStack(processes[i],lumi[i],plot,zScale[i])) + return stacks + +datamubb=[] +datamube=[] +data_mu=[Process(Data2016, normalized=True),Process(Data, normalized=True),Process(Data2018, normalized=True)] +plot_mu_bb = getPlot("massPlotBB") +plot_mu_be = getPlot("massPlotBE") +lumi_mu = [36.3*1000,42.135*1000,61.608*1000] +zScaleFac_mu = [zScale2016["muons"],zScale["muons"],zScale2018["muons"]] +for i in range(3): + datamubb.append(data_mu[i].loadHistogram(plot_mu_bb,lumi_mu[i],zScaleFac_mu[i])) + datamube.append(data_mu[i].loadHistogram(plot_mu_be,lumi_mu[i],zScaleFac_mu[i])) + +bng=[50, 120,150,200,300,400,500,690,900,1250,1610, 2000, 3000,4000, 6070] +#bng=[60,120,400,600,900,1300,1800,6000] +bng=numpy.asarray(bng,dtype=numpy.float64) +datamubb2016=ROOT.TH1D("datamubb2016","datamubb2016",len(bng)-1,bng) +datamubb2017=ROOT.TH1D("datamubb2017","datamubb2017",len(bng)-1,bng) +datamubb2018=ROOT.TH1D("datamubb2018","datamubb2018",len(bng)-1,bng) +datamube2016=ROOT.TH1D("datamube2016","datamube2016",len(bng)-1,bng) +datamube2017=ROOT.TH1D("datamube2017","datamube2017",len(bng)-1,bng) +datamube2018=ROOT.TH1D("datamube2018","datamube2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val=datamubb[0].GetBinContent(i) + err=datamubb[0].GetBinError(i) + datamubb2016.SetBinContent(i,val) + datamubb2016.SetBinError(i,err) + val=datamubb[1].GetBinContent(i) + err=datamubb[1].GetBinError(i) + datamubb2017.SetBinContent(i,val) + datamubb2017.SetBinError(i,err) + val=datamubb[2].GetBinContent(i) + err=datamubb[2].GetBinError(i) + datamubb2018.SetBinContent(i,val) + datamubb2018.SetBinError(i,err) + val=datamube[0].GetBinContent(i) + err=datamube[0].GetBinError(i) + datamube2016.SetBinContent(i,val) + datamube2016.SetBinError(i,err) + val=datamube[1].GetBinContent(i) + err=datamube[1].GetBinError(i) + datamube2017.SetBinContent(i,val) + datamube2017.SetBinError(i,err) + val=datamube[2].GetBinContent(i) + err=datamube[2].GetBinError(i) + datamube2018.SetBinContent(i,val) + datamube2018.SetBinError(i,err) + +dataebb=[] +dataebe=[] +data_e=[Process(Data2016, normalized=True),Process(Data, normalized=True),Process(Data2018, normalized=True)] +plot_e_bb = getPlot("massPlotEleBB") +plot_e_be = getPlot("massPlotEleBE") +lumi_e = [35.9*1000,41.529*1000,59.97*1000] +zScaleFac_e = [zScale2016["electrons"],zScale["electrons"],zScale2018["electrons"]] +for i in range(3): + dataebb.append(data_e[i].loadHistogram(plot_e_bb,lumi_e[i],zScaleFac_e[i][0])) + dataebe.append(data_e[i].loadHistogram(plot_e_be,lumi_e[i],zScaleFac_e[i][0])) + +dataebb2016=ROOT.TH1D("dataebb2016","dataebb2016",len(bng)-1,bng) +dataebb2017=ROOT.TH1D("dataebb2017","dataebb2017",len(bng)-1,bng) +dataebb2018=ROOT.TH1D("dataebb2018","dataebb2018",len(bng)-1,bng) +dataebe2016=ROOT.TH1D("dataebe2016","dataebe2016",len(bng)-1,bng) +dataebe2017=ROOT.TH1D("dataebe2017","dataebe2017",len(bng)-1,bng) +dataebe2018=ROOT.TH1D("dataebe2018","dataebe2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val=dataebb[0].GetBinContent(i) + err=dataebb[0].GetBinError(i) + dataebb2016.SetBinContent(i,val) + dataebb2016.SetBinError(i,err) + val=dataebb[1].GetBinContent(i) + err=dataebb[1].GetBinError(i) + dataebb2017.SetBinContent(i,val) + dataebb2017.SetBinError(i,err) + val=dataebb[2].GetBinContent(i) + err=dataebb[2].GetBinError(i) + dataebb2018.SetBinContent(i,val) + dataebb2018.SetBinError(i,err) + val=dataebe[0].GetBinContent(i) + err=dataebe[0].GetBinError(i) + dataebe2016.SetBinContent(i,val) + dataebe2016.SetBinError(i,err) + val=dataebe[1].GetBinContent(i) + err=dataebe[1].GetBinError(i) + dataebe2017.SetBinContent(i,val) + dataebe2017.SetBinError(i,err) + val=dataebe[2].GetBinContent(i) + err=dataebe[2].GetBinError(i) + dataebe2018.SetBinContent(i,val) + dataebe2018.SetBinError(i,err) + +backgrounds=["Wjets","Other"] + +f=ROOT.TFile("unfoldingData.root","RECREATE") +datamubb2016.Write() +datamubb2017.Write() +datamubb2018.Write() +datamube2016.Write() +datamube2017.Write() +datamube2018.Write() +dataebb2016.Write() +dataebb2017.Write() +dataebb2018.Write() +dataebe2016.Write() +dataebe2017.Write() +dataebe2018.Write() +f.Close() + +print plot_mu_be.useJets +print plot_e_be.useJets +print plot_e_bb.fileName +print plot_e_be.fileName +eventCounts_e = totalNumberOfGeneratedEvents(path,False) +eventCounts_mu = totalNumberOfGeneratedEvents(path,True) +print eventCounts_e +negWeights_mu = negWeightFractions(path,True) +negWeights_e = negWeightFractions(path,False) +processes_mu2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2016=[Process(getattr(Backgrounds2016,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +processes_mu2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2017=[Process(getattr(Backgrounds,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +processes_mu2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_mu,negWeights_mu,normalized=True)] +processes_e2018=[Process(getattr(Backgrounds2018,"Jets"),eventCounts_e,negWeights_e,normalized=True)] +stackmu_bb2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_bb,zScaleFac_mu[0]) +stackmu_bb2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_bb,zScaleFac_mu[1]) +stackmu_bb2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_bb,zScaleFac_mu[2]) +stackmu_be2016 = TheStack(processes_mu2016,lumi_mu[0],plot_mu_be,zScaleFac_mu[0]) +stackmu_be2017 = TheStack(processes_mu2017,lumi_mu[1],plot_mu_be,zScaleFac_mu[1]) +stackmu_be2018 = TheStack(processes_mu2018,lumi_mu[2],plot_mu_be,zScaleFac_mu[2]) +stacke_bb2016 = TheStack(processes_e2016,lumi_e[0],plot_e_bb,zScaleFac_e[0][1]) +stacke_bb2017 = TheStack(processes_e2017,lumi_e[1],plot_e_bb,zScaleFac_e[1][1]) +stacke_bb2018 = TheStack(processes_e2018,lumi_e[2],plot_e_bb,zScaleFac_e[2][1]) +stacke_be2016 = TheStack(processes_e2016,lumi_e[0],plot_e_be,zScaleFac_e[0][2]) +stacke_be2017 = TheStack(processes_e2017,lumi_e[1],plot_e_be,zScaleFac_e[1][2]) +stacke_be2018 = TheStack(processes_e2018,lumi_e[2],plot_e_be,zScaleFac_e[2][2]) +bkgmubb=[stackmu_bb2016.theHistogram,stackmu_bb2017.theHistogram,stackmu_bb2018.theHistogram] +bkgmube=[stackmu_be2016.theHistogram,stackmu_be2017.theHistogram,stackmu_be2018.theHistogram] +bkgebb=[stacke_bb2016.theHistogram,stacke_bb2017.theHistogram,stacke_bb2018.theHistogram] +bkgebe=[stacke_be2016.theHistogram,stacke_be2017.theHistogram,stacke_be2018.theHistogram] +bkgebb2016=ROOT.TH1D("bkgebb2016","bkgebb2016",len(bng)-1,bng) +bkgebb2017=ROOT.TH1D("bkgebb2017","bkgebb2017",len(bng)-1,bng) +bkgebb2018=ROOT.TH1D("bkgebb2018","bkgebb2018",len(bng)-1,bng) +bkgebe2016=ROOT.TH1D("bkgebe2016","bkgebe2016",len(bng)-1,bng) +bkgebe2017=ROOT.TH1D("bkgebe2017","bkgebe2017",len(bng)-1,bng) +bkgebe2018=ROOT.TH1D("bkgebe2018","bkgebe2018",len(bng)-1,bng) +bkgmubb2016=ROOT.TH1D("bkgmubb2016","bkgmubb2016",len(bng)-1,bng) +bkgmubb2017=ROOT.TH1D("bkgmubb2017","bkgmubb2017",len(bng)-1,bng) +bkgmubb2018=ROOT.TH1D("bkgmubb2018","bkgmubb2018",len(bng)-1,bng) +bkgmube2016=ROOT.TH1D("bkgmube2016","bkgmube2016",len(bng)-1,bng) +bkgmube2017=ROOT.TH1D("bkgmube2017","bkgmube2017",len(bng)-1,bng) +bkgmube2018=ROOT.TH1D("bkgmube2018","bkgmube2018",len(bng)-1,bng) +for i in range(len(bng)+1): + val1=0 + val2=0 + val=bkgmubb[0].GetBinContent(i) + err=bkgmubb[0].GetBinError(i) + bkgmubb2016.SetBinContent(i,val) + bkgmubb2016.SetBinError(i,err) + bkgmubb2016.Scale(0.5) + val1+=val + val=bkgmubb[1].GetBinContent(i) + err=bkgmubb[1].GetBinError(i) + bkgmubb2017.SetBinContent(i,val) + bkgmubb2017.SetBinError(i,err) + bkgmubb2017.Scale(0.5) + val1+=val + val=bkgmubb[2].GetBinContent(i) + err=bkgmubb[2].GetBinError(i) + bkgmubb2018.SetBinContent(i,val) + bkgmubb2018.SetBinError(i,err) + bkgmubb2018.Scale(0.5) + val1+=val + val=bkgmube[0].GetBinContent(i) + err=bkgmube[0].GetBinError(i) + bkgmube2016.SetBinContent(i,val) + bkgmube2016.SetBinError(i,err) + bkgmube2016.Scale(0.5) + val2+=val + val=bkgmube[1].GetBinContent(i) + err=bkgmube[1].GetBinError(i) + bkgmube2017.SetBinContent(i,val) + bkgmube2017.SetBinError(i,err) + bkgmube2017.Scale(0.5) + val2+=val + val=bkgmube[2].GetBinContent(i) + err=bkgmube[2].GetBinError(i) + bkgmube2018.SetBinContent(i,val) + bkgmube2018.SetBinError(i,err) + bkgmube2018.Scale(0.5) + val2+=val + val=bkgebb[0].GetBinContent(i) + err=bkgebb[0].GetBinError(i) + bkgebb2016.SetBinContent(i,val) + bkgebb2016.SetBinError(i,err) + val=bkgebb[1].GetBinContent(i) + err=bkgebb[1].GetBinError(i) + bkgebb2017.SetBinContent(i,val) + bkgebb2017.SetBinError(i,err) + val=bkgebb[2].GetBinContent(i) + err=bkgebb[2].GetBinError(i) + bkgebb2018.SetBinContent(i,val) + bkgebb2018.SetBinError(i,err) + val=bkgebe[0].GetBinContent(i) + err=bkgebe[0].GetBinError(i) + bkgebe2016.SetBinContent(i,val) + bkgebe2016.SetBinError(i,err) + val=bkgebe[1].GetBinContent(i) + err=bkgebe[1].GetBinError(i) + bkgebe2017.SetBinContent(i,val) + bkgebe2017.SetBinError(i,err) + val=bkgebe[2].GetBinContent(i) + err=bkgebe[2].GetBinError(i) + bkgebe2018.SetBinContent(i,val) + bkgebe2018.SetBinError(i,err) + print val1 + print val2 +f=ROOT.TFile("unfoldingMC_V2.root","RECREATE") +bkgmubb2016.Write() +bkgmubb2017.Write() +bkgmubb2018.Write() +bkgmube2016.Write() +bkgmube2017.Write() +bkgmube2018.Write() +#bkgebb2016.Write() +#bkgebb2017.Write() +#bkgebb2018.Write() +#bkgebe2016.Write() +#bkgebe2017.Write() +#bkgebe2018.Write() +f.Close() + + + + + From 1799e93f20b71d44eab756ab31f9aeb475baeba8 Mon Sep 17 00:00:00 2001 From: "[Minxi" <[yhzz598@gmail.com]> Date: Thu, 10 Jun 2021 15:14:40 -0400 Subject: [PATCH 4/4] code produce the input for the datacard --- produceCollection.py | 469 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 469 insertions(+) create mode 100644 produceCollection.py diff --git a/produceCollection.py b/produceCollection.py new file mode 100644 index 0000000..34df184 --- /dev/null +++ b/produceCollection.py @@ -0,0 +1,469 @@ +from defs import getPlot, Backgrounds, Backgrounds2016, Backgrounds2018, Signals, Signals2016, Signals2016ADD, Data, Data2016, Data2016el, Data2018, path, plotList, zScale, zScale2016, zScale2018 +import ROOT +from helpers import * + + + + +#savePath="/depot/cms/users/minxi/lepratio/LimitSetting/data/" +#f=ROOT.TFile.Open(savePath+"NewCollection.root","RECREATE") + +def getTag(flavor, year, uncertainty,category): + + print(flavor+"_"+str(year)+"_"+uncertainty+"_"+category) + flavors={"mu":1000,"el":2000} + years={"2016":100,"2017":200,"2018":300} + categories={"bb":10,"be":20} + uncertainties={"scaleUp":1,"scaleDown":1,"reso":2,"ID":3,"PUUp":4,"PUDown":4,"prefireUp":5,"prefireDown":5} + if flavor=="mu": + if uncertainty=="scaleUp" or uncertainty=="scaleDown": + category="bb" + if year==2018: year=2017 + elif uncertainty=="ID" or uncertainty=="reso": + if category=="bb": + if year==2018: + year=2017 + + elif flavor=="el": + if uncertainty.find("PU")!=-1 or uncertainty.find("prefire")!=-1: + year=2016 + category="bb" + + tag=flavors[flavor]+years[str(year)]+uncertainties[uncertainty]+categories[category] + print(str(tag)) + return str(tag) + + + + + +if __name__ == "__main__": + + savePath="/depot/cms/users/minxi/lepratio/LimitSetting/data/" + #savePath="/depot/cms/users/minxi/study_novfit/mass_spectrum/" + bins=binning() + lumi_el = [35.9*1000,41.529*1000,59.97*1000] + lumi_mu = [36.3*1000,42.135*1000,61.608*1000] + muplots_bb={"default":"massPlotBB","scaleUp":"massPlotBBScaleUpNoLog","scaleDown":"massPlotBBScaleDownNoLog","reso":"massPlotBBSmearNoLog","ID":"massPlotBBMuonIDNoLog"} + muplots_be={"default":"massPlotBE","scaleUp":"massPlotBEScaleUpNoLog","scaleDown":"massPlotBEScaleDownNoLog","reso":"massPlotBESmearNoLog","ID":"massPlotBEMuonIDNoLog"} + eleplots_bb={"default":"massPlotEleBB","scaleUp":"massPlotEleBBScaleUpNoLog","scaleDown":"massPlotEleBBScaleDownNoLog","PUUp":"massPlotEleBBPUScaleUpNoLog","PUDown":"massPlotEleBBPUScaleDownNoLog","prefireUp":"massPlotEleBBPrefireUpNoLog","prefireDown":"massPlotEleBBPrefireDownNoLog"} + eleplots_be={"default":"massPlotEleBE","scaleUp":"massPlotEleBEScaleUpNoLog","scaleDown":"massPlotEleBEScaleDownNoLog","PUUp":"massPlotEleBEPUScaleUpNoLog","PUDown":"massPlotEleBEPUScaleDownNoLog","prefireUp":"massPlotEleBEPrefireUpNoLog","prefireDown":"massPlotEleBEPrefireDownNoLog"} + plot_mu_bb={} + plot_el_bb={} + plot_mu_be={} + plot_el_be={} + plotGenbb_mu = getPlot("massPlotBBGenNoLog") + plotGenbe_mu = getPlot("massPlotBEGenNoLog") + plotGenbb_el = getPlot("massPlotEleBBGenNoLog") + plotGenbe_el = getPlot("massPlotEleBEGenNoLog") + plotGenee_el = getPlot("massPlotEleEEGenNoLog") + for key in muplots_bb.keys(): + plot_mu_bb[key] = getPlot(muplots_bb[key]) + plot_mu_bb[key].logX=True + for key in eleplots_bb.keys(): + plot_el_bb[key] = getPlot(eleplots_bb[key]) + plot_el_bb[key].logX=True + for key in muplots_be.keys(): + plot_mu_be[key] = getPlot(muplots_be[key]) + plot_mu_be[key].logX=True + for key in eleplots_be.keys(): + plot_el_be[key] = getPlot(eleplots_be[key]) + plot_el_be[key].logX=True + eventCounts_mu = totalNumberOfGeneratedEvents(path,plot_mu_bb["default"].muon) + eventCounts_el = totalNumberOfGeneratedEvents(path,plot_el_bb["default"].muon) + negWeights_mu = negWeightFractions(path,plot_mu_bb["default"].muon) + negWeights_el = negWeightFractions(path,plot_el_bb["default"].muon) + f=ROOT.TFile.Open(savePath+"NewCollection.root","RECREATE") + #f=ROOT.TFile.Open(savePath+"vfitCollection.root","RECREATE") + #f_jet=ROOT.TFile.Open(".root","RECREATE") + yieldFile=open(savePath+"processYield.txt",'w') + #for year in [2016,2017,2018]: + for year in [2016,2017,2018]: + print year + if year==2016: + processes_other_el=Process(getattr(Backgrounds2016,"OtherEle"),eventCounts_el,negWeights_el) + processes_other_mu=Process(getattr(Backgrounds2016,"Other"),eventCounts_mu,negWeights_mu) + processes_dy_el=Process(getattr(Backgrounds2016,"DrellYan"),eventCounts_el,negWeights_el) + processes_dy_mu=Process(getattr(Backgrounds2016,"DrellYan"),eventCounts_mu,negWeights_mu) + data = Process(Data2016, normalized=True) + data_el = Process(Data2016el, normalized=True) + i=0 + Znorm=zScale2016 + elif year==2018: + processes_other_el=Process(getattr(Backgrounds2018,"Other"),eventCounts_el,negWeights_el) + processes_other_mu=Process(getattr(Backgrounds2018,"Other"),eventCounts_mu,negWeights_mu) + processes_dy_el=Process(getattr(Backgrounds2018,"DrellYan"),eventCounts_el,negWeights_el) + processes_dy_mu=Process(getattr(Backgrounds2018,"DrellYan"),eventCounts_mu,negWeights_mu) + data = Process(Data2018, normalized=True) + i=2 + Znorm=zScale2018 + else: + processes_other_el=Process(getattr(Backgrounds,"Other"),eventCounts_el,negWeights_el) + processes_other_mu=Process(getattr(Backgrounds,"Other"),eventCounts_mu,negWeights_mu) + processes_dy_el=Process(getattr(Backgrounds,"DrellYan"),eventCounts_el,negWeights_el) + processes_dy_mu=Process(getattr(Backgrounds,"DrellYan"),eventCounts_mu,negWeights_mu) + data = Process(Data, normalized=True) + i=1 + Znorm=zScale + + + + hist=data.loadHistogram(plot_mu_bb["default"],lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_bb_data_obs"%str(year)) + Ls=[] + for j in range(hist.GetNbinsX()+1): + Ls.append(hist.GetBinContent(j)) + print(Ls) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=data.loadHistogram(plot_mu_be["default"],lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_be_data_obs"%str(year)) + Ls=[] + for j in range(hist.GetNbinsX()+1): + Ls.append(hist.GetBinContent(j)) + print(Ls) + + #for j in range(hist.GetNbinsX()+1): + # print(hist.GetBinContent(j)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if year==2016: + hist=data_el.loadHistogram(plot_el_bb["default"],lumi_mu[i],Znorm["muons"]) + else: + hist=data.loadHistogram(plot_el_bb["default"],lumi_mu[i],Znorm["muons"]) + hist.SetName("electron_%s_bb_data_obs"%str(year)) + Ls=[] + for j in range(hist.GetNbinsX()+1): + Ls.append(hist.GetBinContent(j)) + print(Ls) + + #for j in range(hist.GetNbinsX()+1): + # print(hist.GetBinContent(j)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if year==2016: + hist=data_el.loadHistogram(plot_el_be["default"],lumi_mu[i],Znorm["muons"]) + else: + hist=data.loadHistogram(plot_el_be["default"],lumi_mu[i],Znorm["muons"]) + hist.SetName("electron_%s_be_data_obs"%str(year)) + Ls=[] + for j in range(hist.GetNbinsX()+1): + Ls.append(hist.GetBinContent(j)) + print(Ls) + + #for j in range(hist.GetNbinsX()+1): + # print(hist.GetBinContent(j)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_dy_mu.loadHistogram(plotGenbb_mu,lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_bb_genMass_DrellYan"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_dy_mu.loadHistogram(plotGenbe_mu,lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_be_genMass_DrellYan"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_other_mu.loadHistogram(plotGenbb_mu,lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_bb_genMass_Other"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_other_mu.loadHistogram(plotGenbe_mu,lumi_mu[i],Znorm["muons"]) + hist.SetName("muon_%s_be_genMass_Other"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_dy_el.loadHistogram(plotGenbb_el,lumi_el[i],Znorm["electrons"][0]) + hist.SetName("electron_%s_bb_genMass_DrellYan"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_dy_el.loadHistogram(plotGenbe_el,lumi_el[i],Znorm["electrons"][1]) + hist.SetName("electron_%s_be_genMass_DrellYan"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + hist=processes_other_el.loadHistogram(plotGenbb_el,lumi_el[i],Znorm["electrons"][0]) + hist.SetName("electron_%s_bb_genMass_Other"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + + hist=processes_other_el.loadHistogram(plotGenbe_el,lumi_el[i],Znorm["electrons"][1]) + hist.SetName("electron_%s_be_genMass_Other"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + + hist=processes_other_el.loadHistogram(plotGenee_el,lumi_el[i],Znorm["electrons"][2]) + hist.SetName("electron_%s_ee_genMass_Other"%str(year)) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + Yield=hist.Integral(bins[0],bins[-1]) + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + norm_el_bb=processes_dy_el.loadHistogram(plot_el_bb["default"],lumi_el[i],Znorm["electrons"][0]).Integral() + norm_el_be=processes_dy_el.loadHistogram(plot_el_be["default"],lumi_el[i],Znorm["electrons"][1]).Integral() + norm_mu_bb=processes_dy_mu.loadHistogram(plot_mu_bb["default"],lumi_mu[i],Znorm["muons"]).Integral() + norm_mu_be=processes_dy_mu.loadHistogram(plot_mu_be["default"],lumi_mu[i],Znorm["muons"]).Integral() + print((norm_el_bb,norm_el_be,norm_mu_bb,norm_mu_be)) + for key in eleplots_bb.keys(): + hist=processes_dy_el.loadHistogram(plot_el_bb[key],lumi_el[i],Znorm["electrons"][0]) + if key == "default": + hist.SetName("electron_%s_bb_recoMass_DrellYan1"%str(year)) + #norm_el_bb=hist.Integral() + else: + tag = getTag("el", year, key, "bb") + hist.SetName(("electron_%s_bb_recoMass_DrellYan1_"%str(year))+tag+key) + + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + hist=processes_dy_el.loadHistogram(plot_el_be[key],lumi_el[i],Znorm["electrons"][1]) + if key == "default": + hist.SetName("electron_%s_be_recoMass_DrellYan1"%str(year)) + #norm_el_be=hist.Integral() + else: + tag = getTag("el", year, key, "be") + hist.SetName(("electron_%s_be_recoMass_DrellYan1_"%str(year))+tag+key) + + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + hist=processes_other_el.loadHistogram(plot_el_bb[key],lumi_el[i],Znorm["electrons"][0]) + if key == "default": + hist.SetName("electron_%s_bb_recoMass_Other"%str(year)) + ee_bb_list=[] + print("ee bb mc") + for j in range(hist.GetNbinsX()+1): + ee_bb_list.append(hist.GetBinContent(j)) + print(ee_bb_list) + else: + tag = getTag("el", year, key, "bb") + hist.SetName(("electron_%s_bb_recoMass_Other_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + hist=processes_other_el.loadHistogram(plot_el_be[key],lumi_el[i],Znorm["electrons"][1]) + if key == "default": + hist.SetName("electron_%s_be_recoMass_Other"%str(year)) + ee_be_list=[] + print("ee be mc") + for j in range(hist.GetNbinsX()+1): + ee_be_list.append(hist.GetBinContent(j)) + print(ee_be_list) + + else: + tag = getTag("el", year, key, "be") + hist.SetName(("electron_%s_be_recoMass_Other_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + + + for key in muplots_bb.keys(): + hist=processes_dy_mu.loadHistogram(plot_mu_bb[key],lumi_mu[i],Znorm["muons"]) + if key == "default": + hist.SetName("muon_%s_bb_recoMass_DrellYan"%str(year)) + #norm_mu_bb=hist.Integral() + else: + tag = getTag("mu", year, key, "bb") + hist.SetName(("muon_%s_bb_recoMass_DrellYan_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if key=="ID" or key=="reso": + tag = getTag("mu", year, key, "bb") + hist.SetName(("muon_%s_bb_recoMass_DrellYan_"%str(year))+tag+key+"Up") + default=processes_dy_mu.loadHistogram(plot_mu_bb["default"],lumi_mu[i],Znorm["muons"]) + default.Scale(2) + default.Add(hist,-1) + default.SetName(("muon_%s_bb_recoMass_DrellYan_"%str(year))+tag+key+"Down") + hist.Scale(norm_el_bb/norm_mu_bb) + default.Scale(norm_el_bb/norm_mu_bb) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + if default.GetBinContent(j)<0: + default.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + f.WriteObject(default,default.GetName()) + else: + print(hist.Integral()) + hist.Scale(norm_el_bb/norm_mu_bb) + print(hist.Integral()) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + + hist=processes_dy_mu.loadHistogram(plot_mu_be[key],lumi_mu[i],Znorm["muons"]) + if key == "default": + hist.SetName("muon_%s_be_recoMass_DrellYan"%str(year)) + #norm_mu_be=hist.Integral() + else: + tag = getTag("mu", year, key, "be") + hist.SetName(("muon_%s_be_recoMass_DrellYan_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if key=="ID" or key=="reso": + tag = getTag("mu", year, key, "be") + hist.SetName(("muon_%s_be_recoMass_DrellYan_"%str(year))+tag+key+"Up") + default=processes_dy_mu.loadHistogram(plot_mu_be["default"],lumi_mu[i],Znorm["muons"]) + default.Scale(2) + default.Add(hist,-1) + default.SetName(("muon_%s_be_recoMass_DrellYan_"%str(year))+tag+key+"Down") + hist.Scale(norm_el_be/norm_mu_be) + default.Scale(norm_el_be/norm_mu_be) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + if default.GetBinContent(j)<0: + default.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + f.WriteObject(default,default.GetName()) + else: + hist.Scale(norm_el_be/norm_mu_be) + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + + hist=processes_other_mu.loadHistogram(plot_mu_bb[key],lumi_mu[i],Znorm["muons"]) + + if key == "default": + hist.SetName("muon_%s_bb_recoMass_Other"%str(year)) + print("mu bb mc") + mu_bb_list=[] + for j in range(hist.GetNbinsX()+1): + mu_bb_list.append(hist.GetBinContent(j)) + print(mu_bb_list) + else: + tag = getTag("mu", year, key, "bb") + hist.SetName(("muon_%s_bb_recoMass_Other_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if key=="ID" or key=="reso": + tag = getTag("mu", year, key, "bb") + hist.SetName(("muon_%s_bb_recoMass_Other_"%str(year))+tag+key+"Up") + default=processes_other_mu.loadHistogram(plot_mu_bb["default"],lumi_mu[i],Znorm["muons"]) + default.Scale(2) + default.Add(hist,-1) + default.SetName(("muon_%s_bb_recoMass_Other_"%str(year))+tag+key+"Down") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + if default.GetBinContent(j)<0: + default.SetBinContent(j,0) + + f.WriteObject(hist,hist.GetName()) + f.WriteObject(default,default.GetName()) + else: + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + + hist=processes_other_mu.loadHistogram(plot_mu_be[key],lumi_mu[i],Znorm["muons"]) + if key == "default": + hist.SetName("muon_%s_be_recoMass_Other"%str(year)) + print("mu be mc") + mu_be_list=[] + for j in range(hist.GetNbinsX()+1): + mu_be_list.append(hist.GetBinContent(j)) + print(mu_be_list) + else: + tag = getTag("mu", year, key, "be") + hist.SetName(("muon_%s_be_recoMass_Other_"%str(year))+tag+key) + Yield=hist.Integral() + yieldFile.write(hist.GetName()+": "+str(Yield)+"\n") + if key=="ID" or key=="reso": + tag = getTag("mu", year, key, "be") + hist.SetName(("muon_%s_be_recoMass_Other_"%str(year))+tag+key+"Up") + default=processes_other_mu.loadHistogram(plot_mu_be["default"],lumi_mu[i],Znorm["muons"]) + default.Scale(2) + default.Add(hist,-1) + default.SetName(("muon_%s_be_recoMass_Other_"%str(year))+tag+key+"Down") + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + if default.GetBinContent(j)<0: + default.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + f.WriteObject(default,default.GetName()) + else: + for j in range(hist.GetNbinsX()+1): + if hist.GetBinContent(j)<0: + hist.SetBinContent(j,0) + f.WriteObject(hist,hist.GetName()) + f.Save() + #print(processes2016_other_el_bb) + +