ATLAS Offline Software
MuonValidation_CreateResolutionProfiles.py
Go to the documentation of this file.
1 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
2 
3 #Dan Mori <dmori@cern.ch>
4 #This macro creates muon resolution profiles from 2D histograms
5 #Usage: python CreateResolutionProfiles.py filename [doAverage doBinned]
6 #Adding the flag doAverage adds a TF1 showing the overall pT scale/resolution to the profiles
7 #Adding the flag doBinned creates resolution plots for each pt, eta, phi bin for prompt muons
8 
9 import sys
10 import os
11 from array import array
12 import itertools
13 from math import sqrt
14 import ROOT
15 from ROOT import *
16 
17 #---------------------------------------------------------------------------
18 
19 def SetYrange( hist ):
20  #adjust y range so negative bins will show
21  ymax = hist.GetBinContent( hist.GetMaximumBin() ) + hist.GetBinError( hist.GetMaximumBin() )
22  ymin = hist.GetBinContent( hist.GetMinimumBin() ) - hist.GetBinError( hist.GetMinimumBin() )
23  #c = 0.2*( ymax - ymin )
24  #ymax += c
25  #ymin -= c
26  ymax += 0.2*ymax
27  ymin -= 0.2*ymax
28  hist.GetYaxis().SetRangeUser( ymin, ymax )
29 
30 #---------------------------------------------------------------------------
31 def GetProfilesFromTH2( hist, doAverage ):
32  binedges = [ hist.GetXaxis().GetBinLowEdge(1) ]
33  binedges += [ hist.GetXaxis().GetBinUpEdge(i) for i in range( 1, hist.GetNbinsX()+1 ) ]
34  if hist.GetXaxis().GetBinUpEdge(hist.GetNbinsX()) == 1000 and hist.GetBinContent(hist.GetNbinsX()) < 3: #remove last bin (200-1000) if nothing there
35  binedges = binedges[:-1]
36  nBins = len(binedges) - 1
37 
38  prof_ave = ROOT.TH1D( 'ave', 'ave', nBins, array( 'f', binedges ) )
39  prof_std = ROOT.TH1D( 'std', 'std', nBins, array( 'f', binedges ) )
40  binnedRes = {}
41  for i in range( 1, nBins + 1 ):
42  h = hist.ProjectionY( hist.GetName()+'_p_y', i, i )
43  n = h.GetEntries()
44  if n<20:
45  continue
46 
47  xmin = h.GetMean()-3*h.GetRMS()
48  xmax = h.GetMean()+3*h.GetRMS()
49  h.GetXaxis().SetRangeUser(xmin,xmax)
50  print '>>> ' , h.GetName()
51  mean, meanErr, sigma, sigmaErr, frame = fit(h, -0.5, 0.5)
52  prof_ave.SetBinContent( i, mean )
53  prof_ave.SetBinError( i, meanErr )
54  prof_std.SetBinContent( i, sigma )
55  prof_std.SetBinError( i, sigmaErr )
56  name = '{0}_{1}'.format( binedges[i-1], binedges[i] )
57  binnedRes[name] = frame.Clone(name)
58  #del h, fit
59  SetYrange( prof_ave )
60  SetYrange( prof_std )
61  return prof_ave,prof_std,binnedRes
62 #---------------------------------------------------------------------------
63 #---------------------------------------------------------------------------
64 
65 # Simple fit method ................................................
66 # F = f1 * Gauss(mean,sigma) + (1-f1) Gauss(mean,sigma) x Exp(tau)
67 # ..................................................................
68 
69 def fit(h, emin,emax):
70  x = RooRealVar('x','x',emin,emax);
71  if (h.GetEntries()<20):
72  return 0,0,0,0,x.frame()
73  mean = RooRealVar('mean','mean',-0.001,-0.01,0.01);
74  sigma = RooRealVar('sigma','#sigma',0.02, 0.01, 0.3);
75  Gauss = RooGaussian('Gauss','Gauss',x,mean,sigma);
76  f1 = RooRealVar('f1','f1',0.8,0.,1.);
77 
78  GaussModel = RooGaussModel('GaussModel','Gauss',x,mean,sigma)
79  tau = RooRealVar('tau','tau',0.02,0.,1.)
80  tailmodel = RooDecay('tailmodel','ExpxGaus',x,tau,GaussModel,RooDecay.DoubleSided);
81  model = RooAddPdf('model','G1+E2',Gauss,tailmodel,f1)
82 
83  hdata = RooDataHist('hdata','hist', RooArgList(x), h);
84  model.fitTo(hdata,RooFit.PrintLevel(-1),RooFit.Verbose(kFALSE));
85  frame = x.frame();
86 
87  hdata.plotOn(frame)
88  model.plotOn(frame,RooFit.LineColor(ROOT.kRed))
89  model.plotOn(frame,RooFit.Components('tailmodel'),RooFit.DrawOption('F'),RooFit.FillColor(ROOT.kOrange),RooFit.MoveToBack())
90  return mean.getVal(), mean.getError(), sigma.getVal(), sigma.getError(), frame
91 #---------------------------------------------------------------------------
92 #---------------------------------------------------------------------------
93 
94 #create pt res and scale profiles vs pt, eta, phi from 2D histograms
95 def CreateProfile( infile, HistDir, HistName, Var, MuonType, doAverage = False, doBinned = False ):
96  hist = infile.Get(HistDir).Get(HistName)
97  if hist.IsA().InheritsFrom(ROOT.TProfile.Class()):
98  return
99  xtitle = MuonType + ' Muon '
100  if 'pT' in Var:
101  xtitle += 'pT [GeV]'
102  else:
103  xtitle += Var
104 
105  prof_ave, prof_std, binnedRes = GetProfilesFromTH2( hist, doAverage = doAverage )
106 
107  prof_ave.SetName( HistName.replace( 'pT_vs', 'PtScale_vs' ) )
108  prof_ave.SetTitle( 'pT Scale vs ' + Var )
109  prof_ave.SetXTitle( xtitle )
110  prof_ave.SetYTitle( 'pT scale' )
111 
112  prof_std.SetName( HistName.replace( 'pT_vs', 'PtResol_vs' ) )
113  prof_std.SetTitle( 'pT Resolution vs ' + Var )
114  prof_std.SetXTitle( xtitle )
115  prof_std.SetYTitle( 'pT resolution' )
116 
117  if not infile.Get( HistDir ).WriteTObject( prof_ave, prof_ave.GetName(), "Overwrite" ):
118  # print( 'INFO Writing histogram to file: ' + HistDir + '/' + prof_ave.GetName() )
119  print( 'WARNING failed to write histogram to file: ' + HistDir + '/' + prof_ave.GetName() )
120  if not infile.Get( HistDir ).WriteTObject( prof_std, prof_std.GetName(), "Overwrite" ):
121  # print( 'INFO Writing histogram to file: ' + HistDir + '/' + prof_std.GetName() )
122  print( 'WARNING failed to write histogram to file: ' + HistDir + '/' + prof_std.GetName() )
123 
124  #create binned res plots
125  bindirname = Var.lower() + 'Bins'
126  bindirpath = HistDir + '/' + bindirname
127  bindir = infile.GetDirectory( bindirpath )
128  if not bindir:
129  bindir = infile.GetDirectory( HistDir ).mkdir( bindirname )
130 
131  if doBinned:
132  canv = ROOT.TCanvas("canv","",800,800);
133  canv.Divide(3,3)
134  ipad = 0
135  icanv = 0
136  storePlots = 'resPlots/'+os.path.split(infile.GetName())[1].replace('.root','')
137  if not os.path.exists(storePlots):
138  os.makedirs(storePlots)
139  print('Creating directory: '+storePlots)
140 
141  for x1,name,plot in sorted([ (float(name.split('_')[0]), name, plot) for name, plot in binnedRes.iteritems() ]):
142  plot.SetName( 'bin_' + name.replace('-','m') )
143  plot.SetTitle( 'pT Resolution {0} < {1} < {2}'.format( name.split('_')[0], Var, name.split('_')[1] ) )
144  plot.SetYTitle( 'Entries' )
145  t=ROOT.TCanvas()
146 # if not bindir.WriteTObject( plot, plot.GetName(), "Overwrite" ):
147 # print('WARNING failed to write histogram to file: ' + bindirpath + '/' + plot.GetName() )
148  ipad+=1
149  if ipad>9:
150  canv.cd()
151  canv.SaveAs(storePlots+'/'+HistDir.replace('/','_')+'_PtResFits_{0}_bins_{1}.pdf'.format(Var,icanv))
152  canv.Clear()
153  canv.Divide(3,3)
154  icanv+=1
155  ipad=1
156  canv.cd(ipad)
157  t = ROOT.TLatex()
158  t.SetNDC(); t.SetTextColor(1);
159  tit = name.replace('m','-').replace('_','<'+Var+'<')
160  mu = mu_err = sigma = sigma_err = 0
161  canv.cd(ipad)
162  plot.Draw()
163  t.DrawLatex(0.2,0.96,plot.GetTitle())
164  resultMeanLab ='#mu = {0:0.2g} #pm {1:0.2g}'.format( mu, mu_err )
165  resultSigmaLab ='#sigma = {0:0.2g} #pm {1:0.2g}'.format( sigma, sigma_err )
166  t.DrawLatex(0.2,0.85,'#splitline{'+resultMeanLab+'}{'+resultSigmaLab+'}')
167 
168  canv.cd()
169  canv.SaveAs(storePlots+'/'+HistDir.replace('/','_')+'_PtResFits_{0}_bins_{1}.pdf'.format(Var,icanv))
170  canv.Close()
171 
172  del prof_ave, prof_std, binnedRes
173 #---------------------------------------------------------------------------
174 
175 def main( args ):
176  if len( args ) > 1:
177  filename = args[1]
178  else:
179  print( 'Usage: python {0} filename [doAverage doBinned]'.format( args[0] ) )
180  return
181 
182  if not os.path.exists( filename ):
183  print ( 'File not found: ' + filename )
184  return
185 
186  doAverage = False
187  doBinned = True
188  if len(args) > 2:
189  doAverage = bool( 'doAverage' in args[2:] )
190  doBinned = bool( 'doBinned' in args[2:] )
191 
192  print( 'Opening file: ' + filename )
193  infile = ROOT.TFile.Open( filename, 'update' )
194 
195  MuonTypes = [ 'Prompt' ]
196  BinVars = [ 'pT','lowpT','highpT','eta', 'phi' ]
197 
198  for MuonType in MuonTypes:
199  if not infile.Get( 'Muons/' + MuonType ):
200  print( 'INFO TDirectory not found: Muons/' + MuonType )
201  continue
202  AuthDir = infile.Get( 'Muons/{0}/matched'.format( MuonType ) )
203  Authors = [ i.GetName() for i in AuthDir.GetListOfKeys() if AuthDir.Get( i.GetName() ).InheritsFrom( 'TDirectory' ) ]
204  for Author in Authors:
205  if Author=='MSTrackParticles' or Author=='CaloTag' or Author=='AllMuons' or Author=='Tight' or Author=='Loose' or Author=='VeryLoose':
206  continue
207  HistDir = 'Muons/{0}/matched/{1}/resolution'.format( MuonType, Author )
208  FileDir = infile.Get( HistDir )
209  if not FileDir:
210  print( 'INFO TDirectory not found: ' + HistDir )
211  continue
212  for Var in BinVars:
213  resTypes= [ '','ID','MS' ]
214  for resType in resTypes:
215  HistName = '_'.join( HistDir.split('/') ) + '_Res{0}_pT_vs_{1}'.format( resType, Var )
216  if not FileDir.Get( HistName ):
217  #print( 'INFO Histogram not found: ' + HistName )
218  continue
219  CreateProfile( infile, HistDir, HistName, Var, MuonType, doAverage = doAverage, doBinned = doBinned )
220 
221  infile.Close()
222 
223 #---------------------------------------------------------------------------
224 
225 if __name__ == "__main__":
226  ROOT.gROOT.SetBatch()
227  main( sys.argv )
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition: hcg.cxx:307
vtune_athena.format
format
Definition: vtune_athena.py:14
MuonValidation_CreateResolutionProfiles.main
def main(args)
Definition: MuonValidation_CreateResolutionProfiles.py:175
plot_material.mkdir
def mkdir(path, recursive=True)
Definition: plot_material.py:16
MuonValidation_CreateResolutionProfiles.CreateProfile
def CreateProfile(infile, HistDir, HistName, Var, MuonType, doAverage=False, doBinned=False)
Definition: MuonValidation_CreateResolutionProfiles.py:95
MuonValidation_CreateResolutionProfiles.GetProfilesFromTH2
def GetProfilesFromTH2(hist, doAverage)
Definition: MuonValidation_CreateResolutionProfiles.py:31
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
DerivationFramework::TriggerMatchingUtils::sorted
std::vector< typename T::value_type > sorted(T begin, T end)
Helper function to create a sorted vector from an unsorted one.
MuonValidation_CreateResolutionProfiles.fit
def fit(h, emin, emax)
Definition: MuonValidation_CreateResolutionProfiles.py:69
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
MuonValidation_CreateResolutionProfiles.SetYrange
def SetYrange(hist)
Definition: MuonValidation_CreateResolutionProfiles.py:19
xAOD::bool
setBGCode setTAP setLVL2ErrorBits bool
Definition: TrigDecision_v1.cxx:60
readCCLHist.float
float
Definition: readCCLHist.py:83
Generated on Sun Dec 22 2024 21:15:33 for ATLAS Offline Software by   doxygen 1.8.18