From 29530ab7f82bd4a75983ba72ce61251f4c7d5234 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 6 May 2024 08:29:14 -0500 Subject: [PATCH 01/14] Reverse RecoY calculation point to the top, but to 244mm this time --- config/TMS_Default_Config.toml | 2 +- src/TMS_Reco.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 0fc6796e..651eaca1 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -75,7 +75,7 @@ # Tilt angle of scintillator planes for calculation of y (tan(90-angle)) TiltAngle = 19.0811366877 #tan(90-3=87) # Y difference for expectation of reconstruction from UV and from X hit [mm] - YDifference = 300.0 + YDifference = 3000.0 [Recon.AStar] #CostMetric = "Manhattan" diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 2216fccc..a83da735 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -736,7 +736,7 @@ double TMS_TrackFinder::CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit) { // TMS_Geom::GetInstance().Scale(positionU[0], positionU[1], positionU[2]); // TMS_Geom::GetInstance().Scale(positionV[0], positionV[1], positionV[2]); - double UV_y = -2949 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); //-1350 + double UV_y = 244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); //-1350/-2949 double X_y = XHit.GetNotZ(); double compared_y = -999999; @@ -1291,7 +1291,7 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &OtherHit) { // TMS_Geom::GetInstance().Scale(positionOne[0], positionOne[1], positionOne[2]); // TMS_Geom::GetInstance().Scale(positionOther[0], positionOther[1], positionOther[2]); - OneHit.SetRecoY(-2949 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(OneHit.GetNotZ() - OtherHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 + OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(OneHit.GetNotZ() - OtherHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 // USING THIS -2949/1350 ALSO IN CompareY FUNCTION. CHANGE THERE AS WELL IF CHANGED HERE!!!!! return; } From 38c28de842e83be43a332f033907e88d2648f512 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 6 May 2024 08:57:48 -0500 Subject: [PATCH 02/14] Preventing a potential logic problem with RecoY calculation by specifying which hit needs to save RecoY exactly --- src/TMS_Reco.cpp | 36 ++++++++++++++++++------------------ src/TMS_Reco.h | 2 +- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index a83da735..7d4f9383 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -851,7 +851,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { //else if (UTracks.front().GetPlaneNumber() >= 20 && std::abs(UTracks.front().GetZ() - XTracks.front().GetZ()) < 90) stereo_view = false; } if (stereo_view) { - CalculateRecoY(UTracks.front(), VTracks.front()); + CalculateRecoY(UTracks.front(), UTracks.front(), VTracks.front()); aTrack.End[0] = 0.5 * (UTracks.front().GetNotZ() + VTracks.front().GetNotZ()); aTrack.End[1] = UTracks.front().GetRecoY(); aTrack.End[2] = UTracks.front().GetZ(); @@ -888,7 +888,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { //else if (VTracks.front().GetPlaneNumber() >= 20 && std::abs(VTracks.front().GetZ() - XTracks.front().GetZ()) < 90) stereo_view = false; } if (stereo_view) { - CalculateRecoY((VTracks.front()), (UTracks.front())); + CalculateRecoY(VTracks.front(), UTracks.front(), VTracks.front()); aTrack.End[0] = 0.5 * (VTracks.front().GetNotZ() + UTracks.front().GetNotZ()); aTrack.End[1] = VTracks.front().GetRecoY(); aTrack.End[2] = VTracks.front().GetZ(); @@ -929,7 +929,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { //else if (VTracks.back().GetPlaneNumber() >= 20 && std::abs(VTracks.back().GetZ() - XTracks.back().GetZ()) < 90) stereo_view = false; } if (stereo_view) { - CalculateRecoY((VTracks.back()), (UTracks.back())); + CalculateRecoY(VTracks.back(), UTracks.back(), VTracks.back()); aTrack.Start[0] = 0.5 * (VTracks.back().GetNotZ() + UTracks.back().GetNotZ()); aTrack.Start[1] = VTracks.back().GetRecoY(); aTrack.Start[2] = VTracks.back().GetZ(); @@ -966,7 +966,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { //else if (UTracks.back().GetPlaneNumber() >= 20 && std::abs(UTracks.back().GetZ() - XTracks.back().GetZ()) < 90) stereo_view = false; } if (stereo_view) { - CalculateRecoY((UTracks.back()), (VTracks.back())); + CalculateRecoY(UTracks.back(), UTracks.back(), VTracks.back()); aTrack.Start[0] = 0.5 * (UTracks.back().GetNotZ() + VTracks.back().GetNotZ()); aTrack.Start[1] = UTracks.back().GetRecoY(); aTrack.Start[2] = UTracks.back().GetZ(); @@ -1051,8 +1051,8 @@ std::vector TMS_TrackFinder::TrackMatching3D() { if ((UTracks[itU]).GetPlaneNumber() == (VTracks[itV]).GetPlaneNumber()) { // Calculate Y info from bar crossing if (stereo_view) { - CalculateRecoY(UTracks[itU], VTracks[itV]); - CalculateRecoY(VTracks[itV], UTracks[itU]); + CalculateRecoY(UTracks[itU], UTracks[itU], VTracks[itV]); + CalculateRecoY(VTracks[itV], UTracks[itU], VTracks[itV]); CalculateRecoX(UTracks[itU], VTracks[itV], UTracks[itU]); CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]); #ifdef DEBUG @@ -1078,8 +1078,8 @@ std::vector TMS_TrackFinder::TrackMatching3D() { #endif if (itX > 0) --itX; } else { - CalculateRecoY(UTracks[itU], VTracks[itV]); - CalculateRecoY(VTracks[itV], UTracks[itU]); + CalculateRecoY(UTracks[itU], UTracks[itU], VTracks[itV]); + CalculateRecoY(VTracks[itV], UTracks[itU], VTracks[itV]); CalculateRecoX(UTracks[itU], VTracks[itV], UTracks[itU]); CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]); #ifdef DEBUG @@ -1099,17 +1099,17 @@ std::vector TMS_TrackFinder::TrackMatching3D() { std::cout << "Hit U: " << UTracks[itU].GetNotZ() << " | " << UTracks[itU].GetZ() << " / V: " << VTracks[itV].GetNotZ() << " | " << VTracks[itV].GetZ() << std::endl; #endif if (itU > 0 && itV > 0) { - CalculateRecoY(VTracks[itV], UTracks[itU - 1]); + CalculateRecoY(VTracks[itV], UTracks[itU - 1], VTracks[itV]); CalculateRecoX(UTracks[itU - 1], VTracks[itV], VTracks[itV]); (aTrack.Hits).push_back(VTracks[itV]); --itV; } else if (itU == 0 && itV > 0) { - CalculateRecoY(VTracks[itV], UTracks[itU]); + CalculateRecoY(VTracks[itV], UTracks[itU], VTracks[itU]); CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]); (aTrack.Hits).push_back(VTracks[itV]); --itV; } else if (itU > 0 && itV == 0) { - CalculateRecoY(UTracks[itU], VTracks[itV]); + CalculateRecoY(UTracks[itU], UTracks[itU], VTracks[itV]); CalculateRecoX(UTracks[itU], VTracks[itV], UTracks[itU]); (aTrack.Hits).push_back(UTracks[itU]); --itU; @@ -1150,17 +1150,17 @@ std::vector TMS_TrackFinder::TrackMatching3D() { std::cout << "Hit U: " << UTracks[itU].GetNotZ() << " | " << UTracks[itU].GetZ() << " / V: " << VTracks[itV].GetNotZ() << " | " << VTracks[itV].GetZ() << std::endl; #endif if (itV > 0 && itU > 0) { - CalculateRecoY(UTracks[itU], VTracks[itV - 1]); + CalculateRecoY(UTracks[itU], UTracks[itU], VTracks[itV - 1]); CalculateRecoX(UTracks[itU], VTracks[itV - 1], UTracks[itU]); (aTrack.Hits).push_back(UTracks[itU]); --itU; } else if (itV == 0 && itU > 0) { - CalculateRecoY(UTracks[itU], VTracks[itV]); + CalculateRecoY(UTracks[itU], UTracks[itU], VTracks[itV]); CalculateRecoX(UTracks[itU], VTracks[itV], UTracks[itU]); (aTrack.Hits).push_back(UTracks[itU]); --itU; } else if (itV > 0 && itU == 0) { - CalculateRecoY(VTracks[itV], UTracks[itU]); + CalculateRecoY(VTracks[itV], UTracks[itU], VTracks[itV]); CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]); (aTrack.Hits).push_back(VTracks[itV]); --itV; @@ -1208,7 +1208,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { else if (UTracks.front().GetPlaneNumber() >= 20 && std::abs(UTracks.front().GetZ() - XTracks.front().GetZ() < 90)) stereo_view = false; } if (stereo_view) { - CalculateRecoY(VTracks.front(), UTracks.front()); + CalculateRecoY(VTracks.front(), UTracks.front(), VTracks.front()); aTrack.End[0] = 0.5 * (VTracks.front().GetNotZ() + UTracks.front().GetNotZ()); aTrack.End[1] = VTracks.front().GetRecoY(); aTrack.End[2] = VTracks.front().GetZ(); @@ -1231,7 +1231,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { else if (UTracks.back().GetPlaneNumber() >= 20 && std::abs(UTracks.back().GetZ() - XTracks.back().GetZ() < 90)) stereo_view = false; } if (stereo_view) { - CalculateRecoY(UTracks.back(), VTracks.back()); + CalculateRecoY(UTracks.back(), UTracks.back(), VTracks.back()); aTrack.Start[0] = 0.5 * (VTracks.back().GetNotZ() + UTracks.back().GetNotZ()); aTrack.Start[1] = VTracks.back().GetRecoY(); aTrack.Start[2] = VTracks.back().GetZ(); @@ -1280,7 +1280,7 @@ std::vector TMS_TrackFinder::TrackMatching3D() { return returned; } -void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &OtherHit) { +void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &UHit, TMS_Hit &VHit) { // TGeoManager *geom = TMS_Geom::GetInstance().GetGeometry(); // Double_t localOne[3] = {OneHit.GetNotZ(), 0, OneHit.GetZ()}; // Double_t localOther[3] = {OtherHit.GetNotZ(), 0, OtherHit.GetZ()}; @@ -1291,7 +1291,7 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &OtherHit) { // TMS_Geom::GetInstance().Scale(positionOne[0], positionOne[1], positionOne[2]); // TMS_Geom::GetInstance().Scale(positionOther[0], positionOther[1], positionOther[2]); - OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(OneHit.GetNotZ() - OtherHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 + OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 // USING THIS -2949/1350 ALSO IN CompareY FUNCTION. CHANGE THERE AS WELL IF CHANGED HERE!!!!! return; } diff --git a/src/TMS_Reco.h b/src/TMS_Reco.h index f1480491..de4c4da6 100644 --- a/src/TMS_Reco.h +++ b/src/TMS_Reco.h @@ -263,7 +263,7 @@ class TMS_TrackFinder { std::vector Extrapolation(const std::vector &TrackHits, const std::vector &Hits); std::vector TrackMatching3D(); void FindPseudoXTrack(); - void CalculateRecoY(TMS_Hit &UHit, TMS_Hit &VHit); + void CalculateRecoY(TMS_Hit &XHit, TMS_Hit &UHit, TMS_Hit &VHit); //XHit is here a placeholder for the hit that gets the RecoY calculated void CalculateRecoX(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit); double CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit); From f5345233daa380c34bb05299186ac5216e726cb4 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 7 May 2024 11:51:43 -0500 Subject: [PATCH 03/14] Updating the performance_reco script with exiting vs. stopping evaluation --- config/TMS_Default_Config.toml | 2 +- scripts/Reco/performance_reco.py | 143 ++++++++++++++++++++++++++++--- 2 files changed, 134 insertions(+), 11 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 651eaca1..2dc98279 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -110,7 +110,7 @@ # Draw PDF of "event display". Slows down reco considerably [Applications] DrawPDF = false - MaximumNEvents = -1 + MaximumNEvents = 10000 # Variables that have to do with the geometry or location of the detector [Geometry] diff --git a/scripts/Reco/performance_reco.py b/scripts/Reco/performance_reco.py index eadfb7c5..776dc777 100755 --- a/scripts/Reco/performance_reco.py +++ b/scripts/Reco/performance_reco.py @@ -5,6 +5,7 @@ import argparse import cppyy.ll import matplotlib.cm as cm +from matplotlib_venn import venn2, venn2_circles # plotstyle red_cbf = '#d55e00' @@ -41,7 +42,7 @@ def draw_performance(out_dir, input_filename): if not truth.GetEntries() > 0: print("Didn't get any entries in Truth_Info, are you sure the input_filename is right?\n", input_filename) - max_n_spills = 10000 # TODO (old) add some meta info to output file with n spill info for file + max_n_spills = 121000 # TODO (old) add some meta info to output file with n spill info for file spill_number_cache = dict() n_events = r.GetEntries() @@ -82,8 +83,8 @@ def draw_performance(out_dir, input_filename): EndPos = np.frombuffer(event.EndPos, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackStart = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackStart, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) - Muon_TrueTrackLength = true_event.Muon_TrueTrackLength - Reco_Track_Length = np.frombuffer(event.Length, dtype = np.float32) + #Muon_TrueTrackLength = true_event.Muon_TrueTrackLength + #Reco_Track_Length = np.frombuffer(event.Length, dtype = np.float32) if RecoTrackPrimaryParticleTruePositionTrackStart[0] > -8000. and not StartPos.size == 0: # checking for muon tracks (minimal length for this are 20 planes traversed -> 890 mm in thin area @@ -91,11 +92,13 @@ def draw_performance(out_dir, input_filename): Reco_Start[i, 0] = StartPos[0] Reco_Start[i, 1] = StartPos[1] Reco_Start[i, 2] = StartPos[2] + #if RecoTrackPrimaryParticleTruePositionTrackStart[2] < 11362: + # print('PDG: ', true_event.PDG[np.frombuffer(true_event.RecoTrackPrimaryParticleIndex, dtype = np.uint8)[0]], ' End: ', RecoTrackPrimaryParticleTruePositionTrackEnd[2]) Primary_True_Start[i, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[0] Primary_True_Start[i, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[1] Primary_True_Start[i, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[2] - True_TrackLength[i] = Muon_TrueTrackLength - Reco_TrackLength[i] = Reco_Track_Length + #True_TrackLength[i] = Muon_TrueTrackLength + #Reco_TrackLength[i] = Reco_Track_Length if RecoTrackPrimaryParticleTruePositionTrackEnd[0] > -8000. and not EndPos.size == 0: if (EndPos[2] - StartPos[2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[2] - RecoTrackPrimaryParticleTruePositionTrackStart[2]) > 890.: Reco_End[i, 0] = EndPos[0] @@ -119,6 +122,9 @@ def draw_performance(out_dir, input_filename): boolean_Reco_TrackLength = (Reco_TrackLength != -9999.) Reco_TrackLength = Reco_TrackLength[boolean_Reco_TrackLength] + # total number of events after filtering + print("#events reconstruction: ", len(Reco_Start), "# events truth: ", len(Primary_True_Start)) + # subtract reconstruction from truth for all directions Diff_Start_x = Primary_True_Start[:, 0] - Reco_Start[:, 0] Diff_Start_y = Primary_True_Start[:, 1] - Reco_Start[:, 1]# - 915 @@ -229,12 +235,12 @@ def draw_performance(out_dir, input_filename): mp.title('End z', fontsize = 'x-large') mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_End_Z_hist.png' % out_dir, bbox_inches = 'tight') - mp.close() + mp.close() # create histograms for all directions for truth and reconstruction - z_bins = np.zeros(101, dtype = float) + z_bins = np.zeros(102, dtype = float) z_bins[0] = 11362 - for i in range(1, 101): + for i in range(1, 102): if i < 40: z_bins[i] = z_bins[0] + i * 55. if i >= 40: @@ -342,7 +348,7 @@ def draw_performance(out_dir, input_filename): ### plot difference in dependence of hit position # create 2d histograms - print(min(min(Primary_True_Start[:, 1]), min(Primary_True_End[:, 1])), max(max(Primary_True_Start[:, 1]), max(Primary_True_End[:, 1])), min(min(Reco_Start[:, 1]), min(Reco_End[:, 1])), max(max(Reco_Start[:, 1]), max(Reco_End[:, 1]))) + #print(min(min(Primary_True_Start[:, 1]), min(Primary_True_End[:, 1])), max(max(Primary_True_Start[:, 1]), max(Primary_True_End[:, 1])), min(min(Reco_Start[:, 1]), min(Reco_End[:, 1])), max(max(Reco_Start[:, 1]), max(Reco_End[:, 1]))) dependence_Start_x_hist, dependence_Start_x_binsX, dependence_Start_x_binsY = np.histogram2d(Primary_True_Start[:, 0], Diff_Start_x, bins = [Primary_True_Start_x_bins, Diff_Start_x_bins]) dependence_Start_y_hist, dependence_Start_y_binsX, dependence_Start_y_binsY = np.histogram2d(Primary_True_Start[:, 1], Diff_Start_y, bins = [Primary_True_Start_y_bins, Diff_Start_y_bins]) @@ -430,9 +436,126 @@ def draw_performance(out_dir, input_filename): mp.colorbar(im) mp.savefig('%s/Tracklength.png' % out_dir, bbox_inches = 'tight') mp.close() - + + + ### Stopping vs exiting evaluation + opposite_direction_counter = 0 + stopping_and_exiting_counter = 0 + stopping_counter = 0 + exiting_counter = 0 + exiting_truth = 0 + exiting_reco = 0 + stopping_truth = 0 + stopping_reco = 0 + + for i in range(len(Reco_Start)): + true_slope, true_intercept = line_slope_intercept(Primary_True_Start[i, 1], Primary_True_End[i, 1], Primary_True_Start[i, 2], Primary_True_End[i, 2]) + reco_slope, reco_intercept = line_slope_intercept(Reco_Start[i, 1], Reco_End[i, 1], Reco_Start[i, 2], Reco_End[i, 2]) + + truth_exiting = False + truth_stopping = False + reco_exiting = False + reco_stopping = False + + # check if line crosses detector volume in yz + if (true_slope * 18294 + true_intercept) < -2949 and Primary_True_End[i, 1] < (-2949 + 338.3): + truth_exiting = True + elif Primary_True_End[i, 2] < 17750: + truth_stopping = True + if (reco_slope * 18294 + reco_intercept) < -2949 and Reco_End[i, 1] < (-2949 + 338.3): + reco_exiting = True + elif Reco_End[i, 2] < 17750: + reco_stopping = True + + if truth_exiting: + exiting_truth += 1 + elif truth_stopping: + stopping_truth += 1 + + if reco_exiting: + exiting_reco += 1 + elif reco_stopping: + stopping_reco += 1 + + if truth_exiting and reco_exiting: + stopping_and_exiting_counter += 1 + exiting_counter += 1 + elif truth_stopping and reco_stopping: + stopping_and_exiting_counter += 1 + stopping_counter += 1 + + if (true_slope < 0 and reco_slope > 0) or (true_slope > 0 and reco_slope < 0): + opposite_direction_counter += 1 + + efficiency_stop, purity_stop, accuracy_stop = eff_pur_acc(stopping_counter, stopping_truth, stopping_reco) + efficiency_exit, purity_exit, accuracy_exit = eff_pur_acc(exiting_counter, exiting_truth, exiting_reco) + + + print("Number of opposite directions: ", opposite_direction_counter, " (", opposite_direction_counter / len(Reco_Start) * 100, "%)") + print("Stopping and exiting total: ", stopping_and_exiting_counter, " (", stopping_and_exiting_counter / len(Reco_Start) * 100, "%)") + print("Exiting: ", exiting_counter, " (", exiting_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_stop * 100, " , pur.: ", purity_stop * 100, " , acc.: ", accuracy_stop * 100) + print("Stopping: ", stopping_counter, " (", stopping_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_exit * 100, " , pur.: ", purity_exit * 100, " , acc.: ", accuracy_exit * 100) + print("Truth: exiting: ", exiting_truth, " stopping: ", stopping_truth) + print("Reco: exiting: ", exiting_reco, " stopping: ", stopping_reco) + + v = venn2(subsets = (stopping_truth - stopping_counter, stopping_reco - stopping_counter, stopping_counter), set_labels = ('Truth', 'Reco'), set_colors = ('#03558e', '#e95125', '#f39e92'), alpha = 1) + c = venn2_circles(subsets = (stopping_truth - stopping_counter, stopping_reco - stopping_counter, stopping_counter), color = 'grey', linestyle = 'dashed', linewidth = 1.5) + v.get_label_by_id('10').set_text('%i' % stopping_truth) + v.get_label_by_id('10').set_color('lightgrey') + v.get_label_by_id('10').set_fontsize('x-large') + v.get_label_by_id('A').set_color('#03558e') + v.get_label_by_id('A').set_fontsize('xx-large') + v.get_patch_by_id('10').set_color('#03558e') + v.get_label_by_id('01').set_text('%i' % stopping_reco) + v.get_label_by_id('01').set_color('black') + v.get_label_by_id('01').set_fontsize('x-large') + v.get_label_by_id('B').set_color('#e95125') + v.get_label_by_id('B').set_fontsize('xx-large') + v.get_patch_by_id('01').set_color('#e95125') + v.get_label_by_id('11').set_text('%i' % stopping_counter) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') + mp.title('Stopping') + mp.savefig('%s/Stopping_venn.png' % out_dir, bbox_inches = 'tight') + mp.close() + + v = venn2(subsets = (exiting_truth - exiting_counter, exiting_reco - exiting_counter, exiting_counter), set_labels = ('Truth', 'Reco'), set_colors = ('#03558e', '#e95125', '#f39e92'), alpha = 1) + c = venn2_circles(subsets = (exiting_truth - exiting_counter, exiting_reco - exiting_counter, exiting_counter), color = 'grey', linestyle = 'dashed', linewidth = 1.5) + v.get_label_by_id('10').set_text('%i' % exiting_truth) + v.get_label_by_id('10').set_color('lightgrey') + v.get_label_by_id('10').set_fontsize('x-large') + v.get_label_by_id('A').set_color('#03558e') + v.get_label_by_id('A').set_fontsize('xx-large') + v.get_patch_by_id('10').set_color('#03558e') + v.get_label_by_id('01').set_text('%i' % exiting_reco) + v.get_label_by_id('01').set_color('black') + v.get_label_by_id('01').set_fontsize('x-large') + v.get_label_by_id('B').set_color('#e95125') + v.get_label_by_id('B').set_fontsize('xx-large') + v.get_patch_by_id('01').set_color('#e95125') + v.get_label_by_id('11').set_text('%i' % exiting_counter) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') + mp.title('Exiting') + mp.savefig('%s/Exiting_venn.png' % out_dir, bbox_inches = 'tight') + mp.close() + return +def eff_pur_acc(overlap, truth, reco): + efficiency = overlap / truth + purity = overlap / reco + accuracy = overlap / (truth + reco - overlap) + return efficiency, purity, accuracy + + +def line_slope_intercept(y_start, y_end, z_start, z_end): + slope = (y_end - y_start) / (z_end - z_start) + intercept = y_start - slope * z_start + return slope, intercept + def histogram_arr_handle(bins, edges):#(u bins,v edges) """Function to calculated x- and y-values from np.histogram From 33cab11cc150075fa11302e9d33a3bc0b6ae971b Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 7 May 2024 11:52:11 -0500 Subject: [PATCH 04/14] Taking out changes in config file --- config/TMS_Default_Config.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 2dc98279..651eaca1 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -110,7 +110,7 @@ # Draw PDF of "event display". Slows down reco considerably [Applications] DrawPDF = false - MaximumNEvents = 10000 + MaximumNEvents = -1 # Variables that have to do with the geometry or location of the detector [Geometry] From c293a1a59e54bc91c5f6ffc944069e10718aac7a Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 13 May 2024 07:57:31 -0500 Subject: [PATCH 05/14] Hopefully y calculation this time correct, changed X anchor point to middle --- src/TMS_Reco.cpp | 52 ++++++++++++++++++++++++++++++++++-------- src/TMS_TreeWriter.cpp | 2 +- 2 files changed, 43 insertions(+), 11 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 7d4f9383..5c4d7c2a 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -735,8 +735,22 @@ double TMS_TrackFinder::CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit) { // geom->GetCurrentMatrix()->LocalToMaster(localV, positionV); // TMS_Geom::GetInstance().Scale(positionU[0], positionU[1], positionU[2]); // TMS_Geom::GetInstance().Scale(positionV[0], positionV[1], positionV[2]); + double UV_y = 2000; + bool above = false; + if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { + if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = true; + else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = false; + } + else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = false; + else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = true; + + if (above) { + UV_y = TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); + } else { + UV_y = TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); + } + //double UV_y = 244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); //-1350/-2949 - double UV_y = 244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); //-1350/-2949 double X_y = XHit.GetNotZ(); double compared_y = -999999; @@ -1290,8 +1304,22 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &UHit, TMS_Hit &VH // geom->GetCurrentMatrix()->LocalToMaster(localOther, positionOther); // TMS_Geom::GetInstance().Scale(positionOne[0], positionOne[1], positionOne[2]); // TMS_Geom::GetInstance().Scale(positionOther[0], positionOther[1], positionOther[2]); - - OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 + + bool above = false; + if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { + if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = true; + else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = false; + } + else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = false; + else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = true; + + if (above) { + OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); + } else { + OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); + } + + //OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 // USING THIS -2949/1350 ALSO IN CompareY FUNCTION. CHANGE THERE AS WELL IF CHANGED HERE!!!!! return; } @@ -1308,16 +1336,20 @@ void TMS_TrackFinder::EvaluateTrackFinding(TMS_Event &event) { // Get the true muon std::vector TrueParticles = event.GetTrueParticles(); TMS_TrueParticle muon; - bool FoundMuon = false; + //bool FoundMuon = false; for (auto &particle: TrueParticles) { - if (particle.GetPDG() == 13 && - particle.GetParent() == -1) { + //if (particle.GetPDG() == 13 && + // particle.GetParent() == -1) { muon = particle; - FoundMuon = true; - } + // FoundMuon = true; + //} } - if (!FoundMuon) return; + for (auto &TrueHit: muon.GetPositionPoints()) { + std::cout << "Truth: " << TrueHit.X() << " | " << TrueHit.Y() << " | " << TrueHit.Z() << std::endl; + } + + /*if (!FoundMuon) return; // Now compare this to our tracks and see how many of the muon hits are included // Loop over the candidate tracks that have been built @@ -1371,7 +1403,7 @@ void TMS_TrackFinder::EvaluateTrackFinding(TMS_Event &event) { if (efficiency > 0) { Efficiency->Fill(muonke, efficiency); Total->Fill(muonke); - } + }*/ } // Calculate the total track energy diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 727a764f..bea29bf5 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -1083,7 +1083,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { } for (unsigned int j = 0; j < RecoTrack->Hits.size(); ++j) { if (RecoTrack->Hits[j].GetBar().GetBarType() != TMS_Bar::kXBar) { - RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetRecoX(); + RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetNotZ();//GetRecoX(); RecoTrackHitPos[itTrack][j][1] = RecoTrack->Hits[j].GetRecoY(); } else if (RecoTrack->Hits[j].GetBar().GetBarType() == TMS_Bar::kXBar) { RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetRecoX(); From 25d6d6b96ab6654d65400712b575c9342abbc358 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 13 May 2024 08:00:44 -0500 Subject: [PATCH 06/14] Reinstating GetRecoX for Reco_Tree --- src/TMS_TreeWriter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index bea29bf5..727a764f 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -1083,7 +1083,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { } for (unsigned int j = 0; j < RecoTrack->Hits.size(); ++j) { if (RecoTrack->Hits[j].GetBar().GetBarType() != TMS_Bar::kXBar) { - RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetNotZ();//GetRecoX(); + RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetRecoX(); RecoTrackHitPos[itTrack][j][1] = RecoTrack->Hits[j].GetRecoY(); } else if (RecoTrack->Hits[j].GetBar().GetBarType() == TMS_Bar::kXBar) { RecoTrackHitPos[itTrack][j][0] = RecoTrack->Hits[j].GetRecoX(); From dcdcbd79c47212bfecedcb239c9c6a51b6a8363c Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 13 May 2024 08:12:32 -0500 Subject: [PATCH 07/14] Correcting above-below check --- src/TMS_Reco.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 5c4d7c2a..5c75f267 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -738,11 +738,11 @@ double TMS_TrackFinder::CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit) { double UV_y = 2000; bool above = false; if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { - if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = true; - else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = false; + if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = false; + else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = true; } - else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = false; - else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = true; + else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = true; + else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = false; if (above) { UV_y = TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); @@ -1307,11 +1307,11 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &UHit, TMS_Hit &VH bool above = false; if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { - if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = true; - else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = false; + if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = false; + else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = true; } - else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = false; - else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = true; + else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = true; + else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = false; if (above) { OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); From ad125a921d888fe65568b2faf44730432b30f383 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Mon, 13 May 2024 10:14:45 -0500 Subject: [PATCH 08/14] Changing direction again. Matches better this way. Why though? --- src/TMS_Reco.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 5c75f267..7f27f6fb 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -738,11 +738,11 @@ double TMS_TrackFinder::CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit) { double UV_y = 2000; bool above = false; if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { - if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = false; - else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = true; + if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = true; + else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = false; } - else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = true; - else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = false; + else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = false; + else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = true; if (above) { UV_y = TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); @@ -1307,16 +1307,16 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &UHit, TMS_Hit &VH bool above = false; if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { - if (std::abs(UHit.GetNotZ()) >= std::abs(VHit.GetNotZ())) above = false; - else if (std::abs(UHit.GetNotZ()) < std::abs(VHit.GetNotZ())) above = true; + if (UHit.GetNotZ() >= VHit.GetNotZ()) above = true; + else if (UHit.GetNotZ() < VHit.GetNotZ()) above = false; } else if (UHit.GetNotZ() > 0 && VHit.GetNotZ() < 0) above = true; else if (UHit.GetNotZ() < 0 && VHit.GetNotZ() > 0) above = false; if (above) { - OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); - } else { OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); + } else { + OneHit.SetRecoY(TMS_Manager::GetInstance().Get_Geometry_YMIDDLE() * 1000 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); } //OneHit.SetRecoY(244 - 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ())); //positionOne[1], positionOther[1] //0.5 * (OneHit.GetY() + OtherHit.GetY())//-1350 From ae2a58dd9be77314d4a8aab55aeee1b78bb7a63f Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 04:32:45 -0500 Subject: [PATCH 09/14] Implement direction calculation and add parameter in config file to allow for tuning of distance to calculate --- config/TMS_Default_Config.toml | 2 ++ src/TMS_Manager.cpp | 1 + src/TMS_Manager.h | 2 ++ src/TMS_Reco.cpp | 6 +++--- 4 files changed, 8 insertions(+), 3 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 0fc6796e..80d77932 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -76,6 +76,8 @@ TiltAngle = 19.0811366877 #tan(90-3=87) # Y difference for expectation of reconstruction from UV and from X hit [mm] YDifference = 300.0 + # Distance from start to calculate track direction for [Number of planes]. Track matching done in plane pairs -> do not set to 1 + DirectionDistance = 2 [Recon.AStar] #CostMetric = "Manhattan" diff --git a/src/TMS_Manager.cpp b/src/TMS_Manager.cpp index e9346035..3224dda3 100644 --- a/src/TMS_Manager.cpp +++ b/src/TMS_Manager.cpp @@ -57,6 +57,7 @@ TMS_Manager::TMS_Manager() { _RECO_TRACKMATCH_YAnchor = toml::find(data, "Recon", "TrackMatch3D", "YAnchor"); _RECO_TRACKMATCH_TiltAngle = toml::find(data, "Recon", "TrackMatch3D", "TiltAngle"); _RECO_TRACKMATCH_YDifference = toml::find(data, "Recon", "TrackMatch3D", "YDifference"); + _RECO_TRACKMATCH_DirectionDistance = toml::find(data, "Recon", "TrackMatch3D", "DirectionDistance"); _RECO_ASTAR_IsGreedy = toml::find (data, "Recon", "AStar", "IsGreedy"); _RECO_ASTAR_CostMetric = toml::find (data, "Recon", "AStar", "CostMetric"); diff --git a/src/TMS_Manager.h b/src/TMS_Manager.h index 9b013f6f..1ced3388 100644 --- a/src/TMS_Manager.h +++ b/src/TMS_Manager.h @@ -51,6 +51,7 @@ class TMS_Manager { float Get_Reco_TRACKMATCH_YAnchor() { return _RECO_TRACKMATCH_YAnchor; }; double Get_Reco_TRACKMATCH_TiltAngle() { return _RECO_TRACKMATCH_TiltAngle; }; float Get_Reco_TRACKMATCH_YDifference() { return _RECO_TRACKMATCH_YDifference; }; + int Get_Reco_TRACKMATCH_DirectionDistance() { return _RECO_TRACKMATCH_DirectionDistance; }; bool Get_Reco_ASTAR_IsGreedy() { return _RECO_ASTAR_IsGreedy; }; std::string Get_Reco_ASTAR_CostMetric() { return _RECO_ASTAR_CostMetric; }; @@ -118,6 +119,7 @@ class TMS_Manager { float _RECO_TRACKMATCH_YAnchor; double _RECO_TRACKMATCH_TiltAngle; float _RECO_TRACKMATCH_YDifference; + int _RECO_TRACKMATCH_DirectionDistance; bool _RECO_ASTAR_IsGreedy; std::string _RECO_ASTAR_CostMetric; diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 2216fccc..0259e96d 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -1258,9 +1258,9 @@ std::vector TMS_TrackFinder::TrackMatching3D() { std::cout << "Added TrackEnergyDeposit: " << aTrack.EnergyDeposit << std::endl; #endif // Track Direction - aTrack.Direction[0] = aTrack.End[0] - aTrack.Start[0]; - aTrack.Direction[1] = aTrack.End[1] - aTrack.Start[1]; - aTrack.Direction[2] = aTrack.End[2] - aTrack.Start[2]; + aTrack.Direction[0] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoX() - aTrack.Start[0]; //aTrack.End[0] - aTrack.Start[0]; + aTrack.Direction[1] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoY() - aTrack.Start[1]; //aTrack.End[1] - aTrack.Start[1]; + aTrack.Direction[2] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetNotZ() - aTrack.Start[2]; //aTrack.End[2] - aTrack.Start[2]; #ifdef DEBUG std::cout << "Start: " << aTrack.Start[0] << " | " << aTrack.Start[1] << " | " << aTrack.Start[2] << std::endl; std::cout << "End: " << aTrack.End[0] << " | " << aTrack.End[1] << " | " << aTrack.End[2] << std::endl; From 76069f5a05328d26642ea5d0c853650384ad6586 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 04:47:29 -0500 Subject: [PATCH 10/14] Get actually z for z coordinate --- config/TMS_Default_Config.toml | 2 +- src/TMS_Reco.cpp | 9 ++++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 80d77932..0c137da9 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -77,7 +77,7 @@ # Y difference for expectation of reconstruction from UV and from X hit [mm] YDifference = 300.0 # Distance from start to calculate track direction for [Number of planes]. Track matching done in plane pairs -> do not set to 1 - DirectionDistance = 2 + DirectionDistance = 20 [Recon.AStar] #CostMetric = "Manhattan" diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 0259e96d..f66df15a 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -1247,6 +1247,9 @@ std::vector TMS_TrackFinder::TrackMatching3D() { #endif } } + // Sort track + SpatialPrio(aTrack.Hits); + // Track Length aTrack.Length = CalculateTrackLength3D(aTrack); #ifdef DEBUG @@ -1258,9 +1261,9 @@ std::vector TMS_TrackFinder::TrackMatching3D() { std::cout << "Added TrackEnergyDeposit: " << aTrack.EnergyDeposit << std::endl; #endif // Track Direction - aTrack.Direction[0] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoX() - aTrack.Start[0]; //aTrack.End[0] - aTrack.Start[0]; - aTrack.Direction[1] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoY() - aTrack.Start[1]; //aTrack.End[1] - aTrack.Start[1]; - aTrack.Direction[2] = aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetNotZ() - aTrack.Start[2]; //aTrack.End[2] - aTrack.Start[2]; + aTrack.Direction[0] = aTrack.Start[0] - aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoX(); //aTrack.End[0] - aTrack.Start[0]; + aTrack.Direction[1] = aTrack.Start[1] - aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetRecoY(); //aTrack.End[1] - aTrack.Start[1]; + aTrack.Direction[2] = aTrack.Start[2] - aTrack.Hits[TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_DirectionDistance()].GetZ(); //aTrack.End[2] - aTrack.Start[2]; #ifdef DEBUG std::cout << "Start: " << aTrack.Start[0] << " | " << aTrack.Start[1] << " | " << aTrack.Start[2] << std::endl; std::cout << "End: " << aTrack.End[0] << " | " << aTrack.End[1] << " | " << aTrack.End[2] << std::endl; From 33f54b0ee370071b6502a0e409cd84867e543980 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 04:48:05 -0500 Subject: [PATCH 11/14] Set DirectionDistance parameter to 10 --- config/TMS_Default_Config.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 0c137da9..2eb1bd05 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -77,7 +77,7 @@ # Y difference for expectation of reconstruction from UV and from X hit [mm] YDifference = 300.0 # Distance from start to calculate track direction for [Number of planes]. Track matching done in plane pairs -> do not set to 1 - DirectionDistance = 20 + DirectionDistance = 10 [Recon.AStar] #CostMetric = "Manhattan" From a975342b57f2890ff18ce206fc0c74f1f80ff626 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 08:36:54 -0500 Subject: [PATCH 12/14] Add multiple tracks per spill option in performance plots --- .../Reco/compare_2_data_sets_performance.py | 118 ++++++------ scripts/Reco/performance_reco.py | 175 ++++++++++-------- 2 files changed, 165 insertions(+), 128 deletions(-) diff --git a/scripts/Reco/compare_2_data_sets_performance.py b/scripts/Reco/compare_2_data_sets_performance.py index 48602e21..5a2b1863 100755 --- a/scripts/Reco/compare_2_data_sets_performance.py +++ b/scripts/Reco/compare_2_data_sets_performance.py @@ -61,15 +61,15 @@ def draw_comparison(out_dir, input_filename1, input_filename2): n_events_2 = r2.GetEntries() - Reco_Start_1 = np.ones((n_events_1, 3), dtype = float) * -9999. - Reco_End_1 = np.ones((n_events_1, 3), dtype = float) * -9999. - Primary_True_Start_1 = np.ones((n_events_1, 3), dtype = float) * -9999. - Primary_True_End_1 = np.ones((n_events_1, 3), dtype = float) * -9999. + Reco_Start_1 = np.ones((n_events_1, 5, 3), dtype = float) * -9999. + Reco_End_1 = np.ones((n_events_1, 5, 3), dtype = float) * -9999. + Primary_True_Start_1 = np.ones((n_events_1, 5, 3), dtype = float) * -9999. + Primary_True_End_1 = np.ones((n_events_1, 5, 3), dtype = float) * -9999. - Reco_Start_2 = np.ones((n_events_2, 3), dtype = float) * -9999. - Reco_End_2 = np.ones((n_events_2, 3), dtype = float) * -9999. - Primary_True_Start_2 = np.ones((n_events_2, 3), dtype = float) * -9999. - Primary_True_End_2 = np.ones((n_events_2, 3), dtype = float) * -9999. + Reco_Start_2 = np.ones((n_events_2, 5, 3), dtype = float) * -9999. + Reco_End_2 = np.ones((n_events_2, 5, 3), dtype = float) * -9999. + Primary_True_Start_2 = np.ones((n_events_2, 5, 3), dtype = float) * -9999. + Primary_True_End_2 = np.ones((n_events_2, 5, 3), dtype = float) * -9999. for current_spill_number in range(max_n_spills): for i in range(n_events_1): @@ -92,29 +92,34 @@ def draw_comparison(out_dir, input_filename1, input_filename2): if true_event == None: truth1.GetEntry(i) true_event = truth1 - - nTracks = event.nTracks - if nTracks <= 0: continue StartPos = np.frombuffer(event.StartPos, dtype = np.float32) EndPos = np.frombuffer(event.EndPos, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackStart = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackStart, dtype = np.float32) - RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) + RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) + + nTracks = event.nTracks + if nTracks <= 0: continue + if nTracks > 4: print("Too many tracks in event. Limit to first 5") + for j in range(nTracks): + if j > 4: break - if RecoTrackPrimaryParticleTruePositionTrackStart[0] > -8000. and not StartPos.size == 0: - Reco_Start_1[i, 0] = StartPos[0] - Reco_Start_1[i, 1] = StartPos[1] - Reco_Start_1[i, 2] = StartPos[2] - Primary_True_Start_1[i, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[0] - Primary_True_Start_1[i, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[1] - Primary_True_Start_1[i, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[2] - if RecoTrackPrimaryParticleTruePositionTrackEnd[0] > -8000. and not EndPos.size == 0: - Reco_End_1[i, 0] = EndPos[0] - Reco_End_1[i, 1] = EndPos[1] - Reco_End_1[i, 2] = EndPos[2] - Primary_True_End_1[i, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[0] - Primary_True_End_1[i, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[1] - Primary_True_End_1[i, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[2] + if RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] > -8000. and not StartPos.size == 0: + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_Start_1[i, j, 0] = StartPos[j*3 + 0] + Reco_Start_1[i, j, 1] = StartPos[j*3 + 1] + Reco_Start_1[i, j, 2] = StartPos[j*3 + 2] + Primary_True_Start_1[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] + Primary_True_Start_1[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 1] + Primary_True_Start_1[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2] + if RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] > -8000. and not EndPos.size == 0: + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_End_1[i, j, 0] = EndPos[j*3 + 0] + Reco_End_1[i, j, 1] = EndPos[j*3 + 1] + Reco_End_1[i, j, 2] = EndPos[j*3 + 2] + Primary_True_End_1[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] + Primary_True_End_1[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 1] + Primary_True_End_1[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] for i in range(n_events_2): try: @@ -136,47 +141,52 @@ def draw_comparison(out_dir, input_filename1, input_filename2): if true_event2 == None: truth2.GetEntry(i) true_event2 = truth2 - - nTracks = event2.nTracks - if nTracks <= 0: continue StartPos = np.frombuffer(event2.StartPos, dtype = np.float32) EndPos = np.frombuffer(event2.EndPos, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackStart = np.frombuffer(true_event2.RecoTrackPrimaryParticleTruePositionTrackStart, dtype = np.float32) - RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event2.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) + RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event2.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) + + nTracks = event2.nTracks + if nTracks <= 0: continue + if nTracks > 4: print("Too many tracks in event. Limit to first 5") + for j in range(nTracks): + if j > 4: break - if RecoTrackPrimaryParticleTruePositionTrackStart[0] > -8000. and not StartPos.size == 0: - Reco_Start_2[i, 0] = StartPos[0] - Reco_Start_2[i, 1] = StartPos[1] - Reco_Start_2[i, 2] = StartPos[2] - Primary_True_Start_2[i, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[0] - Primary_True_Start_2[i, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[1] - Primary_True_Start_2[i, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[2] - if RecoTrackPrimaryParticleTruePositionTrackEnd[0] > -8000. and not EndPos.size == 0: - Reco_End_2[i, 0] = EndPos[0] - Reco_End_2[i, 1] = EndPos[1] - Reco_End_2[i, 2] = EndPos[2] - Primary_True_End_2[i, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[0] - Primary_True_End_2[i, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[1] - Primary_True_End_2[i, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[2] + if RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] > -8000. and not StartPos.size == 0: + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_Start_2[i, j, 0] = StartPos[j*3 + 0] + Reco_Start_2[i, j, 1] = StartPos[j*3 + 1] + Reco_Start_2[i, j, 2] = StartPos[j*3 + 2] + Primary_True_Start_2[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] + Primary_True_Start_2[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 1] + Primary_True_Start_2[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2] + if RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] > -8000. and not EndPos.size == 0: + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_End_2[i, j, 0] = EndPos[j*3 + 0] + Reco_End_2[i, j, 1] = EndPos[j*3 + 1] + Reco_End_2[i, j, 2] = EndPos[j*3 + 2] + Primary_True_End_2[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] + Primary_True_End_2[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 1] + Primary_True_End_2[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] # filter out not filled indice - boolean_Reco_Start_1 = (Reco_Start_1[:, 0] != -9999.) + boolean_Reco_Start_1 = (Reco_Start_1[:, :, 0] != -9999.) Reco_Start_1 = Reco_Start_1[boolean_Reco_Start_1] - boolean_Reco_End_1 = (Reco_End_1[:, 0] != -9999.) + boolean_Reco_End_1 = (Reco_End_1[:, :, 0] != -9999.) Reco_End_1 = Reco_End_1[boolean_Reco_End_1] - boolean_Primary_Start_1 = (Primary_True_Start_1[:, 0] != -9999.) + boolean_Primary_Start_1 = (Primary_True_Start_1[:, :, 0] != -9999.) Primary_True_Start_1 = Primary_True_Start_1[boolean_Primary_Start_1] - boolean_Primary_End_1 = (Primary_True_End_1[:, 0] != -9999.) + boolean_Primary_End_1 = (Primary_True_End_1[:, :, 0] != -9999.) Primary_True_End_1 = Primary_True_End_1[boolean_Primary_End_1] - boolean_Reco_Start_2 = (Reco_Start_2[:, 0] != -9999.) + boolean_Reco_Start_2 = (Reco_Start_2[:, :, 0] != -9999.) Reco_Start_2 = Reco_Start_2[boolean_Reco_Start_2] - boolean_Reco_End_2 = (Reco_End_2[:, 0] != -9999.) + boolean_Reco_End_2 = (Reco_End_2[:, :, 0] != -9999.) Reco_End_2 = Reco_End_2[boolean_Reco_End_2] - boolean_Primary_Start_2 = (Primary_True_Start_2[:, 0] != -9999.) + boolean_Primary_Start_2 = (Primary_True_Start_2[:, :, 0] != -9999.) Primary_True_Start_2 = Primary_True_Start_2[boolean_Primary_Start_2] - boolean_Primary_End_2 = (Primary_True_End_2[:, 0] != -9999.) + boolean_Primary_End_2 = (Primary_True_End_2[:, :, 0] != -9999.) Primary_True_End_2 = Primary_True_End_2[boolean_Primary_End_2] # subtract reconstruction from truth for all directions @@ -299,7 +309,7 @@ def draw_comparison(out_dir, input_filename1, input_filename2): # create histograms for all directions for truth and reconstruction - Reco_Start_x_1_hist, Reco_Start_x_1_bins = np.histogram(Reco_Start_1[:, 0], bins = 3, range = (-35.42, 35.42)) #int((max(Reco_Start_1[:, 0]) - min(Reco_Start_1[:, 0])) / 35.42)) #48 + Reco_Start_x_1_hist, Reco_Start_x_1_bins = np.histogram(Reco_Start_1[:, 0], bins = int((max(Reco_Start_1[:, 0]) - min(Reco_Start_1[:, 0])) / 35.42))#3, range = (-35.42, 35.42)) # #48 Reco_Start_y_1_hist, Reco_Start_y_1_bins = np.histogram(Reco_Start_1[:, 1], bins = 40)#, range = (-4110., -910.)) Reco_Start_z_1_hist, Reco_Start_z_1_bins = np.histogram(Reco_Start_1[:, 2], bins = 100)#, range = (11360., 18320.)) Reco_End_x_1_hist, Reco_End_x_1_bins = np.histogram(Reco_End_1[:, 0], bins = int((max(Reco_End_1[:, 0]) - min(Reco_End_1[:, 0])) / 35.42))#48)#, range = (-3400., 3400.)) @@ -312,7 +322,7 @@ def draw_comparison(out_dir, input_filename1, input_filename2): Primary_True_End_y_1_hist, Primary_True_End_y_1_bins = np.histogram(Primary_True_End_1[:, 1], bins = Reco_End_y_1_bins)#20) Primary_True_End_z_1_hist, Primary_True_End_z_1_bins = np.histogram(Primary_True_End_1[:, 2], bins = 100) - Reco_Start_x_2_hist, Reco_Start_x_2_bins = np.histogram(Reco_Start_2[:, 0], bins = 3, range = (-35.42, 35.42))#10)#int((max(Reco_Start_2[:, 0]) - min(Reco_Start_2[:, 0])) / 35.42)) #48)#, range = (-3400., 3400.)) + Reco_Start_x_2_hist, Reco_Start_x_2_bins = np.histogram(Reco_Start_2[:, 0], bins = int((max(Reco_Start_2[:, 0]) - min(Reco_Start_2[:, 0])) / 35.42))#3, range = (-35.42, 35.42))#10)# #48)# Reco_Start_y_2_hist, Reco_Start_y_2_bins = np.histogram(Reco_Start_2[:, 1], bins = 40)#, range = (-4110., -910.)) Reco_Start_z_2_hist, Reco_Start_z_2_bins = np.histogram(Reco_Start_2[:, 2], bins = 100)#, range = (11360., 18320.)) Reco_End_x_2_hist, Reco_End_x_2_bins = np.histogram(Reco_End_2[:, 0], bins = int((max(Reco_End_2[:, 0]) - min(Reco_End_2[:, 0])) / 35.42)) #48)#, range = (-3400., 3400.)) diff --git a/scripts/Reco/performance_reco.py b/scripts/Reco/performance_reco.py index 776dc777..fb0a3503 100755 --- a/scripts/Reco/performance_reco.py +++ b/scripts/Reco/performance_reco.py @@ -47,12 +47,12 @@ def draw_performance(out_dir, input_filename): spill_number_cache = dict() n_events = r.GetEntries() - Reco_Start = np.ones((n_events, 3), dtype = float) * -9999. - Reco_End = np.ones((n_events, 3), dtype = float) * -9999. - Primary_True_Start = np.ones((n_events, 3), dtype = float) * -9999. - Primary_True_End = np.ones((n_events, 3), dtype = float) * -9999. - True_TrackLength = np.ones(n_events, dtype = float) * -9999. - Reco_TrackLength = np.ones(n_events, dtype = float) * -9999. + Reco_Start = np.ones((n_events, 5, 3), dtype = float) * -9999. + Reco_End = np.ones((n_events, 5, 3), dtype = float) * -9999. + Primary_True_Start = np.ones((n_events, 5, 3), dtype = float) * -9999. + Primary_True_End = np.ones((n_events, 5, 3), dtype = float) * -9999. + True_TrackLength = np.ones((n_events, 5), dtype = float) * -9999. + Reco_TrackLength = np.ones((n_events, 5), dtype = float) * -9999. for current_spill_number in range(max_n_spills): for i in range(n_events): @@ -75,63 +75,83 @@ def draw_performance(out_dir, input_filename): if true_event == None: truth.GetEntry(i) true_event = truth - - nTracks = event.nTracks - if nTracks <= 0: continue + StartPos = np.frombuffer(event.StartPos, dtype = np.float32) EndPos = np.frombuffer(event.EndPos, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackStart = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackStart, dtype = np.float32) RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) - #Muon_TrueTrackLength = true_event.Muon_TrueTrackLength - #Reco_Track_Length = np.frombuffer(event.Length, dtype = np.float32) + Muon_TrueTrackLength = true_event.Muon_TrueTrackLength + Reco_Track_Length = np.frombuffer(event.Length, dtype = np.float32) - if RecoTrackPrimaryParticleTruePositionTrackStart[0] > -8000. and not StartPos.size == 0: - # checking for muon tracks (minimal length for this are 20 planes traversed -> 890 mm in thin area - if (EndPos[2] - StartPos[2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[2] - RecoTrackPrimaryParticleTruePositionTrackStart[2]) > 890.: - Reco_Start[i, 0] = StartPos[0] - Reco_Start[i, 1] = StartPos[1] - Reco_Start[i, 2] = StartPos[2] - #if RecoTrackPrimaryParticleTruePositionTrackStart[2] < 11362: - # print('PDG: ', true_event.PDG[np.frombuffer(true_event.RecoTrackPrimaryParticleIndex, dtype = np.uint8)[0]], ' End: ', RecoTrackPrimaryParticleTruePositionTrackEnd[2]) - Primary_True_Start[i, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[0] - Primary_True_Start[i, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[1] - Primary_True_Start[i, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[2] - #True_TrackLength[i] = Muon_TrueTrackLength - #Reco_TrackLength[i] = Reco_Track_Length - if RecoTrackPrimaryParticleTruePositionTrackEnd[0] > -8000. and not EndPos.size == 0: - if (EndPos[2] - StartPos[2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[2] - RecoTrackPrimaryParticleTruePositionTrackStart[2]) > 890.: - Reco_End[i, 0] = EndPos[0] - Reco_End[i, 1] = EndPos[1] - Reco_End[i, 2] = EndPos[2] - Primary_True_End[i, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[0] - Primary_True_End[i, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[1] - Primary_True_End[i, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[2] + nTracks = event.nTracks + if nTracks <= 0: continue + if nTracks > 4: print("Too many tracks in event. Limit to first 5") + for j in range(nTracks): + if j > 4: break + + if RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] > -8000. and not StartPos.size == 0: + # checking for muon tracks (minimal length for this are 20 planes traversed -> 890 mm in thin area + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_Start[i, j, 0] = StartPos[j*3 + 0] + Reco_Start[i, j, 1] = StartPos[j*3 + 1] + Reco_Start[i, j, 2] = StartPos[j*3 + 2] + #if RecoTrackPrimaryParticleTruePositionTrackStart[2] < 11362: + # print('PDG: ', true_event.PDG[np.frombuffer(true_event.RecoTrackPrimaryParticleIndex, dtype = np.uint8)[0]], ' End: ', RecoTrackPrimaryParticleTruePositionTrackEnd[2]) + Primary_True_Start[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] + Primary_True_Start[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 1] + Primary_True_Start[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2] + True_TrackLength[i, j] = Muon_TrueTrackLength + Reco_TrackLength[i, j] = Reco_Track_Length[j] + if RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] > -8000. and not EndPos.size == 0: + if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: + Reco_End[i, j, 0] = EndPos[j*3 + 0] + Reco_End[i, j, 1] = EndPos[j*3 + 1] + Reco_End[i, j, 2] = EndPos[j*3 + 2] + Primary_True_End[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] + Primary_True_End[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 1] + Primary_True_End[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] # filter out not filled indice - boolean_Reco_Start = (Reco_Start[:, 0] != -9999.) + boolean_Reco_Start = (Reco_Start[:, :, 0] != -9999.) Reco_Start = Reco_Start[boolean_Reco_Start] - boolean_Reco_End = (Reco_End[:, 0] != -9999.) + boolean_Reco_End = (Reco_End[:, :, 0] != -9999.) Reco_End = Reco_End[boolean_Reco_End] - boolean_Primary_Start = (Primary_True_Start[:, 0] != -9999.) + boolean_Primary_Start = (Primary_True_Start[:, :, 0] != -9999.) Primary_True_Start = Primary_True_Start[boolean_Primary_Start] - boolean_Primary_End = (Primary_True_End[:, 0] != -9999.) + boolean_Primary_End = (Primary_True_End[:, :, 0] != -9999.) Primary_True_End = Primary_True_End[boolean_Primary_End] boolean_True_TrackLength = (True_TrackLength != -9999.) True_TrackLength = True_TrackLength[boolean_True_TrackLength] boolean_Reco_TrackLength = (Reco_TrackLength != -9999.) Reco_TrackLength = Reco_TrackLength[boolean_Reco_TrackLength] + # flatten arrays + Reco_Start_x = Reco_Start[:, 0]#.flatten() + Reco_Start_y = Reco_Start[:, 1]#.flatten() + Reco_Start_z = Reco_Start[:, 2]#.flatten() + Reco_End_x = Reco_End[:, 0]#.flatten() + Reco_End_y = Reco_End[:, 1]#.flatten() + Reco_End_z = Reco_End[:, 2]#.flatten() + Primary_True_Start_x = Primary_True_Start[:, 0]#.flatten() + Primary_True_Start_y = Primary_True_Start[:, 1]#.flatten() + Primary_True_Start_z = Primary_True_Start[:, 2]#.flatten() + Primary_True_End_x = Primary_True_End[:, 0]#.flatten() + Primary_True_End_y = Primary_True_End[:, 1]#.flatten() + Primary_True_End_z = Primary_True_End[:, 2]#.flatten() + #True_TrackLength.flatten() + #Reco_TrackLength.flatten() + # total number of events after filtering print("#events reconstruction: ", len(Reco_Start), "# events truth: ", len(Primary_True_Start)) # subtract reconstruction from truth for all directions - Diff_Start_x = Primary_True_Start[:, 0] - Reco_Start[:, 0] - Diff_Start_y = Primary_True_Start[:, 1] - Reco_Start[:, 1]# - 915 - Diff_Start_z = Primary_True_Start[:, 2] - Reco_Start[:, 2] - Diff_End_x = Primary_True_End[:, 0] - Reco_End[:, 0] - Diff_End_y = Primary_True_End[:, 1] - Reco_End[:, 1]# - 915 - Diff_End_z = Primary_True_End[:, 2] - Reco_End[:, 2] + Diff_Start_x = Primary_True_Start_x - Reco_Start_x + Diff_Start_y = Primary_True_Start_y - Reco_Start_y + Diff_Start_z = Primary_True_Start_z - Reco_Start_z + Diff_End_x = Primary_True_End_x - Reco_End_x + Diff_End_y = Primary_True_End_y - Reco_End_y + Diff_End_z = Primary_True_End_z - Reco_End_z # create histograms for the differences Diff_Start_x_hist, Diff_Start_x_bins = np.histogram(Diff_Start_x, bins = 50) @@ -246,18 +266,18 @@ def draw_performance(out_dir, input_filename): if i >= 40: z_bins[i] = 13507. + (i - 40) * 80. - Reco_Start_x_hist, Reco_Start_x_bins = np.histogram(Reco_Start[:, 0], bins = 48, range = (-3520., 3520.)) - Reco_Start_y_hist, Reco_Start_y_bins = np.histogram(Reco_Start[:, 1], bins = 20, range = (-2949., 244.)) #+915 - Reco_Start_z_hist, Reco_Start_z_bins = np.histogram(Reco_Start[:, 2], bins = z_bins) - Reco_End_x_hist, Reco_End_x_bins = np.histogram(Reco_End[:, 0], bins = Reco_Start_x_bins) - Reco_End_y_hist, Reco_End_y_bins = np.histogram(Reco_End[:, 1], bins = Reco_Start_y_bins) #+915 - Reco_End_z_hist, Reco_End_z_bins = np.histogram(Reco_End[:, 2], bins = Reco_Start_z_bins) - Primary_True_Start_x_hist, Primary_True_Start_x_bins = np.histogram(Primary_True_Start[:, 0], bins = Reco_Start_x_bins) - Primary_True_Start_y_hist, Primary_True_Start_y_bins = np.histogram(Primary_True_Start[:, 1], bins = Reco_Start_y_bins) - Primary_True_Start_z_hist, Primary_True_Start_z_bins = np.histogram(Primary_True_Start[:, 2], bins = Reco_Start_z_bins) - Primary_True_End_x_hist, Primary_True_End_x_bins = np.histogram(Primary_True_End[:, 0], bins = Reco_Start_x_bins) - Primary_True_End_y_hist, Primary_True_End_y_bins = np.histogram(Primary_True_End[:, 1], bins = Reco_Start_y_bins) - Primary_True_End_z_hist, Primary_True_End_z_bins = np.histogram(Primary_True_End[:, 2], bins = Reco_Start_z_bins) + Reco_Start_x_hist, Reco_Start_x_bins = np.histogram(Reco_Start_x, bins = 48, range = (-3520., 3520.)) + Reco_Start_y_hist, Reco_Start_y_bins = np.histogram(Reco_Start_y, bins = 20, range = (-2949., 244.)) + Reco_Start_z_hist, Reco_Start_z_bins = np.histogram(Reco_Start_z, bins = z_bins) + Reco_End_x_hist, Reco_End_x_bins = np.histogram(Reco_End_x, bins = Reco_Start_x_bins) + Reco_End_y_hist, Reco_End_y_bins = np.histogram(Reco_End_y, bins = Reco_Start_y_bins) + Reco_End_z_hist, Reco_End_z_bins = np.histogram(Reco_End_z, bins = Reco_Start_z_bins) + Primary_True_Start_x_hist, Primary_True_Start_x_bins = np.histogram(Primary_True_Start_x, bins = Reco_Start_x_bins) + Primary_True_Start_y_hist, Primary_True_Start_y_bins = np.histogram(Primary_True_Start_y, bins = Reco_Start_y_bins) + Primary_True_Start_z_hist, Primary_True_Start_z_bins = np.histogram(Primary_True_Start_z, bins = Reco_Start_z_bins) + Primary_True_End_x_hist, Primary_True_End_x_bins = np.histogram(Primary_True_End_x, bins = Reco_Start_x_bins) + Primary_True_End_y_hist, Primary_True_End_y_bins = np.histogram(Primary_True_End_y, bins = Reco_Start_y_bins) + Primary_True_End_z_hist, Primary_True_End_z_bins = np.histogram(Primary_True_End_z, bins = Reco_Start_z_bins) # make the histograms usable Reco_Start_x_histX, Reco_Start_x_histY = histogram_arr_handle(Reco_Start_x_hist, Reco_Start_x_bins) @@ -350,12 +370,12 @@ def draw_performance(out_dir, input_filename): # create 2d histograms #print(min(min(Primary_True_Start[:, 1]), min(Primary_True_End[:, 1])), max(max(Primary_True_Start[:, 1]), max(Primary_True_End[:, 1])), min(min(Reco_Start[:, 1]), min(Reco_End[:, 1])), max(max(Reco_Start[:, 1]), max(Reco_End[:, 1]))) - dependence_Start_x_hist, dependence_Start_x_binsX, dependence_Start_x_binsY = np.histogram2d(Primary_True_Start[:, 0], Diff_Start_x, bins = [Primary_True_Start_x_bins, Diff_Start_x_bins]) - dependence_Start_y_hist, dependence_Start_y_binsX, dependence_Start_y_binsY = np.histogram2d(Primary_True_Start[:, 1], Diff_Start_y, bins = [Primary_True_Start_y_bins, Diff_Start_y_bins]) - dependence_Start_z_hist, dependence_Start_z_binsX, dependence_Start_z_binsY = np.histogram2d(Primary_True_Start[:, 2], Diff_Start_z, bins = [Primary_True_Start_z_bins, Diff_Start_z_bins]) - dependence_End_x_hist, dependence_End_x_binsX, dependence_End_x_binsY = np.histogram2d(Primary_True_End[:, 0], Diff_End_x, bins = [Primary_True_End_x_bins, Diff_End_x_bins]) - dependence_End_y_hist, dependence_End_y_binsX, dependence_End_y_binsY = np.histogram2d(Primary_True_End[:, 1], Diff_End_y, bins = [Primary_True_End_y_bins, Diff_End_y_bins]) - dependence_End_z_hist, dependence_End_z_binsX, dependence_End_z_binsY = np.histogram2d(Primary_True_End[:, 2], Diff_End_z, bins = [Primary_True_End_z_bins, Diff_End_z_bins]) + dependence_Start_x_hist, dependence_Start_x_binsX, dependence_Start_x_binsY = np.histogram2d(Primary_True_Start_x, Diff_Start_x, bins = [Primary_True_Start_x_bins, Diff_Start_x_bins]) + dependence_Start_y_hist, dependence_Start_y_binsX, dependence_Start_y_binsY = np.histogram2d(Primary_True_Start_y, Diff_Start_y, bins = [Primary_True_Start_y_bins, Diff_Start_y_bins]) + dependence_Start_z_hist, dependence_Start_z_binsX, dependence_Start_z_binsY = np.histogram2d(Primary_True_Start_z, Diff_Start_z, bins = [Primary_True_Start_z_bins, Diff_Start_z_bins]) + dependence_End_x_hist, dependence_End_x_binsX, dependence_End_x_binsY = np.histogram2d(Primary_True_End_x, Diff_End_x, bins = [Primary_True_End_x_bins, Diff_End_x_bins]) + dependence_End_y_hist, dependence_End_y_binsX, dependence_End_y_binsY = np.histogram2d(Primary_True_End_y, Diff_End_y, bins = [Primary_True_End_y_bins, Diff_End_y_bins]) + dependence_End_z_hist, dependence_End_z_binsX, dependence_End_z_binsY = np.histogram2d(Primary_True_End_z, Diff_End_z, bins = [Primary_True_End_z_bins, Diff_End_z_bins]) cmap = cm.get_cmap('cividis'); @@ -411,7 +431,7 @@ def draw_performance(out_dir, input_filename): # histogram # single histogram as well - reco_tracklength_hist, reco_tracklength_bins = np.histogram(Reco_TrackLength / 10, bins = 40, range = (min(min(Reco_TrackLength / 10), min(True_TrackLength)), max(max(Reco_TrackLength / 10), max(True_TrackLength)))) + reco_tracklength_hist, reco_tracklength_bins = np.histogram(Reco_TrackLength / 10, bins = 40, range = (50, 500))#(min(min(Reco_TrackLength / 10), min(True_TrackLength)), max(max(Reco_TrackLength / 10), max(True_TrackLength)))) true_tracklength_hist, true_tracklength_bins = np.histogram(True_TrackLength, bins = reco_tracklength_bins) reco_tracklength_histX, reco_tracklength_histY = histogram_arr_handle(reco_tracklength_hist, reco_tracklength_bins) @@ -493,8 +513,8 @@ def draw_performance(out_dir, input_filename): print("Number of opposite directions: ", opposite_direction_counter, " (", opposite_direction_counter / len(Reco_Start) * 100, "%)") print("Stopping and exiting total: ", stopping_and_exiting_counter, " (", stopping_and_exiting_counter / len(Reco_Start) * 100, "%)") - print("Exiting: ", exiting_counter, " (", exiting_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_stop * 100, " , pur.: ", purity_stop * 100, " , acc.: ", accuracy_stop * 100) - print("Stopping: ", stopping_counter, " (", stopping_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_exit * 100, " , pur.: ", purity_exit * 100, " , acc.: ", accuracy_exit * 100) + print("Exiting: ", exiting_counter, " (", exiting_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_exit * 100, " , pur.: ", purity_exit * 100, " , acc.: ", accuracy_exit * 100) + print("Stopping: ", stopping_counter, " (", stopping_counter / len(Reco_Start) * 100, "%) | eff.: ", efficiency_stop * 100, " , pur.: ", purity_stop * 100, " , acc.: ", accuracy_stop * 100) print("Truth: exiting: ", exiting_truth, " stopping: ", stopping_truth) print("Reco: exiting: ", exiting_reco, " stopping: ", stopping_reco) @@ -512,10 +532,11 @@ def draw_performance(out_dir, input_filename): v.get_label_by_id('B').set_color('#e95125') v.get_label_by_id('B').set_fontsize('xx-large') v.get_patch_by_id('01').set_color('#e95125') - v.get_label_by_id('11').set_text('%i' % stopping_counter) - v.get_label_by_id('11').set_color('black') - v.get_label_by_id('11').set_fontsize('x-large') - v.get_patch_by_id('11').set_color('#f39e92') + if stopping_counter > 0: + v.get_label_by_id('11').set_text('%i' % stopping_counter) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') mp.title('Stopping') mp.savefig('%s/Stopping_venn.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -534,10 +555,11 @@ def draw_performance(out_dir, input_filename): v.get_label_by_id('B').set_color('#e95125') v.get_label_by_id('B').set_fontsize('xx-large') v.get_patch_by_id('01').set_color('#e95125') - v.get_label_by_id('11').set_text('%i' % exiting_counter) - v.get_label_by_id('11').set_color('black') - v.get_label_by_id('11').set_fontsize('x-large') - v.get_patch_by_id('11').set_color('#f39e92') + if exiting_counter > 0: + v.get_label_by_id('11').set_text('%i' % exiting_counter) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') mp.title('Exiting') mp.savefig('%s/Exiting_venn.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -545,9 +567,14 @@ def draw_performance(out_dir, input_filename): return def eff_pur_acc(overlap, truth, reco): - efficiency = overlap / truth - purity = overlap / reco - accuracy = overlap / (truth + reco - overlap) + if truth == 0 or reco == 0: + efficiency = 0 + purity = 0 + accuracy = 0 + else: + efficiency = overlap / truth + purity = overlap / reco + accuracy = overlap / (truth + reco - overlap) return efficiency, purity, accuracy From 37c10b74edd9e62ac0fa064ad629aeb0dbe64e26 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 08:55:33 -0500 Subject: [PATCH 13/14] Making CalculateRecoY explanation clearer --- src/TMS_Reco.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TMS_Reco.h b/src/TMS_Reco.h index de4c4da6..f3c93e1a 100644 --- a/src/TMS_Reco.h +++ b/src/TMS_Reco.h @@ -263,7 +263,7 @@ class TMS_TrackFinder { std::vector Extrapolation(const std::vector &TrackHits, const std::vector &Hits); std::vector TrackMatching3D(); void FindPseudoXTrack(); - void CalculateRecoY(TMS_Hit &XHit, TMS_Hit &UHit, TMS_Hit &VHit); //XHit is here a placeholder for the hit that gets the RecoY calculated + void CalculateRecoY(TMS_Hit &XHit, TMS_Hit &UHit, TMS_Hit &VHit); //XHit is here a placeholder for the hit that gets the RecoY set. This is then to be chosen as either UHit or VHit void CalculateRecoX(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit); double CompareY(TMS_Hit &UHit, TMS_Hit &VHit, TMS_Hit &XHit); From 7c8a679f0c1749603cb0d76f932358cdb2af7c76 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 14 May 2024 08:58:03 -0500 Subject: [PATCH 14/14] Restoring EvaluateTrackFinding --- src/TMS_Reco.cpp | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 063e6871..f488a240 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -1339,20 +1339,16 @@ void TMS_TrackFinder::EvaluateTrackFinding(TMS_Event &event) { // Get the true muon std::vector TrueParticles = event.GetTrueParticles(); TMS_TrueParticle muon; - //bool FoundMuon = false; + bool FoundMuon = false; for (auto &particle: TrueParticles) { - //if (particle.GetPDG() == 13 && - // particle.GetParent() == -1) { + if (particle.GetPDG() == 13 && + particle.GetParent() == -1) { muon = particle; - // FoundMuon = true; - //} - } - - for (auto &TrueHit: muon.GetPositionPoints()) { - std::cout << "Truth: " << TrueHit.X() << " | " << TrueHit.Y() << " | " << TrueHit.Z() << std::endl; + FoundMuon = true; + } } - /*if (!FoundMuon) return; + if (!FoundMuon) return; // Now compare this to our tracks and see how many of the muon hits are included // Loop over the candidate tracks that have been built @@ -1406,7 +1402,7 @@ void TMS_TrackFinder::EvaluateTrackFinding(TMS_Event &event) { if (efficiency > 0) { Efficiency->Fill(muonke, efficiency); Total->Fill(muonke); - }*/ + } } // Calculate the total track energy