Skip to content

Commit

Permalink
Merge pull request #41 from HealthyPear/feature-add_missing_DL1_param…
Browse files Browse the repository at this point in the history 
…eters

Add missing dl1 parameters
  • Loading branch information
@HealthyPear
HealthyPear committed Jun 4, 2020
2 parents c46571c + db806f0 commit a328138
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 10 deletions.
46 changes: 37 additions & 9 deletions protopipe/pipeline/event_preparer.py
View file
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
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

0 comments on commit a328138

Please sign in to comment.