133 def findTracks(self) -> int:
134 hitPosLists = {}
135 hit_detector_ids = {}
136 stationCrossed: dict[int, dict[int, int]] = {}
137 listOfIndices: dict[int, list[int]] = {}
138 trackParams: dict[int, dict] = {}
139 self.fGenFitArray.clear()
140 self.fTrackletsArray.clear()
141 self.fitTrack2MC.clear()
142
143
144 if global_variables.withT0:
145 self.SmearedHits = self.strawtubes.withT0Estimate()
146
147 else:
148 self.SmearedHits = self.strawtubes.smearHits(global_variables.withNoStrawSmearing)
149
150 trackCandidates = []
151
152 if global_variables.realPR:
153
154 track_hits =
shipPatRec.execute(self.SmearedHits, global_variables.ShipGeo, global_variables.realPR)
155 logger.debug("PatRec returned %d track candidates", len(track_hits))
156
157 for i_track in track_hits:
158 atrack = track_hits[i_track]
159 atrack_y12 = atrack["y12"]
160 atrack_stereo12 = atrack["stereo12"]
161 atrack_y34 = atrack["y34"]
162 atrack_stereo34 = atrack["stereo34"]
163 atrack_smeared_hits = (
164 list(atrack_y12) + list(atrack_stereo12) + list(atrack_y34) + list(atrack_stereo34)
165 )
166
167 trackParams[i_track] = {
168 "k_y12": atrack.get("k_y12"),
169 "b_y12": atrack.get("b_y12"),
170 "k_y34": atrack.get("k_y34"),
171 "b_y34": atrack.get("b_y34"),
172 }
173 for sm in atrack_smeared_hits:
174 detID = sm["detID"]
175 station = self.strawtubes.det[sm["digiHit"]].GetStationNumber()
176 trID = i_track
177
178 if trID not in hitPosLists:
179 hitPosLists[trID] = ROOT.std.vector("TVectorD")()
180 listOfIndices[trID] = []
181 stationCrossed[trID] = {}
182 hit_detector_ids[trID] = ROOT.std.vector("int")()
183 hit_detector_ids[trID].push_back(detID)
184 m = array("d", [sm["xtop"], sm["ytop"], sm["z"], sm["xbot"], sm["ybot"], sm["z"], sm["dist"]])
185 hitPosLists[trID].push_back(ROOT.TVectorD(7, m))
186 listOfIndices[trID].append(sm["digiHit"])
187 if station not in stationCrossed[trID]:
188 stationCrossed[trID][station] = 0
189 stationCrossed[trID][station] += 1
190 else:
191 for sm in self.SmearedHits:
192 detID = self.strawtubes.det[sm["digiHit"]].GetDetectorID()
193 station = self.strawtubes.det[sm["digiHit"]].GetStationNumber()
194 trID = self.sTree.strawtubesPoint[sm["digiHit"]].GetTrackID()
195 if trID not in hitPosLists:
196 hitPosLists[trID] = ROOT.std.vector("TVectorD")()
197 listOfIndices[trID] = []
198 stationCrossed[trID] = {}
199 hit_detector_ids[trID] = ROOT.std.vector("int")()
200 hit_detector_ids[trID].push_back(detID)
201 m = array("d", [sm["xtop"], sm["ytop"], sm["z"], sm["xbot"], sm["ybot"], sm["z"], sm["dist"]])
202 hitPosLists[trID].push_back(ROOT.TVectorD(7, m))
203 listOfIndices[trID].append(sm["digiHit"])
204 if station not in stationCrossed[trID]:
205 stationCrossed[trID][station] = 0
206 stationCrossed[trID][station] += 1
207
208 n_too_few_hits = 0
209 n_too_few_stations = 0
210 n_prefit_fail = 0
211 n_fit_fail = 0
212 n_postfit_fail = 0
213 n_no_state = 0
214 n_no_ndf = 0
215
216 for atrack in hitPosLists:
217 if atrack < 0:
218 continue
219 pdg = 13
220 meas = hitPosLists[atrack]
221 detIDs = hit_detector_ids[atrack]
222 nM = len(meas)
223 if nM < 13:
224 n_too_few_hits += 1
225 continue
226 if len(stationCrossed[atrack]) < 3:
227 n_too_few_stations += 1
228 continue
229 if global_variables.debug:
230 self.sTree.MCTrack[atrack]
231
232
233 posM, momM = self._compute_seed_state(atrack, meas, trackParams)
234
235
236 covM = ROOT.TMatrixDSym(6)
237 resolution = self.sigma_spatial
238 if global_variables.withT0:
239 resolution *= 1.4
240 for i in range(3):
241 covM[i][i] = resolution * resolution
242 covM[0][0] = resolution * resolution * 100.0
243 for i in range(3, 6):
244 covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2)
245
246 rep = ROOT.genfit.RKTrackRep(pdg)
247
248 stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep)
249 rep.setPosMomCov(stateSmeared, posM, momM, covM)
250
251 seedState = ROOT.TVectorD(6)
252 seedCov = ROOT.TMatrixDSym(6)
253 rep.get6DStateCov(stateSmeared, seedState, seedCov)
254 theTrack = ROOT.genfit.Track(rep, seedState, seedCov)
255 hitCov = ROOT.TMatrixDSym(7)
256 hitCov[6][6] = resolution * resolution
257 hitID = 0
258 for m, detID in zip(meas, detIDs):
259 tp = ROOT.genfit.TrackPoint(theTrack)
260 measurement = ROOT.genfit.WireMeasurement(
261 m, hitCov, detID, hitID, tp
262 )
263 measurement.setMaxDistance(
264 global_variables.ShipGeo.strawtubes_geo.outer_straw_diameter / 2.0
265 - global_variables.ShipGeo.strawtubes_geo.wall_thickness
266 )
267 tp.addRawMeasurement(measurement)
268 theTrack.insertPoint(tp)
269 hitID += 1
270 trackCandidates.append([theTrack, atrack])
271
272 for entry in trackCandidates:
273
274 atrack = entry[1]
275 theTrack = entry[0]
276 try:
277 theTrack.checkConsistency()
278 except ROOT.genfit.Exception as e:
279 n_prefit_fail += 1
280 logger.warning("Problem with track before fit, not consistent %s %s", atrack, theTrack)
281 logger.warning(e.what())
282 ut.reportError(e)
283
284 try:
285 self.fitter.processTrack(theTrack)
286 except Exception:
287 n_fit_fail += 1
288 if global_variables.debug:
289 print("genfit failed to fit track")
290 error = "genfit failed to fit track"
291 ut.reportError(error)
292 continue
293
294 try:
295 theTrack.checkConsistency()
296 except ROOT.genfit.Exception as e:
297 n_postfit_fail += 1
298 if global_variables.debug:
299 print("Problem with track after fit, not consistent", atrack, theTrack)
300 print(e.what())
301 error = "Problem with track after fit, not consistent"
302 ut.reportError(error)
303 try:
304 fittedState = theTrack.getFittedState()
305 fittedState.getMomMag()
306 except Exception:
307 n_no_state += 1
308 error = "Problem with fittedstate"
309 ut.reportError(error)
310 continue
311 fitStatus = theTrack.getFitStatus()
312 try:
313 fitStatus.isFitConverged()
314 except ROOT.genfit.Exception:
315 error = "Fit not converged"
316 ut.reportError(error)
317 nmeas = fitStatus.getNdf()
318 global_variables.h["nmeas"].Fill(nmeas)
319 if nmeas <= 0:
320 n_no_ndf += 1
321 continue
322 chi2 = fitStatus.getChi2() / nmeas
323 global_variables.h["chi2"].Fill(chi2)
324
325
326 trackCopy = ROOT.genfit.Track(theTrack)
327 ROOT.SetOwnership(trackCopy, False)
328 self.fGenFitArray.push_back(trackCopy)
329 if global_variables.debug:
330 print("save track", theTrack, chi2, nmeas, fitStatus.isFitConverged())
331
332 track_ids = []
333 for index in listOfIndices[atrack]:
334 ahit = self.sTree.strawtubesPoint[index]
335 track_ids += [ahit.GetTrackID()]
336 _frac, tmax = self.fracMCsame(track_ids)
337 self.fitTrack2MC.push_back(tmax)
338
339 indices = ROOT.std.vector("unsigned int")()
340 for index in listOfIndices[atrack]:
341 indices.push_back(index)
342 aTracklet = ROOT.Tracklet(1, indices)
343 self.fTrackletsArray.push_back(aTracklet)
344
345 logger.debug(
346 "findTracks: %d candidates, %d too few hits, %d too few stations, "
347 "%d prefit fail, %d fit fail, %d postfit fail, %d no state, "
348 "%d no NDF, %d fitted tracks saved",
349 len(hitPosLists),
350 n_too_few_hits,
351 n_too_few_stations,
352 n_prefit_fail,
353 n_fit_fail,
354 n_postfit_fail,
355 n_no_state,
356 n_no_ndf,
357 len(self.fGenFitArray),
358 )
359
360
361 if global_variables.debug:
362 print("save tracklets:")
363 for x in self.recoTree.Tracklets:
364 print(x.getType(), len(x.getList()))
365 return len(self.fGenFitArray)
366
def execute(smeared_hits, ship_geo, str method="")