RecoParticleFlow/test/compare.py

0001 #!/usr/bin/env python3
0002 import sys
0003 import os
0004 import ROOT
0005 import argparse
0006
0007 # This is an example of plotting the standard tracking validation
0008 # plots from an explicit set of DQM root files.
0009
0010 from Validation.RecoTrack.plotting.validation import SimpleValidation, SimpleSample
0011
0012 from Validation.RecoTrack.plotting.plotting import Subtract, FakeDuplicate, CutEfficiency, Transform, AggregateBins, ROC, Plot, PlotEmpty, PlotGroup, PlotOnSideGroup, PlotFolder, Plotter
0013 from Validation.RecoTrack.plotting.html import PlotPurpose
0014
0015 from Validation.RecoParticleFlow.defaults_cfi import ptbins, etabins, response_distribution_name, muLowOffset, muHighOffset, npvLowOffset, npvHighOffset, candidateType, offset_name
0016 from offsetStack import offsetStack
0017
0018 def parse_sample_string(ss):
0019     spl = ss.split(":")
0020     if not (len(spl) >= 3):
0021         raise Exception("Sample must be in the format name:DQMfile1.root:DQMfile2.root:...")
0022
0023     name = spl[0]
0024     files = spl[1:]
0025
0026     #check that all supplied files are actually ROOT files
0027     for fi in files:
0028         print("Trying to open DQM file {0} for sample {1}".format(fi, name))
0029         if not os.path.isfile(fi):
0030             raise Exception("Could not read DQM file {0}, it does not exist".format(fi))
0031         tf = ROOT.TFile(fi)
0032         if not tf:
0033             raise Exception("Could not read DQM file {0}, it's not a ROOT file".format(fi))
0034         #KH d = tf.Get("DQMData/Run 1/Physics/Run summary")
0035         d = tf.Get("DQMData/Run 1/ParticleFlow/Run summary")
0036         if not d:
0037             raise Exception("Could not read DQM file {0}, it's does not seem to be a harvested DQM file".format(fi))
0038     return name, files
0039
0040 def parse_plot_string(ss):
0041     spl = ss.split(":")
0042     if not (len(spl) >= 3):
0043         raise Exception("Plot must be in the format folder:name:hist1:hist2:...")
0044
0045     folder = spl[0]
0046     name = spl[1]
0047     histograms = spl[2:]
0048
0049     return folder, name, histograms
0050
0051 def parse_args():
0052     parser = argparse.ArgumentParser()
0053     parser.add_argument("-s", "--sample", type=str, action='append',
0054         required=True,
0055         help="DQM files to compare for a single sample, in the format 'name:file1.root:file2.root:...:fileN.root'",
0056         #default=[
0057         #    "QCD:DQM_V0001_R000000001__Global__CMSSW_X_Y_Z__RECO.root:DQM_V0001_R000000001__Global__CMSSW_X_Y_Z__RECO.root"
0058         #]
0059     )
0060     parser.add_argument("-p", "--plots",
0061         type=str, action='append',
0062         required=False,
0063         help="Plots to put on a single canvas, in the format 'folder:name:plot1:plot2:...:plotN'",
0064         default = [],
0065         #default=[
0066         #    "JetResponse:reso_dist_10_24:reso_dist_10_24_eta05:reso_dist_10_24_eta13"
0067         #]
0068     )
0069     parser.add_argument("--doResponsePlots",
0070         action='store_true',
0071         required=False,
0072         help="If enabled, do all jet response plots"
0073     )
0074     parser.add_argument("--doOffsetPlots",
0075         action='store_true',
0076         required=False,
0077         help="If enabled, do all offset plots"
0078     )
0079     parser.add_argument("--doMETPlots",
0080         action='store_true',
0081         required=False,
0082         help="If enabled, do all JetMET plots"
0083     )
0084     parser.add_argument("--doPFCandPlots",
0085         action='store_true',
0086         required=False,
0087         help="If enabled, do all PFCandidate plots"
0088     )
0089
0090     parser.add_argument( "--offsetVar", type=str,   action='store', default="npv", help="variable to bin offset eT" )
0091     parser.add_argument( "--offsetDR",  type=float, action='store', default=0.4,   help="offset deltaR value" )
0092     args = parser.parse_args()
0093
0094     #collect all the SimpleSample objects
0095     samples = []
0096     plots = []
0097
0098     sample_strings = args.sample
0099     for ss in sample_strings:
0100         name, files = parse_sample_string(ss)
0101         samp = SimpleSample(name, name, [(fn, fn.split('/')[-2]) for fn in files])
0102         samples += [samp]
0103
0104     for ss in args.plots:
0105         folder, name, histograms = parse_plot_string(ss)
0106         plots += [(folder, name, histograms)]
0107
0108
0109     # This needs to be also changed whenever changing binning
0110     if args.doResponsePlots:
0111         # Needs to add extra folders here if the DQM files have other folders of histograms
0112         JetFolderDirs = ["JetResponse/slimmedJets/JEC", "JetResponse/slimmedJets/noJEC", "JetResponse/slimmedJetsPuppi/JEC", "JetResponse/slimmedJetsPuppi/noJEC"]
0113
0114         for JetFolderDir in JetFolderDirs:
0115             plots += [(JetFolderDir, "reso_pt", ["preso_eta05", "preso_eta13",
0116                                                   "preso_eta21","preso_eta25","preso_eta30","preso_eta50"])]
0117             plots += [(JetFolderDir, "reso_pt_rms", ["preso_eta05_rms",
0118                                                       "preso_eta13_rms","preso_eta21_rms","preso_eta25_rms","preso_eta30_rms",
0119                                                       "preso_eta50_rms"])]
0120             plots += [(JetFolderDir, "response_pt", ["presponse_eta05",
0121                                                       "presponse_eta13", "presponse_eta21", "presponse_eta25", "presponse_eta30",
0122                                                       "presponse_eta50"])]
0123             for iptbin in range(len(ptbins)-1):
0124                 pthistograms = []
0125                 for ietabin in range(len(etabins)-1):
0126                     pthistograms += [response_distribution_name(iptbin, ietabin)]
0127                 plots += [(JetFolderDir, "response_{0:.0f}_{1:.0f}".format(ptbins[iptbin], ptbins[iptbin+1]), pthistograms)]
0128
0129     if args.doOffsetPlots:
0130         if args.offsetVar == "npv" :
0131             varHigh, varLow = npvHighOffset, npvLowOffset
0132         else :
0133             varHigh, varLow = muHighOffset, muLowOffset
0134         for ivar in range( varLow, varHigh ) :
0135             offsetHists = []
0136             for itype in candidateType :
0137                 offsetHists += [ offset_name( args.offsetVar, ivar, itype ) ]
0138             plots += [("Offset/{0}Plots/{0}{1}".format(args.offsetVar, ivar), "{0}{1}".format(args.offsetVar, ivar), offsetHists)]
0139
0140     if args.doMETPlots:
0141         doMETPlots(files, plots)
0142
0143
0144     if args.doPFCandPlots:
0145         doPFCandPlots(files, plots)
0146
0147     return samples, plots, args.doOffsetPlots, args.offsetVar, args.offsetDR, args.doPFCandPlots
0148
0149 # function that does METValidation from JetMET
0150 def doMETPlots(files, plots):
0151     #get the names of the histograms
0152     #MetValidation
0153     METHistograms = []
0154     f = ROOT.TFile(files[0])
0155     d = f.Get("DQMData/Run 1/JetMET/Run summary/METValidation/slimmedMETsPuppi")
0156     for i in d.GetListOfKeys():
0157         METHistograms.append([i.GetName()])
0158     # append plots
0159     METFolderDirs = ["METValidation/slimmedMETs","METValidation/slimmedMETsPuppi"]
0160     for METFolderDir in METFolderDirs:
0161         for  METHistogram in  METHistograms:
0162             plots += [(METFolderDir, "", METHistogram)]
0163
0164     #JetValidation
0165     JetHistograms = []
0166     d = f.Get("DQMData/Run 1/JetMET/Run summary/JetValidation/slimmedJets")
0167     for i in d.GetListOfKeys():
0168         JetHistograms.append([i.GetName()])
0169     JetValFolderDirs = ["JetValidation/slimmedJets", "JetValidation/slimmedJetsAK8", "JetValidation/slimmedJetsPuppi"]
0170     for JetValFolderDir in JetValFolderDirs:
0171         for JetHistogram in JetHistograms:
0172             plots += [(JetValFolderDir, "", JetHistogram)]
0173
0174 # does PFCandidate Plots
0175 def doPFCandPlots(files, plots):
0176     #we are going to hard code the end part of the histogram names because there's only 4
0177     hist_list = ["Charge", "Eta", "Phi", "Log10Pt", "PtLow","PtMid", "PtHigh"]
0178     f = ROOT.TFile(files[0])
0179     d = f.Get("DQMData/Run 1/ParticleFlow/Run summary/PackedCandidates")
0180     #get the name of the folders, which can use to complete plot name as well probably
0181     PFFolderNames = []
0182
0183     for i in d.GetListOfKeys():
0184         PFFolderNames.append(i.GetName())
0185
0186     for PFFolderName in PFFolderNames:
0187         for hist in hist_list:
0188             plots += [(PFFolderName, "", [PFFolderName + hist])]
0189
0190
0191 def addPlots(plotter, folder, name, section, histograms, opts, Offset=False):
0192     folders = [folder]
0193     #plots = [PlotGroup(name, [Plot(h, **opts) for h in histograms])]
0194     #KH print plots
0195     if Offset :
0196         plots = [PlotGroup(name, [Plot(h, **opts) for h in histograms])]
0197         plotter.append("Offset", folders, PlotFolder(*plots, loopSubFolders=False, page="offset", section=section))
0198     elif "JetResponse" in folder :
0199         plots = [PlotGroup(name, [Plot(h, **opts) for h in histograms])]
0200         plotter.append("ParticleFlow/" + section, folders, PlotFolder(*plots, loopSubFolders=False, page="pf", section=section))
0201         for plot in plots:
0202             plot.setProperties(ncols=3)
0203             plot.setProperties(legendDw=-0.68)
0204             plot.setProperties(legendDh=0.005)
0205             plot.setProperties(legendDy=0.24)
0206             plot.setProperties(legendDx=0.05)
0207     elif "JetMET" in folder:
0208         for h in histograms:
0209             plots = [PlotGroup(h, [Plot(h, **opts)])]
0210         for plot in plots:
0211             plot.setProperties(legendDw=-0.5)
0212             plot.setProperties(legendDh=0.01)
0213             plot.setProperties(legendDy=0.24)
0214             plot.setProperties(legendDx=0.05)
0215         plotter.append("JetMET" + section, folders, PlotFolder(*plots, loopSubFolders=False, page="JetMET", section=section))
0216     if "PackedCandidates" in folder:
0217         for h in histograms:
0218             if ("PtMid" in h or "PtHigh" in h):
0219                 plots = [PlotGroup(h, [Plot(h, ymin = pow(10,-1), ylog = True)])]
0220             else:
0221                 plots = [PlotGroup(h, [Plot(h, **opts)])]
0222         for plot in plots:
0223             plot.setProperties(legendDw=-0.5)
0224             plot.setProperties(legendDh=0.01)
0225             plot.setProperties(legendDy=0.24)
0226             plot.setProperties(legendDx=0.05)
0227         plotter.append("ParticleFlow/PackedCandidates/" + section, folders, PlotFolder(*plots, loopSubFolders=False, page="PackedCandidates", section= section))
0228
0229
0230 def main():
0231
0232     # plot-dependent style options
0233     # style options can be found from Validation/RecoTrack/python/plotting/plotting.py
0234     styledict_resolution = {"xlog": True, "xgrid":False, "ygrid":False,
0235         "xtitle":"GenJet pT (GeV)", "ytitle":"Jet pT resolution",
0236         "xtitleoffset":7.7,"ytitleoffset":3.8,"adjustMarginLeft":0.00}
0237
0238     styledict_response = {"xlog": True, "xgrid":False, "ygrid":False,
0239         "xtitle":"GenJet pT (GeV)", "ytitle":"Jet response",
0240         "xtitleoffset":7.7,"ytitleoffset":3.8,"adjustMarginLeft":0.00}
0241     plot_opts = {
0242         "reso_pt": styledict_resolution,
0243         "reso_pt_rms": styledict_resolution,
0244         "response_pt": styledict_response
0245     }
0246     for iptbin in range(len(ptbins)-1):
0247         plot_opts["response_{0:.0f}_{1:.0f}".format(ptbins[iptbin], ptbins[iptbin+1])] = {"stat": True}
0248
0249     samples, plots, doOffsetPlots, offsetVar, offsetDR, doPFCandPlots = parse_args()
0250
0251     plotter = Plotter()
0252
0253
0254     for folder, name, histograms in plots:
0255         opts = plot_opts.get(name, {})
0256
0257         #fullfolder =  "DQMData/Run 1/Physics/Run summary/{0}".format(folder)
0258         #fullfolder =  "DQMData/Run 1/ParticleFlow/Run summary/{0}".format(folder)
0259         fullJetFolder = "DQMData/Run 1/ParticleFlow/Run summary/{0}".format(folder)
0260         fullMETFolder = "DQMData/Run 1/JetMET/Run summary/{0}".format(folder)
0261         fullPFCandFolder = "DQMData/Run 1/ParticleFlow/Run summary/PackedCandidates/{0}".format(folder)
0262         print("Booking histogram group {0}={1} from folder {2}".format(name, histograms, folder))
0263         if "Offset/" in folder:
0264             opts = {'xtitle':'Default', 'ytitle':'Default'}
0265             addPlots(plotter, fullJetFolder, name, folder, histograms, opts, True)
0266         if "JetResponse" in folder:
0267             addPlots(plotter, fullJetFolder, name, folder, histograms, opts)
0268         if "METValidation" in folder or "JetValidation" in folder:
0269             addPlots(plotter, fullMETFolder, name, folder, histograms, opts)
0270         if doPFCandPlots:
0271             addPlots(plotter, fullPFCandFolder, name, folder, histograms, opts)
0272
0273     outputDir = "plots" # Plot output directory
0274     description = "Simple ParticleFlow comparison"
0275
0276     plotterDrawArgs = dict(
0277         separate=False, # Set to true if you want each plot in it's own canvas
0278     #    ratio=False,   # Uncomment to disable ratio pad
0279     )
0280
0281
0282     val = SimpleValidation(samples, outputDir)
0283     report = val.createHtmlReport(validationName=description)
0284     val.doPlots([plotter], plotterDrawArgs=plotterDrawArgs)
0285
0286     report.write()
0287
0288     #add tdr-style stack plots to offset html file
0289     if doOffsetPlots :
0290         offsetDir = "OffsetStacks"
0291         fullOffsetDir = os.path.join( outputDir, offsetDir )
0292         os.makedirs( fullOffsetDir )
0293
0294         for s in samples :
0295             offFile = open( outputDir + "/" + s.label() + "_offset.html", "r")
0296             lines = offFile.readlines()
0297             offFile.close()
0298             for f in s.files() :
0299                 fname = f.split('/')[-2]
0300                 outName = offsetStack( [(fname,f)], offsetVar, offsetDR, fullOffsetDir )
0301                 outName = outName.replace("plots/", "") #KH: This "plots" look redundant and causes trouble for .html. Stripping it off.
0302                 addLine( outName, lines )
0303
0304                 for f2 in s.files() :
0305                     if f == f2 : continue
0306                     fname2 = f2.split('/')[-2]
0307                     outName = offsetStack( [(fname,f), (fname2,f2)], offsetVar, offsetDR, fullOffsetDir )
0308                     outName = outName.replace("plots/", "") #KH: This "plots" look redundant and causes trouble for .html. Stripping it off.
0309                     addLine( outName, lines )
0310
0311             offFile = open( outputDir + "/" + s.label() + "_offset.html", "w")
0312             lines = "".join(lines)
0313             offFile.write(lines)
0314             offFile.close()
0315
0316 def addLine(name, oldLines) :
0317     newLines = [
0318         '   <td><a href="{0}">{0}</a></td>\n'.format(name),
0319         '  <br/>\n',
0320         '  <br/>\n'
0321     ]
0322     oldLines[8:len(newLines)] = newLines
0323
0324 if __name__ == "__main__":
0325     main()