File indexing completed on 2024-04-06 12:20:53
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
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
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()
