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