From ea43cfecceb05607bbfea39cda8703bcb6f6a8ef Mon Sep 17 00:00:00 2001 From: Michael Dolce Date: Tue, 27 Aug 2024 11:15:50 -0400 Subject: [PATCH 1/5] Getting 'RecoTrackNHits' doesn't work with the TMS Reco file (from April), and this is the latest file made. Need to wait for a new TMS Reco file to be made, then it will have this info. Then you can un-comment this block. --- scripts/Reco/draw_spill_3D_projections.py | 26 +++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/scripts/Reco/draw_spill_3D_projections.py b/scripts/Reco/draw_spill_3D_projections.py index ea56bcf8..ce9d6838 100644 --- a/scripts/Reco/draw_spill_3D_projections.py +++ b/scripts/Reco/draw_spill_3D_projections.py @@ -264,7 +264,10 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ #print("number of hits: ", nHits) TrackHitPos = np.frombuffer(event.TrackHitPos, dtype = np.float32) - + # I want the number of true hits. + # nTrueHits = true_event.RecoTrackNHits + # True_Hits = np.frombuffer(true_event.RecoTrackTrueHitPosition, dtype=np.float32) + StartPos = np.frombuffer(event.StartPos, dtype = np.float32) EndPos = np.frombuffer(event.EndPos, dtype = np.float32) #Direction = np.frombuffer(event.Direction, dtype = np.float32) @@ -273,7 +276,26 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) for i in range(nTracks): - + + # This isn't working with current file, Aug 2024. File doesn't have true_event.RecoTrackNHits + # # loop through those true hits in the TMS + # for hit in range(nTrueHits[i]): + # thit_x = True_Hits[i * 600 + hit * 4 + 0] + # thit_y = True_Hits[i * 600 + hit * 4 + 1] + # thit_z = True_Hits[i * 600 + hit * 4 + 2] + # + # # print("THit", thit_x, thit_y, thit_z) + # + # if thit_z < 11000.: continue + # + # x_z.scatter(thit_z / 1000.0, thit_x / 1000.0, color=black_cbf, marker='.', alpha=0.5) + # z_y.scatter(thit_z / 1000.0, thit_y / 1000.0, color=black_cbf, marker='.', alpha=0.5) + # x_y.scatter(thit_x / 1000.0, thit_y / 1000.0, color=black_cbf, marker='.', alpha=0.5) + # + # output_filename_thits = os.path.join(out_dir, f"{name}_truth_{current_spill_number:03d}") + # mp.savefig(output_filename_thits + ".png", bbox_inches='tight') + # mp.close() + ### Hit positions of the hits in the track for hit in range(nHits[i]): From b5755e11e6d3958a90f32778ca91344ae7a142c9 Mon Sep 17 00:00:00 2001 From: Michael Dolce Date: Wed, 28 Aug 2024 17:32:03 -0400 Subject: [PATCH 2/5] the index is wrong, need one over truth.PDG and not hits, and calc KE. --- scripts/Reco/draw_spill_3D_projections.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/scripts/Reco/draw_spill_3D_projections.py b/scripts/Reco/draw_spill_3D_projections.py index ce9d6838..26faa644 100644 --- a/scripts/Reco/draw_spill_3D_projections.py +++ b/scripts/Reco/draw_spill_3D_projections.py @@ -4,6 +4,7 @@ import os import argparse import cppyy.ll +from math import sqrt # plotstyle red_cbf = '#d55e00' @@ -28,6 +29,8 @@ delta_y = 0.3383 # uncertainty from +/-3 degree tilted bars delta_z = 0.02 # space of scintilattor with air gap +MUON_MASS = 105.7 # MeV/c^2 + ### Function for upper limit of tilted bar 'hit' def upper_limit(hit_x, hit_y, x, orientation_bar): if orientation_bar == 'kVBar': # assumption VBar is tilted in positive way and UBar then in negative @@ -367,6 +370,19 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ x_y.scatter(RecoTrackPrimaryParticleTruePositionTrackEnd[i*4 + 0] / 1000.0, RecoTrackPrimaryParticleTruePositionTrackEnd[i*4 + 1] / 1000.0, c = green_cbf, marker = '1', alpha = 0.5) + # Write the True Muon KE to each spill plot. + if only_true_tms_muons: + for idx, pdg in enumerate(true_event.PDG): + print('index: ', idx, 'pdg: ', pdg) + if pdg != abs(13): continue + + muon_ke_lar = true_event.Muon_TrueKE + p_tms_start = ROOT.TVector3(truth.MomentumTMSStart[4 * idx], truth.MomentumTMSStart[4 * idx + 1], truth.MomentumTMSStart[4 * idx + 2]) + muon_ke_tms_start = sqrt(p_tms_start.Mag2() + MUON_MASS ** 2) - MUON_MASS + + x_z.text(3.6, -2.5, f'Muon KE at LAr: {muon_ke_lar}', rotation = 'vertical', fontsize = 12, fontweight = 'bold', color = orange_cbf) + x_z.text(3.6, -2.75, f'Muon KE entering TMS: {muon_ke_tms_start}', rotation= 'vertical', fontsize = 12, fontweight = 'bold', color = orange_cbf) + output_filename = os.path.join(out_dir, f"{name}_{current_spill_number:03d}") mp.savefig(output_filename + ".png", bbox_inches = 'tight') mp.close() From 8b20095be7cc356a70d70b14150e0377d48944c9 Mon Sep 17 00:00:00 2001 From: Michael Dolce Date: Wed, 28 Aug 2024 17:54:09 -0400 Subject: [PATCH 3/5] Let's add print statement so we know which events are above 5 GeV and record their event and spill No. (Might be helpful to redirect output to a text file). And adjust positioning of the KE truth text, rotate and divide to get GeV, two decimal places. And had to change the style, b/c of error that seaborn-poster wasn't available --- scripts/Reco/draw_spill_3D_projections.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/scripts/Reco/draw_spill_3D_projections.py b/scripts/Reco/draw_spill_3D_projections.py index 26faa644..45bf9ae2 100644 --- a/scripts/Reco/draw_spill_3D_projections.py +++ b/scripts/Reco/draw_spill_3D_projections.py @@ -13,7 +13,7 @@ magenta_cbf = '#cc79a7' black_cbf = '#000000' green_cbf = '#009e73' -mp.style.use('seaborn-poster') +mp.style.use('Solarize_Light2') #mp.style.use('seaborn-poster') # wasn't working. mp.rc('axes', labelsize = 12) # fontsize of the x and y labels mp.rc('xtick', labelsize = 12) # fontsize of the tick labels @@ -376,12 +376,15 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ print('index: ', idx, 'pdg: ', pdg) if pdg != abs(13): continue - muon_ke_lar = true_event.Muon_TrueKE + muon_ke_lar = true_event.Muon_TrueKE / 1000.0 p_tms_start = ROOT.TVector3(truth.MomentumTMSStart[4 * idx], truth.MomentumTMSStart[4 * idx + 1], truth.MomentumTMSStart[4 * idx + 2]) muon_ke_tms_start = sqrt(p_tms_start.Mag2() + MUON_MASS ** 2) - MUON_MASS + muon_ke_tms_start /= 1000.0 + x_z.text(11, 4, f'Muon KE at birth (LAr): {muon_ke_lar:.2f} GeV', fontsize = 12, fontweight = 'bold', color = orange_cbf) + x_z.text(11, 5, f'Muon KE entering TMS: {muon_ke_tms_start:.2f} GeV', fontsize = 12, fontweight = 'bold', color = orange_cbf) - x_z.text(3.6, -2.5, f'Muon KE at LAr: {muon_ke_lar}', rotation = 'vertical', fontsize = 12, fontweight = 'bold', color = orange_cbf) - x_z.text(3.6, -2.75, f'Muon KE entering TMS: {muon_ke_tms_start}', rotation= 'vertical', fontsize = 12, fontweight = 'bold', color = orange_cbf) + if muon_ke_tms_start or muon_ke_lar > 5.0: # GeV + print(f'Event: {i}, Spill {spill_number_cache[i]}, Muon KE at birth (LAr): {muon_ke_lar}, Muon KE entering TMS: {muon_ke_tms_start}, GeV.') output_filename = os.path.join(out_dir, f"{name}_{current_spill_number:03d}") mp.savefig(output_filename + ".png", bbox_inches = 'tight') From 5d74b2715617c3115276237be52fe1c7a4a97107 Mon Sep 17 00:00:00 2001 From: Michael Dolce Date: Wed, 28 Aug 2024 18:35:02 -0400 Subject: [PATCH 4/5] remove the prints, too much info. Rename the arg so it makes (some) sense to user. And had wrong index for spill number. --- scripts/Reco/draw_spill_3D_projections.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/Reco/draw_spill_3D_projections.py b/scripts/Reco/draw_spill_3D_projections.py index 45bf9ae2..3fea0d80 100644 --- a/scripts/Reco/draw_spill_3D_projections.py +++ b/scripts/Reco/draw_spill_3D_projections.py @@ -146,7 +146,7 @@ def hit_size(hit_x, hit_y, orientation, hit_z): return np.array(size_array[0]), size_array[1, 0], size_array[1, 1] ### Actual function that loops through the spills -def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, only_true_tms_muons = False): +def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, add_ke_truth_to_plot = False): if not os.path.exists(input_filename): raise ValueError(f"Cannor find input_filename {input_filename}") if readout_filename != "" and not os.path.exists(readout_filename): raise ValueError(f"Cannot find readout_filename {readout_filename}") if spill_number < -1: raise ValueError(f"Got spill_number = {spill_number}") @@ -371,9 +371,8 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ # Write the True Muon KE to each spill plot. - if only_true_tms_muons: + if add_ke_truth_to_plot: for idx, pdg in enumerate(true_event.PDG): - print('index: ', idx, 'pdg: ', pdg) if pdg != abs(13): continue muon_ke_lar = true_event.Muon_TrueKE / 1000.0 @@ -383,8 +382,8 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ x_z.text(11, 4, f'Muon KE at birth (LAr): {muon_ke_lar:.2f} GeV', fontsize = 12, fontweight = 'bold', color = orange_cbf) x_z.text(11, 5, f'Muon KE entering TMS: {muon_ke_tms_start:.2f} GeV', fontsize = 12, fontweight = 'bold', color = orange_cbf) - if muon_ke_tms_start or muon_ke_lar > 5.0: # GeV - print(f'Event: {i}, Spill {spill_number_cache[i]}, Muon KE at birth (LAr): {muon_ke_lar}, Muon KE entering TMS: {muon_ke_tms_start}, GeV.') + if muon_ke_tms_start > 5.0 or muon_ke_lar > 5.0: # GeV + print(f'Event: {i}, Spill {spill_number}, Muon KE at birth (LAr): {muon_ke_lar}, Muon KE entering TMS: {muon_ke_tms_start}, GeV.') output_filename = os.path.join(out_dir, f"{name}_{current_spill_number:03d}") mp.savefig(output_filename + ".png", bbox_inches = 'tight') @@ -441,7 +440,7 @@ def calculate_layers(Xlayers): parser.add_argument('--spillnum', "-s", type = int, help = "The spill to draw. -1 for all", default = -1) parser.add_argument('--timeslice', "-t", type = int, help = "The time slice to draw. -1 for all", default = -1) parser.add_argument('--readout_filename', "-r", type = str, help = "(optional) A file with the raw readout.", default = "") - parser.add_argument('--only_true_tms_muons', help = "Only draw events with true muons inside the TMS", action = argparse.BooleanOptionalAction) + parser.add_argument('--add_ke_truth_to_plot', help = "Add the true KE of muon to plot.", action = argparse.BooleanOptionalAction) parser.add_argument('--Xlayers', "-X", help = "Does the geometry use X (90 degree orientated) scintillator layers? Yes -> --Xlayers, No -> --no-Xlayers", action = argparse.BooleanOptionalAction) args = parser.parse_args() @@ -452,11 +451,11 @@ def calculate_layers(Xlayers): spill_number = args.spillnum time_slice = args.timeslice readout_filename = args.readout_filename - only_true_tms_muons = args.only_true_tms_muons + add_ke_truth_to_plot = args.add_ke_truth_to_plot Xlayers = args.Xlayers print(Xlayers) calculate_layers(Xlayers) print(layer_dict) - draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, only_true_tms_muons) + draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, add_ke_truth_to_plot) From 6457c640eca9e2c314afbaa53bf52238307b8cdf Mon Sep 17 00:00:00 2001 From: Michael Dolce Date: Wed, 18 Sep 2024 21:02:30 -0500 Subject: [PATCH 5/5] addressing Asa's comments: fix style to seaborn, and simplify the arg --- scripts/Reco/draw_spill_3D_projections.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/scripts/Reco/draw_spill_3D_projections.py b/scripts/Reco/draw_spill_3D_projections.py index 3fea0d80..9a8d0a4f 100644 --- a/scripts/Reco/draw_spill_3D_projections.py +++ b/scripts/Reco/draw_spill_3D_projections.py @@ -13,7 +13,7 @@ magenta_cbf = '#cc79a7' black_cbf = '#000000' green_cbf = '#009e73' -mp.style.use('Solarize_Light2') #mp.style.use('seaborn-poster') # wasn't working. +mp.style.use('seaborn-v0_8-poster') mp.rc('axes', labelsize = 12) # fontsize of the x and y labels mp.rc('xtick', labelsize = 12) # fontsize of the tick labels @@ -146,7 +146,7 @@ def hit_size(hit_x, hit_y, orientation, hit_z): return np.array(size_array[0]), size_array[1, 0], size_array[1, 1] ### Actual function that loops through the spills -def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, add_ke_truth_to_plot = False): +def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, report_true_ke = False): if not os.path.exists(input_filename): raise ValueError(f"Cannor find input_filename {input_filename}") if readout_filename != "" and not os.path.exists(readout_filename): raise ValueError(f"Cannot find readout_filename {readout_filename}") if spill_number < -1: raise ValueError(f"Got spill_number = {spill_number}") @@ -371,7 +371,7 @@ def draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_ # Write the True Muon KE to each spill plot. - if add_ke_truth_to_plot: + if report_true_ke: for idx, pdg in enumerate(true_event.PDG): if pdg != abs(13): continue @@ -440,7 +440,7 @@ def calculate_layers(Xlayers): parser.add_argument('--spillnum', "-s", type = int, help = "The spill to draw. -1 for all", default = -1) parser.add_argument('--timeslice', "-t", type = int, help = "The time slice to draw. -1 for all", default = -1) parser.add_argument('--readout_filename', "-r", type = str, help = "(optional) A file with the raw readout.", default = "") - parser.add_argument('--add_ke_truth_to_plot', help = "Add the true KE of muon to plot.", action = argparse.BooleanOptionalAction) + parser.add_argument('--report_true_ke', help = "Add the true KE of muon to plot.", action = argparse.BooleanOptionalAction) parser.add_argument('--Xlayers', "-X", help = "Does the geometry use X (90 degree orientated) scintillator layers? Yes -> --Xlayers, No -> --no-Xlayers", action = argparse.BooleanOptionalAction) args = parser.parse_args() @@ -451,11 +451,11 @@ def calculate_layers(Xlayers): spill_number = args.spillnum time_slice = args.timeslice readout_filename = args.readout_filename - add_ke_truth_to_plot = args.add_ke_truth_to_plot + report_true_ke = args.report_true_ke Xlayers = args.Xlayers print(Xlayers) calculate_layers(Xlayers) print(layer_dict) - draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, add_ke_truth_to_plot) + draw_spill(out_dir, name, input_filename, spill_number, time_slice, readout_filename, report_true_ke)