Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 import ROOT
0002 import itertools
0003 import math
0004 from DataFormats.FWLite import Events, Handle
0005 from array import array
0006 
0007 import numpy
0008 def median(lst):
0009     return numpy.median(numpy.array(lst))
0010 
0011 
0012 
0013 
0014 def fetchSegmentsPhi(event,ontime=True,twinMux=True):
0015     phiSeg    = Handle  ('L1MuDTChambPhContainer')
0016     if twinMux:
0017         event.getByLabel('simTwinMuxDigis',phiSeg)
0018     else:
0019         event.getByLabel('simDtTriggerPrimitiveDigis',phiSeg)
0020     if ontime:
0021         filtered=filter(lambda x: x.bxNum()==0, phiSeg.product().getContainer())
0022         return filtered
0023     else:
0024         return phiSeg.product().getContainer()
0025 
0026 def fetchSegmentsEta(event,ontime=True):
0027     thetaSeg  = Handle  ('L1MuDTChambThContainer')
0028     event.getByLabel('dtTriggerPrimitiveDigis',thetaSeg)
0029     if ontime:
0030         filtered=filter(lambda x: x.bxNum()==0, thetaSeg.product().getContainer())
0031         return filtered
0032     else:
0033         return thetaSeg.product().getContainer()    
0034 
0035 def fetchGEANT(event):
0036     geantH  = Handle  ('vector<PSimHit>')
0037     event.getByLabel('g4SimHits:MuonDTHits',geantH)
0038     geant=filter(lambda x: x.pabs()>0.5 and abs(x.particleType())==13,geantH.product())
0039     return geant
0040 
0041 def fetchGEN(event,etaMax=1.2):
0042     genH  = Handle  ('vector<reco::GenParticle>')
0043     event.getByLabel('genParticles',genH)
0044     genMuons=filter(lambda x: abs(x.pdgId())==13 and x.status()==1 and abs(x.eta())<etaMax,genH.product())
0045     return genMuons
0046 
0047 def segINT(seg,f1=1,f2=1):
0048     return seg.phi()*f1,seg.phiB()*f2
0049 
0050 
0051 def qPTInt(qPT,bits=14):
0052     lsb = lsBIT(bits)
0053     floatbinary = int(math.floor(abs(qPT)/lsb))
0054     return int((qPT/abs(qPT))*floatbinary)
0055 
0056 def lsBIT(bits=14):
0057     maximum=1.25
0058     lsb = 1.25/pow(2,bits-1)
0059     return lsb
0060 
0061 
0062 def getTrueCurvature(muon,geant,segments):
0063     thisMuonGEANT = filter(lambda x: (muon.charge()>0 and x.particleType()==13) or ((muon.charge()<0) and x.particleType()==-13),geant)
0064     energyInfo={1:[], 2:[],3:[],4:[]}
0065     qInfo={1:0.0, 2:0.0,3:0.0,4:0.0}
0066     qInfoINT={1:0, 2:0,3:0,4:0}
0067     for p in thisMuonGEANT:
0068         detid=ROOT.DTChamberId(p.detUnitId())
0069         station = detid.station()
0070         for s in segments:
0071             if s.stNum()==detid.station() and s.whNum()==detid.wheel() and s.scNum()==detid.sector()-1:
0072                 energyInfo[station].append(p.pabs()*muon.pt()/muon.energy())
0073                 break;
0074 
0075             
0076     for s in [1,2,3,4]:
0077         if len(energyInfo[s])==0:
0078             continue
0079         p = median(energyInfo[s])
0080         qInfo[s]=muon.charge()/p 
0081         qInfoINT[s] = qPTInt(qInfo[s])
0082     return qInfo,qInfoINT    
0083 
0084 def matchTrack(muon,segments,geant):
0085     thisMuonGEANT = filter(lambda x: (muon.charge()>0 and x.particleType()==13) or ((muon.charge()<0) and x.particleType()==-13),geant)
0086     chambers=[]
0087     for p in thisMuonGEANT:        
0088         detid=ROOT.DTChamberId(p.detUnitId())
0089         chambers.append(p.detUnitId())
0090         
0091     chambers=list(set(chambers))
0092 
0093 
0094     assocSeg=[]   
0095     for s in segments:
0096         for c in chambers:
0097             detid=ROOT.DTChamberId(c)
0098             if s.whNum()==detid.wheel() and s.stNum()==detid.station() and s.scNum()==detid.sector()-1:
0099                 if not (s in assocSeg):
0100                     assocSeg.append(s)
0101 
0102     return assocSeg
0103 
0104 
0105 
0106 
0107 events = Events([
0108     'file:singleMuonOfficial.root',
0109 ]
0110 )
0111 
0112 
0113 
0114 
0115 stations=[1,2,3,4]
0116 PHISCALE=pow(2,11)
0117 PHIBSCALE=pow(2,9)
0118 PHIFACTOR = 1
0119 PHIBFACTOR =8
0120 RELFACTOR = 1
0121 
0122 
0123 DROR = {4:0.147*RELFACTOR,3:0.173*RELFACTOR,2:0.154*RELFACTOR}
0124 DRORB = {4:(1+0.147),3:(1+0.173),2:(1+0.154)}
0125 alpha = {4:-0.0523,3:-0.0793,2:-0.0619}
0126 beta = {4:0.069,3:0.079,2:0.055}
0127 
0128 
0129 DRORCHI = {4: (726.-433.)/433. ,
0130            3: (619.-433.)/433. ,
0131            2: (512.-433.)/433.}
0132 
0133 
0134 binsk = 200
0135 maxk=8192
0136 
0137 
0138 histos={}
0139 
0140 offset={1:0.156,2:0.138,3:0.775,4:0.0}
0141 offsetINV={1:0.207,2:0.,3:0.,4:0.0}
0142 
0143 histos['phiProp']={}
0144 histos['phiPropChi']={}
0145 histos['phiBProp']={}
0146 histos['curvFromPhiB']={}
0147 histos['curvFromDPhi']={}
0148 histos['phiBFromCurv']={}
0149 histos['phiPropChiV']={}
0150 histos['deltaPhiVsPhiB']={}
0151 
0152 
0153 for i,j in itertools.permutations([1,2,3,4],2):
0154     if not (i in histos['deltaPhiVsPhiB'].keys()):
0155         histos['deltaPhiVsPhiB'][i]={}
0156     histos['deltaPhiVsPhiB'][i][j]=ROOT.TH2D("deltaPhiVsPhiB_"+str(i)+"_"+str(j),"",256,-511,512,2048,-2047,2048)
0157 
0158     if not (i in histos['curvFromDPhi'].keys()):
0159         histos['curvFromDPhi'][i]={}
0160     histos['curvFromDPhi'][i][j]=ROOT.TH2D("curvFromDPhi_"+str(i)+"_"+str(j),"",512,-2047,2048,1024,-8192,8192)
0161     
0162     
0163 for s in [1,2,3,4]:
0164     histos['curvFromPhiB'][s]=ROOT.TH2D("curvFromPhiB_"+str(s),"",1024,-512,511,4*400,-8*400,8*400)
0165     histos['phiBFromCurv'][s]=ROOT.TH2D("phiBFromCurv_"+str(s),"",256,-512,511,256,-511,512)
0166     histos['phiProp'][s]=ROOT.TH2D("phiProp_"+str(s),"",binsk,-maxk,maxk,50,-200,200)
0167     histos['phiPropChiV'][s]=ROOT.TH2D("phiPropChiV_"+str(s),"",binsk,-maxk,maxk,50,-200,200)
0168     histos['phiPropChi'][s]=ROOT.TH2D("phiPropChi_"+str(s),"",binsk,-3000,3000,50,-200,200)
0169     histos['phiBProp'][s]=ROOT.TH2D("phiBProp_"+str(s),"",binsk,-maxk,maxk,100,-2000,2000)
0170 
0171 
0172     
0173 N=0
0174 for event in events:
0175     N=N+1
0176     if N==100000:
0177         break;
0178     genMuons=fetchGEN(event)
0179     segments=fetchSegmentsPhi(event)
0180     segmentsTheta=fetchSegmentsEta(event)
0181     geant=fetchGEANT(event)
0182     segmentsTheta=sorted(segmentsTheta,key=lambda x: x.stNum())
0183 
0184     
0185     for g in genMuons:
0186         trueK,trueKINT = getTrueCurvature(g,geant,segments)
0187         cotTheta = int(g.eta()/0.010875)
0188         segTheta=matchTrack(g,segmentsTheta,geant)
0189         seg=matchTrack(g,segments,geant)
0190 
0191 
0192         for s in seg:
0193             phi,phiB=segINT(s,PHIFACTOR,PHIBFACTOR)
0194             histos['curvFromPhiB'][s.stNum()].Fill(s.phiB(),trueKINT[s.stNum()])
0195 #            histos['phiBFromCurv'][s.stNum()].Fill(trueKINT[s.stNum()]>>4,phiB)
0196             histos['phiBFromCurv'][s.stNum()].Fill(qPTInt(g.charge()/g.pt())>>4,s.phiB())
0197                 
0198         for s1,s2 in itertools.permutations(seg,2):
0199             phi1,phiB1=segINT(s1,PHIFACTOR,PHIBFACTOR)
0200             phi2,phiB2 = segINT(s2,PHIFACTOR,PHIBFACTOR)
0201 
0202             if (s2.scNum()==s1.scNum()+1) or (s1.scNum()==11 and s2.scNum()==0) :
0203                 phi2=phi2+2144
0204             if (s2.scNum()==s1.scNum()-1) or (s1.scNum()==0 and s2.scNum()==11) :
0205                 phi2=phi2-2144
0206                 
0207             
0208             
0209             if s1.code()>4 and (s1.stNum()!=s2.stNum()):
0210                 histos['deltaPhiVsPhiB'][s1.stNum()][s2.stNum()].Fill(s1.phiB(),phi2-phi1)
0211                 histos['curvFromDPhi'][s1.stNum()][s2.stNum()].Fill(phi2-phi1,qPTInt(g.charge()/g.pt()))
0212 
0213             if s1.stNum()+1==s2.stNum():                
0214                 if s2.scNum()==s1.scNum()+1 or (s2.scNum()==0 and s1.scNum()==11):
0215                     phi2=phi2+2144
0216                 if s1.scNum()==s2.scNum()+1 or (s2.scNum()==11 and s1.scNum()==0):
0217                     phi2=phi2-2144
0218                 st=s2.stNum()    
0219                 qPT=trueKINT[st]
0220                 propPhi = phi2-phiB2*DROR[st]+alpha[st]*qPT    
0221                 propPhiB =DRORB[st]*phiB2+beta[st]*qPT 
0222                 
0223                 histos['phiProp'][s2.stNum()].Fill(trueKINT[s2.stNum()],(phi1-phi2)+DROR[s2.stNum()]*phiB2)
0224                 histos['phiBProp'][s2.stNum()].Fill(trueKINT[s2.stNum()],phiB1-DRORB[s2.stNum()]*phiB2)
0225 
0226 
0227                 # for chi 2 lookmonly from station 1 -> 2,3,4
0228             if s1.stNum()==1 and s2.stNum()!=1: 
0229                 histos['phiPropChi'][s2.stNum()].Fill(trueKINT[s1.stNum()],(phi2-phi1)+DRORCHI[s2.stNum()]*phiB1)
0230                 histos['phiPropChiV'][s2.stNum()].Fill(qPTInt(g.charge()/g.pt()),(phi2-phi1))
0231                 
0232                 
0233 f=ROOT.TFile("calibrationConstants.root","RECREATE")
0234 for s in [4,3,2,1]:
0235     histos['phiProp'][s].Write()
0236     histos['phiPropChi'][s].Write()
0237     histos['phiPropChiV'][s].Write()
0238 
0239     histos['phiBProp'][s].Write()
0240     histos['curvFromPhiB'][s].Write()
0241     histos['phiBFromCurv'][s].Write()
0242 
0243 
0244 for i,j in itertools.permutations([1,2,3,4],2):
0245     histos['deltaPhiVsPhiB'][i][j].Write()
0246     histos['curvFromDPhi'][i][j].Write()
0247 
0248     
0249 f.Close()