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