test/kalmanTools/kmtfAnalysis.py

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 #verbose=False
0011 #tag='singleMuonOfficial'
0012 #isData=False
0013
0014 tag='zerobias'
0015 #tag='zskim'
0016 isData=True
0017
0018
0019
0020
0021
0022 ##A class to keep BMTF data
0023 class BMTFMuon:
0024     def __init__(self,mu,pt,eta,phi):
0025         self.muon=mu
0026         self.p4 = ROOT.reco.Candidate.PolarLorentzVector(pt,eta,phi,0.105)
0027
0028     def quality(self):
0029         return self.muon.hwQual()
0030
0031     def phiINT(self):
0032         return self.muon.hwPhi()
0033
0034     def processor(self):
0035         return self.muon.processor()
0036
0037     def hasFineEta(self):
0038         return self.muon.hwHF()
0039
0040     def ptUnconstrained(self):
0041         return self.muon.hwPtUnconstrained()
0042
0043     def dxy(self):
0044         return self.muon.hwDXY()
0045
0046     def charge(self):
0047         if self.muon.hwSign()>0:
0048             return -1
0049         else:
0050             return +1
0051
0052     def __getattr__(self, name):
0053         return getattr(self.p4,name)
0054
0055
0056 ###Common methods############
0057
0058 def getQual(track):
0059     q=0
0060     for stub in track.stubs():
0061         q+=stub.quality()
0062     return q;
0063
0064
0065 def fetchTP(event,etaMax=0.83):
0066     tH  = Handle  ('trigger::TriggerEvent')
0067     mH  = Handle  ('std::vector<reco::Muon>')
0068
0069
0070     event.getByLabel('hltTriggerSummaryAOD','','HLT',tH)
0071     event.getByLabel('muons',mH)
0072
0073     muons=filter( lambda x: x.passed(ROOT.reco.Muon.CutBasedIdMediumPrompt) and x.numberOfMatches()>1 and x.pt()>10.0 and abs(x.eta())<2.4 and x.isolationR03().sumPt/x.pt()<0.15,mH.product())
0074     if len(muons)<2:
0075         return []
0076     trigger=tH.product()
0077 #    for f in range(0,trigger.sizeFilters()):
0078 #        print(f,trigger.filterLabel(f))
0079 #    import pdb;pdb.set_trace()
0080     obj = trigger.getObjects()
0081     index = trigger.filterIndex(ROOT.edm.InputTag("hltL3fL1sMu22Or25L1f0L2f10QL3Filtered27Q::HLT"))
0082     if index==trigger.sizeFilters():
0083         return []
0084     keys = trigger.filterKeys(index)
0085     hlt=[]
0086     for key in keys:
0087         hlt.append(obj[key])
0088     if len(hlt)==0:
0089         return []
0090
0091     tags =[]
0092     for x in muons:
0093         for triggered in hlt:
0094             if deltaR(x.eta(),x.phi(),triggered.eta(),triggered.phi())<0.3 and x.pt()>25:
0095                 tags.append(x)
0096                 break
0097     if len(tags)==0:
0098         return []
0099
0100
0101     probes=[]
0102     for mu in muons:
0103         isProbe=False
0104         for tag in tags:
0105             if deltaR(mu.eta(),mu.phi(),tag.eta(),tag.phi())>1.0:
0106                 if abs(mu.eta())<etaMax:
0107                     isProbe=True
0108                     break
0109         if isProbe:
0110             probes.append(mu)
0111
0112
0113     return probes
0114
0115
0116
0117 def fetchStubs(event,ontime=True):
0118     phiSeg2    = Handle  ('std::vector<L1MuKBMTCombinedStub>')
0119     event.getByLabel('simKBmtfStubs',phiSeg2)
0120     if ontime:
0121         filtered=filter(lambda x: x.bxNum()==0, phiSeg2.product())
0122         return filtered
0123     else:
0124         return phiSeg2.product()
0125
0126 def fetchStubsOLD(event,ontime=False,isData=True):
0127     phiSeg    = Handle  ('L1MuDTChambPhContainer')
0128     if not isData:
0129         event.getByLabel('simTwinMuxDigis',phiSeg)
0130     else:
0131         event.getByLabel('BMTFStage2Digis',phiSeg)
0132     if ontime:
0133         filtered=filter(lambda x: x.bxNum()==0, phiSeg.product().getContainer())
0134         return filtered
0135     else:
0136         return phiSeg.product().getContainer()
0137
0138 def fetchStubsOLDTheta(event,isData=True):
0139     phiSeg    = Handle  ('L1MuDTChambThContainer')
0140     if not isData:
0141         event.getByLabel('dtTriggerPrimitiveDigis',phiSeg)
0142     else:
0143         event.getByLabel('BMTFStage2Digis',phiSeg)
0144     return phiSeg.product().getContainer()
0145
0146
0147 def fetchGEN(event,etaMax=1.2):
0148     genH  = Handle  ('vector<reco::GenParticle>')
0149     event.getByLabel('genParticles',genH)
0150     genMuons=filter(lambda x: abs(x.pdgId())==13 and x.status()==1 and abs(x.eta())<etaMax,genH.product())
0151     return genMuons
0152
0153
0154 def fetchRECO(event,etaMax=0.8):
0155     genH  = Handle  ('vector<reco::Muon>')
0156     event.getByLabel('muons',genH)
0157     genMuons=filter(lambda x: x.passed(ROOT.reco.Muon.CutBasedIdMediumPrompt) and x.pt()>10.0 and abs(x.eta())<2.4 and x.isolationR03().sumPt/x.pt()<0.15 and abs(x.eta())<etaMax,genH.product())
0158     return genMuons
0159
0160
0161 def globalBMTFPhi(muon,calib=None):
0162     temp=muon.processor()*48+muon.hwPhi()
0163     temp=temp*2*math.pi/576.0-math.pi*15.0/180.0;
0164     if temp>math.pi:
0165         temp=temp-2*math.pi;
0166
0167     if calib=='KMTF':
0168         K=1.0/muon.hwPt()
0169         if muon.hwSign()>0:
0170             K=-1.0/muon.hwPt()
0171 #        return temp+3.752*K/(1.0-2.405*K*K/abs(K))
0172         return temp+5.740*K
0173
0174     if calib=='BMTF':
0175         K=1.0/muon.hwPt()
0176         if muon.hwSign()>0:
0177             K=-1.0/muon.hwPt()
0178         return temp+5.740*K
0179
0180     return temp;
0181
0182
0183 def rawPhi(muon):
0184     temp=muon.processor()*48+muon.hwPhi()
0185     temp=temp*2*math.pi/576.0-math.pi*15.0/180.0;
0186     if temp>math.pi:
0187         temp=temp-2*math.pi;
0188     return temp;
0189
0190
0191 def fetchBMTF(event,isData,etaMax=1.2):
0192     bmtfH  = Handle  ('BXVector<l1t::RegionalMuonCand>')
0193     if isData:
0194         event.getByLabel('BMTFStage2Digis:BMTF',bmtfH)
0195     else:
0196         event.getByLabel('simBmtfDigis:BMTF',bmtfH)
0197
0198     bmtf=bmtfH.product()
0199     bmtfMuons=[]
0200     for bx in [0]:
0201         for j in range(0,bmtf.size(bx)):
0202             mu = bmtf.at(bx,j)
0203             pt = mu.hwPt()*0.5
0204             K=1.0/pt
0205             K=1.181*K/(1+0.4867*K)
0206             pt=1.0/K
0207             ####
0208             phi=globalBMTFPhi(mu,'BMTF')
0209             rawP = rawPhi(mu)
0210             eta = mu.hwEta()*0.010875
0211             if abs(eta)<=etaMax:
0212                 b = BMTFMuon(mu,pt,eta,phi)
0213                 b.rawPhi=rawP
0214                 bmtfMuons.append(b)
0215     return sorted(bmtfMuons,key=lambda x: x.pt(),reverse=True)
0216
0217
0218 def fetchKMTFNew(event,etaMax=1.2,saturate=True):
0219     kbmtfH  = Handle  ('BXVector<l1t::RegionalMuonCand>')
0220     event.getByLabel('simKBmtfDigis:BMTF',kbmtfH)
0221     kbmtf=kbmtfH.product()
0222     kbmtfMuons=[]
0223     for bx in [0]:
0224         for j in range(0,kbmtf.size(bx)):
0225             mu = kbmtf.at(bx,j)
0226             pt =mu.hwPt()*0.5
0227             K=1.0/pt
0228             K=1.14*K
0229             pt=1.0/K
0230             if pt>140.0 and saturate:
0231                 pt=140.0
0232             phi=globalBMTFPhi(mu,'KMTF')
0233             rawP = rawPhi(mu)
0234             eta = mu.hwEta()*0.010875
0235             if abs(eta)<=etaMax:
0236                 b = BMTFMuon(mu,pt,eta,phi)
0237                 b.rawPhi=rawP
0238                 kbmtfMuons.append(b)
0239     return sorted(kbmtfMuons,key=lambda x: x.pt(),reverse=True)
0240
0241
0242 def qPTInt(qPT,bits=14):
0243     lsb = lsBIT(bits)
0244     floatbinary = int(math.floor(abs(qPT)/lsb))
0245     return int((qPT/abs(qPT))*floatbinary)
0246
0247 def lsBIT(bits=14):
0248     maximum=1.25
0249     lsb = 1.25/pow(2,bits-1)
0250     return lsb
0251
0252 #import pdb;pdb.set_trace()
0253
0254 def fetchKMTF(event,etaMax=0.83,patterns=[],comps=[],chis=[],pts=[]):
0255     kmtfH  = Handle('BXVector<L1MuKBMTrack>')
0256     event.getByLabel('simKBmtfDigis',kmtfH)
0257     kmtf=kmtfH.product()
0258     out=[]
0259     for bx in [0]:
0260         for j in range(0,kmtf.size(bx)):
0261             mu =  kmtf.at(bx,j)
0262             if abs(mu.eta())<etaMax:
0263                 veto=False
0264                 for pattern,comp,chi,pt in zip(patterns,comps,chis,pts):
0265                     if mu.hitPattern()==p and mu.pt()>pt and (mu.trackCompatibility()>comp or mu.approxChi2()>chi):
0266                         veto=True
0267                         break;
0268                 if not veto:
0269                     out.append(mu)
0270     return sorted(out,key=lambda x: x.pt(),reverse=True)
0271
0272 def curvResidual(a,b,factor=1.0):
0273     return (a.charge()/a.pt()-factor*b.charge()/b.pt())*b.pt()/(factor*b.charge())
0274
0275 def ptResidual(a,b):
0276     return (a.pt()-b.pt())/b.pt()
0277
0278 def curvResidualSTA(a,b):
0279     return (a.charge()/a.ptUnconstrained()-b.charge()/b.pt())*b.pt()/b.charge()
0280
0281
0282
0283
0284 def deltaPhi( p1, p2):
0285     '''Computes delta phi, handling periodic limit conditions.'''
0286     res = p1 - p2
0287     while res > math.pi:
0288         res -= 2*math.pi
0289     while res < -math.pi:
0290         res += 2*math.pi
0291     return res
0292
0293 def deltaR( *args ):
0294     return math.sqrt( deltaR2(*args) )
0295
0296 def deltaR2( e1, p1, e2, p2):
0297     de = e1 - e2
0298     dp = deltaPhi(p1, p2)
0299     return de*de + dp*dp
0300
0301
0302 def log(counter,mystubs,gen,kmtfFull,kmtf,bmtf):
0303     print("--------EVENT"+str(counter)+"------------")
0304     print("-----------------------------")
0305     print("-----------------------------")
0306     print('Stubs:')
0307     for stub in mystubs:
0308         print('wheel={w} sector={sc} station={st} high/low={ts} phi={phi} phiB={phiB} qual={qual} BX={BX} eta1={eta1} eta2={eta2}'.format(w=stub.whNum(),sc=stub.scNum(),st=stub.stNum(),ts=stub.tag(),phi=stub.phi(),phiB=stub.phiB(),qual=stub.quality(),BX=stub.bxNum(),eta1=stub.eta1(),eta2=stub.eta2()))
0309     print('Gen muons:')
0310     for g in gen:
0311         print("Generated muon charge={q} pt={pt} eta={eta} phi={phi}".format(q=g.charge(),pt=g.pt(),eta=g.eta(),phi=g.phi()))
0312     print('BMTF:')
0313     for g in bmtf :
0314         print("BMTF sector={sector} charge={q} pt={pt} eta={eta} phi={phi} qual={qual} dxy={dxy} pt2={pt2} hasFineEta={HF} rawPhi={hwPHI}".format(sector=g.processor(),q=g.charge(),pt=g.pt(),eta=g.eta(),phi=g.phi(),qual=g.quality(),dxy=g.dxy(),pt2=g.ptUnconstrained(),HF=g.hasFineEta(),hwPHI=g.phiINT()))
0315     print('KMTF:')
0316     for g in kmtf :
0317         print("KMTF sector={sector} charge={q} pt={pt} eta={eta} phi={phi} qual={qual} dxy={dxy} pt2={pt2} hasFineEta={HF} rawPhi={hwPHI}".format(sector=g.processor(),q=g.charge(),pt=g.pt(),eta=g.eta(),phi=g.phi(),qual=g.quality(),dxy=g.dxy(),pt2=g.ptUnconstrained(),HF=g.hasFineEta(),hwPHI=g.phiINT()))
0318     print('KMTF Full:')
0319     for g in kmtfFull :
0320         print("KMTF charge={q} pt={pt} ptSTA={PTSTA} eta={eta} phi={phi} pattern={pattern} chi={chi1} comp={comp}".format(q=g.charge(),pt=g.pt(),PTSTA=g.ptUnconstrained(),eta=g.eta(),phi=g.phi(),pattern=g.hitPattern(),chi1=g.approxChi2(),comp=g.trackCompatibility() ))
0321
0322
0323
0324
0325     print("-----------------------------")
0326     print("-----------------------------")
0327     print("c + enter to continue")
0328     import pdb;pdb.set_trace()
0329
0330 #########Histograms#############
0331 resKMTF = ROOT.TH2D("resKMTF","resKF",50,3,103,60,-2,2)
0332
0333
0334 resKMTFTrack={}
0335 for i in [0,3,5,6,7,9,10,11,12,13,14,15]:
0336     resKMTFTrack[i]=ROOT.TH2D("resKMTFTrack_"+str(i),"resKF",70,3,143,60,-2,2)
0337
0338
0339
0340 chiMatched={}
0341 trackComp={}
0342 trackCompAll={}
0343 chiAll={}
0344
0345 for i in [0,3,5,6,7,9,10,11,12,13,14,15]:
0346     chiMatched[i] = ROOT.TH2D("chiMatched_"+str(i),"resKF",32,0,1024,64,0,256)
0347     trackComp[i] = ROOT.TH2D("trackComp_"+str(i),"resKF",32,0,1024,100,0,100)
0348     trackCompAll[i] = ROOT.TH2D("trackCompAll_"+str(i),"resKF",20,0,1024,100,0,100)
0349     chiAll[i] = ROOT.TH2D("chiAll_"+str(i),"resKF",32,0,1024,64,0,256)
0350
0351
0352
0353 quality = ROOT.TH1D("quality","resKF",24,0,24)
0354
0355 resKMTFEta = ROOT.TH2D("resKMTFEta","resKF",8,0,0.8,60,-2,2)
0356
0357 resSTAKMTF = ROOT.TH2D("resSTAKMTF","resKF",100,0,100,100,-8,8)
0358 resBMTF = ROOT.TH2D("resBMTF","resKF",50,3,103,60,-2,2)
0359 resBMTFEta = ROOT.TH2D("resBMTFEta","resKF",8,0,0.8,60,-2,2)
0360
0361 resPTKMTF = ROOT.TH2D("resPTKMTF","resKF",100,0,100,60,-2,2)
0362 resPTBMTF = ROOT.TH2D("resPTBMTF","resKF",100,0,100,60,-2,2)
0363 resEtaKMTF = ROOT.TH2D("resEtaKMTF","resKF",5,-1.2,1.2,50,-20.5*0.010875,20.5*0.010875)
0364 resEtaBMTF = ROOT.TH2D("resEtaBMTF","resKF",5,-1.2,1.2,50,-20.5*0.010875,20.5*0.010875)
0365
0366
0367 resPhiKMTF = ROOT.TH2D("resPhiKMTF","resKF",50,3,103,250,-0.5,0.5)
0368 phiCalibKMTF = ROOT.TH2D("phiCalibKMTF","resKF",100,-1.0/3.2,1.0/3.2,101,-0.5,0.5)
0369 phiCalibBMTF = ROOT.TH2D("phiCalibBMTF","resKF",100,-1.0/3.2,1.0/3.2,101,-0.5,0.5)
0370
0371 resPhiBMTF = ROOT.TH2D("resPhiBMTF","resKF",50,3,103,250,-0.5,0.5)
0372 resRBMTF = ROOT.TH1D("resRBMTF","resKF",250,0,8)
0373 resRKMTF = ROOT.TH1D("resRKMTF","resKF",250,0,8)
0374
0375 fineEtaKMTF = ROOT.TH1D("fineEtaKMTF","fineEtaKMTF",2,0,2)
0376 fineEtaBMTF = ROOT.TH1D("fineEtaBMTF","fineEtaBMTF",2,0,2)
0377
0378
0379
0380 genPt=ROOT.TH1F("genPt","genPt",50,0,100)
0381
0382 PTThresholds=[0,5,7,10,15,25,30]
0383 genPtKMTF={}
0384 genPtBMTF={}
0385 genEtaKMTF={}
0386 genEtaBMTF={}
0387
0388 etaArr=[]
0389 for i in range(0,21):
0390     etaArr.append(-0.833+ 2*0.833*i/20.0)
0391
0392
0393 for p in PTThresholds:
0394     genPtKMTF[p]=ROOT.TH1F("genPtKMTF_"+str(p),"genPt",50,0,100)
0395     genPtBMTF[p]=ROOT.TH1F("genPtBMTF_"+str(p),"genPt",50,0,100)
0396     genEtaKMTF[p]=ROOT.TH1F("genEtaKMTF"+str(p),"genPt",len(etaArr)-1,array('f',etaArr))
0397     genEtaBMTF[p]=ROOT.TH1F("genEtaBMTF"+str(p),"genEta",len(etaArr)-1,array('f',etaArr))
0398
0399
0400
0401 kfCalibPlus={}
0402 kfCalibMinus={}
0403 ratePerTrack={}
0404 for track in [3,5,6,7,9,10,11,12,13,14,15]:
0405     kfCalibPlus[track] = ROOT.TH2D("kfCalibPlus_"+str(track),"resKF",560,3.2,143.2,100,0,10)
0406     kfCalibMinus[track] = ROOT.TH2D("kfCalibMinus_"+str(track),"resKF",560,3.2,143.2,100,0,10)
0407     ratePerTrack[track] = ROOT.TH1F("ratePerTrack_"+str(track),"rateKMTF",20,2.5,102.5)
0408
0409 kfCalib = ROOT.TH2D("kfCalib","resKF",560,3.2,143.2,100,0,10)
0410 bmtfCalib = ROOT.TH2D("bmtfCalib","resKF",280,0.,140,100,0,10)
0411
0412
0413
0414
0415
0416 #etaArr = [-0.8,-0.5,-0.3,-0.15,0.0,0.15,0.3,0.5,0.8]
0417
0418
0419
0420 genEta=ROOT.TH1F("genEta","genEta",len(etaArr)-1,array('f',etaArr))
0421
0422 genPhi=ROOT.TH1F("genPhi","genEta",50,-math.pi,math.pi)
0423 genPhiKMTF=ROOT.TH1F("genPhiKMTF","genPt",50,-math.pi,math.pi)
0424 genPhiBMTF=ROOT.TH1F("genPhiBMTF","genEta",50,-math.pi,math.pi)
0425
0426
0427 qualityKMTF = ROOT.TH1F("qualityKMTF","quality",16,0,16)
0428 qualityBMTF = ROOT.TH1F("qualityBMTF","quality",16,0,16)
0429
0430 dxyKMTF = ROOT.TH1F("dxyKMTF","chiBest",4,0,4)
0431
0432 etaBMTF = ROOT.TH1F("etaBMTF","rateBMTF",24,-1.2,1.2)
0433 etaKMTF = ROOT.TH1F("etaKMTF","rateKMTF",24,-1.2,1.2)
0434
0435 rateBMTF = ROOT.TH1F("rateBMTF","rateBMTF",50,0,50)
0436 rateKMTF = ROOT.TH1F("rateKMTF","rateKMTF",50,0,50)
0437
0438
0439
0440 rateBMTFp7 = ROOT.TH1F("rateBMTFp7","rateBMTF",20,2.5,102.5)
0441 rateKMTFp7 = ROOT.TH1F("rateKMTFp7","rateKMTF",20,2.5,102.5)
0442
0443 ##############################
0444
0445
0446 events=Events([tag+'.root'])
0447
0448
0449
0450
0451 counter=-1
0452 for event in events:
0453     counter=counter+1
0454     #fetch stubs
0455     stubs=[]
0456     kmtf=[]
0457     bmtf=[]
0458     stubsOLD=[]
0459
0460
0461     stubs = fetchStubs(event)
0462 #    stubsOLD = fetchStubsOLD(event,True,isData)
0463 #    stubsOLDTheta = fetchStubsOLDTheta(event,isData)
0464 #    stubs = []
0465     stubsOLD = []
0466 #    stubsOLDTheta = []
0467
0468
0469 #    if counter==1000:
0470 #        break;
0471
0472 #    print('OLD stubs')
0473 #    for s in stubsOLD:
0474 #        print(s.bxNum(),s.whNum(),s.scNum(),s.stNum(),s.phi(),s.phiB(),s.code())
0475
0476
0477 #    reco=fetchRECO(event,2.4)
0478     reco=[]
0479
0480     #fetch gen
0481     if isData:
0482         if tag=="zerobias" or tag=="output":
0483             gen=[]
0484         else:
0485             gen=fetchTP(event,0.83)
0486
0487     else:
0488         gen  = fetchGEN(event,0.83)
0489 #        bmtf = []
0490 #        bmtf = fetchBMTF(event,isData,1.5)
0491
0492     #fetch kalman (prompt)
0493     kmtf = fetchKMTFNew(event,1.5)
0494     bmtf = fetchBMTF(event,isData,1.5)
0495 #    bmtf=[]
0496
0497     #fetch detailed kalman
0498     kmtfFull = fetchKMTF(event,1.5)
0499
0500
0501 #    for g in gen:
0502 #        print('GEN', g.pt(),g.eta(),g.phi() )
0503
0504 #    for k in kmtf:
0505 #        print('L1', k.pt(),k.eta(),k.phi() )
0506
0507 #    log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0508
0509
0510     for track in kmtfFull:
0511         chiAll[track.hitPattern()].Fill(abs(track.curvatureAtVertex()),track.approxChi2())
0512         trackCompAll[track.hitPattern()].Fill(abs(track.curvatureAtVertex()),min(track.trackCompatibility(),99.5));
0513
0514         quality.Fill(getQual(track))
0515
0516     for track in kmtf:
0517         dxyKMTF.Fill(abs(track.dxy()))
0518         qualityKMTF.Fill(track.quality())
0519         etaKMTF.Fill(track.eta())
0520         fineEtaKMTF.Fill(track.hasFineEta())
0521
0522     if len(kmtf)>0:
0523         PT=kmtf[0].pt()
0524         if kmtf[0].pt()>49.99:
0525             PT=49.99
0526
0527         if abs(kmtf[0].eta())<0.7:
0528             rateKMTFp7.Fill(PT)
0529
0530         rateKMTF.Fill(PT)
0531
0532     if len(kmtfFull)>0:
0533
0534         PT=kmtfFull[0].pt()
0535         if kmtfFull[0].pt()>49.99:
0536             PT=49.99
0537         ratePerTrack[kmtfFull[0].hitPattern()].Fill(PT)
0538
0539
0540     for track in bmtf:
0541         qualityBMTF.Fill(track.quality())
0542         etaBMTF.Fill(track.eta())
0543         fineEtaBMTF.Fill(track.hasFineEta())
0544
0545     if (len(bmtf)>0):
0546         PT=bmtf[0].pt()
0547         if bmtf[0].pt()>49.99:
0548             PT=49.99
0549         if abs(bmtf[0].eta())<0.7:
0550             rateBMTFp7.Fill(PT)
0551         rateBMTF.Fill(PT)
0552
0553 #    if ( len(kmtfFull)>0)  and (kmtfFull[0].pt()>50):
0554 #        log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0555
0556
0557 #    if (len(kmtf)>0 and kmtf[0].pt()>20) and  (len(bmtf)==0 or (len(bmtf)>0 and bmtf[0].pt()<10)):
0558 #        log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0559
0560
0561 #    log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0562
0563     ##loop on gen and fill resolutuion and efficiencies
0564     for g in gen:
0565         if abs(g.eta())>0.83:
0566             continue
0567         gK=g.charge()/g.pt()
0568         genPhiAt2 = g.phi()-2.675*gK;
0569
0570         genPt.Fill(g.pt())
0571         ##the eta efficiency we want at the plateau to see strucuture
0572         if g.pt()>40.0:
0573             genEta.Fill(g.eta())
0574             genPhi.Fill(g.phi())
0575
0576         #match *(loosely because we still use coarse eta)
0577         matchedBMTF = filter(lambda x: deltaR(g.eta(),g.phi(),x.eta(),x.phi())<0.3,bmtf)
0578         matchedKMTF = filter(lambda x: deltaR(g.eta(),g.phi(),x.eta(),x.phi())<0.3,kmtf)
0579         matchedKMTFFull = filter(lambda x: deltaR(g.eta(),g.phi(),x.eta(),x.phi())<0.3,kmtfFull)
0580
0581 #        matchedBMTF = filter(lambda x: abs(deltaPhi(g.phi(),x.phi()))<2.5,bmtf)
0582 #        matchedKMTF = filter(lambda x: abs(deltaPhi(g.phi(),x.phi()))<2.5,kmtf)
0583 #        matchedKMTFFull = filter(lambda x: abs(deltaPhi(g.phi(),x.phi()))<2.5,kmtfFull)
0584
0585
0586
0587
0588         bestBMTF=None
0589         if len(matchedBMTF)>0:
0590             bestBMTF = max(matchedBMTF,key = lambda x:  x.quality()*1000+x.pt())
0591             resBMTF.Fill(g.pt(),curvResidual(bestBMTF,g))
0592             resBMTFEta.Fill(abs(g.eta()),curvResidual(bestBMTF,g))
0593
0594             resPTBMTF.Fill(g.pt(),ptResidual(bestBMTF,g))
0595             resEtaBMTF.Fill(g.eta(),bestBMTF.eta()-g.eta())
0596             resPhiBMTF.Fill(g.pt(),bestBMTF.rawPhi-genPhiAt2)
0597             if g.pt()<140:
0598                 phiCalibBMTF.Fill(g.charge()/g.pt(),bestBMTF.rawPhi-g.phi())
0599             resRBMTF.Fill(deltaR(g.eta(),g.phi(),bestBMTF.eta(),bestBMTF.phi()))
0600             bmtfCalib.Fill(bestBMTF.pt(),bestBMTF.pt()/g.pt())
0601
0602             # for the turn on , first cut on pt and then match
0603             for threshold in PTThresholds:
0604                 filteredBMTF = filter(lambda x: x.pt()>=float(threshold),matchedBMTF)
0605                 if len(filteredBMTF)>0:
0606                     genPtBMTF[threshold].Fill(g.pt())
0607                     if g.pt()>40:
0608                         genEtaBMTF[threshold].Fill(g.eta())
0609
0610
0611         bestKMTF=None
0612         if len(matchedKMTF)>0:
0613             bestKMTF = max(matchedKMTF,key = lambda x:  1000*x.quality()+x.pt())
0614
0615
0616             resKMTF.Fill(g.pt(),curvResidual(bestKMTF,g))
0617             resKMTFEta.Fill(abs(g.eta()),curvResidual(bestKMTF,g))
0618             resSTAKMTF.Fill(g.pt(),curvResidualSTA(bestKMTF,g))
0619             resPTKMTF.Fill(g.pt(),ptResidual(bestKMTF,g))
0620
0621
0622
0623             resEtaKMTF.Fill(g.eta(),bestKMTF.eta()-g.eta())
0624             resPhiKMTF.Fill(g.pt(),bestKMTF.rawPhi-genPhiAt2)
0625             if g.pt()<140:
0626                 phiCalibKMTF.Fill(g.charge()/g.pt(),bestKMTF.rawPhi-g.phi())
0627
0628             resRKMTF.Fill(deltaR(g.eta(),g.phi(),bestKMTF.eta(),bestKMTF.phi()))
0629             K = bestKMTF.charge()/bestKMTF.pt()
0630             if K==0:
0631                 K=1;
0632
0633 #        if len(matchedKMTF)>0 and len(matchedBMTF)>0:
0634 #            if bestBMTF.hasFineEta() and (not bestKMTF.hasFineEta()):
0635
0636 #                for s in stubsOLDTheta:
0637 #                    d=[]
0638 #                    for i in range(0,7):
0639 #                        d.append(s.position(i))
0640 #                    print(s.bxNum(),s.scNum(), s.whNum(), s.stNum(),d    )
0641
0642
0643
0644
0645 #            if abs(curvResidual(bestKMTF,g))>2.:
0646
0647
0648             for threshold in PTThresholds:
0649                 filteredKMTF = filter(lambda x: x.pt()>=float(threshold),matchedKMTF)
0650                 if len(filteredKMTF)>0:
0651                     genPtKMTF[threshold].Fill(g.pt())
0652                 if len(filteredKMTF)>0 and g.pt()>40:
0653                     genEtaKMTF[threshold].Fill(g.eta())
0654
0655 #        if (bestKMTF==None or (bestKMTF!=None and bestKMTF.pt()<15))  and bestBMTF!=None and  g.pt()>30 and abs(g.eta())>0.15 and abs(g.eta())<0.32:
0656 #            log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0657
0658 #        if bestKMTF!=None  and bestBMTF!=None and g.pt()>30 and abs(g.eta())>0.15 and abs(g.eta())<0.3 and bestKMTF.pt()<25 and bestBMTF.pt()>25:
0659 #            print('Residual Kalman=',abs(genPhiAt2-bestKMTF.rawPhi),'raw=',bestKMTF.rawPhi,'int=',bestKMTF.phiINT())
0660 #            print('Residual BMTF=',abs(genPhiAt2-bestBMTF.rawPhi),'raw=',bestBMTF.rawPhi,'int=',bestBMTF.phiINT())
0661 #            log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0662
0663 #        if bestKMTF!=None  and g.pt()<5 and  curvResidual(bestKMTF,g)>0.2:
0664 #            log(counter,stubs,gen,kmtfFull,kmtf,bmtf)
0665
0666
0667 #        if bestKMTF!=None and bestBMTF!=None and abs(bestKMTF.eta()-g.eta())>abs(bestBMTF.eta()-g.eta()):
0668 #            print('Input Theta stubs')
0669 #            for s in stubsOLDTheta:
0670 #                if s.bxNum()!=0:
0671 #                    continue
0672 #                a=[]
0673 #                for i in range(0,7):
0674 #                    a.append(s.position(i))
0675 #                print(s.whNum(),s.scNum(),s.stNum(),':',a)
0676 #            print('Combined stubs')
0677 #
0678 #            for s in stubs:
0679 #                print(s.whNum(),s.scNum(),s.stNum(),s.eta1(),s.eta2(),s.qeta1(),s.qeta2())
0680 #            import pdb;pdb.set_trace()
0681
0682         if len(matchedKMTFFull)>0:
0683             bestKMTFFull = max(matchedKMTFFull,key = lambda x: x.rank()*1000+x.pt())
0684             chiMatched[bestKMTFFull.hitPattern()].Fill(abs(bestKMTFFull.curvatureAtVertex()),bestKMTFFull.approxChi2());
0685             trackComp[bestKMTFFull.hitPattern()].Fill(abs(bestKMTFFull.curvatureAtVertex()),min(bestKMTFFull.trackCompatibility(),99.5));
0686
0687
0688             resKMTFTrack[bestKMTFFull.hitPattern()].Fill(g.pt(),curvResidual(bestKMTFFull,g))
0689             resKMTFTrack[0].Fill(g.pt(),curvResidual(bestKMTFFull,g))
0690             if bestKMTFFull.charge()>0:
0691                 kfCalibPlus[bestKMTFFull.hitPattern()].Fill(bestKMTFFull.pt(),bestKMTFFull.pt()/g.pt())
0692             else:
0693                 kfCalibMinus[bestKMTFFull.hitPattern()].Fill(bestKMTFFull.pt(),bestKMTFFull.pt()/g.pt())
0694             kfCalib.Fill(bestKMTFFull.pt(),bestKMTFFull.pt()/g.pt())
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709 f=ROOT.TFile("results_"+tag+".root","RECREATE")
0710
0711
0712
0713
0714
0715 quality.Write()
0716
0717
0718 resKMTF.Write()
0719 for n,t in resKMTFTrack.iteritems():
0720     t.Write()
0721
0722 resKMTFEta.Write()
0723 resSTAKMTF.Write()
0724 resPTKMTF.Write()
0725 resBMTF.Write()
0726 resBMTFEta.Write()
0727 resPTBMTF.Write()
0728 resEtaKMTF.Write()
0729 resEtaBMTF.Write()
0730 resPhiKMTF.Write()
0731 resRKMTF.Write()
0732 resPhiBMTF.Write()
0733
0734 phiCalibBMTF.Write()
0735 phiCalibKMTF.Write()
0736
0737
0738 resRBMTF.Write()
0739 #bmtfCalib.Write()
0740 #kfCalib.Write()
0741 genPt.Write()
0742
0743 for p in PTThresholds:
0744     genPtKMTF[p].Write()
0745     genPtBMTF[p].Write()
0746     genEtaKMTF[p].Write()
0747     genEtaBMTF[p].Write()
0748
0749     kmtfEff = ROOT.TGraphAsymmErrors(genPtKMTF[p],genPt)
0750     kmtfEff.Write("efficiencyVsPtKMTF"+str(p))
0751     bmtfEff = ROOT.TGraphAsymmErrors(genPtBMTF[p],genPt)
0752     bmtfEff.Write("efficiencyVsPtBMTF"+str(p))
0753
0754     kmtfEff = ROOT.TGraphAsymmErrors(genEtaKMTF[p],genEta)
0755     kmtfEff.Write("efficiencyVsEtaKMTF"+str(p))
0756     bmtfEff = ROOT.TGraphAsymmErrors(genEtaBMTF[p],genEta)
0757     bmtfEff.Write("efficiencyVsEtaBMTF"+str(p))
0758
0759
0760
0761     genPtKMTF[p].Add(genPtBMTF[p],-1)
0762     genPtKMTF[p].Divide(genPt)
0763     genPtKMTF[p].Write("efficiencyDiffVsPt"+str(p))
0764
0765     genEtaKMTF[p].Add(genEtaBMTF[p],-1)
0766     genEtaKMTF[p].Divide(genEta)
0767     genEtaKMTF[p].Write("efficiencyDiffVsEta"+str(p))
0768
0769
0770
0771
0772
0773 etaKMTF.Write()
0774 etaBMTF.Write()
0775
0776
0777 dxyKMTF.Write()
0778 qualityKMTF.Write()
0779 qualityBMTF.Write()
0780
0781
0782
0783 rateBMTF.Write()
0784 c = rateBMTF.GetCumulative(False)
0785 c.SetLineWidth(3)
0786 c.SetLineColor(ROOT.kBlack)
0787 c.Write("normRateBMTF")
0788
0789 for track in kfCalibPlus.keys():
0790     kfCalibPlus[track].Write()
0791     kfCalibMinus[track].Write()
0792     ratePerTrack[track].Write()
0793     chiMatched[track].Write()
0794     chiAll[track].Write()
0795     trackComp[track].Write()
0796     trackCompAll[track].Write()
0797
0798 kfCalib.Write()
0799
0800 bmtfCalib.Write()
0801
0802 rateKMTF.Write()
0803 c = rateKMTF.GetCumulative(False)
0804 c.Write("normRateKMTF")
0805
0806
0807 d = rateBMTF.GetCumulative(False)
0808 d.Divide(c)
0809 d.Write("rateRatioBMTFoverKMTF")
0810
0811
0812 rateBMTFp7.Write()
0813 c = rateBMTFp7.GetCumulative(False)
0814 c.Write("normRateBMTFEtaP7")
0815
0816 rateKMTFp7.Write()
0817 c = rateKMTFp7.GetCumulative(False)
0818 c.Write("normRateKMTFEtaP7")
0819
0820
0821
0822 fineEtaKMTF.Write()
0823 fineEtaBMTF.Write()
0824
0825 f.Close()