File indexing completed on 2024-04-06 12:18:55
0001
0002
0003
0004
0005
0006
0007
0008
0009 from __future__ import print_function
0010 import sys
0011 import optparse
0012 import os
0013 import re
0014
0015
0016 if '-h' not in sys.argv:
0017 sys.argv.append('-b')
0018 import ROOT
0019 ROOT.gROOT.Macro('rootlogon.C')
0020 sys.argv.remove('-b')
0021 ROOT.gStyle.SetOptStat(0)
0022 ROOT.gStyle.SetOptFit(0)
0023 ROOT.gStyle.SetErrorX(0.5)
0024 ROOT.gErrorIgnoreLevel = ROOT.kWarning
0025 c1 = ROOT.TCanvas("c1","c1",800,600)
0026 c1.GetFrame().SetBorderSize(6)
0027 c1.GetFrame().SetBorderMode(-1)
0028
0029
0030
0031 class RootFile:
0032 def __init__(self, file_name):
0033 self.name = file_name[0:file_name.find(".root")]
0034 self.file = ROOT.TFile(file_name, "read")
0035 if self.file.IsZombie():
0036 print("Error opening %s, exiting..." % file_name)
0037 sys.exit(1)
0038 def Get(self, object_name):
0039 return self.file.Get(object_name)
0040
0041
0042
0043 def main():
0044
0045 usage="usage: %prog [options] file1.root file2.root file3.root ..."
0046 parser = optparse.OptionParser(usage=usage)
0047 parser.add_option('-g', '--gen', action="store_true", default=False,
0048 help="skip plots which match to reconstructed muons")
0049 parser.add_option('-r', '--rec', action="store_true", default=False,
0050 help="skip plots which match to generated muons")
0051 parser.add_option('-s', '--stats', action="store_true", default=False,
0052 help="print a table rather than plotting")
0053 parser.add_option('--stats-format', default="twiki",
0054 help="choose 'twiki' (default) or 'html' table format")
0055 parser.add_option('-p', '--path', default="HLT_Mu9",
0056 help="specify which HLT path to focus on")
0057 parser.add_option('--file-type', default="pdf", metavar="EXT",
0058 help="choose an output format other than pdf")
0059 options, arguments = parser.parse_args()
0060 path = options.path
0061 source_types = ["gen", "rec"]
0062 if options.gen: source_types.remove("rec")
0063 if options.rec: source_types.remove("gen")
0064 if len(arguments) == 0:
0065 print("Please provide at least one file as an argument")
0066 return
0067
0068
0069 global file_type, path_muon, path_dist, cross_channel_format, colors
0070 file_type = options.file_type
0071 path_muon = "/DQMData/Run 1/HLT/Run summary/Muon"
0072 path_dist = "%s/Distributions" % path_muon
0073 cross_channel_format = re.compile("HLT_[^_]*_[^_]*$")
0074 colors = [ROOT.kBlue, ROOT.kRed, ROOT.kGreen+2, ROOT.kMagenta+2,
0075 ROOT.kYellow+2, ROOT.kCyan+2, ROOT.kBlue+3, ROOT.kRed+3]
0076
0077
0078 files, path_names = get_files_and_path_names(arguments)
0079 efficiencies = get_global_efficiencies(files, path)
0080 filter_names, me_values = get_filters_and_MEs(files, path)
0081 if options.stats:
0082 for source_type in source_types:
0083 make_efficiency_table(files, source_type, filter_names,
0084 efficiencies, options.stats_format)
0085 sys.exit(0)
0086 print("Generating plots for %s path" % path)
0087 make_efficiency_summary_plot(files, me_values)
0088 if len(files) == 1:
0089 for source_type in source_types:
0090 for filter in filter_names:
0091 make_overlaid_turnon_plot(files, source_type, filter, path,
0092 efficiencies, me_values)
0093 plot_vars = ["TurnOn1", "EffEta", "EffPhi"]
0094 if "SingleMu" in files[0].name or len(files) == 1:
0095 plot_vars.remove("TurnOn1")
0096 for source_type in source_types:
0097 for plot_var in plot_vars:
0098 for filter in filter_names:
0099 plot_name = "%s%s_%s" % (source_type, plot_var, filter)
0100 make_efficiency_plot(files, plot_name, path,
0101 efficiencies, me_values)
0102 if "SingleMu" not in files[0].name:
0103 for source_type in source_types[0:1]:
0104 for path in path_names:
0105 if "HighMultiplicity" in path: continue
0106 filter_names, me_values = get_filters_and_MEs(files, path)
0107 for filter in filter_names:
0108 plot_name = "%s%s_%s" % ("gen", "TurnOn1", filter)
0109 make_efficiency_plot(files, plot_name, path,
0110 efficiencies, me_values)
0111 if file_type == "pdf":
0112 merge_pdf_output(files)
0113
0114
0115
0116 def counter_generator():
0117 k = 0
0118 while True:
0119 k += 1
0120 yield k
0121 next_counter = counter_generator().next
0122
0123
0124
0125 def get_files_and_path_names(arguments):
0126 files = [RootFile(filename) for filename in arguments]
0127 path_names = []
0128 keys = files[0].file.GetDirectory(path_dist).GetListOfKeys()
0129 obj = keys.First()
0130 while obj:
0131 if obj.IsFolder():
0132 path_names.append(obj.GetName())
0133 obj = keys.After(obj)
0134 return files, path_names
0135
0136
0137
0138 def get_global_efficiencies(files, path):
0139 efficiencies = []
0140 for file_index, file in enumerate(files):
0141 try:
0142 eff_hist = ROOT.TH1F(file.Get("%s/%s/globalEfficiencies" %
0143 (path_dist, path)))
0144 except:
0145 print("No global efficiency in %s" % file.name)
0146 sys.exit(1)
0147 eff_hist.LabelsDeflate()
0148 efficiencies.append({})
0149 for i in range(eff_hist.GetNbinsX()):
0150 label = eff_hist.GetXaxis().GetBinLabel(i + 1)
0151 eff = 100. * eff_hist.GetBinContent(i + 1)
0152 err = 100. * eff_hist.GetBinError(i + 1)
0153 efficiencies[file_index][label] = [eff, err]
0154 return efficiencies
0155
0156
0157
0158 def get_filters_and_MEs(files, path):
0159 filter_names = []
0160 me_values = {}
0161 regexp = ROOT.TPRegexp("^<(.*)>(i|f|s)=(.*)</\\1>$")
0162 full_path = "%s/%s" % (path_dist, path)
0163 keys = files[0].file.GetDirectory(full_path).GetListOfKeys()
0164 obj = keys.First()
0165 while obj:
0166 name = obj.GetName()
0167 filter_name = name[name.find("_") + 1:len(name)]
0168 if "genPassPhi" in name:
0169 if "L1" in name or "L2" in name or "L3" in name:
0170 filter_names.append(filter_name)
0171 if regexp.Match(name):
0172 array = ROOT.TObjArray(regexp.MatchS(name))
0173 me_values[array.At(1).GetName()] = float(array.At(3).GetName())
0174 obj = keys.After(obj)
0175 return filter_names, me_values
0176
0177
0178
0179 def make_efficiency_table(files, source_type, filter_names,
0180 efficiencies, stats_format):
0181 if "twiki" in stats_format:
0182 format = "| %s | %s |"
0183 separator = " | "
0184 if "html" in stats_format:
0185 format = "<tr><th>%s<td>%s"
0186 separator = "<td>"
0187 print("<table>")
0188 sample = files[0].name[0:files[0].name.find("_")]
0189 print((format % (sample, separator.join(filter_names))).replace('td', 'th'))
0190 for file_index, file in enumerate(files):
0191 entries = []
0192 for filter in filter_names:
0193 label = "%sEffPhi_%s" % (source_type, filter)
0194 eff, err = efficiencies[file_index][label]
0195 entries.append("%.1f ± %.1f" % (eff, err))
0196 version = file.name
0197 version = version[version.find("_")+1:version.find("-")]
0198 print(format % (version, separator.join(entries)))
0199 if "html" in stats_format:
0200 print("</table>")
0201
0202
0203 def make_efficiency_summary_plot(files, me_values):
0204 class Bin:
0205 def __init__(self, label, content, error):
0206 self.label = label
0207 self.content = float(content)
0208 self.error = float(error)
0209 hists = ROOT.TList()
0210 hist_path = "%s/Summary/Efficiency_Summary" % path_muon
0211 for file_index, file in enumerate(files):
0212 hist = file.Get(hist_path)
0213 if not hist:
0214 hist = file.Get(hist_path + "_Muon")
0215 if not hist:
0216 print("No efficiency summary plot found in %s" % file.name)
0217 return
0218 n_bins = hist.GetNbinsX()
0219 bins = [Bin(hist.GetXaxis().GetBinLabel(i + 1),
0220 hist.GetBinContent(i + 1),
0221
0222 hist.GetBinContent(i + 1)/1000)
0223 for i in range(n_bins)]
0224 n_cross = i = 0
0225 while i != (len(bins) - n_cross):
0226 if cross_channel_format.match(bins[i].label):
0227 bins.append(bins.pop(i))
0228 n_cross += 1
0229 else:
0230 i += 1
0231 new_hist = ROOT.TH1F("%s_clone%i" % (hist.GetName(), file_index),
0232 file.name, n_bins,
0233 hist.GetXaxis().GetBinLowEdge(1),
0234 hist.GetXaxis().GetBinUpEdge(n_bins))
0235 for i, bin in enumerate(bins):
0236 new_hist.GetXaxis().SetBinLabel(i + 1, bin.label)
0237 new_hist.SetBinContent(i + 1, bin.content)
0238 new_hist.SetBinError(i + 1, bin.error)
0239 hists.Add(new_hist)
0240 if hists.At(0):
0241 plot(hists, "Summary", "Efficiency Summary")
0242
0243
0244
0245 def make_efficiency_plot(files, plot_name, path, efficiencies, me_values):
0246 hists = ROOT.TList()
0247 hist_path = "%s/%s/%s" % (path_dist, path, plot_name)
0248 hist_title = files[0].Get(hist_path).GetTitle()
0249 for file in files:
0250 try: hist = ROOT.TH1F(file.Get(hist_path))
0251 except:
0252 try: hist = ROOT.TProfile(file.Get(hist_path))
0253 except: print("Failed to find %s!!" % plot_name); return
0254 new_hist = hist.Clone()
0255 new_hist.SetTitle(file.name)
0256 hists.Add(new_hist)
0257 if hists.At(0):
0258 plot(hists, plot_name, hist_title, efficiencies, me_values, path)
0259
0260
0261
0262 def make_overlaid_turnon_plot(files, source_type, filter, path,
0263 efficiencies, me_values):
0264 hists = ROOT.TList()
0265 for plot_name in ["%sTurnOn%s_%s" % (source_type, number, filter)
0266 for number in [1, 2]]:
0267 hist_path = "%s/%s/%s" % (path_dist, path, plot_name)
0268 try: hist = ROOT.TH1F(files[0].Get(hist_path))
0269 except: hist = ROOT.TProfile(files[0].Get(hist_path))
0270 new_hist = hist.Clone()
0271 x_title = hist.GetXaxis().GetTitle()
0272 x_title = x_title.replace("Leading", "Leading/Next-to-Leading")
0273 new_hist.SetTitle("%s;%s" % (new_hist.GetTitle(), x_title))
0274 hists.Add(new_hist)
0275 if hists.At(0):
0276 plot(hists, "%sTurnOn", "Turn-On Curves for %s" % filter,
0277 efficiencies, me_values, path)
0278
0279
0280
0281 def plot(hists, hist_type, hist_title,
0282 efficiencies=None, me_values=None, path=None):
0283 plot_num = next_counter()
0284 first = hists.First()
0285 last = hists.Last()
0286 hist = last
0287 x_title = first.GetXaxis().GetTitle()
0288 if "Gen" in x_title: source_type = "gen"
0289 elif "Rec" in x_title: source_type = "reco"
0290 expression = ("0.5 * [2] * (" +
0291 "TMath::Erf((x / [0] + 1.) / (TMath::Sqrt(2.) * [1])) + " +
0292 "TMath::Erf((x / [0] - 1.) / (TMath::Sqrt(2.) * [1])))")
0293 function_turn_on = ROOT.TF1("turnOn", expression, 10, 40)
0294 while hist:
0295 hist_index = int(hists.IndexOf(hist))
0296 hist.Draw()
0297 ROOT.gPad.SetLogx(False)
0298 ROOT.gPad.SetTickx(1)
0299 title = hist.GetTitle()
0300 hist.SetLineWidth(2)
0301 hist.SetLineColor(colors[hist_index % len(colors)])
0302 if hist.GetMaximum() <= 1.0:
0303 hist.Scale(100.)
0304 hist.GetYaxis().SetRangeUser(0., 100.)
0305 if "Summary" in hist_type:
0306 hist.GetXaxis().LabelsOption("u")
0307 c1.SetTickx(0)
0308 if "Eff" in hist_type or "TurnOn" in hist_type:
0309 yTitle = hist.GetYaxis().GetTitle()
0310 slashIndex = yTitle.find("/")
0311 yTitle = "#frac{%s}{%s} (%%)" % (yTitle[0:slashIndex - 1],
0312 yTitle[slashIndex + 2:100])
0313 hist.GetYaxis().SetTitle(yTitle)
0314 hist.GetYaxis().SetTitleOffset(1.5)
0315 hist.GetYaxis().SetTitleSize(0.04)
0316 if "Eff" in hist_type:
0317 eff, err = efficiencies[hist_index][hist_type]
0318 hist.SetTitle("%s (%.1f#pm%.1f%%)" % (title, eff, err))
0319 if "TurnOn" in hist_type:
0320 ROOT.gPad.SetLogx(True)
0321 hist.GetXaxis().SetRangeUser(2., 300.)
0322 function_turn_on.SetParameters(1,20,100)
0323 function_turn_on.SetLineColor(hist.GetLineColor())
0324 hist.Fit("turnOn","q")
0325 hist.Draw()
0326 eff = function_turn_on.GetParameter(2)
0327 if eff < 100 and eff > 0:
0328 hist.SetTitle("%s (%.1f%%)" % (title, eff))
0329 hist = hists.Before(hist)
0330 last.Draw()
0331 hist = hists.Before(last)
0332 while hist:
0333 hist.Draw("same")
0334 hist = hists.Before(hist)
0335 lower_bound_y = 0.15
0336 upper_bound_y = lower_bound_y + (0.055 * hists.Capacity())
0337 if "Summary" in hist_type:
0338 cuts = None
0339 elif "TurnOn" in hist_type:
0340 cuts = "|#eta| < %.1f" % me_values["CutMaxEta"]
0341 else:
0342 cuts = "p_{T} #geq %.0f, |#eta| < %.1f" % (me_values["CutMinPt"],
0343 me_values["CutMaxEta"])
0344 legend = ROOT.gPad.BuildLegend(0.22, lower_bound_y, 0.90, upper_bound_y)
0345 legend.SetFillColor(0)
0346 legend.SetFillStyle(1001)
0347 if "Summary" not in hist_type:
0348 latex = ROOT.TLatex()
0349 latex.SetNDC()
0350 latex.SetTextAlign(31)
0351 latex.DrawLatex(0.93, upper_bound_y + 0.015,
0352 "Cuts on #mu_{%s}: %s" % (source_type, cuts))
0353 hist_title = "%.3i: %s" % (plot_num, hist_title)
0354 if path: hist_title += " step in %s" % path
0355 last.SetTitle(hist_title)
0356 ROOT.gPad.Update()
0357 c1.SaveAs("%.3i.%s" % (plot_num, file_type))
0358
0359
0360
0361 def merge_pdf_output(files):
0362 os.system("gs -q -dBATCH -dNOPAUSE -sDEVICE=pdfwrite -dAutoRotatePages=/All -sOutputFile=merged.pdf [0-9][0-9][0-9].pdf")
0363 if len(files) == 1:
0364 pdfname = "%s.pdf" % files[0].name
0365 elif len(files) == 2:
0366 pdfname = "%s_vs_%s.pdf" % (files[0].name, files[1].name)
0367 else:
0368 pdfname = "%s_vs_Many.pdf" % files[0].name
0369 os.system("cp merged.pdf %s" % pdfname)
0370 os.system("rm [0-9]*.pdf")
0371 print("Wrote %i plots to %s" % (next_counter() - 1, pdfname))
0372
0373
0374
0375 if __name__ == "__main__":
0376 sys.exit(main())
0377
