diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 0fc6796e..571ad10c 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -75,7 +75,9 @@ # 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 + # Distance from start to calculate track direction for [Number of planes]. Track matching done in plane pairs -> do not set to 1 + DirectionDistance = 10 [Recon.AStar] #CostMetric = "Manhattan" 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 eadfb7c5..fb0a3503 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,17 +42,17 @@ 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() - 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): @@ -74,10 +75,8 @@ 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) @@ -85,47 +84,74 @@ def draw_performance(out_dir, input_filename): 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] - 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) @@ -229,29 +255,29 @@ 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: 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) @@ -342,14 +368,14 @@ 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]) - 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'); @@ -405,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) @@ -430,9 +456,133 @@ 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_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) + + 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') + 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() + + 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') + 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() + return +def eff_pur_acc(overlap, truth, reco): + 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 + + +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 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..f488a240 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 = -2949 + 0.5 * TMS_Manager::GetInstance().Get_Reco_TRACKMATCH_TiltAngle() * std::abs(UHit.GetNotZ() - VHit.GetNotZ()); //-1350 double X_y = XHit.GetNotZ(); double compared_y = -999999; @@ -851,7 +865,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 +902,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 +943,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 +980,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 +1065,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 +1092,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 +1113,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 +1164,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 +1222,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 +1245,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(); @@ -1247,6 +1261,9 @@ std::vector TMS_TrackFinder::TrackMatching3D() { #endif } } + // Sort track + SpatialPrio(aTrack.Hits); + // Track Length aTrack.Length = CalculateTrackLength3D(aTrack); #ifdef DEBUG @@ -1258,9 +1275,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.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; @@ -1280,7 +1297,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()}; @@ -1290,8 +1307,22 @@ void TMS_TrackFinder::CalculateRecoY(TMS_Hit &OneHit, TMS_Hit &OtherHit) { // 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(-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 + + bool above = false; + if ((UHit.GetNotZ() > 0 && VHit.GetNotZ() > 0) || (UHit.GetNotZ() < 0 && VHit.GetNotZ() < 0)) { + 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())); + } + + //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..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 &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 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);