CMS 3D CMS Logo

compare.py
Go to the documentation of this file.
1 from officialStyle import officialStyle
2 from array import array
3 from ROOT import gROOT, gStyle, TH1F, TH1D, TF1, TFile, TCanvas, TH2F, TLegend, TGraphAsymmErrors, Double, TLatex
4 import os, copy, sys
5 import six
6 
7 gROOT.SetBatch(True)
8 officialStyle(gStyle)
9 gStyle.SetOptTitle(0)
10 #set_palette("color")
11 #gStyle.SetPaintTextFormat("2.0f")
12 
13 
14 argvs = sys.argv
15 argc = len(argvs)
16 
17 if argc != 2:
18  print 'Please specify the runtype : python tauPOGplot.py <ZTT, ZEE, ZMM, QCD>'
19  sys.exit(0)
20 
21 runtype = argvs[1]
22 print 'You selected', runtype
23 
24 
25 tlabel = 'Z #rightarrow #tau#tau'
26 xlabel = 'gen. tau p_{T}^{vis} (GeV)'
27 
28 if runtype == 'QCD':
29  tlabel = 'QCD'
30  xlabel = 'jet p_{T} (GeV)'
31 elif runtype == 'ZEE':
32  tlabel = 'Z #rightarrow ee'
33  xlabel = 'electron p_{T} (GeV)'
34 elif runtype == 'ZMM':
35  tlabel = 'Z #rightarrow #mu#mu'
36  xlabel = 'muon p_{T} (GeV)'
37 
38 
39 
40 
41 def ensureDir(directory):
42  if not os.path.exists(directory):
43  os.makedirs(directory)
44 
45 def save(canvas, name):
46  ensureDir('compare_' + runtype)
47  canvas.SaveAs(name.replace(' ','').replace('&&','')+'.pdf')
48  canvas.SaveAs(name.replace(' ','').replace('&&','')+'.gif')
49 
50 
51 def LegendSettings(leg, ncolumn):
52  leg.SetNColumns(ncolumn)
53  leg.SetBorderSize(0)
54  leg.SetFillColor(10)
55  leg.SetLineColor(0)
56  leg.SetFillStyle(0)
57  leg.SetTextSize(0.03)
58  leg.SetTextFont(42)
59 
60 
61 
62 def makeCompareVars(tree, var, sel, leglist, nbin, xmin, xmax, xtitle, ytitle, scale, name):
63 
64 # print leglist
65 
66  c = TCanvas()
67 
68  hists = []
69  col = [1,2,4,8,6]
70  ymax = 0
71 
72  for ii, isel in enumerate(sel):
73  hist = TH1F('h_' + str(ii), 'h_' + str(ii), nbin, xmin, xmax)
74  hist.GetXaxis().SetTitle(xtitle)
75  hist.GetYaxis().SetTitle(ytitle)
76  hist.GetYaxis().SetNdivisions(507)
77  hist.SetLineColor(col[ii])
78  hist.SetLineWidth(len(sel)+1-ii)
79  hist.SetLineStyle(1)
80  hist.SetMarkerSize(0)
81  hist.SetMinimum(0)
82  hist.Sumw2()
83 
84 # print hist.GetName(), var, isel
85  tree.Project(hist.GetName(), var, isel)
86  hist.Scale(1./hist.GetEntries())
87 
88  if ymax < hist.GetMaximum():
89  ymax = hist.GetMaximum()
90 
91  hists.append(hist)
92 
93 
94  leg = TLegend(0.6,0.65,0.91,0.9)
95  LegendSettings(leg,1)
96 
97  for ii, ihist in enumerate(hists):
98  ihist.SetMaximum(ymax*1.2)
99  ihist.SetMinimum(0.)
100 
101  if ii ==0:
102  ihist.Draw('h')
103  else:
104  ihist.Draw('hsame')
105 
106  if leglist[ii] != 'None':
107  leg.AddEntry(ihist, leglist[ii], "l")
108 
109 
110  if leglist[0]!='None':
111  leg.Draw()
112 
113 
114 # save(c, 'compare_' + runtype + '/compare_' + name)
115 
116 
117 
118 
119 
120 def overlay(hists, ytitle, header, addon):
121 
122  print 'number of histograms = ', len(hists)
123 
124  canvas = TCanvas()
125  leg = TLegend(0.2,0.7,0.5,0.9)
126  LegendSettings(leg, 1)
127 
128  col = [1,2,4,6,8,9,12]
129 
130  ymax = -1
131  ymin = 100
132 
133  for ii, hist in enumerate(hists):
134  hist.GetYaxis().SetTitle('efficiency')
135  hist.SetLineColor(col[ii])
136  hist.SetMarkerColor(col[ii])
137  hist.SetLineWidth(2)
138  hist.SetMarkerSize(1)
139 
140 
141  for ip in range(hist.GetN()):
142  x = Double(-1)
143  y = Double(-1)
144  hist.GetPoint(ip, x, y)
145 
146  if ymin > y:
147  ymin = y
148  if ymax < y:
149  ymax = y
150 
151 #
152 # if ymax < hist.GetMaximum():
153 # ymax = hist.GetMaximum()
154 # if ymin > hist.GetMinimum():
155 # ymin = hist.GetMinimum()
156 
157  if ii==0:
158  hist.Draw("Azp")
159  else:
160 
161  hist.Draw("pzsame")
162 
163 # print hist.GetName(), hist.GetTitle()
164  legname = hist.GetName()
165 
166  leg.AddEntry(hist, legname, 'lep')
167 
168 
169  for hist in hists:
170  hist.SetMaximum(ymax*2)
171 # hist.SetMinimum(ymin*0.5)
172 
173  leg.Draw()
174 
175  tex = TLatex( hists[-1].GetXaxis().GetXmin() + 0.01*(hists[-1].GetXaxis().GetXmax() - hists[-1].GetXaxis().GetXmin()), ymax*2.1, addon.replace('tau_',''))
176 
177  tex.SetTextFont(42)
178  tex.SetTextSize(0.03)
179  tex.Draw()
180 
181  tex2 = TLatex( hists[-1].GetXaxis().GetXmin() + 0.87*(hists[-1].GetXaxis().GetXmax() - hists[-1].GetXaxis().GetXmin()), ymax*2.1, tlabel)
182 
183  tex2.SetTextFont(42)
184  tex2.SetTextSize(0.03)
185  tex2.Draw()
186 
187 
188 
189  save(canvas, 'compare_' + runtype + '/' + header)
190 
191 
192 
193 def hoverlay(hists, xtitle, ytitle, name):
194 
195  c = TCanvas()
196 
197  ymax = 0
198  for hist in hists:
199  if ymax < hist.GetMaximum():
200  ymax = hist.GetMaximum()
201 
202 
203  leg = TLegend(0.6,0.65,0.91,0.9)
204  LegendSettings(leg,1)
205 
206  for ii, ihist in enumerate(hists):
207  ihist.SetMaximum(ymax*1.2)
208  ihist.SetMinimum(0.)
209  ihist.SetMarkerSize(0.)
210  ihist.GetXaxis().SetTitle(xtitle)
211  ihist.GetYaxis().SetTitle(ytitle)
212 
213  if ii ==0:
214  ihist.Draw('h')
215  else:
216  ihist.Draw('hsame')
217 
218  leg.AddEntry(ihist, ihist.GetName(), "l")
219 
220 
221  leg.Draw()
222 
223 
224  save(c, 'compare_' + runtype + '/hist_' + name)
225 
226 
227 
228 
229 def makeEffPlotsVars(tree, varx, vary, sel, nbinx, xmin, xmax, nbiny, ymin, ymax, xtitle, ytitle, leglabel = None, header='', addon='', option='pt', marker=20):
230 
231  binning = [20,30,40,50,60,70,80,100,150,200]
232 
233  c = TCanvas()
234 
235  if option=='pt':
236  _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, len(binning)-1, array('d',binning))
237  _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, len(binning)-1, array('d',binning))
238  elif option=='eta':
239  _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, nbinx, xmin, xmax)
240  _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, nbinx, xmin, xmax)
241  elif option=='nvtx':
242  _hist_ = TH1F('h_effp_' + addon, 'h_effp' + addon, len(vbinning)-1, array('d',vbinning))
243  _ahist_ = TH1F('ah_effp_' + addon, 'ah_effp' + addon, len(vbinning)-1, array('d',vbinning))
244 
245 
246  tree.Draw(varx + ' >> ' + _hist_.GetName(), sel)
247  tree.Draw(varx + ' >> ' + _ahist_.GetName(), sel + ' && ' + vary)
248 
249  g_efficiency = TGraphAsymmErrors()
250  g_efficiency.BayesDivide(_ahist_, _hist_)
251  g_efficiency.GetXaxis().SetTitle(xtitle)
252  g_efficiency.GetYaxis().SetTitle('efficiency')
253  g_efficiency.GetYaxis().SetNdivisions(507)
254  g_efficiency.SetLineWidth(3)
255  g_efficiency.SetName(header)
256  g_efficiency.SetMinimum(0.)
257  g_efficiency.GetYaxis().SetTitleOffset(1.3)
258  g_efficiency.SetMarkerStyle(marker)
259  g_efficiency.SetMarkerSize(1)
260  g_efficiency.Draw('ap')
261 
262 # save(c, 'plots/' + addon)
263 
264  return copy.deepcopy(g_efficiency)
265 
266 
267 
268 
269 
270 if __name__ == '__main__':
271 
272 
273  vardict = {
274  'againstMuonLoose3':{'var':'tau_againstMuonLoose3 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstMuonLoose3'},
275  'againstMuonTight3':{'var':'tau_againstMuonTight3 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstMuonTight3'},
276 
277  'againstElectronVLooseMVA5':{'var':'tau_againstElectronVLooseMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronVLooseMVA5'},
278  'againstElectronLooseMVA5':{'var':'tau_againstElectronLooseMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronLooseMVA5'},
279  'againstElectronMediumMVA5':{'var':'tau_againstElectronMediumMVA5 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronMediumMVA5'},
280 
281  'againstElectronVLooseMVA6':{'var':'tau_againstElectronVLooseMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronVLooseMVA6'},
282  'againstElectronLooseMVA6':{'var':'tau_againstElectronLooseMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronLooseMVA6'},
283  'againstElectronMediumMVA6':{'var':'tau_againstElectronMediumMVA6 > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'againstElectronMediumMVA6'},
284 
285  'byLoosePileupWeightedIsolation3Hits':{'var':'tau_byLoosePileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLoosePileupWeightedIsolation3Hits'},
286  'byMediumPileupWeightedIsolation3Hits':{'var':'tau_byMediumPileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumPileupWeightedIsolation3Hits'},
287  'byTightPileupWeightedIsolation3Hits':{'var':'tau_byTightPileupWeightedIsolation3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightPileupWeightedIsolation3Hits'},
288 
289  'byLooseCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byLooseCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseCombinedIsolationDeltaBetaCorr3Hits'},
290  'byMediumCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byMediumCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumCombinedIsolationDeltaBetaCorr3Hits'},
291  'byTightCombinedIsolationDeltaBetaCorr3Hits':{'var':'tau_byTightCombinedIsolationDeltaBetaCorr3Hits > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightCombinedIsolationDeltaBetaCorr3Hits'},
292 
293  'byLooseIsolationMVA3oldDMwLT':{'var':'tau_byLooseIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseIsolationMVA3oldDMwLT'},
294  'byMediumIsolationMVA3oldDMwLT':{'var':'tau_byMediumIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumIsolationMVA3oldDMwLT'},
295  'byTightIsolationMVA3oldDMwLT':{'var':'tau_byTightIsolationMVA3oldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightIsolationMVA3oldDMwLT'},
296 
297  'byLooseIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byLooseIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byLooseIsolationMVArun2v1DBoldDMwLT'},
298  'byMediumIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byMediumIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byMediumIsolationMVArun2v1DBoldDMwLT'},
299  'byTightIsolationMVArun2v1DBoldDMwLT':{'var':'tau_byTightIsolationMVArun2v1DBoldDMwLT > 0.5 && tau_decayModeFindingOldDMs > 0.5', 'nbin':2, 'min':-0.5, 'max':1.5, 'title':'byTightIsolationMVArun2v1DBoldDMwLT'},
300 
301  }
302 
303 
304  reco_cut = 'tau_pt > 20 && abs(tau_eta) < 2.3'
305  loose_id = 'tau_decayModeFindingOldDMs > 0.5 && tau_byLooseCombinedIsolationDeltaBetaCorr3Hits > 0.5'
306 
307  sampledict = {
308 # '7_6_0_pre7':{'file':'Myroot_7_6_0_pre7_' + runtype + '.root', 'col':2, 'marker':20, 'width':3},
309 # '7_6_0':{'file':'Myroot_7_6_0_' + runtype + '.root', 'col':1, 'marker':21, 'width':4},
310 # '7_6_1':{'file':'Myroot_7_6_1_' + runtype + '.root', 'col':4, 'marker':22, 'width':2},
311  '7_6_1_v3':{'file':'Myroot_7_6_1_v3_' + runtype + '.root', 'col':3, 'marker':23, 'width':1},
312  }
313 
314  for hname, hdict in sorted(six.iteritems(vardict)):
315 
316  hists = []
317 
318  for rel, rdict in sorted(six.iteritems(sampledict)):
319 
320  if rel.find('7_6_1')==-1 and (hname.find('MVA6')!=-1 or hname.find('MVArun2')!=-1): continue
321 
322  tfile = TFile(rdict['file'])
323  tree = tfile.Get('per_tau')
324 
325  num_sel = reco_cut
326  den_sel = '1'
327 
328  if hname.find('against')!=-1:
329  num_sel = '1'
330  den_sel = reco_cut + ' && ' + loose_id
331 
332  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']))
333 
334 
335 # if rel=='7_6_1' and (hname.find('MVA5')!=-1 or hname.find('IsolationMVA3')!=-1):
336 # xvar = hdict['var'].replace('IsolationMVA3', 'IsolationMVArun2v1DB').replace('MVA5','MVA6')
337 # print 'adding', xvar
338 #
339 # 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']))
340 
341 
342 
343  overlay(hists, hname, hname, hdict['title'])
344 
345 
346 
347 
348  hvardict = {
349  'tau_dm':{'var':'tau_dm', 'nbin':12, 'min':0., 'max':12, 'title':'decay Mode', 'sel':'1'},
350  'tau_mass_1prong':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 1prong', 'sel':'tau_dm==0'},
351  'tau_mass_1prongp0':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 1prong+#pi^{0}', 'sel':'tau_dm==1'},
352  'tau_mass_2prong':{'var':'tau_mass', 'nbin':30, 'min':0., 'max':2.5, 'title':'Tau mass, 2prong', 'sel':'(tau_dm==5 || tau_dm==6)'},
353  '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)'},
354  'pt_resolution_1prong':{'var':'(tau_genpt-tau_pt)/(tau_genpt)', 'nbin':30, 'min':-1., 'max':1., 'title':'pT resolution, 1prong', 'sel':'tau_dm==0'},
355  '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'},
356  '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)'},
357  '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)'},
358  }
359 
360 
361  for hname, hdict in sorted(six.iteritems(hvardict)):
362 
363  hists = []
364 
365  if runtype != 'ZTT' and hname.find('pt_resolution')!=-1: continue
366 
367 
368  for rel, rdict in sorted(six.iteritems(sampledict)):
369 
370  tfile = TFile(rdict['file'])
371  tree = tfile.Get('per_tau')
372 
373  hist = TH1F('h_' + hname + '_' + rel, 'h_' + hname + '_' + rel, hdict['nbin'], hdict['min'], hdict['max'])
374  hist.GetYaxis().SetNdivisions(507)
375  hist.SetLineColor(rdict['col'])
376  hist.SetLineWidth(rdict['width'])
377  hist.SetMinimum(0)
378  hist.SetName(rel)
379  hist.Sumw2()
380 
381  tree.Project(hist.GetName(), hdict['var'], hdict['sel'])
382  hist.Scale(1./hist.GetEntries())
383 
384  hists.append(copy.deepcopy(hist))
385 
386  hoverlay(hists, hdict['title'], 'a.u.', hname)
def replace(string, replacements)
def save(canvas, name)
Definition: compare.py:45
def LegendSettings(leg, ncolumn)
Definition: compare.py:51
def overlay(hists, ytitle, header, addon)
Definition: compare.py:120
def hoverlay(hists, xtitle, ytitle, name)
Definition: compare.py:193
def makeEffPlotsVars(tree, varx, vary, sel, nbinx, xmin, xmax, nbiny, ymin, ymax, xtitle, ytitle, leglabel=None, header='', addon='', option='pt', marker=20)
Definition: compare.py:229
def ensureDir(directory)
Definition: compare.py:41
#define str(s)
def makeCompareVars(tree, var, sel, leglist, nbin, xmin, xmax, xtitle, ytitle, scale, name)
Definition: compare.py:62