RecoTau/python/compare.py

0001 from __future__ import print_function
0002 from __future__ import absolute_import
0003 from builtins import range
0004 from .officialStyle import officialStyle
0005 from array import array
0006 from ROOT import gROOT, gStyle, TH1F, TH1D, TF1, TFile, TCanvas, TH2F, TLegend, TGraphAsymmErrors, Double, TLatex
0007 import os, copy, sys
0008
0009 gROOT.SetBatch(True)
0010 officialStyle(gStyle)
0011 gStyle.SetOptTitle(0)
0012 #set_palette("color")
0013 #gStyle.SetPaintTextFormat("2.0f")
0014
0015
0016 argvs = sys.argv
0017 argc = len(argvs)
0018
0019 if argc != 2:
0020     print('Please specify the runtype : python3 tauPOGplot.py <ZTT, ZEE, ZMM, QCD>')
0021     sys.exit(0)
0022
0023 runtype = argvs[1]
0024 print('You selected', runtype)
0025
0026
0027 tlabel = 'Z #rightarrow #tau#tau'
0028 xlabel = 'gen. tau p_{T}^{vis} (GeV)'
0029
0030 if runtype == 'QCD':
0031     tlabel = 'QCD'
0032     xlabel = 'jet p_{T} (GeV)'
0033 elif runtype == 'ZEE':
0034     tlabel = 'Z #rightarrow ee'
0035     xlabel = 'electron p_{T} (GeV)'
0036 elif runtype == 'ZMM':
0037     tlabel = 'Z #rightarrow #mu#mu'
0038     xlabel = 'muon p_{T} (GeV)'
0039
0040
0041
0042
0043 def ensureDir(directory):
0044     if not os.path.exists(directory):
0045         os.makedirs(directory)
0046
0047 def save(canvas, name):
0048     ensureDir('compare_' + runtype)
0049     canvas.SaveAs(name.replace(' ','').replace('&&','')+'.pdf')
0050     canvas.SaveAs(name.replace(' ','').replace('&&','')+'.gif')
0051
0052
0053 def LegendSettings(leg, ncolumn):
0054     leg.SetNColumns(ncolumn)
0055     leg.SetBorderSize(0)
0056     leg.SetFillColor(10)
0057     leg.SetLineColor(0)
0058     leg.SetFillStyle(0)
0059     leg.SetTextSize(0.03)
0060     leg.SetTextFont(42)
0061
0062
0063
0064 def makeCompareVars(tree, var, sel, leglist, nbin, xmin, xmax, xtitle, ytitle, scale, name):
0065
0066 #    print leglist
0067
0068     c = TCanvas()
0069
0070     hists = []
0071     col = [1,2,4,8,6]
0072     ymax = 0
0073
0074     for ii, isel in enumerate(sel):
0075         hist = TH1F('h_' + str(ii), 'h_' + str(ii), nbin, xmin, xmax)
0076         hist.GetXaxis().SetTitle(xtitle)
0077         hist.GetYaxis().SetTitle(ytitle)
0078         hist.GetYaxis().SetNdivisions(507)
0079         hist.SetLineColor(col[ii])
0080         hist.SetLineWidth(len(sel)+1-ii)
0081         hist.SetLineStyle(1)
0082         hist.SetMarkerSize(0)
0083         hist.SetMinimum(0)
0084         hist.Sumw2()
0085
0086 #        print hist.GetName(), var, isel
0087         tree.Project(hist.GetName(), var, isel)
0088         hist.Scale(1./hist.GetEntries())
0089
0090         if ymax < hist.GetMaximum():
0091             ymax = hist.GetMaximum()
0092
0093         hists.append(hist)
0094
0095
0096     leg = TLegend(0.6,0.65,0.91,0.9)
0097     LegendSettings(leg,1)
0098
0099     for ii, ihist in enumerate(hists):
0100         ihist.SetMaximum(ymax*1.2)
0101         ihist.SetMinimum(0.)
0102
0103         if ii ==0:
0104             ihist.Draw('h')
0105         else:
0106             ihist.Draw('hsame')
0107
0108         if leglist[ii] != 'None':
0109             leg.AddEntry(ihist, leglist[ii], "l")
0110
0111
0112     if leglist[0]!='None':
0113         leg.Draw()
0114
0115
0116 #    save(c, 'compare_' + runtype + '/compare_' + name)
0117
0118
0119
0120
0121
0122 def overlay(hists, ytitle, header, addon):
0123
0124     print('number of histograms = ', len(hists))
0125
0126     canvas = TCanvas()
0127     leg = TLegend(0.2,0.7,0.5,0.9)
0128     LegendSettings(leg, 1)
0129
0130     col = [1,2,4,6,8,9,12]
0131
0132     ymax = -1
0133     ymin = 100
0134
0135     for ii, hist in enumerate(hists):
0136         hist.GetYaxis().SetTitle('efficiency')
0137         hist.SetLineColor(col[ii])
0138         hist.SetMarkerColor(col[ii])
0139         hist.SetLineWidth(2)
0140         hist.SetMarkerSize(1)
0141
0142
0143         for ip in range(hist.GetN()):
0144             x = Double(-1)
0145             y = Double(-1)
0146             hist.GetPoint(ip, x, y)
0147
0148             if ymin > y:
0149                 ymin = y
0150             if ymax < y:
0151                 ymax = y
0152
0153 #
0154 #        if ymax < hist.GetMaximum():
0155 #            ymax = hist.GetMaximum()
0156 #        if ymin > hist.GetMinimum():
0157 #            ymin = hist.GetMinimum()
0158
0159         if ii==0:
0160             hist.Draw("Azp")
0161         else:
0162
0163             hist.Draw("pzsame")
0164
0165 #        print hist.GetName(), hist.GetTitle()
0166         legname = hist.GetName()
0167
0168         leg.AddEntry(hist, legname, 'lep')
0169
0170
0171     for hist in hists:
0172         hist.SetMaximum(ymax*2)
0173 #        hist.SetMinimum(ymin*0.5)
0174
0175     leg.Draw()
0176
0177     tex = TLatex( hists[-1].GetXaxis().GetXmin() + 0.01*(hists[-1].GetXaxis().GetXmax() - hists[-1].GetXaxis().GetXmin()), ymax*2.1, addon.replace('tau_',''))
0178
0179     tex.SetTextFont(42)
0180     tex.SetTextSize(0.03)
0181     tex.Draw()
0182
0183     tex2 = TLatex( hists[-1].GetXaxis().GetXmin() + 0.87*(hists[-1].GetXaxis().GetXmax() - hists[-1].GetXaxis().GetXmin()), ymax*2.1, tlabel)
0184
0185     tex2.SetTextFont(42)
0186     tex2.SetTextSize(0.03)
0187     tex2.Draw()
0188
0189
0190
0191     save(canvas, 'compare_' + runtype + '/' + header)
0192
0193
0194
0195 def hoverlay(hists, xtitle, ytitle, name):
0196
0197     c = TCanvas()
0198
0199     ymax = 0
0200     for hist in hists:
0201         if ymax < hist.GetMaximum():
0202             ymax = hist.GetMaximum()
0203
0204
0205     leg = TLegend(0.6,0.65,0.91,0.9)
0206     LegendSettings(leg,1)
0207
0208     for ii, ihist in enumerate(hists):
0209         ihist.SetMaximum(ymax*1.2)
0210         ihist.SetMinimum(0.)
0211         ihist.SetMarkerSize(0.)
0212         ihist.GetXaxis().SetTitle(xtitle)
0213         ihist.GetYaxis().SetTitle(ytitle)
0214
0215         if ii ==0:
0216             ihist.Draw('h')
0217         else:
0218             ihist.Draw('hsame')
0219
0220         leg.AddEntry(ihist, ihist.GetName(), "l")
0221
0222
0223     leg.Draw()
0224
0225
0226     save(c, 'compare_' + runtype + '/hist_' + name)
0227
0228
0229
0230
0231 def makeEffPlotsVars(tree, varx, vary, sel, nbinx, xmin, xmax, nbiny, ymin, ymax, xtitle, ytitle, leglabel = None, header='', addon='', option='pt', marker=20):
0232
0233     binning = [20,30,40,50,60,70,80,100,150,200]
0234
0235     c = TCanvas()
0236
0237     if option=='pt':
0238         _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, len(binning)-1, array('d',binning))
0239         _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, len(binning)-1, array('d',binning))
0240     elif option=='eta':
0241         _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, nbinx, xmin, xmax)
0242         _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, nbinx, xmin, xmax)
0243     elif option=='nvtx':
0244         _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, len(vbinning)-1, array('d',vbinning))
0245         _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, len(vbinning)-1, array('d',vbinning))
0246
0247
0248     tree.Draw(varx + ' >> ' + _hist_.GetName(), sel)
0249     tree.Draw(varx + ' >> ' + _ahist_.GetName(), sel + ' && ' + vary)
0250
0251     g_efficiency = TGraphAsymmErrors()
0252     g_efficiency.BayesDivide(_ahist_, _hist_)
0253     g_efficiency.GetXaxis().SetTitle(xtitle)
0254     g_efficiency.GetYaxis().SetTitle('efficiency')
0255     g_efficiency.GetYaxis().SetNdivisions(507)
0256     g_efficiency.SetLineWidth(3)
0257     g_efficiency.SetName(header)
0258     g_efficiency.SetMinimum(0.)
0259     g_efficiency.GetYaxis().SetTitleOffset(1.3)
0260     g_efficiency.SetMarkerStyle(marker)
0261     g_efficiency.SetMarkerSize(1)
0262     g_efficiency.Draw('ap')
0263
0264 #    save(c, 'plots/' + addon)
0265
0266     return copy.deepcopy(g_efficiency)
0267
0268
0269
0270
0271
0272 if __name__ == '__main__':
0273
0274
0275     vardict = {
0276         'againstMuonLoose3':{'var':'tau_againstMuonLoose3 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstMuonLoose3'},
0277         'againstMuonTight3':{'var':'tau_againstMuonTight3 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstMuonTight3'},
0278
0279         'againstElectronVLooseMVA5':{'var':'tau_againstElectronVLooseMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronVLooseMVA5'},
0280         'againstElectronLooseMVA5':{'var':'tau_againstElectronLooseMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronLooseMVA5'},
0281         'againstElectronMediumMVA5':{'var':'tau_againstElectronMediumMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronMediumMVA5'},
0282
0283         'againstElectronVLooseMVA6':{'var':'tau_againstElectronVLooseMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronVLooseMVA6'},
0284         'againstElectronLooseMVA6':{'var':'tau_againstElectronLooseMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronLooseMVA6'},
0285         'againstElectronMediumMVA6':{'var':'tau_againstElectronMediumMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronMediumMVA6'},
0286
0287         'byLoosePileupWeightedIsolation3Hits':{'var':'tau_byLoosePileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLoosePileupWeightedIsolation3Hits'},
0288         'byMediumPileupWeightedIsolation3Hits':{'var':'tau_byMediumPileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumPileupWeightedIsolation3Hits'},
0289         'byTightPileupWeightedIsolation3Hits':{'var':'tau_byTightPileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightPileupWeightedIsolation3Hits'},
0290
0291         'byLooseCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byLooseCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseCombinedIsolationDeltaBetaCorr3Hits'},
0292         'byMediumCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byMediumCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumCombinedIsolationDeltaBetaCorr3Hits'},
0293         'byTightCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byTightCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightCombinedIsolationDeltaBetaCorr3Hits'},
0294
0295         'byLooseIsolationMVA3oldDMwLT':{'var':'tau_byLooseIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseIsolationMVA3oldDMwLT'},
0296         'byMediumIsolationMVA3oldDMwLT':{'var':'tau_byMediumIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumIsolationMVA3oldDMwLT'},
0297         'byTightIsolationMVA3oldDMwLT':{'var':'tau_byTightIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightIsolationMVA3oldDMwLT'},
0298
0299         'byLooseIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byLooseIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseIsolationMVArun2v1DBoldDMwLT'},
0300         'byMediumIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byMediumIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumIsolationMVArun2v1DBoldDMwLT'},
0301         'byTightIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byTightIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightIsolationMVArun2v1DBoldDMwLT'},
0302
0303         }
0304
0305
0306     reco_cut = 'tau_pt > 20 && abs(tau_eta) < 2.3'
0307     loose_id = 'tau_decayModeFindingOldDMs > 0.5 && tau_byLooseCombinedIsolationDeltaBetaCorr3Hits > 0.5'
0308
0309     sampledict = {
0310 #        '7_6_0_pre7':{'file':'Myroot_7_6_0_pre7_' + runtype + '.root', 'col':2, 'marker':20, 'width':3},
0311 #        '7_6_0':{'file':'Myroot_7_6_0_' + runtype + '.root', 'col':1, 'marker':21, 'width':4},
0312 #        '7_6_1':{'file':'Myroot_7_6_1_' + runtype + '.root', 'col':4, 'marker':22, 'width':2},
0313         '7_6_1_v3':{'file':'Myroot_7_6_1_v3_' + runtype + '.root', 'col':3, 'marker':23, 'width':1},
0314         }
0315
0316     for hname, hdict in sorted(vardict.items()):
0317
0318         hists = []
0319
0320         for rel, rdict in sorted(sampledict.items()):
0321
0322             if rel.find('7_6_1')==-1 and (hname.find('MVA6')!=-1 or hname.find('MVArun2')!=-1): continue
0323
0324             tfile = TFile(rdict['file'])
0325             tree = tfile.Get('per_tau')
0326
0327             num_sel = reco_cut
0328             den_sel = '1'
0329
0330             if hname.find('against')!=-1:
0331                 num_sel = '1'
0332                 den_sel = reco_cut + ' && ' + loose_id
0333
0334             hists.append(makeEffPlotsVars(tree, 'tau_genpt', num_sel + '&&' + hdict['var'], den_sel, 30, 0, 300, hdict['nbin'], hdict['min'], hdict['max'], xlabel, hdict['title'], '', rel, rel, 'pt', rdict['marker']))
0335
0336
0337 #            if rel=='7_6_1' and (hname.find('MVA5')!=-1 or hname.find('IsolationMVA3')!=-1):
0338 #                xvar = hdict['var'].replace('IsolationMVA3', 'IsolationMVArun2v1DB').replace('MVA5','MVA6')
0339 #                print 'adding', xvar
0340 #
0341 #                hists.append(makeEffPlotsVars(tree, 'tau_genpt', num_sel + '&&' + xvar, den_sel, 30, 0, 300, hdict['nbin'], hdict['min'], hdict['max'], xlabel, hdict['title'], '', rel + '(' + xvar.replace('tau_','').replace('> 0.5','').replace(' && decayModeFindingOldDMs ','') + ')', rel, 'pt', rdict['marker']))
0342
0343
0344
0345         overlay(hists, hname, hname, hdict['title'])
0346
0347
0348
0349
0350     hvardict = {
0351         'tau_dm':{'var':'tau_dm', 'nbin':12, 'min':0., 'max':12, 'title':'decay Mode', 'sel':'1'},
0352         'tau_mass_1prong':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 1prong', 'sel':'tau_dm==0'},
0353         'tau_mass_1prongp0':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 1prong+#pi^{0}', 'sel':'tau_dm==1'},
0354         'tau_mass_2prong':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 2prong', 'sel':'(tau_dm==5 || tau_dm==6)'},
0355         'tau_mass_3prong':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 3prong (+#pi^{0})', 'sel':'(tau_dm==10 || tau_dm==11)'},
0356         'pt_resolution_1prong':{'var':'(tau_genpt-tau_pt)/(tau_genpt)', 'nbin':30, 'min':-1., 'max':1., 'title':'pT resolution, 1prong', 'sel':'tau_dm==0'},
0357         'pt_resolution_1prongp0':{'var':'(tau_genpt-tau_pt)/(tau_genpt)', 'nbin':30, 'min':-1., 'max':1., 'title':'pT resolution, 1prong+#pi^{0}', 'sel':'tau_dm==1'},
0358         'pt_resolution_2prong':{'var':'(tau_genpt-tau_pt)/(tau_genpt)', 'nbin':30, 'min':-1., 'max':1., 'title':'pT resolution, 2prong', 'sel':'(tau_dm==5 || tau_dm==6)'},
0359         'pt_resolution_3prong':{'var':'(tau_genpt-tau_pt)/(tau_genpt)', 'nbin':30, 'min':-1., 'max':1., 'title':'pT resolution, 3prong (+#pi^{0})', 'sel':'(tau_dm==10 || tau_dm==11)'},
0360         }
0361
0362
0363     for hname, hdict in sorted(hvardict.items()):
0364
0365         hists = []
0366
0367         if runtype != 'ZTT' and hname.find('pt_resolution')!=-1: continue
0368
0369
0370         for rel, rdict in sorted(sampledict.items()):
0371
0372             tfile = TFile(rdict['file'])
0373             tree = tfile.Get('per_tau')
0374
0375             hist = TH1F('h_' + hname + '_' + rel, 'h_' + hname + '_' + rel, hdict['nbin'], hdict['min'], hdict['max'])
0376             hist.GetYaxis().SetNdivisions(507)
0377             hist.SetLineColor(rdict['col'])
0378             hist.SetLineWidth(rdict['width'])
0379             hist.SetMinimum(0)
0380             hist.SetName(rel)
0381             hist.Sumw2()
0382
0383             tree.Project(hist.GetName(), hdict['var'], hdict['sel'])
0384             hist.Scale(1./hist.GetEntries())
0385
0386             hists.append(copy.deepcopy(hist))
0387
0388         hoverlay(hists, hdict['title'], 'a.u.', hname)