Loading...
Searching...
No Matches
SndlhcTracking.Tracking Class Reference
Inheritance diagram for SndlhcTracking.Tracking:
Collaboration diagram for SndlhcTracking.Tracking:
Public Member Functions | |
| Init (self, online=False) | |
| SetTrackClassType (self, genfitTrack) | |
| FinishEvent (self) | |
| ExecuteTask (self, option='ScifiDS') | |
| DStrack (self) | |
| Scifi_track (self) | |
| scifiCluster (self) | |
| dsCluster (self) | |
| patternReco (self) | |
| fitTrack (self, hitlist) | |
| trackDir (self, theTrack) | |
| multipleTrackCandidates (self, planesWithClusters=10, nMaxCl=8, dGap=0.2, dMax=0.8, dMax3=0.8, ovMax=1, doublet=True, debug=False) | |
Detailed Description
Definition at line 24 of file SndlhcTracking.py.
Member Function Documentation
◆ dsCluster()
| SndlhcTracking.Tracking.dsCluster | ( | self | ) |
Definition at line 382 of file SndlhcTracking.py.
382 def dsCluster(self):
383 clusters = []
384 hitDict = {}
385 for k in range(self.event.Digi_MuFilterHits.GetEntries()):
386 d = self.event.Digi_MuFilterHits[k]
387 if (d.GetDetectorID()//10000)<3 or (not d.isValid()): continue
388 hitDict[d.GetDetectorID()] = k
389 hitList = list(hitDict.keys())
390 if len(hitList)>0:
391 hitList.sort()
392 tmp = [ hitList[0] ]
393 cprev = hitList[0]
394 ncl = 0
395 last = len(hitList)-1
396 hitvector = ROOT.std.vector("MuFilterHit*")()
397 for i in range(len(hitList)):
398 if i==0 and len(hitList)>1: continue
399 c=hitList[i]
400 neighbour = False
401 if (c-cprev)==1 or (c-cprev)==2: # allow for one missing channel
402 neighbour = True
403 tmp.append(c)
404 if not neighbour or c==hitList[last] or c%1000==59:
405 first = tmp[0]
406 N = len(tmp)
407 hitvector.clear()
408 for aHit in tmp: hitvector.push_back( self.event.Digi_MuFilterHits[hitDict[aHit]])
409 aCluster = ROOT.sndCluster(first,N,hitvector,self.mufiDet,False)
410 clusters.append(aCluster)
411 if c!=hitList[last]:
412 ncl+=1
413 tmp = [c]
414 elif not neighbour : # save last channel
415 hitvector.clear()
416 hitvector.push_back( self.event.Digi_MuFilterHits[hitDict[c]])
417 aCluster = ROOT.sndCluster(c,1,hitvector,self.mufiDet,False)
418 clusters.append(aCluster)
419 cprev = c
420 self.clusMufi.Delete()
421 for c in clusters: self.clusMufi.Add(c)
422
◆ DStrack()
| SndlhcTracking.Tracking.DStrack | ( | self | ) |
Definition at line 138 of file SndlhcTracking.py.
138 def DStrack(self):
139 event = self.event
140 trackCandidates = []
141 clusters = self.clusMufi
142
143 stations = {}
144 clusPerStation = {}
145 planesPerProjection = {0:0,1:0}
146 s = 3
147 for p in range(self.systemAndPlanes[s]+1):
148 stations[s*10+p] = {}
149 clusPerStation[s*10+p] = 0
150 k=-1
151 for aCl in clusters:
152 k+=1
153 detID = aCl.GetFirst()
154 if detID//10000 < 3: continue
155 p = (detID//1000)%10
156 bar = detID%1000
157 plane = s*10+p
158 if bar<60:
159 plane = s*10+2*p # even numbers horizontal planes, odd numbers vertical planes
160 else:
161 plane = s*10+2*p+1
162 stations[plane][k] = aCl
163 clusPerStation[plane] +=1
164 for p in clusPerStation:
165 if clusPerStation[p]>0:
166 planesPerProjection[p%2]+=1
167
168 failed = False
169 if planesPerProjection[1]<self.DSnPlanes or planesPerProjection[0]<self.DSnPlanes: return trackCandidates
170 if self.DSnPlanes==2 and (planesPerProjection[1]==2 or planesPerProjection[0]==2):
171 # require max 1 cluster per plane if only 2 planes hit per projection
172 for p in clusPerStation:
173 if clusPerStation[p]>1:
174 failed = True
175 break
176 if failed: return trackCandidates
177 hitlist = {}
178 for p in stations:
179 for k in stations[p]: hitlist[k] = stations[p][k]
180 trackCandidates.append(hitlist)
181 return trackCandidates
182
183 for p in clusPerStation:
184 if clusPerStation[p]>self.DSnHits: return trackCandidates
185
186# require one plane with 1 cluster as seed
187
188# proj = 0, horizontal, max 3 planes
189# proj = 1, vertex, max 4 planes
190 seed = -1
191 combinations = {}
192 hitlist = {}
193 for proj in range(2):
194 for plane in range(self.systemAndPlanes[s]+1):
195 if not plane%2==proj: continue
196 if clusPerStation[s*10+plane]==1:
197 seed = s*10+plane
198 break
199 if seed < 0: return trackCandidates
200 combinations[proj] = {}
201 for kA in stations[seed]:
202 clA = stations[seed][kA]
203 clA.GetPosition(A,B)
204 posA = (A+B)/2.
205 clInPlane = {}
206 for p2 in range(self.systemAndPlanes[s]+1):
207 if not p2%2==proj: continue
208 planeB = s*10+p2
209 if planeB == seed: continue
210 l = len(stations[planeB])
211 if l>0: clInPlane[planeB] = l
212 srt = sorted(clInPlane)
213 if len(stations[srt[0]])+len(stations[srt[1]]) > self.DSnHits+1: return trackCandidates
214 planeB = srt[0]
215 for kB in stations[planeB]:
216 clB = stations[planeB][kB]
217 clB.GetPosition(A,B)
218 posB = (A+B)/2.
219 delBA = posB-posA
220 if proj==0:
221 lam = delBA[1]/delBA[2]
222 b = posA[1]-lam*posA[2]
223 else:
224 lam = delBA[0]/delBA[2]
225 b = posA[0]-lam*posA[2]
226 for p3 in range(self.systemAndPlanes[s]+1):
227 if not p3%2==proj: continue
228 planeC = s*10+p3
229 if planeC == seed or planeC == planeB: continue
230 for kC in stations[planeC]:
231 clC = stations[planeC][kC]
232 clC.GetPosition(A,B)
233 posC = (A+B)/2.
234 eX = posC[2]*lam+b
235 if proj==0: res = abs(eX-posC[1])
236 else: res = abs(eX-posC[0])
237 if proj==0: combinations[proj][res] = [[seed,kA],[planeB,kB],[planeC,kC]]
238 else:
239 for p4 in range(self.systemAndPlanes[s]+1):
240 if not p4%2==proj: continue
241 planeD = s*10+p4
242 if planeD == seed or planeD == planeB or planeD == planeC: continue
243 if len(stations[planeD])==0:
244 combinations[proj][res] = [[seed,kA],[planeB,kB],[planeC,kC]]
245 for kD in stations[planeD]:
246 clD = stations[planeD][kD]
247 clD.GetPosition(A,B)
248 posD = (A+B)/2.
249 eX = posD[2]*lam+b
250 if proj==0: res+= abs(eX-posD[1])
251 else: res+= abs(eX-posD[0])
252 combinations[proj][res] = [[seed,kA],[planeB,kB],[planeC,kC],[planeD,kD]]
253 # find combination with smallest residual
254 srt = sorted(combinations[proj])[0]
255 for x in combinations[proj][srt]:
256 hitlist[x[1]] = stations[x[0]][x[1]]
257 trackCandidates.append(hitlist)
258 return trackCandidates
259
◆ ExecuteTask()
| SndlhcTracking.Tracking.ExecuteTask | ( | self, | |
option = 'ScifiDS' |
|||
| ) |
Definition at line 94 of file SndlhcTracking.py.
94 def ExecuteTask(self,option='ScifiDS'):
95 self.trackCandidates = {}
96 if not option.find('DS')<0:
97 self.clusMufi.Delete()
98 self.dsCluster()
99 self.trackCandidates['DS'] = self.DStrack()
100 if not option.find('Scifi')<0:
101 self.clusScifi.Delete()
102 self.scifiCluster()
103 self.trackCandidates['Scifi'] = self.Scifi_track()
104 i_muon = -1
105 for x in self.trackCandidates:
106 for aTrack in self.trackCandidates[x]:
107 rc = self.fitTrack(aTrack)
108 if type(rc)==type(1):
109 print('trackfit failed',rc,aTrack)
110 else:
111 i_muon += 1
112 if x=='DS': rc.SetUniqueID(3)
113 if x=='Scifi': rc.SetUniqueID(1)
114 # add the tracks
115 if self.genfitTrack:
116 self.fittedTracks.Add(rc)
117 else:
118 # Load items into snd track class object
119 #if not rc.getFitStatus().isFitConverged(): continue
120 this_track = ROOT.sndRecoTrack(rc)
121 pointTimes = ROOT.std.vector(ROOT.std.vector('float'))()
122 if x=='DS':
123 for pnt in rc.getPointsWithMeasurement():
124 hitID = pnt.getRawMeasurement().getHitId()
125 aCl = self.clusMufi[hitID]
126 pointTimes.push_back([aCl.GetTime()])
127 if x=='Scifi':
128 for pnt in rc.getPointsWithMeasurement():
129 hitID = pnt.getRawMeasurement().getHitId()
130 aCl = self.clusScifi[hitID]
131 pointTimes.push_back([aCl.GetTime()])
132 this_track.setRawMeasTimes(pointTimes)
133 this_track.setTrackType(rc.GetUniqueID())
134 # Store the track in sndRecoTrack format
135 self.fittedTracks[i_muon] = this_track
136
137
◆ FinishEvent()
| SndlhcTracking.Tracking.FinishEvent | ( | self | ) |
Definition at line 91 of file SndlhcTracking.py.
91 def FinishEvent(self):
92 pass
93
◆ fitTrack()
| SndlhcTracking.Tracking.fitTrack | ( | self, | |
| hitlist | |||
| ) |
Definition at line 453 of file SndlhcTracking.py.
453 def fitTrack(self,hitlist):
454# hitlist: clusterID: [A,B] endpoints of scifiCluster
455 hitPosLists={}
456 trID = 0
457
458 posM = ROOT.TVector3(0, 0, 0.)
459 momM = ROOT.TVector3(0,0,100.) # default track with high momentum
460
461# approximate covariance
462 covM = ROOT.TMatrixDSym(6)
463 for k in hitlist:
464 aCl = hitlist[k]
465 if hasattr(aCl,"GetFirst"):
466 detID = aCl.GetFirst()
467 else:
468 detID = aCl.GetDetectorID()
469 if detID<40000: res = self.sigmaMufiDS_spatial
470 else: res = self.sigmaScifi_spatial
471 break
472
473 for i in range(3): covM[i][i] = res*res
474 for i in range(3,6): covM[i][i] = ROOT.TMath.Power(res / (4.*2.) / ROOT.TMath.Sqrt(3), 2)
475 rep = ROOT.genfit.RKTrackRep(13)
476
477# start state
478 state = ROOT.genfit.MeasuredStateOnPlane(rep)
479 rep.setPosMomCov(state, posM, momM, covM)
480
481# create track
482 seedState = ROOT.TVectorD(6)
483 seedCov = ROOT.TMatrixDSym(6)
484 rep.get6DStateCov(state, seedState, seedCov)
485 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
486
487# make measurements sorted in z
488 unSortedList = {}
489 tmpList = {}
490 A,B = ROOT.TVector3(),ROOT.TVector3()
491 for k in hitlist:
492 aCl = hitlist[k]
493 if hasattr(aCl,"GetFirst"):
494 detID = aCl.GetFirst()
495 aCl.GetPosition(A,B)
496 detSys = 1
497 if detID<40000: detSys=3
498 else:
499 detID = aCl.GetDetectorID()
500 detSys = 1
501 if detID<40000: detSys=3
502 if detSys==3: self.mufiDet.GetPosition(detID,A,B)
503 if detSys==1: self.scifiDet.GetSiPMPosition(detID,A,B)
504 distance = 0
505 tmp = array('d',[A[0],A[1],A[2],B[0],B[1],B[2],distance])
506 unSortedList[A[2]] = [ROOT.TVectorD(7,tmp),detID,k,detSys]
507 sorted_z=list(unSortedList.keys())
508 sorted_z.sort()
509 for z in sorted_z:
510 tp = ROOT.genfit.TrackPoint() # note how the point is told which track it belongs to
511 hitCov = ROOT.TMatrixDSym(7)
512 detSys = unSortedList[z][3]
513 if detSys==3:
514 res = self.sigmaMufiDS_spatial
515 maxDis = 1.0
516 elif detSys==2:
517 res = self.sigmaMufiUS_spatial
518 maxDis = 5.0
519 else:
520 res = self.sigmaScifi_spatial
521 maxDis = 0.1
522 hitCov[6][6] = res*res
523 measurement = ROOT.genfit.WireMeasurement(unSortedList[z][0],hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to
524 measurement.setMaxDistance(maxDis)
525 measurement.setDetId(unSortedList[z][1])
526 measurement.setHitId(unSortedList[z][2])
527 tp.addRawMeasurement(measurement) # package measurement in the TrackPoint
528 theTrack.insertPoint(tp) # add point to Track
529 if not theTrack.checkConsistency():
530 print("track not consistent")
531 theTrack.Delete()
532 return -2
533# do the fit
534 self.fitter.processTrack(theTrack) # processTrackWithRep(theTrack,rep,True)
535 fitStatus = theTrack.getFitStatus()
536 if self.Debug: print("Fit result: converged chi2 Ndf",fitStatus.isFitConverged(),fitStatus.getChi2(),fitStatus.getNdf())
537 if not fitStatus.isFitConverged() and 0>1:
538 theTrack.Delete()
539 return -1
540 if self.Debug:
541 chi2 = fitStatus.getChi2()/(fitStatus.getNdf()+1E-15)
542 fittedState = theTrack.getFittedState()
543 P = fittedState.getMomMag()
544 print("track fitted Ndf #Meas P",fitStatus.getNdf(), theTrack.getNumPointsWithMeasurement(),P)
545 for p in theTrack.getPointsWithMeasurement():
546 rawM = p.getRawMeasurement()
547 info = p.getFitterInfo()
548 if not info: continue
549 detID = rawM.getDetId()
550 print(detID,"weights",info.getWeights()[0],info.getWeights()[1],fitStatus.getNdf())
551 return theTrack
552
◆ Init()
| SndlhcTracking.Tracking.Init | ( | self, | |
online = False |
|||
| ) |
Definition at line 26 of file SndlhcTracking.py.
26 def Init(self,online=False):
27 geoMat = ROOT.genfit.TGeoMaterialInterface()
28 bfield = ROOT.genfit.ConstField(0,0,0) # constant field of zero
29 fM = ROOT.genfit.FieldManager.getInstance()
30 fM.init(bfield)
31 ROOT.genfit.MaterialEffects.getInstance().init(geoMat)
32 ROOT.genfit.MaterialEffects.getInstance().setNoEffects()
33 lsOfGlobals = ROOT.gROOT.GetListOfGlobals()
34 self.scifiDet = lsOfGlobals.FindObject('Scifi')
35 self.mufiDet = lsOfGlobals.FindObject('MuFilter')
36
37 # internal storage of clusters
38 self.clusScifi = ROOT.TObjArray(100)
39 self.DetID2Key = {}
40 self.clusMufi = ROOT.TObjArray(100)
41
42 self.fitter = ROOT.genfit.KalmanFitter()
43 self.fitter.setMaxIterations(50)
44 #internal storage of fitted tracks
45 # output is in genfit::track or sndRecoTrack format
46 # default is genfit::Track
47 if hasattr(self, "genfitTrack"): pass
48 else: self.genfitTrack = True
49
50 if self.genfitTrack:
51 self.fittedTracks = ROOT.TObjArray(10)
52 else:
53 self.fittedTracks = ROOT.TClonesArray("sndRecoTrack", 10)
54 self.sigmaScifi_spatial = 2*150.*u.um
55 self.sigmaMufiUS_spatial = 2.*u.cm
56 self.sigmaMufiDS_spatial = 0.3*u.cm
57 self.scifi_vsignal = 15.*u.cm/u.ns
58 self.firstScifi_z = 300*u.cm
59 self.dZ = 13. # distance between scifi stations = 13cm
60 self.multipleTrackStore = {}
61 self.Debug = False
62 self.ioman = ROOT.FairRootManager.Instance()
63 self.sink = self.ioman.GetSink()
64 # online mode: raw data in, converted data in output
65 # offline read only: converted data in, no output
66 # offline read/write: converted data in, converted data out
67
68 # for scifi tracking
69 self.nPlanes = 3
70 self.nClusters = 5
71 self.sigma=150*u.um
72 self.maxRes=50
73 self.maskPlane=-1
74 # for DS tracking
75 self.DSnPlanes = 3
76 self.DSnHits = 5
77 self.nDSPlanesVert = self.mufiDet.GetConfParI("MuFilter/NDownstreamPlanes")
78 self.nDSPlanesHor = self.nDSPlanesVert-1
79
80 if online:
81 self.event = self.sink.GetOutTree()
82 else:
83 self.event = self.ioman.GetInChain() # should contain all digis, but not necessarily the tracks and scifi clusters
84
85 self.systemAndPlanes = {1:2,2:5,3:7}
86 return 0
87
◆ multipleTrackCandidates()
| SndlhcTracking.Tracking.multipleTrackCandidates | ( | self, | |
planesWithClusters = 10, |
|||
nMaxCl = 8, |
|||
dGap = 0.2, |
|||
dMax = 0.8, |
|||
dMax3 = 0.8, |
|||
ovMax = 1, |
|||
doublet = True, |
|||
debug = False |
|||
| ) |
Definition at line 590 of file SndlhcTracking.py.
590 def multipleTrackCandidates(self,planesWithClusters=10,nMaxCl=8,dGap=0.2,dMax=0.8,dMax3=0.8,ovMax = 1,doublet=True,debug=False):
591 A,B = ROOT.TVector3(),ROOT.TVector3()
592 h = self.multipleTrackStore
593 h['trackCand'] = {0:{},1:{}}
594 h['sortedClusters'] = {}
595 sortedClusters = h['sortedClusters']
596 self.scifiCluster()
597 clusters = self.clusScifi
598 for aCl in clusters:
599 so = aCl.GetFirst()//100000
600 if not so in sortedClusters:
601 sortedClusters[so]=[]
602 aCl.GetPosition(A,B)
603 if aCl.GetFirst()//100000%10 == 0: pos = A[1]
604 else: pos = A[0]
605 sortedClusters[so].append([pos,A[2],aCl])
606# select events with clusters in each plane
607 if len(sortedClusters)<planesWithClusters: return
608#
609 h['mergedClusters'] = {}
610 mergedClusters = h['mergedClusters']
611 for so in sortedClusters:
612 pos={}
613 mergedClusters[so]=[]
614 for k in range(len(sortedClusters[so])):
615 pos[k] = sortedClusters[so][k][0]
616 sorted_pos = sorted(pos.items(), key=lambda x: x[1])
617 merged = -1
618 for i in range(len(sorted_pos)):
619 x = sorted_pos[i]
620 aClobj = sortedClusters[so][x[0]]
621 if merged < 0:
622 merged+=1
623 mergedClusters[so].append( [aClobj[0],aClobj[1],[aClobj[2]]] )
624 else:
625 prevPos = mergedClusters[so][merged][0]
626 pos = aClobj[0]
627 if pos-prevPos > dGap:
628 merged+=1
629 mergedClusters[so].append( [aClobj[0],aClobj[1],[aClobj[2]]] )
630 else:
631 N = len(mergedClusters[so][merged][2])
632 newPos = (prevPos*N+pos)/(N+1)
633 mergedClusters[so][merged][0]=newPos
634 mergedClusters[so][merged][2].append(aClobj[2])
635# not more than nMaxCl in a plane and at least 2
636 for so in mergedClusters:
637 if len(mergedClusters[so]) > nMaxCl: return
638 if len(mergedClusters[so]) < 2: return
639 for p in range(2):
640 if debug: print('-------- p=',p)
641 if doublet:
642# method using doublets
643 h['doublet'] = {}
644 for j1 in range(1,5):
645 h['doublet'][j1] = {}
646 for k1 in range(len(mergedClusters[j1*10+p])):
647 aCl1 = mergedClusters[j1*10+p][k1]
648 h['doublet'][j1][k1] = []
649 j2 = j1+1
650 for k2 in range(len(mergedClusters[j2*10+p])):
651 aCl2 = mergedClusters[j2*10+p][k2]
652 D = aCl2[0]-aCl1[0]
653 if debug: print(j1,j2,'x1',aCl1[0],'x2',aCl2[0],'D',D)
654 if abs(D ) > dMax: continue
655 h['doublet'][j1][k1].append(10*j2+k2)
656 if len(h['doublet'][j1][k1]) == 0 and j2<5: # allow one missing plane
657 j2 = j1+2
658 for k2 in range(len(mergedClusters[j2*10+p])):
659 aCl2 = mergedClusters[j2*10+p][k2]
660 D = aCl2[0]-aCl1[0]
661 if debug: print(j1,j2,'x1',aCl1[0],'x2',aCl2[0],'D',D)
662 if abs(D ) > dMax: continue
663 h['doublet'][j1][k1].append(10*j2+k2)
664 h['trackCand'][p] = {}
665 for j1 in [1,2]:
666 for k1 in h['doublet'][j1]:
667 if j1==2:
668 alreadyUsed = False
669 for k0 in h['doublet'][1]:
670 for sk0 in h['doublet'][1][k0]:
671 if k1== (sk0%10) :
672 alreadyUsed = True
673 break
674 if alreadyUsed: continue
675 trackId = (k1+1)*10**(j1-1)
676 if debug: print(j1,k1,trackId)
677 for sk2 in h['doublet'][j1][k1]:
678 j2 = sk2//10
679 k2 = sk2%10
680 trackId = (k1+1)*10**(j1-1) + (k2+1)*10**(j2-1)
681 if debug: print(j2,k2,trackId)
682 for sk3 in h['doublet'][j2][k2]:
683 j3 = sk3//10
684 k3 = sk3%10
685 if j3>4: continue
686 trackId = (k1+1)*10**(j1-1) + (k2+1)*10**(j2-1) + (k3+1)*10**(j3-1)
687 if debug: print(j3,k3,trackId)
688 for sk4 in h['doublet'][j3][k3]:
689 j4 = sk4//10
690 k4 = sk4%10
691 trackId = (k1+1)*10**(j1-1) + (k2+1)*10**(j2-1) + (k3+1)*10**(j3-1) + (k4+1)*10**(j4-1)
692 if debug: print(j4,k4,trackId)
693 found = False
694 if j4<5:
695 for sk5 in h['doublet'][j4][k4]:
696 found = True
697 j5 = sk5//10
698 k5 = sk5%10
699 trackId = (k1+1)*10**(j1-1) + (k2+1)*10**(j2-1) + (k3+1)*10**(j3-1) + (k4+1)*10**(j4-1) + (k5+1)*10**(j5-1)
700 h['trackCand'][p][trackId] = ROOT.TGraph()
701 if not found:
702 h['trackCand'][p][trackId] = ROOT.TGraph()
703 for trackId in h['trackCand'][p]:
704 N = -1
705 for s in range(1,6):
706 k = (trackId//10**(s-1))%10-1
707 if not k<0:
708 N+=1
709 aCl = mergedClusters[s*10+p][k]
710 h['trackCand'][p][trackId].SetPoint(N,aCl[1],aCl[0])
711# 5 plane combinatoric
712 else:
713 for k1 in range(len(mergedClusters[10+p])):
714 aCl1 = mergedClusters[10+p][k1]
715 for k2 in range(len(mergedClusters[20+p])):
716 aCl2 = mergedClusters[20+p][k2]
717 D = aCl2[0]-aCl1[0]
718 if debug: print(2,'x1',aCl1[0],'x2',aCl2[0],'D',D)
719 if abs(D) > dMax: continue
720 for k3 in range(len(mergedClusters[30+p])):
721 aCl3 = mergedClusters[30+p][k3]
722 D3 = 2*D+aCl1[0]-aCl3[0]
723 if debug: print(3,'x1',aCl1[0],'x2',aCl2[0],'x3',aCl3[0],'D',D,'D3',D3)
724 if abs(D3) > dMax3: continue
725 trackId = 1000+100*k3+10*k2+k1
726 trackId3 = trackId
727 h['trackCand'][p][trackId]=ROOT.TGraph()
728 h['trackCand'][p][trackId].SetPoint(0,aCl1[1],aCl1[0])
729 h['trackCand'][p][trackId].SetPoint(1,aCl2[1],aCl2[0])
730 h['trackCand'][p][trackId].SetPoint(2,aCl3[1],aCl3[0])
731 extr4 = h['trackCand'][p][trackId].Eval(aCl3[1]+self.dZ)
732 for k4 in range(len(mergedClusters[40+p])):
733 aCl4 = mergedClusters[40+p][k4]
734 D4 = aCl4[0]-extr4
735 if debug: print(4,'D4',D4,'x4',aCl4[0],'extr',extr4)
736 if abs(D4) > dMax3: continue
737 trackId = 10000+1000*k4+100*k3+10*k2+k1
738 trackId4 = trackId
739 if debug: print(4,trackId)
740 h['trackCand'][p][trackId] = h['trackCand'][p][trackId3].Clone()
741 h['trackCand'][p][trackId].SetPoint(3,aCl4[1],aCl4[0])
742 extr5 = h['trackCand'][p][trackId].Eval(aCl4[1]+self.dZ)
743 for k5 in range(len(mergedClusters[50+p])):
744 aCl5 = mergedClusters[50+p][k5]
745 D5 = aCl5[0]-extr5
746 if debug: print(5,'D5',D5,'x5',aCl5[0],'extr',extr5)
747 if abs(D5) > dMax3: continue
748 trackId = 100000+10000*k5+1000*k4+100*k3+10*k2+k1
749 h['trackCand'][p][trackId] = h['trackCand'][p][trackId4].Clone()
750 h['trackCand'][p][trackId].SetPoint(4,aCl5[1],aCl5[0])
751 if debug: print('final',trackId)
752 '''
753clonekiller
754if more than 2 clusters are shared by a track, take the track with better chi2
755 '''
756 h['cloneCand'] = {0:[],1:[]}
757 discarded = h['cloneCand']
758 for p in range(2):
759 keys = list(h['trackCand'][p])
760 for k1 in range( len(keys)-1):
761 if keys[k1] < 100000 and not doublet: continue
762 if keys[k1] in discarded[p]: continue
763 if doublet: test1 = str(keys[k1]).zfill(5)
764 else: test1 = str(keys[k1]%100000).zfill(5)
765 for k2 in range(k1+1,len(keys) ):
766 if keys[k2] < 100000 and not doublet: continue
767 if keys[k2] in discarded[p]: continue
768 if doublet: test2 = str(keys[k2]).zfill(5)
769 else: test2 = str(keys[k2]%100000).zfill(5)
770 ov = 0
771 for i in range(5):
772 if test1[i]==test2[i]: ov+=1
773 if ov>ovMax:
774 if doublet and test1.find('0')<0 and not test2.find('0')<0: discarded[p].append(keys[k2])
775 elif doublet and test2.find('0')<0 and not test1.find('0')<0: discarded[p].append(keys[k1])
776 else:
777 rc1 = h['trackCand'][p][keys[k1]].Fit('pol1','QS')
778 rc2 = h['trackCand'][p][keys[k2]].Fit('pol1','QS')
779 if rc1.Get().Chi2() < rc2.Get().Chi2(): discarded[p].append(keys[k2])
780 else: discarded[p].append(keys[k1])
◆ patternReco()
| SndlhcTracking.Tracking.patternReco | ( | self | ) |
Definition at line 423 of file SndlhcTracking.py.
423 def patternReco(self):
424# very simple for the moment, take all scifi clusters
425 trackCandidates = []
426 hitlist = {}
427 ScifiStations = {}
428 for k in range(len(self.event.Cluster_Scifi)):
429 hitlist[k] = self.event.Cluster_Scifi[k]
430 ScifiStations[hitlist[k].GetFirst()//1000000] = 1
431# take fired muonFilter bars if more than 2 SiPMs have fired
432 nMin = 1
433 MuFiPlanes = {}
434 for k in range(self.event.Digi_MuFilterHits.GetEntries()):
435 aHit = self.event.Digi_MuFilterHits[k]
436 if not aHit.isValid(): continue
437 detID = aHit.GetDetectorID()
438 sy = detID//10000
439 l = (detID%10000)//1000 # plane number
440 bar = (detID%1000)
441 nSiPMs = aHit.GetnSiPMs()
442 nSides = aHit.GetnSides()
443 nFired = 0
444 for i in range(nSides*nSiPMs):
445 if aHit.GetSignal(i) > 0: nFired+=1
446 if nMin > nFired: continue
447 hitlist[k*1000] = self.event.Digi_MuFilterHits[k]
448 MuFiPlanes[sy*100+l] = 1
449 if (len(ScifiStations) == 5 or len(MuFiPlanes)>4) and len(hitlist)<20:
450 trackCandidates.append(hitlist)
451 return trackCandidates
452
◆ Scifi_track()
| SndlhcTracking.Tracking.Scifi_track | ( | self | ) |
Definition at line 260 of file SndlhcTracking.py.
260 def Scifi_track(self):
261# check for low occupancy and enough hits in Scifi
262 event = self.event
263 trackCandidates = []
264 clusters = self.clusScifi
265 stations = {}
266 projClusters = {0:{},1:{}}
267 for s in range(1,6):
268 for o in range(2):
269 stations[s*10+o] = []
270 k=0
271 for cl in clusters:
272 detID = cl.GetFirst()
273 s = detID//1000000
274 o = (detID//100000)%10
275 if self.maskPlane != s:
276 stations[s*10+o].append(detID)
277 projClusters[o][detID] = [cl,k]
278 k+=1
279 nclusters = 0
280 check = {}
281 ignore = []
282 for o in range(2):
283 check[o]={}
284 for s in range(1,6):
285 if len(stations[s*10+o]) > self.nClusters:
286 ignore.append(s*10+o)
287 elif len(stations[s*10+o]) > 0 :
288 check[o][s] = 1
289 nclusters+=len(stations[s*10+o])
290
291 if len(check[0])<self.nPlanes or len(check[1])<self.nPlanes: return trackCandidates
292# build trackCandidate
293# PR logic, fit straight line in x/y projection, remove outliers. Ignore tilt.
294 hitlist = {}
295 sortedClusters = {}
296 check[0]=0
297 check[1]=0
298 for o in range(2):
299 sortedClusters[o]=sorted(projClusters[o])
300 g = ROOT.TGraph()
301 n = 0
302 mean = {}
303 points = {}
304 for detID in sortedClusters[o]:
305 s = detID//1000000
306 if (s*10+o) in ignore: continue
307 if not s in mean:
308 mean[s]=[0,0,0]
309 projClusters[o][detID][0].GetPosition(A,B)
310 z = (A[2]+B[2])/2.
311 if o==0: y = (A[1]+B[1])/2.
312 else: y = (A[0]+B[0])/2.
313 points[detID] = [z,y]
314 mean[s][0]+=z
315 mean[s][1]+=y
316 mean[s][2]+=1
317 for s in mean:
318 for n in range(2):
319 mean[s][n]=mean[s][n]/mean[s][2]
320 g.AddPoint(mean[s][0],mean[s][1])
321 rc = g.Fit('pol1','SQ')
322 fun = g.GetFunction('pol1')
323 for detID in points:
324 z = points[detID][0]
325 y = points[detID][1]
326 res = abs(y-fun.Eval(z))/self.sigma
327 if res < self.maxRes:
328 k = projClusters[o][detID][1]
329 hitlist[k] = projClusters[o][detID][0]
330 check[o]+=1
331 if check[0]>2 and check[1]>2:
332 trackCandidates.append(hitlist)
333 return trackCandidates
334
◆ scifiCluster()
| SndlhcTracking.Tracking.scifiCluster | ( | self | ) |
Definition at line 335 of file SndlhcTracking.py.
335 def scifiCluster(self):
336 clusters = []
337 self.DetID2Key.clear()
338 hitDict = {}
339 for k in range(self.event.Digi_ScifiHits.GetEntries()):
340 d = self.event.Digi_ScifiHits[k]
341 if not d.isValid(): continue
342 hitDict[d.GetDetectorID()] = k
343 hitList = list(hitDict.keys())
344 if len(hitList)>0:
345 hitList.sort()
346 tmp = [ hitList[0] ]
347 cprev = hitList[0]
348 ncl = 0
349 last = len(hitList)-1
350 hitvector = ROOT.std.vector("sndScifiHit*")()
351 for i in range(len(hitList)):
352 if i==0 and len(hitList)>1: continue
353 c=hitList[i]
354 neighbour = False
355 if (c-cprev)==1: # does not account for neighbours across sipms
356 neighbour = True
357 tmp.append(c)
358 if not neighbour or c==hitList[last]:
359 first = tmp[0]
360 N = len(tmp)
361 hitvector.clear()
362 for aHit in tmp: hitvector.push_back( self.event.Digi_ScifiHits[hitDict[aHit]])
363 aCluster = ROOT.sndCluster(first,N,hitvector,self.scifiDet,False)
364 clusters.append(aCluster)
365 if c!=hitList[last]:
366 ncl+=1
367 tmp = [c]
368 elif not neighbour : # save last channel
369 hitvector.clear()
370 hitvector.push_back( self.event.Digi_ScifiHits[hitDict[c]])
371 aCluster = ROOT.sndCluster(c,1,hitvector,self.scifiDet,False)
372 clusters.append(aCluster)
373 cprev = c
374 self.clusScifi.Delete()
375 for c in clusters:
376 self.clusScifi.Add(c)
377# map clusters to hit keys
378 for nHit in range(self.event.Digi_ScifiHits.GetEntries()):
379 if self.event.Digi_ScifiHits[nHit].GetDetectorID()==c.GetFirst():
380 self.DetID2Key[c.GetFirst()] = nHit
381
◆ SetTrackClassType()
| SndlhcTracking.Tracking.SetTrackClassType | ( | self, | |
| genfitTrack | |||
| ) |
Definition at line 88 of file SndlhcTracking.py.
88 def SetTrackClassType(self,genfitTrack):
89 self.genfitTrack = genfitTrack
90
◆ trackDir()
| SndlhcTracking.Tracking.trackDir | ( | self, | |
| theTrack | |||
| ) |
Definition at line 553 of file SndlhcTracking.py.
553 def trackDir(self,theTrack):
554 if theTrack.GetUniqueID()>1: return False # for the moment, only the scifi is time calibrated
555 fitStatus = theTrack.getFitStatus()
556 if not fitStatus.isFitConverged() : return [100,-100]
557 state = ROOT.getFittedState(theTrack,0)
558 pos = state.getPos()
559 mom = state.getMom()
560 lam = (self.firstScifi_z-pos.z())/mom.z()
561 # nominal first position
562 pos1 = ROOT.TVector3(pos.x()+lam*mom.x(),pos.y()+lam*mom.y(),self.firstScifi_z)
563
564 self.Tline = ROOT.TGraph()
565 meanT = 0
566 for nM in range(theTrack.getNumPointsWithMeasurement()):
567 state = ROOT.getFittedState(theTrack,nM)
568 if not state: continue
569 posM = state.getPos()
570 M = theTrack.getPointWithMeasurement(nM)
571 W = M.getRawMeasurement()
572 detID = W.getDetId()
573 clkey = W.getHitId()
574 aCl = self.clusScifi[clkey]
575 aHit = self.event.Digi_ScifiHits[ self.DetID2Key[detID] ]
576 self.scifiDet.GetSiPMPosition(detID,A,B)
577 if aHit.isVertical(): X = B-posM
578 else: X = A-posM
579 L = X.Mag()/self.scifi_vsignal
580 # need to correct for signal propagation along fibre
581 corTime = self.scifiDet.GetCorrectedTime(detID, aCl.GetTime(), 0)
582 trajLength = (posM-pos1).Mag()
583 T = corTime - L - trajLength/u.speedOfLight
584 self.Tline.AddPoint(trajLength,T)
585 rc = self.Tline.Fit('pol1','SQ')
586 fitResult = rc.Get()
587 slope = fitResult.Parameter(1)
588 return [slope,slope/(fitResult.ParError(1)+1E-13),fitResult.Parameter(0)]
589
Member Data Documentation
◆ clusMufi
| SndlhcTracking.Tracking.clusMufi |
Definition at line 40 of file SndlhcTracking.py.
◆ clusScifi
| SndlhcTracking.Tracking.clusScifi |
Definition at line 38 of file SndlhcTracking.py.
◆ Debug
| SndlhcTracking.Tracking.Debug |
Definition at line 61 of file SndlhcTracking.py.
◆ DetID2Key
| SndlhcTracking.Tracking.DetID2Key |
Definition at line 39 of file SndlhcTracking.py.
◆ DSnHits
| SndlhcTracking.Tracking.DSnHits |
Definition at line 76 of file SndlhcTracking.py.
◆ DSnPlanes
| SndlhcTracking.Tracking.DSnPlanes |
Definition at line 75 of file SndlhcTracking.py.
◆ dZ
| SndlhcTracking.Tracking.dZ |
Definition at line 59 of file SndlhcTracking.py.
◆ event
| SndlhcTracking.Tracking.event |
Definition at line 81 of file SndlhcTracking.py.
◆ firstScifi_z
| SndlhcTracking.Tracking.firstScifi_z |
Definition at line 58 of file SndlhcTracking.py.
◆ fittedTracks
| SndlhcTracking.Tracking.fittedTracks |
Definition at line 51 of file SndlhcTracking.py.
◆ fitter
| SndlhcTracking.Tracking.fitter |
Definition at line 42 of file SndlhcTracking.py.
◆ genfitTrack
| SndlhcTracking.Tracking.genfitTrack |
Definition at line 48 of file SndlhcTracking.py.
◆ ioman
| SndlhcTracking.Tracking.ioman |
Definition at line 62 of file SndlhcTracking.py.
◆ maskPlane
| SndlhcTracking.Tracking.maskPlane |
Definition at line 73 of file SndlhcTracking.py.
◆ maxRes
| SndlhcTracking.Tracking.maxRes |
Definition at line 72 of file SndlhcTracking.py.
◆ mufiDet
| SndlhcTracking.Tracking.mufiDet |
Definition at line 35 of file SndlhcTracking.py.
◆ multipleTrackStore
| SndlhcTracking.Tracking.multipleTrackStore |
Definition at line 60 of file SndlhcTracking.py.
◆ nClusters
| SndlhcTracking.Tracking.nClusters |
Definition at line 70 of file SndlhcTracking.py.
◆ nDSPlanesHor
| SndlhcTracking.Tracking.nDSPlanesHor |
Definition at line 78 of file SndlhcTracking.py.
◆ nDSPlanesVert
| SndlhcTracking.Tracking.nDSPlanesVert |
Definition at line 77 of file SndlhcTracking.py.
◆ nPlanes
| SndlhcTracking.Tracking.nPlanes |
Definition at line 69 of file SndlhcTracking.py.
◆ scifi_vsignal
| SndlhcTracking.Tracking.scifi_vsignal |
Definition at line 57 of file SndlhcTracking.py.
◆ scifiDet
| SndlhcTracking.Tracking.scifiDet |
Definition at line 34 of file SndlhcTracking.py.
◆ sigma
| SndlhcTracking.Tracking.sigma |
Definition at line 71 of file SndlhcTracking.py.
◆ sigmaMufiDS_spatial
| SndlhcTracking.Tracking.sigmaMufiDS_spatial |
Definition at line 56 of file SndlhcTracking.py.
◆ sigmaMufiUS_spatial
| SndlhcTracking.Tracking.sigmaMufiUS_spatial |
Definition at line 55 of file SndlhcTracking.py.
◆ sigmaScifi_spatial
| SndlhcTracking.Tracking.sigmaScifi_spatial |
Definition at line 54 of file SndlhcTracking.py.
◆ sink
| SndlhcTracking.Tracking.sink |
Definition at line 63 of file SndlhcTracking.py.
◆ systemAndPlanes
| SndlhcTracking.Tracking.systemAndPlanes |
Definition at line 85 of file SndlhcTracking.py.
◆ Tline
| SndlhcTracking.Tracking.Tline |
Definition at line 564 of file SndlhcTracking.py.
◆ trackCandidates
| SndlhcTracking.Tracking.trackCandidates |
Definition at line 95 of file SndlhcTracking.py.
The documentation for this class was generated from the following file:
- python/SndlhcTracking.py
