HcalRecAlgos/test/plotMahi.py

0001 #!/usr/bin/env python
0002
0003 # Plot MAHI multi-pulse fit results
0004 import ROOT
0005 ROOT.gStyle.SetOptStat(0)
0006 c = ROOT.TCanvas('c','c',500,500)
0007 c.cd()
0008 c.SetRightMargin(0.05)
0009 c.SetLeftMargin(0.12)
0010 c.SetTopMargin(0.1)
0011 c.SetTopMargin(0.05)
0012
0013 color = [ROOT.kAzure+3, ROOT.kAzure+2, ROOT.kAzure+1, ROOT.kYellow-7, ROOT.kRed+1, ROOT.kRed+2, ROOT.kRed+3, ROOT.kRed+4]
0014
0015 import os
0016 plotdir = 'plotMahi'
0017 if not os.path.exists(plotdir):
0018     os.mkdir(plotdir)
0019
0020 # Parameters for Run 3 HBHE
0021 nTS = 10 # number of TS to plot
0022 nOOT = 7 # number of OOT samples
0023 offset = 3 # soi
0024
0025 hists = {}
0026 hists['digi'] = ROOT.TH1D('digi', ';TS;E [GeV]', nTS, -0.5, nTS-0.5)
0027 hists['digi'].SetMarkerStyle(20)
0028 for i in range(nOOT+1):
0029     hists[i] = ROOT.TH1D('bx%i' % (i-3), ';TS;E [GeV]', nTS, -0.5, nTS-0.5)
0030     hists[i].SetFillColor(color[i])
0031     hists[i].SetLineColor(ROOT.kBlack)
0032
0033 tags = ['']#, '_shape206', '_timeslew', '_mt6', '_qcd', '_qcdnopu']
0034
0035 from tqdm import tqdm
0036
0037 for tag in tags:
0038     tree = ROOT.TChain('mahiDebugger/HcalTree')
0039     tree.Add('mahidebugger%s.root' % tag)
0040     count = 0
0041     for rh in tqdm(tree):
0042         if abs(rh.ieta) > 16:
0043             continue
0044         soiEnergy = rh.mahiEnergy*rh.inGain
0045         if soiEnergy < 1:
0046             continue
0047
0048         energy = {}
0049         soi = rh.soi
0050         assert soi == offset
0051
0052         for i in range(nTS-1):
0053             hists['digi'].SetBinContent(i+1, rh.inputTS[i]*rh.inGain)
0054             hists[soi].SetBinContent(i+1, rh.itPulse[i]*rh.mahiEnergy*rh.inGain)
0055             energy[soi] = soiEnergy
0056
0057         for o in range(nOOT):
0058             oh = o+1 if o>=soi else o
0059             ootPulse = []
0060             for i in range(nTS*o, nTS*(o+1)):
0061                 ootPulse.append(rh.ootPulse[i])
0062
0063             for i in range(nTS):
0064                 hists[oh].SetBinContent(i+1, max(0, ootPulse[i]*rh.ootEnergy[o]*rh.inGain))
0065                 energy[oh] = rh.ootEnergy[o]*rh.inGain
0066
0067         stack = ROOT.THStack('stack', '')
0068         for i in range(nOOT+1):
0069             stack.Add(hists[i])
0070
0071         hists['digi'].GetXaxis().SetRangeUser(-0.5, 7.5)
0072         hists['digi'].GetYaxis().SetRangeUser(0, max(hists['digi'].GetBinContent(4), hists['digi'].GetBinContent(5))*1.5)
0073         hists['digi'].Draw('P0')
0074         stack.Draw('same')
0075         hists['digi'].Draw('P0,same')
0076
0077         legend = ROOT.TLegend(0.6,0.5,0.93,0.93)
0078         legend.SetLineWidth(0)
0079         legend.SetFillStyle(0)
0080         legend.AddEntry(hists['digi'], 'Digi', 'P')
0081         for i in range(nOOT+1):
0082             legend.AddEntry(hists[i], 'BX %+i, E = %.1f GeV' % (i-offset, energy[i]), 'F')
0083
0084         tex = ROOT.TLatex()
0085         tex.SetTextSize(0.025)
0086         tex.DrawLatexNDC(0.15, 0.90, 'run:'+str(rh.run)+' evt:'+str(rh.evt))
0087         tex.DrawLatexNDC(0.15, 0.85, 'ieta:'+str(rh.ieta)+' iphi:'+str(rh.iphi)+' depth:'+str(rh.depth))
0088         tex.DrawLatexNDC(0.15, 0.80, '#chi^{2} = %.1f, TDC=%i, mahiTime=%.1f' % (rh.chiSq, rh.inputTDC[3], rh.arrivalTime))
0089
0090         wTime = 0.
0091         sumEn = 0.
0092         for i in range(len(energy)):
0093             thisEnergy = rh.inputTS[i]*rh.inGain
0094             wTime += thisEnergy * i
0095             sumEn += thisEnergy
0096         if sum(energy) > 0:
0097             wTime /= sumEn
0098         tex.DrawLatexNDC(0.15, 0.75, 'Mean time = %.1f TS' % (wTime))
0099
0100         legend.Draw()
0101
0102         c.Print(plotdir+'/'+str(rh.run)+'_'+str(rh.evt)+'_'+str(rh.ieta)+'_'+str(rh.iphi)+'_'+str(rh.depth)+tag+'.pdf')
0103         c.Print(plotdir+'/'+str(rh.run)+'_'+str(rh.evt)+'_'+str(rh.ieta)+'_'+str(rh.iphi)+'_'+str(rh.depth)+tag+'.png')
0104         count += 1