Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:20:59

0001 #!/usr/bin/env python
0002 
0003 ## Utility used to overlay muon HLT plots from several releases
0004 
0005 ## For more information about available options, use the -h option:
0006 ##     ./compareReleases.py -h
0007 
0008 ## Import python libraries
0009 from __future__ import print_function
0010 import sys
0011 import optparse
0012 import os
0013 import re
0014 
0015 ## Set up ROOT in batch mode
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     ## Parse command line options
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     ## Define constants
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     ## Make table/plots
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 # Kludge
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 &plusmn; %.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                     ## hist.GetBinError(i + 1) kludge below!
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