File indexing completed on 2023-03-17 11:13:09
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
0011
0012
0013 idCut=[
0014 36, 56, 64, 64, 32, 32, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0015 28, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 28, 28, 28, 28,
0016 28, 52, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0017 28, 52, 52, 52, 56, 56, 56, 56, 56, 56, 56, 28, 56, 56, 56, 56,
0018 28, 36, 36, 36, 36, 36, 36, 36, 36, 36, 56, 56, 36, 56, 36, 36,
0019 28, 56, 60, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0020 28, 92, 120, 120, 92, 92, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,
0021 28, 120, 152, 128, 100, 124, 80, 80, 64, 64, 64, 64, 64, 64, 64,64,
0022 28, 88, 92, 64, 92, 60, 60, 60, 92, 60, 60, 60, 60, 60, 60, 60,
0023 36, 88, 88, 64, 64, 60, 60, 60, 60, 60, 60, 60, 92, 92, 92,92,
0024 28, 64, 92, 88, 88, 88, 88, 88, 88, 88, 88, 88, 120, 120, 120,120,
0025 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,
0026 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,
0027 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,
0028 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,
0029 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64]
0030
0031 idCut=[
0032 40, 56, 64, 88, 64, 32, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0033 28, 56, 56, 54, 36, 36, 36, 36, 36, 36, 36, 36, 32, 28, 28, 28,
0034 28, 52, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0035 28, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0036 28, 52, 36, 36, 36, 36, 36, 36, 36, 36, 56, 56, 36, 56, 36, 36,
0037 28, 60, 60, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
0038 28, 120, 120, 120, 92, 92, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,
0039 28, 120, 152, 128, 100, 124, 80, 80, 64, 100, 64, 64, 64, 64, 64,64,
0040 28, 88, 92, 64, 92, 60, 60, 60, 92, 60, 60, 60, 60, 60, 60, 60,
0041 36, 88, 88, 64, 64, 60, 60, 60, 60, 60, 60, 60, 92, 92, 92,92,
0042 28, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88,88,
0043 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,
0044 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,
0045 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,
0046 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,
0047 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64]
0048
0049
0050 def ID(eta,pt,quality):
0051 p=pt
0052 if p>4095:
0053 p=4095
0054 p=p/256
0055 e=int(abs(eta))
0056 e=e/256
0057 addr=(e<<4) | p
0058 if (quality-min(64,pt/32))>=(idCut[addr]):
0059 return True
0060 else:
0061 return False
0062
0063
0064
0065
0066
0067
0068 def strAppend(s):
0069 return "root://cmsxrootd.fnal.gov/"+s
0070
0071 from dy200 import *
0072 from jpsi import *
0073 from minbias import *
0074 tag='/uscmst1b_scratch/lpc1/3DayLifetime/bachtis/DY200'
0075
0076
0077 events=Events(map(strAppend,dy200))
0078
0079
0080
0081
0082 def fetchTracks(event):
0083 trackHandle = Handle('vector<TTTrack<edm::Ref<edm::DetSetVector<Phase2TrackerDigi>,Phase2TrackerDigi,edm::refhelper::FindForDetSetVector<Phase2TrackerDigi> > > >')
0084 event.getByLabel("l1tTTTracksFromTrackletEmulation:Level1TTTracks",trackHandle)
0085 return trackHandle.product()
0086
0087
0088 def fetchTPS(event,tag,etamax=3.0):
0089 phiSeg2 = Handle ('vector<l1t::TrackerMuon>')
0090 event.getByLabel(tag,phiSeg2)
0091 return filter(lambda x: abs(x.eta())<etamax,phiSeg2.product())
0092
0093
0094 def fetchStubs(event,tag):
0095 phiSeg2 = Handle ('vector<l1t::MuonStub>')
0096 event.getByLabel(tag,phiSeg2)
0097 return phiSeg2.product()
0098
0099
0100 def fetchGEN(event,etaMax=3.0):
0101 genH = Handle ('vector<reco::GenParticle>')
0102 event.getByLabel('genParticles',genH)
0103 genMuons=filter(lambda x: abs(x.pdgId())==13 and x.status()==1 and abs(x.eta())<etaMax,genH.product())
0104 return genMuons
0105
0106
0107 def deltaPhi( p1, p2):
0108 '''Computes delta phi, handling periodic limit conditions.'''
0109 res = p1 - p2
0110 while res > math.pi:
0111 res -= 2*math.pi
0112 while res < -math.pi:
0113 res += 2*math.pi
0114 return res
0115
0116 def deltaR( *args ):
0117 return math.sqrt( deltaR2(*args) )
0118
0119 def deltaR2( e1, p1, e2, p2):
0120 de = e1 - e2
0121 dp = deltaPhi(p1, p2)
0122 return de*de + dp*dp
0123
0124
0125
0126
0127 quality = ROOT.TH3D("quality","",16,0,4096,16,0,4096,128,0,512)
0128 qualityAll = ROOT.TH3D("qualityAll","",16,0,4096,16,0,4096,128,0,512)
0129 histoData={'effpt':{}, 'effeta':{},'geneta':{},'genphi':{},'effphi':{}}
0130 histoData['genpt'] = ROOT.TH1D("genpt","genpt",50,0,50)
0131
0132 thresholds=[1,3,5,15,20]
0133 for t in thresholds:
0134 histoData['effpt'][t] = ROOT.TH1D("effpt_"+str(t),"effpt_"+str(t),50,0,50)
0135 histoData['effeta'][t] = ROOT.TH1D("effeta_"+str(t),"effeta_"+str(t),48,-2.4,2.4)
0136 histoData['effphi'][t] = ROOT.TH1D("effphi_"+str(t),"effphi_"+str(t),32,-math.pi,math.pi)
0137 histoData['geneta'][t] = ROOT.TH1D("geneta_"+str(t),"effeta_"+str(t),48,-2.4,2.4)
0138 histoData['genphi'][t] = ROOT.TH1D("genphi_"+str(t),"genphi_"+str(t),32,-math.pi,math.pi)
0139
0140
0141
0142
0143 histoData['rate']= ROOT.TH1D("rate","rate",50,0,100)
0144 histoData['rate20eta']= ROOT.TH1D("rate20eta","rate",8,0,8)
0145
0146
0147 BUNCHFACTOR=40000*2760.0/3564.0
0148
0149 QUALITYCUT=0
0150 verbose=0
0151 counter=-1
0152 for event in events:
0153 counter=counter+1
0154 if counter==100000:
0155 break;
0156 gen=fetchGEN(event,2.4)
0157 stubs = fetchStubs(event,'l1tGMTStubs')
0158 tps = filter(lambda x: ID(x.hwEta(),x.hwPt(),x.hwQual()),fetchTPS(event,'l1tGMTMuons'))
0159
0160 tracks=fetchTracks(event)
0161
0162 if verbose:
0163 for t in sorted(tracks,key=lambda x: x.momentum().transverse(),reverse=True):
0164 print("Track ",t.eta(),t.phi(),t.momentum().transverse(),t.getStubRefs().size())
0165 for t in sorted(tps,key=lambda x: x.pt(),reverse=True):
0166 print("Tracker Muon pt={pt} eta={eta} phi={phi} qual={qual}".format(pt=t.pt(),eta=t.eta(),phi=t.phi(),qual=t.hwQual))
0167
0168
0169 for s in stubs:
0170 print(" All Stub depth={depth} eta={eta1},{eta1I},{eta2},{eta2I} phi={phi1},{phi1I},{phi2},{phi2I} qualities={etaQ},{phiQ}".format(depth=s.depthRegion(),eta1=s.offline_eta1(),eta1I=s.eta1(),eta2=s.offline_eta2(),eta2I=s.eta2(),phi1=s.offline_coord1(),phi1I=s.coord1(),phi2=s.offline_coord2(),phi2I=s.coord2(),etaQ=s.etaQuality(),phiQ=s.quality()))
0171
0172
0173 if len(tps)>0:
0174 tpsInEta =filter(lambda x: abs(x.eta())<2.4,tps)
0175 if len(tpsInEta)>0:
0176 best=max(tpsInEta,key=lambda x: x.pt())
0177 qualityAll.Fill(abs(best.hwEta()),best.hwPt(),best.hwQual()-min(63,best.hwPt()/32))
0178 pt=best.pt()
0179 if pt>99.9:
0180 pt=99.9
0181 histoData['rate'].Fill(pt)
0182 if pt>20:
0183 histoData['rate20eta'].Fill(abs(best.hwEta())/512)
0184
0185
0186
0187 if counter %1000==0:
0188 print(counter)
0189 if verbose:
0190 print ("EVENT:",counter)
0191
0192
0193 for g in gen:
0194 histoData['genpt'].Fill(g.pt())
0195 for thres in thresholds:
0196 if g.pt()>(thres+2.0):
0197 histoData['geneta'][thres].Fill(g.eta())
0198 histoData['genphi'][thres].Fill(g.phi())
0199
0200 if verbose:
0201 print("Gen Muon pt={pt} eta={eta} phi={phi} charge={charge}".format(pt=g.pt(),eta=g.eta(),phi=g.phi(),charge=g.charge()))
0202
0203 foundMatch=False
0204 tpsMatched = sorted(filter(lambda x: deltaR(g.eta(),g.phi(),x.eta(),x.phi())<0.3, tps),key=lambda x:x.pt(),reverse=True)
0205
0206
0207
0208 if len(tpsMatched)>0:
0209 foundMatch=True
0210
0211 if verbose:
0212 for t in tpsMatched:
0213 print(" -> matched track ->Tracker Muon pt={pt} eta={eta} phi={phi} qual={qual}".format(pt=t.pt(),eta=t.eta(),phi=t.phi(),qual=t.hwQual()))
0214 best=tpsMatched[0]
0215 quality.Fill(abs(best.hwEta()),best.hwPt(),best.hwQual()-min(63,best.hwPt()/32))
0216 for thres in thresholds:
0217 overThreshold = filter(lambda x: x.pt()>thres,tpsMatched)
0218 if len(overThreshold)>0:
0219 histoData['effpt'][thres].Fill(g.pt())
0220 if g.pt()>(thres+2.0):
0221 histoData['effeta'][thres].Fill(g.eta())
0222 histoData['effphi'][thres].Fill(g.phi())
0223
0224
0225
0226
0227
0228
0229
0230 stubsMatched = filter(lambda x: deltaR(g.eta(),g.phi(),x.offline_eta1(),x.offline_coord1())<0.3, stubs)
0231 if verbose:
0232 for s in stubsMatched:
0233 print(" -> matched stub -> Muon Stub eta={eta1},{eta2} phi={phi1},{phi2} qualities={etaQ},{phiQ}".format(eta1=s.offline_eta1(),eta2=s.offline_eta2(),phi1=s.offline_coord1(),phi2=s.offline_coord2(),etaQ=s.etaQuality(),phiQ=s.quality()))
0234
0235
0236
0237 f=ROOT.TFile("tpsResults_output.root","RECREATE")
0238 f.cd()
0239 quality.Write()
0240 qualityAll.Write()
0241
0242 c = histoData['rate'].GetCumulative(False)
0243 c.Scale(float(BUNCHFACTOR)/float(counter))
0244 c.Write("rate")
0245
0246 for t in thresholds:
0247 g = ROOT.TGraphAsymmErrors(histoData['effpt'][t],histoData['genpt'])
0248 g.Write("eff_"+str(t))
0249 g = ROOT.TGraphAsymmErrors(histoData['effeta'][t],histoData['geneta'][t])
0250 g.Write("effeta_"+str(t))
0251 g = ROOT.TGraphAsymmErrors(histoData['effphi'][t],histoData['genphi'][t])
0252 g.Write("effphi_"+str(t))
0253
0254 f.Close()
0255