Phase2L1GMT/test/gmtStudy.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 #0.005
0008 #idCut=[
0009 #20, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 48, 44, 44, 44, 60, 20, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 48, 20, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 60, 20, 24, 20, 44, 44, 44, 44, 44, 44, 44, 44, 20, 44, 44, 44, 60, 20, 24, 24, 24, 24, 44, 24, 24, 24, 20, 44, 40, 44, 44, 48, 44, 20, 44, 44, 48, 48, 48, 48, 60, 48, 48, 60, 60, 48, 60, 60, 60, 20, 44, 60, 60, 60, 60, 60, 60, 60, 48, 60, 60, 44, 60, 60, 60, 20, 40, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 16, 20, 60, 60, 60, 60, 60, 60, 60, 60, 48, 60, 44, 60, 60, 60, 16, 20, 44, 60, 60, 60, 60, 60, 48, 60, 60, 60, 60, 60, 60, 60, 20, 20, 40, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 48, 60, 20, 16, 20, 44, 48, 44, 44, 44, 60, 60, 60, 60, 48, 48, 48, 60, 20, 20, 20, 60, 44, 60, 60, 60, 60, 44, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60]
0010
0011 #0.006
0012
0013 idCut=[
0014 36, 56, 64, 64, 32, 32, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #0
0015 28, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 28, 28, 28, 28,       #16
0016 28, 52, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #32
0017 28, 52, 52, 52, 56, 56, 56, 56, 56, 56, 56, 28, 56, 56, 56, 56,       #48
0018 28, 36, 36, 36, 36, 36, 36, 36, 36, 36, 56, 56, 36, 56, 36, 36,       #64
0019 28, 56, 60, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #80
0020 28, 92, 120, 120, 92, 92, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,     #96
0021 28, 120, 152, 128, 100, 124, 80, 80, 64, 64, 64, 64, 64, 64, 64,64,   #112
0022 28, 88, 92, 64, 92, 60, 60, 60, 92, 60, 60, 60, 60, 60, 60, 60,       #128
0023 36, 88, 88, 64, 64, 60, 60, 60, 60, 60, 60, 60, 92, 92, 92,92,        #144
0024 28, 64, 92, 88, 88, 88, 88, 88, 88, 88, 88, 88, 120, 120, 120,120,    #160
0025 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,        #176
0026 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,        #192
0027 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,        #208
0028 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,        #224
0029 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64]       #240
0030
0031 idCut=[
0032 40, 56, 64, 88, 64, 32, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #0
0033 28, 56, 56, 54, 36, 36, 36, 36, 36, 36, 36, 36, 32, 28, 28, 28,       #16
0034 28, 52, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #32
0035 28, 56, 52, 52, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #48
0036 28, 52, 36, 36, 36, 36, 36, 36, 36, 36, 56, 56, 36, 56, 36, 36,       #64
0037 28, 60, 60, 60, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,       #80
0038 28, 120, 120, 120, 92, 92, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,     #96
0039 28, 120, 152, 128, 100, 124, 80, 80, 64, 100, 64, 64, 64, 64, 64,64,   #112
0040 28, 88, 92, 64, 92, 60, 60, 60, 92, 60, 60, 60, 60, 60, 60, 60,       #128
0041 36, 88, 88, 64, 64, 60, 60, 60, 60, 60, 60, 60, 92, 92, 92,92,        #144
0042 28, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88, 88,88,    #160
0043 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,        #176
0044 28, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92, 92,92,        #192
0045 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,        #208
0046 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64,64,        #224
0047 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64]       #240
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 #verbose=False
0064 #tag='singleMuonOfficial'
0065 #isData=False
0066 #tag='signal'
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 #events=Events(['/uscmst1b_scratch/lpc1/3DayLifetime/bachtis/MinBias.root'
0076 #          ])
0077 events=Events(map(strAppend,dy200))
0078
0079 #events=Events(minbias)
0080 #events=Events(['reprocess.root'])
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 #    tps = fetchTPS(event,'gmtMuons')
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          #         import pdb;pdb.set_trace()
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 #                 import pdb;pdb.set_trace()
0185
0186
0187     if counter %1000==0:
0188          print(counter)
0189     if verbose:
0190          print ("EVENT:",counter)
0191
0192     #efficiencies -loop on gen
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 #        if len(tpsMatched)==0 and g.pt()>10:
0206 #             print("Not Matched ",g.pt(),g.eta(),len(tps))
0207 #             import pdb;pdb.set_trace()
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