Scifi_monitoring.Scifi_residuals Class Reference

Inheritance diagram for Scifi_monitoring.Scifi_residuals:

Collaboration diagram for Scifi_monitoring.Scifi_residuals:

Public Member Functions
	Init (self, options, monitor)
	ExecuteEvent (self, event)
	Plot (self)

Public Attributes
	unbiased
	M
	projs
	parallelToZ
	OT
	nav
	trackTask
	xing
	zExVeto
	zEx
	res

Detailed Description

Definition at line 165 of file Scifi_monitoring.py.

Member Function Documentation

ExecuteEvent()

Scifi_monitoring.Scifi_residuals.ExecuteEvent	(		self,
			event
	)

Definition at line 219 of file Scifi_monitoring.py.

   def ExecuteEvent(self,event):
       h = self.M.h
       W = self.M.Weight
       nav = self.nav
       if not hasattr(event,"Cluster_Scifi"):
               self.trackTask.scifiCluster()
               clusters = self.trackTask.clusScifi
       else:
               clusters = event.Cluster_Scifi
# overall tracking
       MufiTracks    = []
       ScifiTracks = []
       k = -1
       for theTrack in self.M.Reco_MuonTracks:
          k+=1
          fitStatus = theTrack.getFitStatus()
          if not fitStatus.isFitConverged(): continue
          if theTrack.GetUniqueID()==1:   ScifiTracks.append(k)
          if theTrack.GetUniqueID()==3:   MufiTracks.append(k)
       if len(MufiTracks)==0 or len(ScifiTracks)==0: return
       vetoHits = []
       k = -1
       for aHit in event.Digi_MuFilterHits:
           k+=1
           Minfo = self.M.MuFilter_PlaneBars(aHit.GetDetectorID())
           s,l,bar = Minfo['station'],Minfo['plane'],Minfo['bar']
           if s>1: continue
           X = aHit.GetAllSignals()
           if len(X)<5: continue   # number of fired SiPMs
           vetoHits.append(k)
       xExTag,yExTag = -10000,-10000
       for kMu in MufiTracks:
          theTrack = self.M.Reco_MuonTracks[kMu]
          fstate =  theTrack.getFittedState()
          posT,momT  = fstate.getPos(),fstate.getMom()
          slopeXT = momT.X()/momT.Z()
          slopeYT = momT.Y()/momT.Z()
          if not abs(slopeXT)<0.1 or not abs(slopeYT)<0.1: continue
          lam      = (self.zEx-posT.z())/momT.z()
          yExTag      = posT.y()+lam*momT.y()
          xExTag      = posT.x()+lam*momT.x()
          # eventually require hit in veto to remove ghost tracks
       if xExTag < -9000: return
       ok = False
       for k in vetoHits:
              aHit = event.Digi_MuFilterHits[k]
              self.M.MuFilter.GetPosition(aHit.GetDetectorID(),A,B)
# calculate DOCA
              lam = (self.zExVeto-posT.z())/momT.z()
              xExV,yExV = posT.x()+lam*momT.x(),posT.y()+lam*momT.y()
              pq = A-posT
              uCrossv= (B-A).Cross(momT)
              doca = pq.Dot(uCrossv)/uCrossv.Mag()
              if abs(doca)<self.res: 
                ok=True
                break
       if not ok: return
       # DS track confirmed by Veto hit
           
       theTrack = False
       for theTrack in self.M.Reco_MuonTracks:
          if theTrack.GetUniqueID()==1:
              fitStatus = theTrack.getFitStatus()
              if  fitStatus.isFitConverged():
                 state = theTrack.getFittedState()
                 pos    = state.getPos()
                 mom = state.getMom()
                 slopeX = mom.X()/mom.Z()
                 slopeY = mom.Y()/mom.Z()
                 rc = h[detector+'trackChi2/ndof'].Fill(fitStatus.getChi2()/(fitStatus.getNdf()+1E-10),fitStatus.getNdf())
                 '''
                 self.M.fillHist2(detector+'trackSlopes',slopeX*1000-pos.X()/48.2,slopeY*1000-pos.Y()/48.2)
                 self.M.fillHist2(detector+'trackSlopesXL',slopeX-pos.X()/48200,slopeY-pos.Y()/48200)
                 self.M.fillHist2(detector+'trackPos',pos.X(),pos.Y())
                 if abs(slopeX)<0.1 and abs(slopeY)<0.1:  self.M.fillHist2(detector+'trackPosBeam',pos.X(),pos.Y())
                 '''
 
       if not theTrack: return
 
       sortedClusters={}
       for aCl in clusters:
           so = aCl.GetFirst()//100000
           if not so in sortedClusters: sortedClusters[so]=[]
           sortedClusters[so].append(aCl)
# select events with clusters in each plane for making residuals
       if len(sortedClusters)<10: return
       goodEvent = True
       for s in sortedClusters:
          if len(sortedClusters[s])>3:
             goodEvent=False
             break
       if not goodEvent: return
 
       for s in range(1,6):
            if self.unbiased:
# build trackCandidate
              hitlist = {}
              if self.unbiased or s==1:
                k=0
                Nproj = {0:0,1:0}
                for so in sortedClusters:
                    if so//10 == s and self.unbiased: continue
                    Nproj[so%2]+=1
                    for x in sortedClusters[so]:
                       hitlist[k] = x
                       k+=1
                if Nproj[0]<3 or Nproj[1]<3: continue
                theTrack = self.trackTask.fitTrack(hitlist)
                if not hasattr(theTrack,"getFittedState"): continue
# check residuals
                fitStatus = theTrack.getFitStatus()
                if not fitStatus.isFitConverged(): 
                  theTrack.Delete()
                  continue
# confirm track with DS track
            fstate =  theTrack.getFittedState()
            pos,mom  = fstate.getPos(),fstate.getMom()
            lam      = (self.zExVeto-pos.z())/mom.z()
            yEx      = pos.y()+lam*mom.y()
            xEx      = pos.x()+lam*mom.x()
            dx = xExTag-xEx
            dy = yExTag-yEx
            rc = h['dx'].Fill(dx)
            rc = h['dy'].Fill(dy)
            if abs(dy)>self.res or abs(dx)>self.res:  
               if self.unbiased: continue
               else:             break
 
# test plane
            for o in range(2):
                testPlane = s*10+o
                z = self.M.zPos['Scifi'][testPlane]
                rep     = ROOT.genfit.RKTrackRep(13)
                state  = ROOT.genfit.StateOnPlane(rep)
# find closest track state
                mClose = 0
                mZmin = 999999.
                for m in range(0,theTrack.getNumPointsWithMeasurement()):
                   st   = ROOT.getFittedState(theTrack,m)
                   if not st: break
                   Pos = st.getPos()
                   if abs(z-Pos.z())<mZmin:
                      mZmin = abs(z-Pos.z())
                      mClose = m
                if mZmin>10000:
                    print("something wrong here with measurements",mClose,mZmin,theTrack.getNumPointsWithMeasurement())
                    break
                fstate =  theTrack.getFittedState(mClose)
                pos,mom = fstate.getPos(),fstate.getMom()
                rep.setPosMom(state,pos,mom)
                NewPosition = ROOT.TVector3(0., 0., z)   # assumes that plane in global coordinates is perpendicular to z-axis, which is not true for TI18 geometry.
                rep.extrapolateToPlane(state, NewPosition, parallelToZ )
                pos = state.getPos()
                xEx,yEx = pos.x(),pos.y()
                rc = h['track_Scifi'+str(testPlane)].Fill(xEx,yEx,W)
                for aCl in sortedClusters[testPlane]:
                   aCl.GetPosition(A,B)
                   detID = aCl.GetFirst()
                   channel = detID%1000 + ((detID%10000)//1000)*128 + (detID%100000//10000)*512
# calculate DOCA
                   pq = A-pos
                   uCrossv= (B-A).Cross(mom)
                   doca = pq.Dot(uCrossv)/uCrossv.Mag()
                   rc = h['resC'+self.projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,channel,W)
                   rc = h['res'+self.projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,W)
                   rc = h['resX'+self.projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,xEx,W)
                   rc = h['resY'+self.projs[o]+'_Scifi'+str(testPlane)].Fill(doca/u.um,yEx,W)
 
            if self.unbiased: theTrack.Delete()
 
# analysis and plots 

Init()

Scifi_monitoring.Scifi_residuals.Init	(		self,
			options,
			monitor
	)

Definition at line 167 of file Scifi_monitoring.py.

   def Init(self,options,monitor):
       NbinsRes = options.ScifiNbinsRes
       xmin        = options.Scifixmin
       alignPar   = options.ScifialignPar
       self.unbiased = options.ScifiResUnbiased
 
       self.M = monitor
       h = self.M.h
       self.projs = {1:'V',0:'H'}
       self.parallelToZ = ROOT.TVector3(0., 0., 1.)
       run = ROOT.FairRunAna.Instance()
       ioman = ROOT.FairRootManager.Instance()
       self.OT = ioman.GetSink().GetOutTree()
       self.nav = ROOT.gGeoManager.GetCurrentNavigator()
       self.trackTask = self.M.FairTasks['simpleTracking']
       if not self.trackTask: self.trackTask = run.GetTask('houghTransform')
 
       for s in range(1,6):
          for o in range(2):
             for p in self.projs:
               proj = self.projs[p]
               xmax = -xmin
               ut.bookHist(h,'res'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax)
               ut.bookHist(h,'resX'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,100,-50.,0.)
               ut.bookHist(h,'resY'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,100,10.,60.)
               ut.bookHist(h,'resC'+proj+'_Scifi'+str(s*10+o),'residual '+proj+str(s*10+o)+'; [#mum]',NbinsRes,xmin,xmax,128*4*3,-0.5,128*4*3-0.5)
               ut.bookHist(h,'track_Scifi'+str(s*10+o),'track x/y '+str(s*10+o)+'; x [cm]; y [cm]',80,-70.,10.,80,0.,80.)
       ut.bookHist(h,'dx','DS track X - scifi X; X[cm]',100,-20.,20.)
       ut.bookHist(h,'dy','DS track Y - scifi Y; Y[cm]',100,-20.,20.)
 
       ut.bookHist(h,detector+'trackChi2/ndof','track chi2/ndof vs ndof; #chi^{2}/Ndof; Ndof',100,0,100,20,0,20)
# type of crossing, check for b1only,b2nob1,nobeam
       self.xing = {'':True,'B1only':False,'B2noB1':False,'noBeam':False}
       '''
       for xi in self.xing:
          if not self.M.fsdict and not self.M.hasBunchInfo and xi!='': continue
          ut.bookHist(h,detector+'trackSlopes'+xi,'track slope; x/z [mrad]; y/z [mrad]',1000,-100,100,1000,-100,100)
          ut.bookHist(h,detector+'trackSlopesXL'+xi,'track slope; x/z [rad]; y/z [rad]',2200,-1.1,1.1,2200,-1.1,1.1)
          ut.bookHist(h,detector+'trackPos'+xi,'track pos; x [cm]; y [cm]',100,-90,10.,80,0.,80.)
          ut.bookHist(h,detector+'trackPosBeam'+xi,'beam track pos slopes<0.1rad; x [cm]; y [cm]',100,-90,10.,80,0.,80.)
          for item in h:
            if isinstance(h[item], ROOT.TH2) and (item.find('trackPos')>0 or item.find('trackSlopes')>0): 
              h[item].SetTitleOffset(1.1, 'Y')
       '''
       if alignPar:
            for x in alignPar:
               self.M.Scifi.SetConfPar(x,alignPar[x])
               
       self.zExVeto = self.M.zPos['MuFilter'][10]
       self.zEx = self.M.zPos['Scifi'][10]
       self.res = 10.
 

Plot()

Scifi_monitoring.Scifi_residuals.Plot

(

self

)

Definition at line 390 of file Scifi_monitoring.py.

   def Plot(self):
       h = self.M.h
       P = {'':'','X':'colz','Y':'colz','C':'colz'}
       Par = {'mean':1,'sigma':2}
       h['globalPos']   = {'meanH':ROOT.TGraphErrors(),'sigmaH':ROOT.TGraphErrors(),'meanV':ROOT.TGraphErrors(),'sigmaV':ROOT.TGraphErrors()}
       h['globalPosM'] = {'meanH':ROOT.TGraphErrors(),'sigmaH':ROOT.TGraphErrors(),'meanV':ROOT.TGraphErrors(),'sigmaV':ROOT.TGraphErrors()}
       for x in h['globalPosM']:
            h['globalPos'][x].SetMarkerStyle(21)
            h['globalPos'][x].SetMarkerColor(ROOT.kBlue)
            h['globalPosM'][x].SetMarkerStyle(21)
            h['globalPosM'][x].SetMarkerColor(ROOT.kBlue)
       globalPos = h['globalPos']
       for proj in P:
           ut.bookCanvas(h,'scifiRes'+proj,'',1600,1900,2,5)
           k=1
           j = {0:0,1:0}
           for s in range(1,6):
               for o in range(2):
                  so = s*10+o
                  tc = h['scifiRes'+proj].cd(k)
                  k+=1
                  hname = 'res'+proj+self.projs[o]+'_Scifi'+str(so)
                  h[hname].Draw(P[proj])
                  if proj == '':
                     rc = h[hname].Fit('gaus','SQ')
                     fitResult = rc.Get()
                     if not fitResult: continue
                     for p in Par:
                          globalPos[p+self.projs[o]].SetPoint(s-1,s,fitResult.Parameter(Par[p]))
                          globalPos[p+self.projs[o]].SetPointError(s-1,0.5,fitResult.ParError(1))
                  if proj == 'C':
                       for m in range(3):
                             h[hname+str(m)] = h[hname].ProjectionX(hname+str(m),m*512,m*512+512)
                             rc = h[hname+str(m)].Fit('gaus','SQ0')
                             fitResult = rc.Get()
                             if not fitResult: continue
                             for p in Par:
                                 h['globalPosM'][p+self.projs[o]].SetPoint(j[o], s*10+m,   fitResult.Parameter(Par[p]))
                                 h['globalPosM'][p+self.projs[o]].SetPointError(j[o],0.5,fitResult.ParError(1))
                             j[o]+=1
       
       S  = ctypes.c_double()
       M = ctypes.c_double()
       h['alignPar'] = {}
       alignPar = h['alignPar']
       for p in globalPos:
           ut.bookCanvas(h,p,p,750,750,1,1)
           tc = h[p].cd()
           globalPos[p].SetTitle(p+';station; offset [#mum]')
           globalPos[p].Draw("ALP")
           if p.find('mean')==0:
               for n in range(globalPos[p].GetN()):
                  rc = globalPos[p].GetPoint(n,S,M)
                  print("station %i: offset %s =  %5.2F um"%(S.value,p[4:5],M.value))
                  s = int(S.value*10)
                  if p[4:5] == "V": s+=1
                  alignPar["Scifi/LocD"+str(s)] = M.value
 
       ut.bookCanvas(h,'mean&sigma',"mean and sigma",1200,1200,2,2)
       k=1
       for p in h['globalPosM']:
           ut.bookCanvas(h,p+'M',p,750,750,1,1)
           tc = h[p+'M'].cd()
           h['globalPosM'][p].SetTitle(p+';mat ; offset [#mum]')
           h['globalPosM'][p].Draw("ALP")
           tc = h['mean&sigma'].cd(k)
           h['globalPosM'][p].Draw("ALP")
           k+=1
           if p.find('mean')==0:
              for n in range(h['globalPosM'][p].GetN()):
                 rc = h['globalPosM'][p].GetPoint(n,S,M)
                 print("station %i: offset %s =  %5.2F um"%(S.value,p[4:5],M.value))
                 s = int(S.value*10)
                 if p[4:5] == "V": s+=1
                 alignPar["Scifi/LocM"+str(s)] = M.value
       T = ['mean&sigma']
       for proj in P: T.append('scifiRes'+proj)
       for canvas in T:
           self.M.myPrint(self.M.h[canvas],"Scifi-"+canvas,subdir='scifi/expert')
       '''
       for xi in self.xing:
           if not self.M.fsdict and not self.M.hasBunchInfo and xi!='': continue
           tname = detector+'trackDir'+xi
           ut.bookCanvas(h,tname,"track directions",1600,1800,3,2)
           h[tname].cd(1)
           rc = h[detector+'trackSlopes'+xi].Draw('colz')
           h[tname].cd(2)
           rc = h[detector+'trackSlopes'+xi].ProjectionX("slopeX"+xi)
           rc.Draw()
           rc.SetTitle('track X slope')
           h[tname].cd(3)
           rc = h[detector+'trackSlopes'+xi].ProjectionY("slopeY"+xi)
           rc.Draw()
           rc.SetTitle('track Y slope')
           h[tname].cd(4)
           rc = h[detector+'trackSlopesXL'+xi].Draw('colz')
           h[tname].cd(5)
           rc = h[detector+'trackSlopesXL'+xi].ProjectionX("slopeXL"+xi)
           rc.Draw()
           rc.SetTitle('track X slope')
           h[tname].cd(6)
           rc = h[detector+'trackSlopesXL'+xi].ProjectionY("slopeYL"+xi)
           rc.Draw()
           rc.SetTitle('track Y slope')
           if x=='': self.M.myPrint(self.M.h[tname],tname,subdir='scifi/shifter')
           else:     self.M.myPrint(self.M.h[tname],tname,subdir='scifi/shifter/'+xi)
           tname = detector+'TtrackPos'+xi
           ut.bookCanvas(h,tname,"track position first state",600,1200,1,2)
           h[tname].cd(1)
           rc = h[detector+'trackPosBeam'+xi].Draw('colz')
           h[tname].cd(2)
           rc = h[detector+'trackPos'+xi].Draw('colz')
           if x=='': self.M.myPrint(self.M.h[tname],detector+'trackPos'+xi,subdir='scifi/shifter')
           else:     self.M.myPrint(self.M.h[tname],detector+'trackPos'+xi,subdir='scifi/shifter/'+xi)
       '''
           

Member Data Documentation

M

Scifi_monitoring.Scifi_residuals.M

Definition at line 173 of file Scifi_monitoring.py.

nav

Scifi_monitoring.Scifi_residuals.nav

Definition at line 180 of file Scifi_monitoring.py.

OT

Scifi_monitoring.Scifi_residuals.OT

Definition at line 179 of file Scifi_monitoring.py.

parallelToZ

Scifi_monitoring.Scifi_residuals.parallelToZ

Definition at line 176 of file Scifi_monitoring.py.

projs

Scifi_monitoring.Scifi_residuals.projs

Definition at line 175 of file Scifi_monitoring.py.

res

Scifi_monitoring.Scifi_residuals.res

Definition at line 217 of file Scifi_monitoring.py.

trackTask

Scifi_monitoring.Scifi_residuals.trackTask

Definition at line 181 of file Scifi_monitoring.py.

unbiased

Scifi_monitoring.Scifi_residuals.unbiased

Definition at line 171 of file Scifi_monitoring.py.

xing

Scifi_monitoring.Scifi_residuals.xing

Definition at line 199 of file Scifi_monitoring.py.

zEx

Scifi_monitoring.Scifi_residuals.zEx

Definition at line 216 of file Scifi_monitoring.py.

zExVeto

Scifi_monitoring.Scifi_residuals.zExVeto

Definition at line 215 of file Scifi_monitoring.py.

---

The documentation for this class was generated from the following file:

• shipLHC/scripts/Scifi_monitoring.py