RecoParticleFlow/test/offsetStack.py

0001 #chad harrington 2019
0002 #execute as python3 test/offsetStack.py -f HS_old:tmp/HS_old/DQMTotal.root HS:tmp/HS/DQMTotal.root -o .
0003 import argparse
0004 import ROOT
0005
0006 from Validation.RecoParticleFlow.defaults_cfi import candidateType,muHighOffset,npvHighOffset
0007
0008 ROOT.gROOT.SetBatch(True)
0009
0010 def arg2List( args ) :
0011   list = []
0012   for pair in args :
0013     split = pair.split(":")
0014     l = len(split)
0015     if l < 2 : raise Exception("Must provide a label and file, as in label:file")
0016
0017     label, file = split[0], ":".join(split[1:l]) #join in the case file contains ':'
0018     list += [(label, file)]
0019   return list
0020
0021 def main() :
0022   parser = argparse.ArgumentParser()
0023   parser.add_argument( "-f", "--files", nargs='+', required=True, help="must provide at least one file. Eg, label1:file1 label2:file2" )
0024   parser.add_argument( "-v", "--var",    type=str,   action='store', default="npv", help="variable to bin offset eT" )
0025   parser.add_argument( "-r", "--deltaR", type=float, action='store', default=0.4,   help="deltaR value" )
0026   parser.add_argument( "-o", "--outDir", type=str,   action='store', required=True, help="output directory" )
0027   args = parser.parse_args()
0028
0029   offsetStack( arg2List(args.files), args.var, args.deltaR, args.outDir )
0030
0031 def offsetStack( files, var, r, outdir ) :
0032
0033   file1 = ROOT.TFile.Open( files[0][1] )
0034   if file1 == None : raise Exception( "{} does not exist!".format(files[0][1]) )
0035   label1 = files[0][0]
0036   label = label1
0037
0038   h1 = file1.FindObjectAny( var )
0039   if h1 == None : raise Exception( "Could not find {} histogram in {}!".format(var, files[0][1]) )
0040   avg1 = int( h1.GetMean()+0.5 )
0041
0042   d_hist1 = getHists( file1, var, avg1, r )
0043
0044   setStyle()
0045   c = ROOT.TCanvas("c", "c", 600, 600)
0046
0047   stack1 = ROOT.THStack( "stack1", ";#eta;<Offset Energy_{{T}}(#pi{:.1f}^{{2}})>/<{}> [GeV]".format(r, "N_{PV}" if var=="npv" else "#mu") )
0048   setStack( stack1, d_hist1 )
0049   stack1.Draw("hist")
0050   stack1.SetMaximum( stack1.GetMaximum()*1.4 )
0051   stack1.Draw("hist")
0052
0053   legPF = "F"
0054   leg = ROOT.TLegend(.4,.67,.65,.92)
0055
0056   #file2 included
0057   if len(files) > 1 :
0058     file2 = ROOT.TFile.Open( files[1][1] )
0059     if file2 == None : raise Exception( "{} does not exist!".format(files[1][1]) )
0060     label2 = files[1][0]
0061     label = label1 + "_vs_" + label2
0062
0063     h2 = file2.FindObjectAny( var )
0064     if h2 == None : raise Exception( "Could not find {} histogram in {}!".format(var, files[1][1]) )
0065     avg2 = int( h2.GetMean()+0.5 )
0066
0067     d_hist2 = getHists( file2, var, avg2, r )
0068
0069     stack2 = ROOT.THStack( "stack2", "stack2" )
0070     setStack( stack2, d_hist2 )
0071     stack2.Draw("samepe");
0072
0073     legPF = "PF"
0074     leg.SetHeader( "#bf{{Markers: {}, Histograms: {}}}".format(label2, label1) )
0075
0076     #Draw Markers for EM Deposits and Hadronic Deposits in two separate regions
0077     hfe_clone = d_hist2["hfe"].Clone("hfe_clone")
0078     hfh_clone = d_hist2["hfh"].Clone("hfh_clone")
0079
0080     cloneStack = ROOT.THStack( "cloneStack", "cloneStack" )
0081     cloneStack.Add(d_hist2["ne"])
0082     cloneStack.Add(hfe_clone)
0083     cloneStack.Add(d_hist2["nh"])
0084     cloneStack.Add(hfh_clone)
0085     cloneStack.Draw("samepe")
0086
0087     #
0088     # Disable plotting restrictions for fixing Phase 2 offset plots.
0089     # Better to have busy plots than missing info...
0090     # Some clever trick for skipping zero fractions could be implemented
0091     # for the second PR
0092     # -Juska 24 June 2019
0093     #
0094
0095     #d_hist2["ne"] .SetAxisRange(-2.9,2.9)
0096     #d_hist2["hfe"].SetAxisRange(-5,-2.6)
0097     #d_hist2["nh"] .SetAxisRange(-2.9,2.9)
0098     #d_hist2["hfh"].SetAxisRange(-5,-2.6)
0099     #d_hist2["chu"].SetAxisRange(-2.9,2.9)
0100     #d_hist2["chm"].SetAxisRange(-2.9,2.9)
0101
0102     #hfe_clone     .SetAxisRange(2.6,5)
0103     #hfh_clone     .SetAxisRange(2.6,5)
0104
0105   leg.SetBorderSize(0)
0106   leg.SetFillColor(0)
0107   leg.SetFillStyle(0)
0108   leg.SetTextSize(0.04)
0109   leg.SetTextFont(42)
0110
0111   leg.AddEntry( d_hist1["ne"],  "Photons",                  legPF )
0112   leg.AddEntry( d_hist1["hfe"], "EM Deposits",              legPF )
0113   leg.AddEntry( d_hist1["nh"],  "Neutral Hadrons",          legPF )
0114   leg.AddEntry( d_hist1["hfh"], "Hadronic Deposits",        legPF )
0115   leg.AddEntry( d_hist1["chu"], "Unassoc. Charged Hadrons", legPF )
0116   leg.AddEntry( d_hist1["chm"], "Assoc. Charged Hadrons",   legPF )
0117
0118   leg.Draw()
0119
0120   text = ROOT.TLatex()
0121   text.SetNDC()
0122
0123   text.SetTextSize(0.065)
0124   text.SetTextFont(61)
0125   text.DrawLatex(0.2, 0.87, "CMS")
0126
0127   text.SetTextSize(0.045)
0128   text.SetTextFont(42)
0129   text.DrawLatex(1-len(label)/41., 0.96, label)
0130
0131   outName = outdir + "/stack_" + label + ".pdf"
0132   c.Print( outName )
0133   return outName
0134
0135 def getHists( file, var, var_val, r ) :
0136   dict = {}
0137
0138   var_val_range=var_val;
0139   if var=="mu" and var_val>=muHighOffset: var_val_range=muHighOffset-1
0140   if var=="npv" and var_val>=npvHighOffset: var_val_range=npvHighOffset-1
0141
0142   for pf in candidateType :
0143
0144     name = "p_offset_eta_{}{}_{}".format( var, var_val_range, pf )
0145     p = file.FindObjectAny(name)
0146     if p == None : raise Exception( "Could not find {} profile in {}!".format(name, file.GetName()) )
0147     dict[pf] = p.ProjectionX( pf )
0148     dict[pf].Scale( r*r / 2 / var_val )
0149
0150     xbins = p.GetXaxis().GetXbins().GetArray()
0151     for i in range(1, p.GetNbinsX()+1) :
0152       dict[pf].SetBinContent( i, dict[pf].GetBinContent(i) / (xbins[i]-xbins[i-1]) )
0153       dict[pf].SetBinError( i, dict[pf].GetBinError(i) / (xbins[i]-xbins[i-1]) )
0154
0155   return dict
0156
0157 def setStack( stack, hists ) :
0158
0159   stack.Add( hists["ne"] )
0160   stack.Add( hists["hfe"] )
0161   stack.Add( hists["nh"] )
0162   stack.Add( hists["hfh"] )
0163   stack.Add( hists["chu"] )
0164   stack.Add( hists["chm"] )
0165
0166   hists["ne"] .SetMarkerStyle(ROOT.kMultiply)
0167   hists["hfe"].SetMarkerStyle(ROOT.kOpenStar)
0168   hists["nh"] .SetMarkerStyle(ROOT.kOpenDiamond)
0169   hists["hfh"].SetMarkerStyle(ROOT.kOpenTriangleUp)
0170   hists["chu"].SetMarkerStyle(ROOT.kOpenCircle)
0171   hists["chm"].SetMarkerStyle(ROOT.kOpenCircle)
0172
0173   hists["ne"] .SetFillColor(ROOT.kBlue)
0174   hists["hfe"].SetFillColor(ROOT.kViolet+2)
0175   hists["nh"] .SetFillColor(ROOT.kGreen)
0176   hists["hfh"].SetFillColor(ROOT.kPink+6)
0177   hists["chu"].SetFillColor(ROOT.kRed-9)
0178   hists["chm"].SetFillColor(ROOT.kRed)
0179
0180   hists["ne"] .SetLineColor(ROOT.kBlack)
0181   hists["hfe"].SetLineColor(ROOT.kBlack)
0182   hists["nh"] .SetLineColor(ROOT.kBlack)
0183   hists["hfh"].SetLineColor(ROOT.kBlack)
0184   hists["chu"].SetLineColor(ROOT.kBlack)
0185   hists["chm"].SetLineColor(ROOT.kBlack)
0186
0187 def setStyle() :
0188
0189   ROOT.gStyle.SetPadTopMargin(0.05)
0190   ROOT.gStyle.SetPadBottomMargin(0.1)
0191   ROOT.gStyle.SetPadLeftMargin(0.16)
0192   ROOT.gStyle.SetPadRightMargin(0.02)
0193
0194   ROOT.gStyle.SetOptStat(0)
0195   ROOT.gStyle.SetOptTitle(0)
0196
0197   ROOT.gStyle.SetTitleFont(42, "XYZ")
0198   ROOT.gStyle.SetTitleSize(0.05, "XYZ")
0199   ROOT.gStyle.SetTitleXOffset(0.9)
0200   ROOT.gStyle.SetTitleYOffset(1.4)
0201
0202   ROOT.gStyle.SetLabelFont(42, "XYZ")
0203   ROOT.gStyle.SetLabelOffset(0.007, "XYZ")
0204   ROOT.gStyle.SetLabelSize(0.04, "XYZ")
0205
0206   ROOT.gStyle.SetPadTickX(1)
0207   ROOT.gStyle.SetPadTickY(1)
0208
0209 if __name__ == "__main__":
0210     main()