diff --git a/protopipe/pipeline/event_preparer.py b/protopipe/pipeline/event_preparer.py index 87d62606..9bce6aa3 100644 --- a/protopipe/pipeline/event_preparer.py +++ b/protopipe/pipeline/event_preparer.py @@ -11,7 +11,10 @@ from ctapipe.containers import ReconstructedShowerContainer from ctapipe.calib import CameraCalibrator from ctapipe.image.extractor import TwoPassWindowSum -from ctapipe.image import leakage_parameters, number_of_islands, largest_island +from ctapipe.image import (leakage_parameters, + number_of_islands, + largest_island, + concentration_parameters) from ctapipe.utils.CutFlow import CutFlow from ctapipe.coordinates import GroundFrame, TelescopeFrame, CameraFrame @@ -56,6 +59,7 @@ "hillas_dict", "hillas_dict_reco", "leakage_dict", + "concentration_dict", "n_tels", "max_signals", "n_cluster_dict", @@ -78,6 +82,7 @@ def stub( hillas_dict_reco, n_tels, leakage_dict, + concentration_dict ): """Default container for images that did not survive cleaning.""" return PreparedEvent( @@ -91,6 +96,7 @@ def stub( hillas_dict=hillas_dict, hillas_dict_reco=hillas_dict_reco, leakage_dict=leakage_dict, + concentration_dict=concentration_dict, n_tels=n_tels, max_signals=dict.fromkeys(hillas_dict_reco.keys(), np.nan), # no charge n_cluster_dict=dict.fromkeys(hillas_dict_reco.keys(), 0), # no clusters @@ -360,6 +366,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict_reco = {} # for direction reconstruction hillas_dict = {} # for discrimination leakage_dict = {} + concentration_dict = {} n_tels = { "Triggered": len(event.r1.tel.keys()), "LST_LST_LSTCam": 0, @@ -614,6 +621,13 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False print("Image parameters from all-clusters cleaning:") print(moments) + # Add concentration parameters + concentrations = {} + concentrations_extended = concentration_parameters(camera_extended_tel, image_extended, moments) + concentrations["concentration_cog"] = concentrations_extended["cog"] + concentrations["concentration_core"] = concentrations_extended["core"] + concentrations["concentration_pixel"] = concentrations_extended["pixel"] + # =================================================== # PARAMETRIZED IMAGE SELECTION # =================================================== @@ -675,6 +689,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict_reco[tel_id] = moments_reco n_pixel_dict[tel_id] = len(np.where(image_extended > 0)[0]) leakage_dict[tel_id] = leakages + concentration_dict[tel_id] = concentrations except ( FloatingPointError, @@ -696,6 +711,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False ] = HillasParametersTelescopeFrameContainer() n_pixel_dict[tel_id] = len(np.where(image_extended > 0)[0]) leakage_dict[tel_id] = leakages + concentration_dict[tel_id] = concentrations # END OF THE CYCLE OVER THE TELESCOPES @@ -743,6 +759,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict_reco, n_tels, leakage_dict, + concentration_dict ) continue else: @@ -832,6 +849,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict_reco, n_tels, leakage_dict, + concentration_dict ) else: continue @@ -862,6 +880,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict_reco, n_tels, leakage_dict, + concentration_dict ) else: continue @@ -884,6 +903,7 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False hillas_dict=hillas_dict, hillas_dict_reco=hillas_dict_reco, leakage_dict=leakage_dict, + concentration_dict=concentration_dict, n_tels=n_tels, max_signals=max_signals, n_cluster_dict=n_cluster_dict, diff --git a/protopipe/scripts/data_training.py b/protopipe/scripts/data_training.py index 66a63c43..68b47ec5 100755 --- a/protopipe/scripts/data_training.py +++ b/protopipe/scripts/data_training.py @@ -169,36 +169,39 @@ def main(): leakage_intensity_width_2_reco=tb.Float32Col(dflt=np.nan, pos=29), leakage_intensity_width_1=tb.Float32Col(dflt=np.nan, pos=30), leakage_intensity_width_2=tb.Float32Col(dflt=np.nan, pos=31), + concentration_cog=tb.Float32Col(dflt=np.nan, pos=32), + concentration_core=tb.Float32Col(dflt=np.nan, pos=33), + concentration_pixel=tb.Float32Col(dflt=np.nan, pos=34), # The following are missing from current ctapipe DL1 output # Not sure if it's worth to add them - hillas_ellipticity_reco=tb.FloatCol(dflt=1, pos=32), - hillas_ellipticity=tb.FloatCol(dflt=1, pos=33), - max_signal_cam=tb.Float32Col(dflt=1, pos=34), - pixels=tb.Int16Col(dflt=1, pos=35), - clusters=tb.Int16Col(dflt=-1, pos=36), + hillas_ellipticity_reco=tb.FloatCol(dflt=1, pos=35), + hillas_ellipticity=tb.FloatCol(dflt=1, pos=36), + max_signal_cam=tb.Float32Col(dflt=1, pos=37), + pixels=tb.Int16Col(dflt=1, pos=38), + clusters=tb.Int16Col(dflt=-1, pos=39), # ====================================================================== # DL2 - DIRECTION RECONSTRUCTION - impact_dist=tb.Float32Col(dflt=1, pos=37), - h_max=tb.Float32Col(dflt=1, pos=38), - alt=tb.Float32Col(dflt=np.nan, pos=39), - az=tb.Float32Col(dflt=np.nan, pos=40), - err_est_pos=tb.Float32Col(dflt=1, pos=41), - err_est_dir=tb.Float32Col(dflt=1, pos=42), - xi=tb.Float32Col(dflt=np.nan, pos=43), - offset=tb.Float32Col(dflt=np.nan, pos=44), - mc_core_x=tb.FloatCol(dflt=1, pos=45), - mc_core_y=tb.FloatCol(dflt=1, pos=46), - reco_core_x=tb.FloatCol(dflt=1, pos=47), - reco_core_y=tb.FloatCol(dflt=1, pos=48), - mc_h_first_int=tb.FloatCol(dflt=1, pos=49), - mc_x_max=tb.Float32Col(dflt=np.nan, pos=50), - is_valid=tb.BoolCol(dflt=False, pos=51), - good_image=tb.Int16Col(dflt=1, pos=52), + impact_dist=tb.Float32Col(dflt=1, pos=40), + h_max=tb.Float32Col(dflt=1, pos=41), + alt=tb.Float32Col(dflt=np.nan, pos=42), + az=tb.Float32Col(dflt=np.nan, pos=43), + err_est_pos=tb.Float32Col(dflt=1, pos=44), + err_est_dir=tb.Float32Col(dflt=1, pos=45), + xi=tb.Float32Col(dflt=np.nan, pos=46), + offset=tb.Float32Col(dflt=np.nan, pos=47), + mc_core_x=tb.FloatCol(dflt=1, pos=48), + mc_core_y=tb.FloatCol(dflt=1, pos=49), + reco_core_x=tb.FloatCol(dflt=1, pos=50), + reco_core_y=tb.FloatCol(dflt=1, pos=51), + mc_h_first_int=tb.FloatCol(dflt=1, pos=52), + mc_x_max=tb.Float32Col(dflt=np.nan, pos=53), + is_valid=tb.BoolCol(dflt=False, pos=54), + good_image=tb.Int16Col(dflt=1, pos=55), # ====================================================================== # DL2 - ENERGY ESTIMATION - true_energy=tb.FloatCol(dflt=1, pos=53), - reco_energy=tb.FloatCol(dflt=np.nan, pos=54), - reco_energy_tel=tb.Float32Col(dflt=np.nan, pos=55), + true_energy=tb.FloatCol(dflt=1, pos=56), + reco_energy=tb.FloatCol(dflt=np.nan, pos=57), + reco_energy_tel=tb.Float32Col(dflt=np.nan, pos=58), # ====================================================================== # DL1 IMAGES # this is optional data saved by the user @@ -235,6 +238,7 @@ def main(): hillas_dict, hillas_dict_reco, leakage_dict, + concentration_dict, n_tels, max_signals, n_cluster_dict, @@ -330,6 +334,9 @@ def main(): "leakage_intensity_width_2_reco": [leakage_dict[tel_id]['leak2_reco']], "leakage_intensity_width_1": [leakage_dict[tel_id]['leak1']], "leakage_intensity_width_2": [leakage_dict[tel_id]['leak2']], + "concentration_cog": [concentration_dict[tel_id]['concentration_cog']], + "concentration_core": [concentration_dict[tel_id]['concentration_core']], + "concentration_pixel": [concentration_dict[tel_id]['concentration_pixel']], "az": [reco_result.az.to("deg").value], "alt": [reco_result.alt.to("deg").value], "h_max": [h_max.value], @@ -483,6 +490,15 @@ def main(): outData[cam_id]["leakage_intensity_width_2"] = leakage_dict[tel_id][ "leak2" ] + outData[cam_id]["concentration_cog"] = concentration_dict[tel_id][ + "concentration_cog" + ] + outData[cam_id]["concentration_core"] = concentration_dict[tel_id][ + "concentration_core" + ] + outData[cam_id]["concentration_pixel"] = concentration_dict[tel_id][ + "concentration_pixel" + ] # ======================= # IMAGES INFORMATION diff --git a/protopipe/scripts/write_dl2.py b/protopipe/scripts/write_dl2.py index f30beacb..86f5db29 100755 --- a/protopipe/scripts/write_dl2.py +++ b/protopipe/scripts/write_dl2.py @@ -283,6 +283,7 @@ class RecoEvent(tb.IsDescription): hillas_dict, hillas_dict_reco, leakage_dict, + concentration_dict, n_tels, max_signals, n_cluster_dict, @@ -386,6 +387,9 @@ class RecoEvent(tb.IsDescription): "leakage_intensity_width_2_reco": [leakage_dict[tel_id]['leak2_reco']], "leakage_intensity_width_1": [leakage_dict[tel_id]['leak1']], "leakage_intensity_width_2": [leakage_dict[tel_id]['leak2']], + "concentration_cog": [concentration_dict[tel_id]['concentration_cog']], + "concentration_core": [concentration_dict[tel_id]['concentration_core']], + "concentration_pixel": [concentration_dict[tel_id]['concentration_pixel']], "az": [reco_result.az.to("deg").value], "alt": [reco_result.alt.to("deg").value], "h_max": [h_max.value], @@ -473,6 +477,9 @@ class RecoEvent(tb.IsDescription): "leakage_intensity_width_2_reco": [leakage_dict[tel_id]['leak2_reco']], "leakage_intensity_width_1": [leakage_dict[tel_id]['leak1']], "leakage_intensity_width_2": [leakage_dict[tel_id]['leak2']], + "concentration_cog": [concentration_dict[tel_id]['concentration_cog']], + "concentration_core": [concentration_dict[tel_id]['concentration_core']], + "concentration_pixel": [concentration_dict[tel_id]['concentration_pixel']], "az": [reco_result.az.to("deg").value], "alt": [reco_result.alt.to("deg").value], "h_max": [h_max.value],