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