test/kalmanTools/validation.py

0001 from __future__ import print_function
0002 import ROOT,itertools,math      #
0003 from array import array         #
0004 from DataFormats.FWLite import Events, Handle
0005 ROOT.FWLiteEnabler.enable()
0006 #
0007
0008
0009
0010
0011 tag='output'
0012
0013
0014 ##A class to keep BMTF data
0015
0016
0017
0018
0019 ###Common methods############
0020
0021 def fetchStubsOLD(event,ontime=False,isData=True):
0022     phiSeg    = Handle  ('L1MuDTChambPhContainer')
0023     if not isData:
0024         event.getByLabel('simTwinMuxDigis',phiSeg)
0025     else:
0026         event.getByLabel('bmtfDigis',phiSeg)
0027     if ontime:
0028         filtered=filter(lambda x: x.bxNum()==0, phiSeg.product().getContainer())
0029         return filtered
0030     else:
0031         return phiSeg.product().getContainer()
0032
0033
0034 def fetchStubs(event,ontime=True):
0035     phiSeg2    = Handle  ('std::vector<L1MuKBMTCombinedStub>')
0036     event.getByLabel('simKBmtfStubs',phiSeg2)
0037     if ontime:
0038         filtered=filter(lambda x: x.bxNum()==0, phiSeg2.product())
0039         return filtered
0040     else:
0041         return phiSeg2.product()
0042 def globalBMTFPhi(muon):
0043     temp=muon.processor()*48+muon.hwPhi()
0044     temp=temp*2*math.pi/576.0-math.pi*15.0/180.0;
0045     if temp>math.pi:
0046         temp=temp-2*math.pi;
0047
0048     K=1.0/muon.hwPt()
0049     if muon.hwSign()>0:
0050         K=-1.0/muon.hwPt()
0051     return temp+5.740*K
0052
0053
0054
0055 def fetchKMTF(event,etaMax,collection):
0056     kbmtfH  = Handle  ('BXVector<l1t::RegionalMuonCand>')
0057     event.getByLabel(collection,kbmtfH)
0058     kbmtf=kbmtfH.product()
0059     kbmtfMuons={}
0060     for bx in [-3,-2,-1,0,1,2,3]:
0061         kbmtfMuons[bx]=[]
0062     for bx in range(kbmtf.getFirstBX(),kbmtf.getLastBX()+1):
0063         for j in range(0,kbmtf.size(bx)):
0064             mu = kbmtf.at(bx,j)
0065             kbmtfMuons[bx].append(mu)
0066 #        kbmtfMuons[bx]=sorted(kbmtfMuons[bx],key=lambda x: x.hwPt(),reverse=True)
0067     return kbmtfMuons
0068
0069 def curvResidual(a,b):
0070     return (a.charge()/a.pt()-b.charge()/b.pt())*b.pt()/b.charge()
0071
0072 def ptResidual(a,b):
0073     return (a.pt()-b.pt())/b.pt()
0074
0075 def curvResidualSTA(a,b):
0076     return (a.charge()/a.ptUnconstrained()-b.charge()/b.pt())*b.pt()/b.charge()
0077
0078
0079
0080
0081 def deltaPhi( p1, p2):
0082     '''Computes delta phi, handling periodic limit conditions.'''
0083     res = p1 - p2
0084     while res > math.pi:
0085         res -= 2*math.pi
0086     while res < -math.pi:
0087         res += 2*math.pi
0088     return res
0089
0090 def deltaR( *args ):
0091     return math.sqrt( deltaR2(*args) )
0092
0093 def deltaR2( e1, p1, e2, p2):
0094     de = e1 - e2
0095     dp = deltaPhi(p1, p2)
0096     return de*de + dp*dp
0097
0098
0099 def log(event,counter,mystubs,kmtf,bmtf):
0100     print("--------EVENT"+str(counter)+"------------")
0101     print('RUN={run} LUMI={lumi} EVENT={event}'.format(run=event.eventAuxiliary().id().run(),lumi=event.eventAuxiliary().id().luminosityBlock(),event=event.eventAuxiliary().id().event()))
0102     print("-----------------------------")
0103     print("-----------------------------")
0104     print('Stubs:')
0105     for stub in mystubs:
0106         print('wheel={w} sector={sc} station={st} high/low={ts} phi={phi} phiB={phiB} qual={qual} BX={BX}'.format(w=stub.whNum(),sc=stub.scNum(),st=stub.stNum(),ts=stub.Ts2Tag(),phi=stub.phi(),phiB=stub.phiB(),qual=stub.code(),BX=stub.bxNum()))
0107     print('EMU:')
0108     for g in bmtf :
0109         print("EMU sector={sector} pt={pt} eta={eta} phi={phi} qual={qual} dxy={dxy} pt2={pt2} hasFineEta={HF}".format(sector=g.processor(), pt=g.hwPt(),eta=g.hwEta(),phi=g.hwPhi(),qual=g.hwQual(),dxy=g.hwDXY(),pt2=g.hwPtUnconstrained(),HF=g.hwHF()))
0110     print('DATA:')
0111     for g in kmtf :
0112         print("DATA sector={sector} pt={pt} eta={eta} phi={phi} qual={qual} dxy={dxy} pt2={pt2} hasFineEta={HF}".format(sector=g.processor(),pt=g.hwPt(),eta=g.hwEta(),phi=g.hwPhi(),qual=g.hwQual(),dxy=g.hwDXY(),pt2=g.hwPtUnconstrained(),HF=g.hwHF()))
0113     print("-----------------------------")
0114     print("-----------------------------")
0115     print("c + enter to continue")
0116     import pdb;pdb.set_trace()
0117
0118 ###############################
0119
0120 #########Histograms#############
0121 histos={}
0122 histos['fw']={}
0123 histos['fw']['pt1']=ROOT.TH1D("fw_pt1","HW p_{T}",512,0,511)
0124 histos['fw']['eta1']=ROOT.TH1D("fw_eta1","HW #eta",256,-127,128)
0125 histos['fw']['phi1']=ROOT.TH1D("fw_phi1","HW #phi",256,-127,128)
0126 histos['fw']['HF1']=ROOT.TH1D("fw_HF1","HW HF",256,-127,128)
0127 histos['fw']['qual1']=ROOT.TH1D("fw_qual1","HW qual",16,0,16)
0128 histos['fw']['dxy1']=ROOT.TH1D("fw_dxy1","HW DXY",4,0,4)
0129 histos['fw']['ptSTA1']=ROOT.TH1D("fw_ptSTA1","HW STA PT",256,0,255)
0130
0131 histos['fw']['pt2']=ROOT.TH1D("fw_pt2","HW p_{T}",512,0,511)
0132 histos['fw']['eta2']=ROOT.TH1D("fw_eta2","HW #eta",256,-127,128)
0133 histos['fw']['phi2']=ROOT.TH1D("fw_phi2","HW #phi",256,-127,128)
0134 histos['fw']['HF2']=ROOT.TH1D("fw_HF2","HW HF",256,-127,128)
0135 histos['fw']['qual2']=ROOT.TH1D("fw_qual2","HW qual",16,0,16)
0136 histos['fw']['dxy2']=ROOT.TH1D("fw_dxy2","HW DXY",4,0,4)
0137 histos['fw']['ptSTA2']=ROOT.TH1D("fw_ptSTA2","HW STA PT",256,0,255)
0138
0139 histos['fw']['pt3']=ROOT.TH1D("fw_pt3","HW p_{T}",512,0,511)
0140 histos['fw']['eta3']=ROOT.TH1D("fw_eta3","HW #eta",256,-127,128)
0141 histos['fw']['phi3']=ROOT.TH1D("fw_phi3","HW #phi",256,-127,128)
0142 histos['fw']['HF3']=ROOT.TH1D("fw_HF3","HW HF",256,-127,128)
0143 histos['fw']['qual3']=ROOT.TH1D("fw_qual3","HW qual",16,0,16)
0144 histos['fw']['dxy3']=ROOT.TH1D("fw_dxy3","HW DXY",4,0,4)
0145 histos['fw']['ptSTA3']=ROOT.TH1D("fw_ptSTA3","HW STA PT",256,0,255)
0146
0147
0148
0149 histos['emu']={}
0150
0151 histos['emu']['pt1']=ROOT.TH1D("emu_pt1","HW p_{T}",512,0,511)
0152 histos['emu']['eta1']=ROOT.TH1D("emu_eta1","HW #eta",256,-127,128)
0153 histos['emu']['phi1']=ROOT.TH1D("emu_phi1","HW #phi",256,-127,128)
0154 histos['emu']['HF1']=ROOT.TH1D("emu_HF1","HW HF",256,-127,128)
0155 histos['emu']['qual1']=ROOT.TH1D("emu_qual1","HW qual",16,0,16)
0156 histos['emu']['dxy1']=ROOT.TH1D("emu_dxy1","HW DXY",4,0,4)
0157 histos['emu']['ptSTA1']=ROOT.TH1D("emu_ptSTA1","HW STA PT",256,0,255)
0158
0159 histos['emu']['pt2']=ROOT.TH1D("emu_pt2","HW p_{T}",512,0,511)
0160 histos['emu']['eta2']=ROOT.TH1D("emu_eta2","HW #eta",256,-127,128)
0161 histos['emu']['phi2']=ROOT.TH1D("emu_phi2","HW #phi",256,-127,128)
0162 histos['emu']['HF2']=ROOT.TH1D("emu_HF2","HW HF",256,-127,128)
0163 histos['emu']['qual2']=ROOT.TH1D("emu_qual2","HW qual",16,0,16)
0164 histos['emu']['dxy2']=ROOT.TH1D("emu_dxy2","HW DXY",4,0,4)
0165 histos['emu']['ptSTA2']=ROOT.TH1D("emu_ptSTA2","HW STA PT",256,0,255)
0166
0167 histos['emu']['pt3']=ROOT.TH1D("emu_pt3","HW p_{T}",512,0,511)
0168 histos['emu']['eta3']=ROOT.TH1D("emu_eta3","HW #eta",256,-127,128)
0169 histos['emu']['phi3']=ROOT.TH1D("emu_phi3","HW #phi",256,-127,128)
0170 histos['emu']['HF3']=ROOT.TH1D("emu_HF3","HW HF",256,-127,128)
0171 histos['emu']['qual3']=ROOT.TH1D("emu_qual3","HW qual",16,0,16)
0172 histos['emu']['dxy3']=ROOT.TH1D("emu_dxy3","HW DXY",4,0,4)
0173 histos['emu']['ptSTA3']=ROOT.TH1D("emu_ptSTA3","HW STA PT",256,0,255)
0174
0175
0176 for key,histo in histos['fw'].iteritems():
0177     histo.Sumw2()
0178
0179
0180 def fill(info,mu):
0181     if len(mu)>0:
0182         info['pt1'].Fill(mu[0].hwPt())
0183         info['eta1'].Fill(mu[0].hwEta())
0184         info['phi1'].Fill(mu[0].hwPhi())
0185         info['HF1'].Fill(mu[0].hwHF())
0186         info['qual1'].Fill(mu[0].hwQual())
0187         info['dxy1'].Fill(mu[0].hwDXY())
0188         info['ptSTA1'].Fill(mu[0].hwPtUnconstrained())
0189     else:
0190         info['pt1'].Fill(0)
0191         info['eta1'].Fill(0)
0192         info['phi1'].Fill(0)
0193         info['HF1'].Fill(0)
0194         info['qual1'].Fill(0)
0195         info['dxy1'].Fill(0)
0196         info['ptSTA1'].Fill(0)
0197
0198     if len(mu)>1:
0199         info['pt2'].Fill(mu[1].hwPt())
0200         info['eta2'].Fill(mu[1].hwEta())
0201         info['phi2'].Fill(mu[1].hwPhi())
0202         info['HF2'].Fill(mu[1].hwHF())
0203         info['qual2'].Fill(mu[1].hwQual())
0204         info['dxy2'].Fill(mu[1].hwDXY())
0205         info['ptSTA2'].Fill(mu[1].hwPtUnconstrained())
0206     else:
0207         info['pt2'].Fill(0)
0208         info['eta2'].Fill(0)
0209         info['phi2'].Fill(0)
0210         info['HF2'].Fill(0)
0211         info['qual2'].Fill(0)
0212         info['dxy2'].Fill(0)
0213         info['ptSTA2'].Fill(0)
0214
0215     if len(mu)>2:
0216         info['pt3'].Fill(mu[2].hwPt())
0217         info['eta3'].Fill(mu[2].hwEta())
0218         info['phi3'].Fill(mu[2].hwPhi())
0219         info['HF3'].Fill(mu[2].hwHF())
0220         info['qual3'].Fill(mu[2].hwQual())
0221         info['dxy3'].Fill(mu[2].hwDXY())
0222         info['ptSTA3'].Fill(mu[2].hwPtUnconstrained())
0223     else:
0224         info['pt3'].Fill(0)
0225         info['eta3'].Fill(0)
0226         info['phi3'].Fill(0)
0227         info['HF3'].Fill(0)
0228         info['qual3'].Fill(0)
0229         info['dxy3'].Fill(0)
0230         info['ptSTA3'].Fill(0)
0231
0232
0233
0234
0235
0236
0237
0238 ##############################
0239
0240 BUNCHES=[0]
0241
0242
0243 events=Events([tag+'.root'])
0244 counter=-1
0245 for event in events:
0246     counter=counter+1
0247     #fetch stubs
0248     stubs=fetchStubsOLD(event,True)
0249     unpacker=fetchKMTF(event,100.0,'bmtfDigis:kBMTF')
0250     emulator=fetchKMTF(event,100.0,'simKBmtfDigis:BMTF')
0251
0252
0253     for processor in range(0,12):
0254         for bx in BUNCHES:
0255             emu=filter(lambda x: x.processor()==processor,emulator[bx])
0256             data=filter(lambda x: x.processor()==processor,unpacker[bx])
0257             if (len(emu)+len(data))>0:
0258
0259                 fill(histos['emu'],emu)
0260                 fill(histos['fw'],data)
0261 #                if len(emu)!=0 and len(data)==0:
0262 #                    log(event,counter,stubs,data,emu)
0263 #                    import pdb;pdb.set_trace()
0264
0265 f=ROOT.TFile("validationResults.root","RECREATE")
0266 for key,histo in histos['fw'].iteritems():
0267     histo.SetMarkerStyle(7)
0268     histo.Write()
0269 for key,histo in histos['emu'].iteritems():
0270     histo.SetLineColor(ROOT.kRed)
0271     histo.Write()
0272
0273
0274 #make fancy plots
0275 histonames=['pt1','eta1','phi1','HF1','qual1','dxy1','ptSTA1']
0276
0277 for h in histonames:
0278     c=ROOT.TCanvas(h)
0279     c.cd()
0280     histos['emu'][h].Draw("HIST")
0281     histos['emu'][h].GetXaxis().SetTitle(histos['emu'][h].GetTitle())
0282     histos['emu'][h].GetYaxis().SetTitle("events")
0283     histos['fw'][h].Draw("SAME")
0284     c.SetLogy()
0285     l=ROOT.TLegend(0.6,0.6,0.9,0.8)
0286     l.AddEntry(histos['emu'][h],"emulator","l")
0287     l.AddEntry(histos['fw'][h],"data","p")
0288     l.Draw()
0289     c.Write("plot_"+h)
0290
0291
0292
0293
0294
0295
0296 f.Close()