Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:23:17

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()