Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add missing dl1 parameters #41

Merged
merged 13 commits into from
Jun 4, 2020
Merged

Add missing dl1 parameters #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
a674fc2
Merge branch 'patch-gain_selection_threshold' into feature-add_missin…
Nov 25, 2019
1144d37
Add all COG missing coordinate data.
Nov 25, 2019
b250927
Add missing declaration of recently added COG cooordinates data.
Nov 26, 2019
70dee29
Use correct number of islands in case of extended multi-cluster clean…
Dec 12, 2019
e290813
Added leakage parameter to the DL1 output.
Dec 12, 2019
043da1e
Merged 'master' branch and solved merge-conflicts.
Jan 27, 2020
d93b9ad
Removed stale code.
Jan 27, 2020
264d45d
Added missing COG coordinates for both cleanings.
Jan 29, 2020
548cdbe
Merge branch 'master' into feature-add_missing_DL1_parameters
Feb 7, 2020
37336b5
Names of leakage-related columns in DL1 file more similar to new data…
Feb 7, 2020
3cfb7f3
Re-ordening of some DL1 output.
Feb 7, 2020
abcc851
Forgot to modify column descriptors to account for new names.
Feb 13, 2020
db806f0
Merge branch 'master' into feature-add_missing_DL1_parameters
May 5, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
46 changes: 37 additions & 9 deletions protopipe/pipeline/event_preparer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,7 @@
from ctapipe.io.containers import TimingParametersContainer
from ctapipe.calib import CameraCalibrator
from ctapipe.calib.camera.gainselection import GainSelector

# from ctapipe.image.extractor import LocalPeakWindowSum
from ctapipe.image import hillas
from ctapipe.image import hillas, leakage, concentration
from ctapipe.image.cleaning import tailcuts_clean
from ctapipe.utils.CutFlow import CutFlow
from ctapipe.coordinates import GroundFrame
Expand Down Expand Up @@ -64,6 +62,7 @@
"n_pixel_dict",
"hillas_dict",
"hillas_dict_reco",
"leakage_dict",
"n_tels",
"tot_signal",
"max_signals",
Expand Down Expand Up @@ -698,11 +697,12 @@ def stub(event):
return PreparedEvent(
event=event,
dl1_phe_image=None, # container for the calibrated image in phe
dl1_phe_image_mask_reco =None, # container for the reco cleaning mask
dl1_phe_image_mask_reco=None, # container for the reco cleaning mask
mc_phe_image=None, # container for the simulated image in phe
n_pixel_dict=None,
hillas_dict=None,
hillas_dict_reco=None,
leakage_dict=None,
n_tels=None,
tot_signal=None,
max_signals=None,
Expand Down Expand Up @@ -903,6 +903,7 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
n_pixel_dict = {}
hillas_dict_reco = {} # for direction reconstruction
hillas_dict = {} # for discrimination
leakage_dict = {}
n_tels = {
"tot": len(event.dl0.tels_with_data),
"LST_LST_LSTCam": 0,
Expand Down Expand Up @@ -947,7 +948,7 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
else: # no other image extractor has 2 passes
dl1_phe_image_1stPass[tel_id] = pmt_signal
calibration_status[tel_id] = np.nan

mc_phe_image[tel_id] = event.mc.tel[tel_id].photo_electron_image

if self.cleaner_reco.mode == "tail": # tail uses only ctapipe
Expand All @@ -956,6 +957,19 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
image_biggest, mask_reco = self.cleaner_reco.clean_image(
pmt_signal, camera
)

# calculate the leakage (before filtering)
leakages = {} # this is needed by both cleanings
# The check on SIZE shouldn't be here, but for the moment
# I prefer to sacrifice fancincess
if np.sum(image_biggest[mask_reco]) != 0.0:
leakage_biggest = leakage(camera, image_biggest, mask_reco)
leakages["leak1_reco"] = leakage_biggest["leakage1_intensity"]
leakages["leak2_reco"] = leakage_biggest["leakage2_intensity"]
else:
leakages["leak1_reco"] = 0.0
leakages["leak2_reco"] = 0.0

# find all islands using this cleaning
num_islands, labels = number_of_islands(camera, mask_reco)

Expand All @@ -965,17 +979,17 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
camera_biggest = camera[mask_reco]
image_biggest = image_biggest[mask_reco]
if save_images is True:
dl1_phe_image_mask_reco[tel_id] = mask_reco
dl1_phe_image_mask_reco[tel_id] = mask_reco

elif num_islands > 1: # if more islands survived..
# ...find the biggest one
mask_biggest = largest_island(labels)
# and also reduce dimensions
camera_biggest = camera[mask_biggest]
image_biggest = image_biggest[mask_biggest]
if save_images is True:
dl1_phe_image_mask_reco[tel_id] = mask_biggest
dl1_phe_image_mask_reco[tel_id] = mask_biggest

else: # if no islands survived use old camera and image
camera_biggest = camera

Expand All @@ -984,6 +998,18 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
pmt_signal, camera
)

# calculate the leakage (before filtering)
# this part is not well coded, but for the moment it works
if np.sum(image_extended[mask_extended]) != 0.0:
leakage_extended = leakage(
camera, image_extended, mask_extended
)
leakages["leak1"] = leakage_extended["leakage1_intensity"]
leakages["leak2"] = leakage_extended["leakage2_intensity"]
else:
leakages["leak1"] = 0.0
leakages["leak2"] = 0.0

# find all islands with this cleaning
# we will also register how many have been found
n_cluster_dict[tel_id], labels = number_of_islands(
Expand Down Expand Up @@ -1100,6 +1126,7 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
hillas_dict_reco[tel_id] = moments_reco
n_pixel_dict[tel_id] = len(np.where(image_extended > 0)[0])
tot_signal += moments.intensity
leakage_dict[tel_id] = leakages

n_tels["reco"] = len(hillas_dict_reco)
n_tels["discri"] = len(hillas_dict)
Expand Down Expand Up @@ -1175,6 +1202,7 @@ def prepare_event(self, source, return_stub=False, save_images=False, debug=Fals
n_pixel_dict=n_pixel_dict,
hillas_dict=hillas_dict,
hillas_dict_reco=hillas_dict_reco,
leakage_dict=leakage_dict,
n_tels=n_tels,
tot_signal=tot_signal,
max_signals=max_signals,
Expand Down
32 changes: 31 additions & 1 deletion protopipe/scripts/write_dl1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,16 @@ class EventFeatures(tb.IsDescription):
psi = tb.Float32Col(dflt=1, pos=42)
psi_reco = tb.Float32Col(dflt=1, pos=43)
sum_signal_cam_reco = tb.Float32Col(dflt=1, pos=44)
cog_x = tb.Float32Col(dflt=1, pos=45)
cog_y = tb.Float32Col(dflt=1, pos=46)
phi = tb.Float32Col(dflt=1, pos=47)
cog_x_reco = tb.Float32Col(dflt=1, pos=48)
cog_y_reco = tb.Float32Col(dflt=1, pos=49)
phi_reco = tb.Float32Col(dflt=1, pos=50)
leakage_pixels_width_1_reco = tb.Float32Col(dflt=np.nan, pos=51)
leakage_pixels_width_2_reco = tb.Float32Col(dflt=np.nan, pos=52)
leakage_pixels_width_1 = tb.Float32Col(dflt=np.nan, pos=53)
leakage_pixels_width_2 = tb.Float32Col(dflt=np.nan, pos=54)

feature_outfile = tb.open_file(args.outfile, mode="w")
feature_table = {}
Expand Down Expand Up @@ -202,6 +212,7 @@ class EventFeatures(tb.IsDescription):
n_pixel_dict,
hillas_dict,
hillas_dict_reco,
leakage_dict,
n_tels,
tot_signal,
max_signals,
Expand Down Expand Up @@ -317,7 +328,11 @@ class EventFeatures(tb.IsDescription):
feature_events[cam_id]["err_est_pos"] = np.nan
feature_events[cam_id]["err_est_dir"] = np.nan
feature_events[cam_id]["mc_energy"] = event.mc.energy.to("TeV").value
feature_events[cam_id]["cog_x"] = moments.x.to("m").value
feature_events[cam_id]["cog_y"] = moments.y.to("m").value
feature_events[cam_id]["phi"] = moments.phi.to("deg").value
feature_events[cam_id]["local_distance"] = moments.r.to("m").value

feature_events[cam_id]["n_pixel"] = n_pixel_dict[tel_id]
feature_events[cam_id]["obs_id"] = event.r0.obs_id
feature_events[cam_id]["event_id"] = event.r0.event_id
Expand All @@ -326,7 +341,6 @@ class EventFeatures(tb.IsDescription):
feature_events[cam_id]["reco_energy"] = reco_energy
feature_events[cam_id]["ellipticity"] = ellipticity.value
feature_events[cam_id]["n_cluster"] = n_cluster_dict[tel_id]
feature_events[cam_id]["n_tel_reco"] = n_tels["reco"]
feature_events[cam_id]["n_tel_discri"] = n_tels["discri"]
feature_events[cam_id]["mc_core_x"] = event.mc.core_x.to("m").value
feature_events[cam_id]["mc_core_y"] = event.mc.core_y.to("m").value
Expand All @@ -341,6 +355,10 @@ class EventFeatures(tb.IsDescription):
feature_events[cam_id]["az"] = reco_result.az.to("deg").value
feature_events[cam_id]["reco_energy_tel"] = reco_energy_tel[tel_id]
# Variables from hillas_dist_reco
feature_events[cam_id]["n_tel_reco"] = n_tels["reco"]
feature_events[cam_id]["cog_x_reco"] = moments_reco.x.to("m").value
feature_events[cam_id]["cog_y_reco"] = moments_reco.y.to("m").value
feature_events[cam_id]["phi_reco"] = moments_reco.phi.to("deg").value
feature_events[cam_id]["ellipticity_reco"] = ellipticity_reco.value
feature_events[cam_id]["local_distance_reco"] = moments_reco.r.to(
"m"
Expand All @@ -353,6 +371,18 @@ class EventFeatures(tb.IsDescription):
).value
feature_events[cam_id]["psi_reco"] = moments_reco.psi.to("deg").value
feature_events[cam_id]["sum_signal_cam_reco"] = moments_reco.intensity
feature_events[cam_id]["leakage_pixels_width_1_reco"] = leakage_dict[
tel_id
]["leak1_reco"]
feature_events[cam_id]["leakage_pixels_width_2_reco"] = leakage_dict[
tel_id
]["leak2_reco"]
feature_events[cam_id]["leakage_pixels_width_1"] = leakage_dict[tel_id][
"leak1"
]
feature_events[cam_id]["leakage_pixels_width_2"] = leakage_dict[tel_id][
"leak2"
]

feature_events[cam_id].append()

Expand Down