RecoTau/python/compare.py

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