Phase2L1GMT/test/saStudy.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 def strAppend(s):
0010      return "root://cmsxrootd.fnal.gov/"+s
0011
0012 def fetchKMTF(event,tag,etamax=3.0):
0013      phiSeg2    = Handle  ('vector<l1t::SAMuon>')
0014      event.getByLabel(tag,phiSeg2)
0015      return list(filter(lambda x: abs(x.eta())<etamax,phiSeg2.product()))
0016
0017 def fetchStubs(event,tag):
0018      phiSeg2    = Handle  ('vector<l1t::MuonStub>')
0019      event.getByLabel(tag,phiSeg2)
0020      return phiSeg2.product()
0021
0022
0023 def fetchGEN(event,etaMax=3.0):
0024      genH  = Handle  ('vector<reco::GenParticle>')
0025      event.getByLabel('genParticles',genH)
0026      genMuons=list(filter(lambda x: abs(x.pdgId())==13 and x.status()==1 and abs(x.eta())<etaMax,genH.product()))
0027      return genMuons
0028
0029
0030 def deltaPhi( p1, p2):
0031      '''Computes delta phi, handling periodic limit conditions.'''
0032      res = p1 - p2
0033      while res > math.pi:
0034          res -= 2*math.pi
0035      while res < -math.pi:
0036          res += 2*math.pi
0037      return res
0038
0039 def deltaR( *args ):
0040      return math.sqrt( deltaR2(*args) )
0041
0042 def deltaR2( e1, p1, e2, p2):
0043      de = e1 - e2
0044      dp = deltaPhi(p1, p2)
0045      return de*de + dp*dp
0046
0047
0048 histoData={'effpt':{}, 'effeta':{},'geneta':{},'genphi':{},'effphi':{},'effptB':{},'effptO':{},'effptE':{}}
0049 histoData['genpt'] = ROOT.TH1D("genpt","genpt",20,0,100)
0050 histoData['genptB'] = ROOT.TH1D("genptB","genptB",20,0,100)
0051 histoData['genptO'] = ROOT.TH1D("genptO","genptO",20,0,100)
0052 histoData['genptE'] = ROOT.TH1D("genptE","genptE",20,0,100)
0053 rate = ROOT.TH1D("rate","rate",20,0,100)
0054 rateB = ROOT.TH1D("rateB","rateB",20,0,100)
0055 rateO = ROOT.TH1D("rateO","rateO",20,0,100)
0056 rateE = ROOT.TH1D("rateE","rateE",20,0,100)
0057
0058 thresholds=[0,1,3,5,15,20]
0059 for t in thresholds:
0060      histoData['effpt'][t] = ROOT.TH1D("effpt_"+str(t),"effpt_"+str(t),20,0,100)
0061      histoData['effptB'][t] = ROOT.TH1D("effptB_"+str(t),"effpt_"+str(t),20,0,100)
0062      histoData['effptO'][t] = ROOT.TH1D("effptO_"+str(t),"effpt_"+str(t),20,0,100)
0063      histoData['effptE'][t] = ROOT.TH1D("effptE_"+str(t),"effpt_"+str(t),20,0,100)
0064      histoData['effeta'][t] = ROOT.TH1D("effeta_"+str(t),"effeta_"+str(t),48,-2.4,2.4)
0065      histoData['effphi'][t] = ROOT.TH1D("effphi_"+str(t),"effphi_"+str(t),32,-math.pi,math.pi)
0066      histoData['geneta'][t] = ROOT.TH1D("geneta_"+str(t),"effeta_"+str(t),48,-2.4,2.4)
0067      histoData['genphi'][t] = ROOT.TH1D("genphi_"+str(t),"genphi_"+str(t),32,-math.pi,math.pi)
0068
0069
0070
0071
0072 etaLUT = ROOT.TH2D("etaLUT","etaLUT",256,0,256,128,0.0,1.0)
0073
0074 from samples200 import *
0075 toProcess = jpsi200
0076
0077 files=[]
0078
0079 for p in toProcess:
0080      files.append(strAppend(p))
0081
0082 #events=Events(['file:/uscmst1b_scratch/lpc1/3DayLifetime/bachtis/test.root'])
0083 events=Events(files)
0084
0085
0086
0087 displaced=0
0088
0089
0090
0091 BUNCHFACTOR=40000*2760.0/3564.0
0092
0093 counter=0;
0094 for event in events:
0095
0096     gen=fetchGEN(event,2.4)
0097     sa=fetchKMTF(event,'l1tSAMuonsGmt:prompt',2.5)
0098     if displaced==1:
0099          gen=fetchGEN(event,0.)
0100          sa=fetchKMTF(event,'l1tKMTFMuonsGmt:displaced',2.5)
0101
0102     highq=[]
0103
0104     #reject single stub muons
0105     for s in sa:
0106          if s.hwQual()>0:
0107               highq.append(s)
0108     sa=highq
0109
0110     maxPT = None
0111     maxPTB = None
0112     maxPTO = None
0113     maxPTE = None
0114
0115     for s in sa:
0116      if maxPT==None or s.pt()>maxPT:
0117           maxPT=s.pt()
0118      if (maxPTB==None or s.pt()>maxPTB) and abs(s.eta())<0.83:
0119           maxPTB=s.pt()
0120      if (maxPTO==None or s.pt()>maxPTO) and abs(s.eta())>0.83 and abs(s.eta())<1.2:
0121           maxPTO=s.pt()
0122      if (maxPTE==None or s.pt()>maxPTE) and abs(s.eta())>1.2:
0123           maxPTE=s.pt()
0124     if maxPT!=None:
0125          rate.Fill(maxPT)
0126     if maxPTB!=None:
0127          rateB.Fill(maxPTB)
0128     if maxPTO!=None:
0129          rateO.Fill(maxPTO)
0130     if maxPTE!=None:
0131          rateE.Fill(maxPTE)
0132
0133
0134
0135 #    print('---------------------NEW EVENT---------------------')
0136     if counter %10000 ==0:
0137          print(counter)
0138
0139 #    if counter==150000:
0140 #         break;
0141 #    print(counter)
0142 #    print('Generated muons')
0143     for m in gen:
0144
0145 #         print("gen pt=",m.pt())
0146          histoData['genpt'].Fill(m.pt())
0147          if abs(m.eta())<0.83:
0148               histoData['genptB'].Fill(m.pt())
0149          elif abs(m.eta())>0.83 and abs(m.eta())<1.2:
0150               histoData['genptO'].Fill(m.pt())
0151          else:
0152               histoData['genptE'].Fill(m.pt())
0153          for t in thresholds:
0154               if m.pt()>(float(t)+10.0):
0155                    histoData['geneta'][t].Fill(m.eta())
0156                    histoData['genphi'][t].Fill(m.phi())
0157          #fill the etaLUT
0158          for r in sa:
0159                if abs(deltaPhi(r.phi(),m.phi()))<0.3 and r.pt()>float(t) and r.hwQual()>0:
0160                     code=0;
0161                     for s in r.stubs():
0162                          code = code| ( (int(abs(s.etaRegion())+1))<<(2*(s.depthRegion()-1)))
0163                     etaLUT.Fill(code,abs(m.eta()))
0164
0165          for t in thresholds:
0166 #              print("Threshold {}".format(t))
0167 #              for s in sa:
0168 #                   print("muon pt={}".format(s.pt()))
0169               matched=False
0170               for r in sa:
0171 #                   print("SA over muon pt={}".format(r.pt()))
0172
0173                    if deltaR(r.eta(),r.phi(),m.eta(),m.phi())<0.5 and r.pt()>float(t) and r.hwQual()>0:
0174                         matched=True
0175                         break
0176 #              print("matched={}".format(matched))
0177               if matched:
0178                    histoData['effpt'][t].Fill(m.pt())
0179                    if abs(m.eta())<0.83:
0180                         histoData['effptB'][t].Fill(m.pt())
0181                    elif abs(m.eta())>0.83 and abs(m.eta())<1.2:
0182                         histoData['effptO'][t].Fill(m.pt())
0183                    else:
0184                         histoData['effptE'][t].Fill(m.pt())
0185                    if m.pt()>(t+10):
0186                         histoData['effeta'][t].Fill(m.eta())
0187                         histoData['effphi'][t].Fill(m.phi())
0188
0189     counter=counter+1
0190
0191 f=ROOT.TFile("saStudy_results.root","RECREATE")
0192 f.cd()
0193 #c = histoData['rate'].GetCumulative(False)
0194 #c.Scale(float(BUNCHFACTOR)/float(counter))
0195 #c.Write("rate")
0196 histoData['genpt'].Write("genpt")
0197 histoData['genptB'].Write("genptB")
0198 histoData['genptO'].Write("genptO")
0199 histoData['genptE'].Write("genptE")
0200 etaLUT.Write()
0201 rate.Write()
0202 rateB.Write()
0203 rateE.Write()
0204 rateO.Write()
0205
0206 c = rate.GetCumulative(False)
0207 c.Scale(float(BUNCHFACTOR)/float(counter))
0208 c.Write("rate")
0209 c = rateB.GetCumulative(False)
0210 c.Scale(float(BUNCHFACTOR)/float(counter))
0211 c.Write("rateBarrel")
0212 c = rateO.GetCumulative(False)
0213 c.Scale(float(BUNCHFACTOR)/float(counter))
0214 c.Write("rateOverlap")
0215 c = rateE.GetCumulative(False)
0216 c.Scale(float(BUNCHFACTOR)/float(counter))
0217 c.Write("rateEndcap")
0218
0219
0220 for t in thresholds:
0221     histoData['effpt'][t].Write("numpt_"+str(t))
0222     histoData['effeta'][t].Write("numeta_"+str(t))
0223     histoData['effphi'][t].Write("numphi_"+str(t))
0224     histoData['geneta'][t].Write("geneta_"+str(t))
0225     histoData['genphi'][t].Write("genphi_"+str(t))
0226
0227     g = ROOT.TGraphAsymmErrors(histoData['effpt'][t],histoData['genpt'])
0228     g.Write("eff_"+str(t))
0229     g = ROOT.TGraphAsymmErrors(histoData['effptB'][t],histoData['genptB'])
0230     g.Write("effB_"+str(t))
0231     g = ROOT.TGraphAsymmErrors(histoData['effptO'][t],histoData['genptO'])
0232     g.Write("effO_"+str(t))
0233     g = ROOT.TGraphAsymmErrors(histoData['effptE'][t],histoData['genptE'])
0234     g.Write("effE_"+str(t))
0235     g = ROOT.TGraphAsymmErrors(histoData['effeta'][t],histoData['geneta'][t])
0236     g.Write("effeta_"+str(t))
0237     g = ROOT.TGraphAsymmErrors(histoData['effphi'][t],histoData['genphi'][t])
0238     g.Write("effphi_"+str(t))
0239
0240 f.Close()