Skip to content

Commit

Permalink
Compute concentration for extended image, add it to TRAINING+DL2 outputs
Browse files Browse the repository at this point in the history
  • Loading branch information
@HealthyPear
HealthyPear committed May 10, 2021
1 parent 9149e12 commit f88d9ad
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 25 deletions.
22 changes: 21 additions & 1 deletion protopipe/pipeline/event_preparer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -56,6 +59,7 @@
"hillas_dict",
"hillas_dict_reco",
"leakage_dict",
"concentration_dict",
"n_tels",
"max_signals",
"n_cluster_dict",
Expand All @@ -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(
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
# ===================================================
Expand Down Expand Up @@ -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,
Expand All @@ -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

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand Down
64 changes: 40 additions & 24 deletions protopipe/scripts/data_training.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -235,6 +238,7 @@ def main():
hillas_dict,
hillas_dict_reco,
leakage_dict,
concentration_dict,
n_tels,
max_signals,
n_cluster_dict,
Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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
Expand Down
7 changes: 7 additions & 0 deletions protopipe/scripts/write_dl2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,7 @@ class RecoEvent(tb.IsDescription):
hillas_dict,
hillas_dict_reco,
leakage_dict,
concentration_dict,
n_tels,
max_signals,
n_cluster_dict,
Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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],
Expand Down

0 comments on commit f88d9ad

Please sign in to comment.