30def run_track_pattern_recognition(input_file, geo_file, output_file, method, dy=None):
31 """
32 Runs all steps of track pattern recognition.
33 Parameters
34 ----------
35 input_file : string
36 Path to an input .root file with events.
37 geo_file : string
38 Path to a file with SHiP geometry.
39 output_file : string
40 Path to an output .root file with quality plots.
41 method : string
42 Name of a track pattern recognition method.
43 """
44
45
46
47
48 try:
49 fgeo = ROOT.TFile(geo_file)
50 except Exception:
51 print("An error with opening the ship geo file.")
52 raise
53
54 sGeo = fgeo.Get("FAIRGeom")
55
56
57 if not fgeo.FindKey("ShipGeo"):
58 if dy:
60 else:
62 else:
63 ShipGeo = load_from_root_file(fgeo, "ShipGeo")
64
65
66 global_variables.ShipGeo = ShipGeo
67
68
69
70 run = ROOT.FairRunSim()
71 run.SetName("TGeant4")
72
73
74 run.SetSink(ROOT.FairRootFileSink(ROOT.TMemFile("output", "recreate")))
75 run.SetUserConfig("g4Config_basic.C")
76 run.GetRuntimeDb()
77
79 run.Init()
80
81
82
83
84
85
86 import geomGeant4
87
88 if hasattr(ShipGeo.Bfield, "fieldMap"):
90 else:
91 print("no fieldmap given, geofile too old, not anymore support")
92 exit(-1)
93 sGeo = fgeo.Get("FAIRGeom")
94 geoMat = ROOT.genfit.TGeoMaterialInterface()
95 ROOT.genfit.MaterialEffects.getInstance().init(geoMat)
96 bfield = ROOT.genfit.FairShipFields()
97 bfield.setField(fieldMaker.getGlobalField())
98 fM = ROOT.genfit.FieldManager.getInstance()
99 fM.init(bfield)
100
101
102
103
104 try:
105 fn = ROOT.TFile(input_file, "update")
106 except Exception:
107 print("An error with opening the input data file.")
108 raise
109
110 sTree = fn.Get("cbmsim")
111 sTree.Write()
112
113
114
115 h = init_book_hist()
116
117
118
119
120 metrics = {
121 "n_hits": [],
122 "reconstructible": 0,
123 "passed_y12": 0,
124 "passed_stereo12": 0,
125 "passed_12": 0,
126 "passed_y34": 0,
127 "passed_stereo34": 0,
128 "passed_34": 0,
129 "passed_combined": 0,
130 "reco_passed": 0,
131 "reco_passed_no_clones": 0,
132 "frac_y12": [],
133 "frac_stereo12": [],
134 "frac_12": [],
135 "frac_y34": [],
136 "frac_stereo34": [],
137 "frac_34": [],
138 "reco_frac_tot": [],
139 "reco_mc_p": [],
140 "reco_mc_theta": [],
141 "fitted_p": [],
142 "fitted_pval": [],
143 "fitted_chi": [],
144 "fitted_x": [],
145 "fitted_y": [],
146 "fitted_z": [],
147 "fitted_mass": [],
148 }
149
150
151 nEvents = sTree.GetEntries()
152
153 for iEvent in range(nEvents):
154 if iEvent % 1000 == 0:
155 print("Event ", iEvent)
156
157
158
159 sTree.GetEvent(iEvent)
160
161
162
163 reconstructible_tracks = getReconstructibleTracks(iEvent, sTree, sGeo, ShipGeo)
164
165 metrics["reconstructible"] += len(reconstructible_tracks)
166 for i_reco in reconstructible_tracks:
167 h["TracksPassed"].Fill("Reconstructible tracks", 1)
168 h["TracksPassedU"].Fill("Reconstructible tracks", 1)
169
170 in_y12 = []
171 in_stereo12 = []
172 in_12 = []
173 in_y34 = []
174 in_stereo34 = []
175 in_34 = []
176 in_combo = []
177
178 found_track_ids = []
179 n_tracks = len(reconstructible_tracks)
180 n_recognized = 0
181 n_clones = 0
182 n_ghosts = 0
183 n_others = 0
184
185 min_eff = 0.0
186
187
188
189 nTracklets = sTree.Tracklets.GetEntriesFast()
190
191 for i_track in range(nTracklets):
192 atracklet = sTree.Tracklets[i_track]
193
194 if atracklet.getType() != 1:
195 continue
196
197 atrack = atracklet.getList()
198
199 if len(atrack) == 0:
200 continue
201
202 hits = {
203 "X": [],
204 "Y": [],
205 "Z": [],
206 "DetID": [],
207 "TrackID": [],
208 "Pz": [],
209 "Px": [],
210 "Py": [],
211 "dist2Wire": [],
212 "Pdg": [],
213 }
214
215 for ihit in atrack:
216 ahit = sTree.strawtubesPoint[ihit]
217 hits["X"] += [ahit.GetX()]
218 hits["Y"] += [ahit.GetY()]
219 hits["Z"] += [ahit.GetZ()]
220 hits["DetID"] += [ahit.GetDetectorID()]
221 hits["TrackID"] += [ahit.GetTrackID()]
222 hits["Pz"] += [ahit.GetPz()]
223 hits["Px"] += [ahit.GetPx()]
224 hits["Py"] += [ahit.GetPy()]
225 hits["dist2Wire"] += [ahit.dist2Wire()]
226 hits["Pdg"] += [ahit.PdgCode()]
227
228
229 for key in hits:
230 hits[key] = numpy.array(hits[key])
231
232
233 decode = global_variables.modules["strawtubes"].StrawDecode(hits["DetID"])
234
235 is_stereo = (decode[1] == 1) + (decode[1] == 2)
236 is_y = (decode[1] == 0) + (decode[1] == 3)
237 is_before = (decode[0] == 1) + (decode[0] == 2)
238 is_after = (decode[0] == 3) + (decode[0] == 4)
239
240
241 metrics["n_hits"] += [get_n_hits(hits["TrackID"])]
242
243
244 frac_y12, tmax_y12 = fracMCsame(hits["TrackID"][is_before * is_y])
245 n_hits_y12 = get_n_hits(hits["TrackID"][is_before * is_y])
246 frac_stereo12, tmax_stereo12 = fracMCsame(hits["TrackID"][is_before * is_stereo])
247 n_hits_stereo12 = get_n_hits(hits["TrackID"][is_before * is_stereo])
248 frac_12, tmax_12 = fracMCsame(hits["TrackID"][is_before])
249 n_hits_12 = get_n_hits(hits["TrackID"][is_before])
250
251 frac_y34, tmax_y34 = fracMCsame(hits["TrackID"][is_after * is_y])
252 n_hits_y34 = get_n_hits(hits["TrackID"][is_after * is_y])
253 frac_stereo34, tmax_stereo34 = fracMCsame(hits["TrackID"][is_after * is_stereo])
254 n_hits_stereo34 = get_n_hits(hits["TrackID"][is_after * is_stereo])
255 frac_34, tmax_34 = fracMCsame(hits["TrackID"][is_after])
256 n_hits_34 = get_n_hits(hits["TrackID"][is_after])
257 frac_tot, tmax_tot = fracMCsame(hits["TrackID"])
258 n_hits_tot = get_n_hits(hits["TrackID"])
259
260 if tmax_y12 == tmax_stereo12 and tmax_y12 == tmax_y34 and tmax_y12 == tmax_stereo34:
261 if (
262 frac_y12 >= min_eff
263 and frac_stereo12 >= min_eff
264 and frac_y34 >= min_eff
265 and frac_stereo34 >= min_eff
266 ):
267 if tmax_y12 in reconstructible_tracks and tmax_y12 not in found_track_ids:
268 n_recognized += 1
269 found_track_ids.append(tmax_y12)
270 elif tmax_y12 in reconstructible_tracks and tmax_y12 in found_track_ids:
271 n_clones += 1
272 elif tmax_y12 not in reconstructible_tracks:
273 n_others += 1
274
275 else:
276 n_ghosts += 1
277 else:
278 n_ghosts += 1
279
280 is_reconstructed = 0
281
282 if tmax_y12 in reconstructible_tracks:
283 metrics["passed_y12"] += 1
284 metrics["frac_y12"] += [frac_y12]
285 h["TracksPassed"].Fill("Y view station 1&2", 1)
286 if tmax_y12 not in in_y12:
287 h["TracksPassedU"].Fill("Y view station 1&2", 1)
288 in_y12.append(tmax_y12)
289
290 if tmax_stereo12 == tmax_y12:
291 metrics["passed_stereo12"] += 1
292 metrics["frac_stereo12"] += [frac_stereo12]
293 h["TracksPassed"].Fill("Stereo station 1&2", 1)
294 if tmax_stereo12 not in in_stereo12:
295 h["TracksPassedU"].Fill("Stereo station 1&2", 1)
296 in_stereo12.append(tmax_stereo12)
297
298 if tmax_12 == tmax_y12:
299 metrics["passed_12"] += 1
300 metrics["frac_12"] += [frac_12]
301 h["TracksPassed"].Fill("station 1&2", 1)
302 if tmax_12 not in in_12:
303 h["TracksPassedU"].Fill("station 1&2", 1)
304 in_12.append(tmax_12)
305
306 if tmax_y34 in reconstructible_tracks:
307 metrics["passed_y34"] += 1
308 metrics["frac_y34"] += [frac_y34]
309 h["TracksPassed"].Fill("Y view station 3&4", 1)
310 if tmax_y34 not in in_y34:
311 h["TracksPassedU"].Fill("Y view station 3&4", 1)
312 in_y34.append(tmax_y34)
313
314 if tmax_stereo34 == tmax_y34:
315 metrics["passed_stereo34"] += 1
316 metrics["frac_stereo34"] += [frac_stereo34]
317 h["TracksPassed"].Fill("Stereo station 3&4", 1)
318 if tmax_stereo34 not in in_stereo34:
319 h["TracksPassedU"].Fill("Stereo station 3&4", 1)
320 in_stereo34.append(tmax_stereo34)
321
322 if tmax_34 == tmax_y34:
323 metrics["passed_34"] += 1
324 metrics["frac_34"] += [frac_34]
325 h["TracksPassed"].Fill("station 3&4", 1)
326 if tmax_34 not in in_34:
327 h["TracksPassedU"].Fill("station 3&4", 1)
328 in_34.append(tmax_34)
329
330 if tmax_12 == tmax_34:
331 metrics["passed_combined"] += 1
332 h["TracksPassed"].Fill("Combined stations 1&2/3&4", 1)
333 metrics["reco_passed"] += 1
334 is_reconstructed = 1
335 if tmax_34 not in in_combo:
336 h["TracksPassedU"].Fill("Combined stations 1&2/3&4", 1)
337 metrics["reco_passed_no_clones"] += 1
338 in_combo.append(tmax_34)
339
340
341 if is_reconstructed == 0:
342 continue
343
344 metrics["reco_frac_tot"] += [frac_tot]
345
346
347 hits["Pz"]
348 hits["Px"]
349 hits["Py"]
350
351 p, px, py, pz = getPtruthFirst(sTree, tmax_tot)
352 pt = math.sqrt(px**2 + py**2)
353
354 Z_true = []
355 X_true = []
356 Y_true = []
357 for ahit in sTree.strawtubesPoint:
358 if ahit.GetTrackID() == tmax_tot:
359 az, ax, ay = ahit.GetZ(), ahit.GetX(), ahit.GetY()
360 Z_true.append(az)
361 X_true.append(ax)
362 Y_true.append(ay)
363
364 metrics["reco_mc_p"] += [p]
365 h["TracksPassed_p"].Fill(p, 1)
366
367
368 Z = hits["Z"][(hits["TrackID"] == tmax_tot) * is_before]
369 X = hits["X"][(hits["TrackID"] == tmax_tot) * is_before]
370 Y = hits["Y"][(hits["TrackID"] == tmax_tot) * is_before]
371 Z = Z - Z[0]
372 X = X - X[0]
373 Y = Y - Y[0]
374 R = numpy.sqrt(X**2 + Y**2 + Z**2)
375 Theta = numpy.arccos(Z[1:] / R[1:])
376 theta = numpy.mean(Theta)
377
378 metrics["reco_mc_theta"] += [theta]
379
380 h["n_hits_reco_y12"].Fill(n_hits_y12)
381 h["n_hits_reco_stereo12"].Fill(n_hits_stereo12)
382 h["n_hits_reco_12"].Fill(n_hits_12)
383 h["n_hits_reco_y34"].Fill(n_hits_y34)
384 h["n_hits_reco_stereo34"].Fill(n_hits_stereo34)
385 h["n_hits_reco_34"].Fill(n_hits_34)
386 h["n_hits_reco"].Fill(n_hits_tot)
387
388 h["n_hits_y12"].Fill(p, n_hits_y12)
389 h["n_hits_stereo12"].Fill(p, n_hits_stereo12)
390 h["n_hits_12"].Fill(p, n_hits_12)
391 h["n_hits_y34"].Fill(p, n_hits_y34)
392 h["n_hits_stereo34"].Fill(p, n_hits_stereo34)
393 h["n_hits_34"].Fill(p, n_hits_34)
394 h["n_hits_total"].Fill(p, n_hits_tot)
395
396 h["frac_y12"].Fill(p, frac_y12)
397 h["frac_stereo12"].Fill(p, frac_stereo12)
398 h["frac_12"].Fill(p, frac_12)
399 h["frac_y34"].Fill(p, frac_y34)
400 h["frac_stereo34"].Fill(p, frac_stereo34)
401 h["frac_34"].Fill(p, frac_34)
402 h["frac_total"].Fill(p, frac_tot)
403
404 h["frac_y12_dist"].Fill(frac_y12)
405 h["frac_stereo12_dist"].Fill(frac_stereo12)
406 h["frac_12_dist"].Fill(frac_12)
407 h["frac_y34_dist"].Fill(frac_y34)
408 h["frac_stereo34_dist"].Fill(frac_stereo34)
409 h["frac_34_dist"].Fill(frac_34)
410 h["frac_total_dist"].Fill(frac_tot)
411
412
413
414 thetrack = sTree.FitTracks[i_track]
415
416 fitStatus = thetrack.getFitStatus()
417
418 nmeas = fitStatus.getNdf()
419 pval = fitStatus.getPVal()
420 chi2 = fitStatus.getChi2() / nmeas
421
422 metrics["fitted_pval"] += [pval]
423 metrics["fitted_chi"] += [chi2]
424
425 h["chi2fittedtracks"].Fill(chi2)
426 h["pvalfittedtracks"].Fill(pval)
427
428 try:
429 fittedState = thetrack.getFittedState()
430 fittedMom = fittedState.getMomMag()
431 fittedMom = fittedMom
432
433 px_fit, py_fit, pz_fit = fittedState.getMom().x(), fittedState.getMom().y(), fittedState.getMom().z()
434 p_fit = fittedMom
435 pt_fit = math.sqrt(px_fit**2 + py_fit**2)
436
437 metrics["fitted_p"] += [p_fit]
438 perr = (p - p_fit) / p
439 h["ptrue-p/ptrue"].Fill(perr)
440 h["perr"].Fill(p, perr)
441 h["perr_direction"].Fill(numpy.rad2deg(theta), perr)
442
443 pterr = (pt - pt_fit) / pt
444 h["pttrue-pt/pttrue"].Fill(pterr)
445
446 pxerr = (px - px_fit) / px
447 h["pxtrue-px/pxtrue"].Fill(pxerr)
448
449 pyerr = (py - py_fit) / py
450 h["pytrue-py/pytrue"].Fill(pyerr)
451
452 pzerr = (pz - pz_fit) / pz
453 h["pztrue-pz/pztrue"].Fill(pzerr)
454
455 if math.fabs(p) > 0.0:
456 h["pvspfitted"].Fill(p, fittedMom)
457 fittedtrackDir = fittedState.getDir()
458 fittedx = math.degrees(math.acos(fittedtrackDir[0]))
459 fittedy = math.degrees(math.acos(fittedtrackDir[1]))
460 fittedz = math.degrees(math.acos(fittedtrackDir[2]))
461 fittedmass = fittedState.getMass()
462 h["momentumfittedtracks"].Fill(fittedMom)
463 h["xdirectionfittedtracks"].Fill(fittedx)
464 h["ydirectionfittedtracks"].Fill(fittedy)
465 h["zdirectionfittedtracks"].Fill(fittedz)
466 h["massfittedtracks"].Fill(fittedmass)
467
468 metrics["fitted_x"] += [fittedx]
469 metrics["fitted_y"] += [fittedy]
470 metrics["fitted_z"] += [fittedz]
471 metrics["fitted_mass"] += [fittedmass]
472
473 Z_fit = []
474 X_fit = []
475 Y_fit = []
476 for az in Z_true:
477 _rc, pos, _mom = extrapolateToPlane(thetrack, az)
478 Z_fit.append(pos.Z())
479 X_fit.append(pos.X())
480 Y_fit.append(pos.Y())
481
482 for i in range(len(Z_true)):
483 xerr = abs(X_fit[i] - X_true[i])
484 yerr = abs(Y_fit[i] - Y_true[i])
485 h["abs(x - x-true)"].Fill(xerr)
486 h["abs(y - y-true)"].Fill(yerr)
487
488 rmse_x = numpy.sqrt(numpy.mean((numpy.array(X_fit) - numpy.array(X_true)) ** 2))
489 rmse_y = numpy.sqrt(numpy.mean((numpy.array(Y_fit) - numpy.array(Y_true)) ** 2))
490 h["rmse_x"].Fill(rmse_x)
491 h["rmse_y"].Fill(rmse_y)
492
493 except Exception:
494 print("Problem with fitted state.")
495
496 h["Reco_tracks"].Fill("N total", n_tracks)
497 h["Reco_tracks"].Fill("N recognized tracks", n_recognized)
498 h["Reco_tracks"].Fill("N clones", n_clones)
499 h["Reco_tracks"].Fill("N ghosts", n_ghosts)
500 h["Reco_tracks"].Fill("N others", n_others)
501
502
503
504 save_hists(h, output_file)
505
506 return metrics
507
508
def addVMCFields(shipGeo, str controlFile="", bool verbose=False, bool withVirtualMC=True)
def create_config(str DecayVolumeMedium="helium", float Yheight=6.0, int strawDesign=10, muShieldGeo=None, str shieldName="New_HA_Design", int nuTargetPassive=1, bool SND=True, SND_design=None, TARGET_YAML=None)
def configure(run, ship_geo)
1.9.6