ATLAS Offline Software
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 
10 import ROOT
11 import os
12 import sys
13 import itertools
14 from array import array
15 ROOT.gROOT.SetBatch( True )
16 
17 doProjX=1
18 doProjY=0
19 
20 #--------------------------------------------------------------------------
21 def 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 
48 def 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 #--------------------------------------------------------------------------
193 def 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 
376 if __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 )
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition: hcg.cxx:307
MuonValidation_CreateEffAndRecoFracPlots.SetBinomialError
def SetBinomialError(ratio, den)
Definition: MuonValidation_CreateEffAndRecoFracPlots.py:21
vtune_athena.format
format
Definition: vtune_athena.py:14
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
plot_material.mkdir
def mkdir(path, recursive=True)
Definition: plot_material.py:16
RootHelpers::FindBin
Int_t FindBin(const TAxis *axis, const double x)
Definition: RootHelpers.cxx:14
MuonValidation_CreateEffAndRecoFracPlots.CreateRatioPlot
def CreateRatioPlot(infile, numHist, denHist, var, projXorY=doProjY, x1=None, x2=None, xtitle='', plottype='', doAverage=False)
Definition: MuonValidation_CreateEffAndRecoFracPlots.py:48
MuonValidation_CreateEffAndRecoFracPlots.main
def main(argv)
Definition: MuonValidation_CreateEffAndRecoFracPlots.py:193
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:179
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
calibdata.exit
exit
Definition: calibdata.py:236
print
void print(char *figname, TCanvas *c1)
Definition: TRTCalib_StrawStatusPlots.cxx:25
array
TCS::join
std::string join(const std::vector< std::string > &v, const char c=',')
Definition: Trigger/TrigT1/L1Topo/L1TopoCommon/Root/StringUtils.cxx:10
Trk::split
@ split
Definition: LayerMaterialProperties.h:38
Generated on Sun Dec 22 2024 21:15:33 for ATLAS Offline Software by   doxygen 1.8.18