File indexing completed on 2023-03-17 11:12:38
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
0015
0016
0017
0018
0019
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
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
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
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
0262
0263
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
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()