test/kalmanTools/validation.py

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