ATLAS Offline Software
Loading...
Searching...
No Matches
MuonValidation_CreateEffAndRecoFracPlots.py
Go to the documentation of this file.
1# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
2
3#Written by Dan Mori
4#Create efficiency plots for each author by dividing matched/author pt by truth/all pt
5
6#Usage: python CreateEffAndRecoFracPlots.py filename [doAverage]
7
8#Added flag doAverage to add TF1 of line showing overall efficiency/reco fraction
9
10import ROOT
11import os
12import sys
13import itertools
14from array import array
15ROOT.gROOT.SetBatch( True )
16
17doProjX=1
18doProjY=0
19
20#--------------------------------------------------------------------------
21def SetBinomialError( ratio, den ):
22 if 'TH2' in ratio.IsA().GetName():
23 for j in range(1,ratio.GetNbinsY()+1):
24 for i in range(1,ratio.GetNbinsX()+1):
25 n = den.Integral(i,i,j,j)
26 p = ratio.GetBinContent(i,j)
27 if n <= 0 or p>=1:
28 ratio.SetBinError( i, j, 0 )
29 continue
30 else:
31 ratio.SetBinError( i, j, (p*(1-p)/n)**0.5 )
32 elif 'TH1' in ratio.IsA().GetName():
33 for i in range(1,ratio.GetNbinsX()+1):
34 n = den.Integral(i,i)
35 p = ratio.GetBinContent(i)
36 if n <= 0 or p>=1:
37 ratio.SetBinError( i, 0 )
38 continue
39 else:
40 ratio.SetBinError( i, (p*(1.-p)/n)**0.5 )
41 else:
42 print( 'WARNING ' + ratio.GetName() + 'is not a TH1 or TH2' )
43
44#--------------------------------------------------------------------------
45#create histograms that require dividing one histogram by another
46#works for both efficiency and reco fraction plots, set by plottype variable
47
48def CreateRatioPlot( infile, numHist, denHist, var, projXorY=doProjY, x1=None, x2=None, xtitle = '', plottype = '', doAverage = False ):
49 #require plottype variable for setting plot title, etc
50 #plottype = 'eff' or 'RecoFrac' to create efficiency or reco fraction histogram
51 if plottype == 'eff':
52 typeName = 'efficiency'
53 elif plottype == 'RecoFrac':
54 typeName = 'Reco Fraction'
55 elif plottype == 'purity':
56 typeName = 'purity'
57 else:
58 print('plottype must be one of the following: \'eff\', \'RecoFrac\', \'purity\'' )
59 return
60
61 PlotNamePrefix = numHist.GetName()[:numHist.GetName().find('_kinematics')] + '_' + typeName.replace(' ','')
62
63 #create ratio hist for given eta range
64 if (x1 != None and x2 != None):
65 if '_' in var:
66 if projXorY==doProjX:
67 xvar=var.split('_')[1]
68 yvar=var.split('_')[0]
69 else:
70 xvar=var.split('_')[0]
71 yvar=var.split('_')[1]
72 #PlotName = PlotNamePrefix + '_' + var.split('_')[1] + '_etaRange_{0}_{1}'.format(x1,x2).replace('-','m').replace('.','p')
73 #PlotTitle = var.split('_',1)[1].capitalize() + ' ' + typeName
74 PlotName = PlotNamePrefix + '_' + yvar + '_'+xvar+'Range_{0}_{1}'.format(x1,x2).replace('-','m').replace('.','p')
75 PlotTitle = yvar.capitalize() + ' ' + typeName
76
77 if x1>=0 and x2>0:
78 if projXorY==doProjX:
79 ibin1 = numHist.GetYaxis().FindBin(x1)
80 ibin2 = numHist.GetYaxis().FindBin(x2)
81 num = numHist.ProjectionX('num_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
82 den = denHist.ProjectionX('den_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
83 else:
84 ibin1 = numHist.GetXaxis().FindBin(x1)
85 ibin2 = numHist.GetXaxis().FindBin(x2)
86 num = numHist.ProjectionY('num_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
87 den = denHist.ProjectionY('den_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
88
89 if yvar=='eta':
90 if projXorY==doProjX:
91 ibin1 = numHist.GetYaxis().FindBin(-1*x2)
92 ibin2 = numHist.GetYaxis().FindBin(-1*x1)
93 num2 = numHist.ProjectionX('num2_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
94 den2 = denHist.ProjectionX('den2_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
95 else:
96 ibin1 = numHist.GetXaxis().FindBin(-1*x2)
97 ibin2 = numHist.GetXaxis().FindBin(-1*x1)
98 num2 = numHist.ProjectionY('num2_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
99 den2 = denHist.ProjectionY('den2_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
100 num.Add(num2)
101 den.Add(den2)
102 PlotTitle = PlotTitle + ' ({0}<|eta|<{1})'.format(x1,x2)
103 del num2, den2
104 else:
105 if yvar=='pt':
106 PlotTitle = PlotTitle + ' ({0}<pt<{1})'.format(x1,x2)
107 # else:
108 # if projXorY==doProjX:
109 # ibin1 = numHist.GetYaxis().FindBin(x1)
110 # ibin2 = numHist.GetYaxis().FindBin(x2)
111 # num = numHist.ProjectionX('num_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
112 # den = denHist.ProjectionX('den_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
113 # #PlotTitle = PlotTitle + ' (pt>{0})'.format(eta1)
114 # PlotTitle = PlotTitle + ' ({0}<{2}<{1})'.format(x1,x2,yvar)
115 # PlotName = PlotNamePrefix + '_' + var.split('_')[1] + '_ptMin{0}'.format(x1).replace('-','m').replace('.','p')
116 # else:
117 # ibin1 = numHist.GetXaxis().FindBin(x1)
118 # ibin2 = numHist.GetXaxis().FindBin(x2)
119 # num = numHist.ProjectionY('num_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
120 # den = denHist.ProjectionY('den_'+xvar+"_"+yvar+"Range",ibin1,ibin2)
121 # PlotTitle = PlotTitle + ' ({0}<{2}<{1})'.format(x1,x2,yvar)
122 else:
123 PlotTitle = var.capitalize() + ' ' + typeName
124 PlotName = PlotNamePrefix + '_' + var
125 num = numHist.Clone()
126 den = denHist.Clone()
127
128 #create variable bin histogram for pt to group high bins together
129 if ('_pt' in var and projXorY==doProjY) or var=='pt':
130 #rebin the initial pt hist that is in bins of 1 GeV
131 ptBins = [ 0, 5, 10, 15, 20, 25, 30, 40, 50, 65, 80, 100, 200 ]
132 #ptBins = [ 0, 2,3,4,5,6, 8, 15 ]
133 numRebinned = ROOT.TH1D( 'num_'+var, '', len(ptBins)-1, array( 'f', ptBins ) )
134 denRebinned = ROOT.TH1D( 'den_'+var, '', len(ptBins)-1, array( 'f', ptBins ) )
135 for i in range( 1, numHist.GetNbinsX() + 1 ):
136 numRebinned.Fill( num.GetXaxis().GetBinCenter(i), num.GetBinContent(i) )
137 denRebinned.Fill( den.GetXaxis().GetBinCenter(i), den.GetBinContent(i) )
138 #rebin eta and phi
139 else:
140 if 'TH1' in den.IsA().GetName():
141 numRebinned = num.Clone()
142 denRebinned = den.Clone()
143 numRebinned.Rebin(2)
144 denRebinned.Rebin(2)
145 else:
146 return
147
148
151
152 ratio = numRebinned.Clone( PlotName )
153 ratio.Divide( numRebinned, denRebinned, 1, 1 )
154 ratio.SetTitle( PlotTitle )
155 SetBinomialError( ratio, denRebinned ) #root binomial error is different - custom done
156
157 # Define x and y-axis name from histogram name (var)
158 if 'TH2' in ratio.IsA().GetName():
159 if '_' in var:
160 ratio.GetXaxis().SetTitle( var.split('_')[0] )
161 ratio.GetYaxis().SetTitle( var.split('_')[1] )
162 ratio.GetZaxis().SetTitle( typeName.lower() )
163 else:
164 print 'WARNING Could not get axes name from histogram name'
165 elif 'TH1' in ratio.IsA().GetName():
166 ymax = ratio.GetBinContent( ratio.GetMaximumBin() ) + ratio.GetBinError( ratio.GetMaximumBin() )
167 ratio.GetYaxis().SetRangeUser(0, ymax*1.2 )
168 ratio.GetYaxis().SetTitle( typeName.lower() )
169 if xtitle == '':
170 xtitle = var
171 ratio.GetXaxis().SetTitle( xtitle )
172
173 if doAverage and den.Integral() > 0 and num.Integral() > 0:
174 tot_eff = ROOT.TF1( 'aveline', "[0]", ratio.GetXaxis().GetBinLowEdge(1), ratio.GetXaxis().GetBinUpEdge( ratio.GetNbinsX() ) )
175 ratio.GetListOfFunctions().Add( tot_eff )
176 f1 = ratio.GetFunction( "aveline" )
177 f1.SetParameter( 0, num.Integral()/den.Integral() )
178 f1.SetLineColor( ROOT.kRed )
179 #f1.SetBit(ROOT.TF1.kNotDraw)
180
181 PlotDirName = PlotNamePrefix.replace('__','_').split('_')
182 if not infile.GetDirectory( '/'.join( PlotDirName ) ):
183 PlotDir = infile.Get( '/'.join( PlotDirName[:-1] ) ).mkdir( PlotDirName[-1] )
184 else:
185 PlotDir = infile.Get( '/'.join( PlotDirName ) )
186 PlotDirName = '/'.join( PlotDirName )
187
188 if not infile.GetDirectory(PlotDirName).WriteTObject( ratio, PlotName, "Overwrite" ):
189 print('WARNING failed to write histogram to file: ' + PlotDirName + '/' + PlotName )
190 del ratio, num, den
191
192#--------------------------------------------------------------------------
193def main( argv ):
194 """
195 Main function to be executed when starting the code.
196 """
197
198 if len( argv ) < 2:
199 print( 'No filename given' )
200 print( 'Usage: python '+argv[0]+' physval_filename [doAverage]' )
201 exit(1)
202
203 filename = argv[1]
204 if not os.path.exists( filename ):
205 print ( 'File not found: ' + filename )
206 exit(1)
207
208 if len(argv) > 2 and argv[2] == 'doAverage':
209 doAverage = True
210 else:
211 doAverage = False
212
213 infile = ROOT.TFile.Open( filename, 'update' )
214
215 muonTypesEff = [ 'All', 'Prompt', 'InFlight', 'NonIsolated' ]
216 muonTypesReco = [ 'Prompt', 'InFlight', 'NonIsolated', 'Rest' ]
217 Variables = [ 'pt', 'eta', 'phi', 'eta_phi', 'eta_pt' ]
218
219 Xtitles = {
220 'pt' : 'p_{T} [GeV]',
221 'eta' : '|#eta|',
222 'phi' : '#phi',
223 'eta_phi' : '#phi',
224 'eta_pt' : 'p_{T} [GeV]' }
225
226 #Efficiency plots
227 Authors = []
228 for muType in muonTypesEff:
229 if not infile.Get( 'Muons/' + muType ):
230 print( 'INFO TDirectory not found: Muons/' + muType )
231 continue
232 #get list of authors from matched dir
233 AuthDir = infile.Get( 'Muons/{0}/matched'.format( muType ) )
234 if Authors == []:
235 Authors = [ i.GetName() for i in AuthDir.GetListOfKeys() if AuthDir.Get( i.GetName() ).InheritsFrom( 'TDirectory' ) ]
236 for author in Authors:
237 truthDirName = 'Muons/{0}/truth/all'.format( muType )
238 matchDirName = 'Muons/{0}/matched/{1}/kinematics'.format( muType, author )
239 truthDir = infile.GetDirectory( truthDirName )
240 matchDir = infile.GetDirectory( matchDirName )
241 if not truthDir:
242 print( 'WARNING Directory not found: '+truthDirName )
243 continue
244 if not matchDir:
245 print( 'WARNING Directory not found: '+matchDirName )
246 continue
247 for var in Variables:
248 truthHistName = truthDirName.replace('/','_') + '_' + var
249 truthHist = truthDir.Get( truthHistName )
250 matchHistName = matchDirName.replace('/','_') + '_' + var
251 matchHist = matchDir.Get( matchHistName )
252 if not truthHist:
253 print( 'WARNING histogram not found: '+truthHistName )
254 continue
255 if not matchHist:
256 print( 'WARNING histogram not found: '+matchHistName )
257 continue
258 CreateRatioPlot( infile, matchHist, truthHist, var, xtitle = muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
259 if var == 'eta_phi' or var == 'eta_pt':
260 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 0, 2.5, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
261 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 0, 0.1, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
262 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 0.1, 1.05, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
263 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 1.05, 2.0, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
264 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 2.0, 2.5, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
265 CreateRatioPlot( infile, matchHist, truthHist, var, doProjY, 2.5, 2.7, muType+' Muon '+Xtitles[var], plottype = 'eff', doAverage = doAverage )
266 if var == 'eta_pt' and author == 'CaloTag' :
267 CreateRatioPlot( infile, matchHist, truthHist, var, doProjX, 10, 1000, muType+' Muon '+Xtitles['eta'], plottype = 'eff', doAverage = doAverage )
268 CreateRatioPlot( infile, matchHist, truthHist, var, doProjX, 15, 1000, muType+' Muon '+Xtitles['eta'], plottype = 'eff', doAverage = doAverage )
269 CreateRatioPlot( infile, matchHist, truthHist, var, doProjX, 20, 1000, muType+' Muon '+Xtitles['eta'], plottype = 'eff', doAverage = doAverage )
270 CreateRatioPlot( infile, matchHist, truthHist, var, doProjX, 25, 1000, muType+' Muon '+Xtitles['eta'], plottype = 'eff', doAverage = doAverage )
271 #Reco Fraction plots
272 for muType in muonTypesReco:
273 if not infile.Get( 'Muons/' + muType ):
274 print( 'INFO TDirectory not found: Muons/' + muType )
275 continue
276 #get list of authors from matched dir
277 AuthDir = infile.Get( 'Muons/{0}/matched'.format( muType ) )
278 if Authors == []:
279 Authors = [ i.GetName() for i in AuthDir.GetListOfKeys() if AuthDir.Get( i.GetName() ).InheritsFrom( 'TDirectory' ) ]
280 for author in Authors:
281 typedir = 'Muons/{0}/reco/{1}/kinematics'.format( muType, author )
282 alldir = 'Muons/All/reco/{0}/kinematics'.format( author )
283 typeRecoDir = infile.Get( typedir )
284 allRecoDir = infile.Get( alldir )
285 if not typeRecoDir:
286 print( 'INFO TDirectory not found: '+typedir )
287 continue
288 if not allRecoDir:
289 print( 'INFO TDirectory not found: '+alldir )
290 continue
291 for var in Variables:
292 typeplot = typedir.replace('/','_') + '_' + var
293 allplot = alldir.replace('/','_') + '_' + var
294 typeRecoHist = typeRecoDir.Get( typeplot )
295 allRecoHist = allRecoDir.Get( allplot )
296 if not typeRecoHist:
297 print( 'WARNING plot not found: ' + typeplot )
298 continue
299 if not allRecoHist:
300 print( 'WARNING plot not found: ' + allplot )
301 continue
302 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, xtitle = muType + ' Muon ' + Xtitles[var], plottype = 'RecoFrac', doAverage = doAverage )
303
304 #Purity plots (matched/reco)
305
306 muType = 'All'
307 if not infile.Get( 'Muons/' + muType ):
308 print( 'INFO TDirectory not found: Muons/' + muType )
309 else:
310 #get list of authors from matched dir
311 AuthDir = infile.Get( 'Muons/{0}/matched'.format( muType ) )
312 if Authors == []:
313 Authors = [ i.GetName() for i in AuthDir.GetListOfKeys() if AuthDir.Get( i.GetName() ).InheritsFrom( 'TDirectory' ) ]
314 for author in Authors:
315 typedir = 'Muons/{0}/matched/{1}/kinematicsReco'.format( muType, author )
316 alldir = 'Muons/{0}/reco/{1}/kinematics'.format( muType, author )
317 typeRecoDir = infile.Get( typedir )
318 allRecoDir = infile.Get( alldir )
319 if not typeRecoDir:
320 print( 'INFO TDirectory not found: '+typedir )
321 continue
322 if not allRecoDir:
323 print( 'INFO TDirectory not found: '+alldir )
324 continue
325 for var in Variables:
326 typeplot = typedir.replace('/','_') + '_' + var
327 allplot = alldir.replace('/','_') + '_' + var
328 typeRecoHist = typeRecoDir.Get( typeplot )
329 allRecoHist = allRecoDir.Get( allplot )
330 if not typeRecoHist:
331 print( 'WARNING plot not found: ' + typeplot )
332 continue
333 if not allRecoHist:
334 print( 'WARNING plot not found: ' + allplot )
335 continue
336 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, xtitle = muType + ' Muon ' + Xtitles[var], plottype = 'purity', doAverage = doAverage )
337 if var == 'eta_pt':
338 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjX, 4, 500, muType+' Muon '+Xtitles['pt'], plottype = 'purity', doAverage = doAverage )
339 if var == 'eta_pt' and author == 'CaloTag' :
340 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjX, 10, 500, muType+' Muon '+Xtitles['pt'], plottype = 'purity', doAverage = doAverage )
341 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjX, 15, 500, muType+' Muon '+Xtitles['pt'], plottype = 'purity', doAverage = doAverage )
342 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjX, 20, 500, muType+' Muon '+Xtitles['pt'], plottype = 'purity', doAverage = doAverage )
343 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjX, 25, 500, muType+' Muon '+Xtitles['pt'], plottype = 'purity', doAverage = doAverage )
344 if (var == 'eta_phi' or var == 'eta_pt') and author is not 'CaloTag' :
345 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 0, 0.1, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
346 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 0.1, 1.05, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
347 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 1.05, 2.0, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
348 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 2.0, 2.5, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
349 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 0, 2.5, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
350 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, doProjY, 2.5, 2.7, muType+' Muon '+Xtitles[var], plottype = 'purity', doAverage = doAverage )
351
352 #unmatched muon reco fraction
353 muType = 'UnmatchedRecoMuons'
354 typedir = 'Muons/{0}/kinematics'.format( muType )
355 allnames = [ i for i in Authors if i == 'AllMuons' or i == 'AllAuthors' ]
356 if len(allnames) == 0:
357 return
358 alldir = 'Muons/All/reco/{0}/kinematics'.format( allnames[0] )
359 if not infile.GetDirectory( typedir ):
360 print( 'INFO directory not found: ' + typedir )
361 elif not infile.GetDirectory( alldir ):
362 print( 'INFO directory not found: ' + alldir )
363 else:
364 for var in Variables:
365 typeplot = 'Muons_{0}__kinematics_{1}'.format( muType, var )
366 allplot = alldir.replace('/','_') + '_{0}'.format( var )
367 #print('Working on Muons_{0}__{1}'.format(muType,var))
368 typeRecoHist = infile.GetDirectory( typedir ).Get( typeplot )
369 allRecoHist = infile.GetDirectory( alldir ).Get( allplot )
370 if typeRecoHist and allRecoHist:
371 CreateRatioPlot( infile, typeRecoHist, allRecoHist, var, xtitle = 'Unmatched Reco Muon '+Xtitles[var], plottype = 'RecoFrac', doAverage = doAverage )
372 infile.Close()
373
374#===============================================================================
375
376if __name__ == "__main__":
377 """
378 Here the code should appear that is executed when running the plotter directly
379 (and not import it in another python file via 'import Plotter')
380 """
381 # start main program
382 main( sys.argv )
print
void print(char *figname, TCanvas *c1)
Definition TRTCalib_StrawStatusPlots.cxx:26
array
STL class.
Get
T * Get(TFile &f, const std::string &n, const std::string &dir="", const chainmap_t *chainmap=0, std::vector< std::string > *saved=0)
get a histogram given a path, and an optional initial directory if histogram is not found,...
Definition comparitor.cxx:181
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition hcg.cxx:310
find
std::string find(const std::string &s)
return a remapped string
Definition hcg.cxx:138
split
std::vector< std::string > split(const std::string &s, const std::string &t=":")
Definition hcg.cxx:177
main
int main()
Definition hello.cxx:18
MuonValidation_CreateEffAndRecoFracPlots.SetBinomialError
SetBinomialError(ratio, den)
Definition MuonValidation_CreateEffAndRecoFracPlots.py:21
MuonValidation_CreateEffAndRecoFracPlots.CreateRatioPlot
CreateRatioPlot(infile, numHist, denHist, var, projXorY=doProjY, x1=None, x2=None, xtitle='', plottype='', doAverage=False)
Definition MuonValidation_CreateEffAndRecoFracPlots.py:48
Generated on for ATLAS Offline Software by doxygen 1.14.0