Skip to content

Commit

Permalink
Calibrate (MAGIC & LST)
Browse files Browse the repository at this point in the history
  • Loading branch information
@Elisa-Visentin
Elisa-Visentin committed Oct 4, 2023
1 parent cd9b396 commit e0756c5
Show file tree
Hide file tree
Showing 3 changed files with 14 additions and 146 deletions.
9 changes: 3 additions & 6 deletions magicctapipe/image/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,17 +10,14 @@
)

from .calib import (
Calibrate_LST,
Calibrate_MAGIC
calibrate,
)

__all__ = [
"MAGICClean",
"PixelTreatment",
"get_num_islands_MAGIC",
"calibrate",
"clean_image_params",
"get_leakage",
"Calibrate_LST",
"Calibrate_MAGIC",
"Calibrate"
"get_leakage",
]
143 changes: 7 additions & 136 deletions magicctapipe/image/calib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,143 +13,11 @@


__all__ = [
"Calibrate_LST", "Calibrate_MAGIC", "Calibrate"
"Calibrate"
]

def Calibrate_LST(event, tel_id, obs_id, config_lst, camera_geoms, calibrator_lst):

"""
This function computes and returns some information for a single event of a telescope of LST type
Parameters
----------
event: event
From an EventSource
tel_id: int
Telescope ID
obs_id: int
Observation ID
config_lst: dictionary
Parameters for image extraction and calibration
camera_geoms: telescope.camera.geometry
Camera geometry
calibrator_lst: CameraCalibrator (ctapipe.calib)
ctapipe object needed to calibrate the camera
Returns
-------
signal_pixels: Mask of the pixels selected by the cleaning
image: Array of number of p.e. in the camera pixels
peak_time: Array of the signal peak time in the camera pixels
"""

calibrator_lst._calibrate_dl0(event, tel_id)
calibrator_lst._calibrate_dl1(event, tel_id)

image = event.dl1.tel[tel_id].image.astype(np.float64)
peak_time = event.dl1.tel[tel_id].peak_time.astype(np.float64)
increase_nsb = config_lst["increase_nsb"].pop("use")
increase_psf = config_lst["increase_psf"]["use"]
use_time_delta_cleaning = config_lst["time_delta_cleaning"].pop("use")
use_dynamic_cleaning = config_lst["dynamic_cleaning"].pop("use")
use_only_main_island = config_lst["use_only_main_island"]

if increase_nsb:
rng = np.random.default_rng(obs_id)
# Add extra noise in pixels
image = add_noise_in_pixels(rng, image, **config_lst["increase_nsb"])

if increase_psf:
set_numba_seed(obs_id)
# Smear the image
image = random_psf_smearer(
image=image,
fraction=config_lst["increase_psf"]["fraction"],
indices=camera_geoms[tel_id].neighbor_matrix_sparse.indices,
indptr=camera_geoms[tel_id].neighbor_matrix_sparse.indptr,
)

# Apply the image cleaning
signal_pixels = tailcuts_clean(
camera_geoms[tel_id], image, **config_lst["tailcuts_clean"]
)

if use_time_delta_cleaning:
signal_pixels = apply_time_delta_cleaning(
geom=camera_geoms[tel_id],
mask=signal_pixels,
arrival_times=peak_time,
**config_lst["time_delta_cleaning"],
)

if use_dynamic_cleaning:
signal_pixels = apply_dynamic_cleaning(
image, signal_pixels, **config_lst["dynamic_cleaning"]
)

if use_only_main_island:
_, island_labels = number_of_islands(camera_geoms[tel_id], signal_pixels)
n_pixels_on_island = np.bincount(island_labels.astype(np.int64))

# The first index means the pixels not surviving
# the cleaning, so should not be considered
n_pixels_on_island[0] = 0
max_island_label = np.argmax(n_pixels_on_island)
signal_pixels[island_labels != max_island_label] = False

return signal_pixels, image, peak_time


def Calibrate_MAGIC(event, tel_id, config_magic, magic_clean, calibrator_magic):

"""
This function computes and returns some information for a single event of a telescope of MAGIC type
Parameters
----------
event: event
From an EventSource
tel_id: int
telescope ID
config_magic: dictionary
Parameters for image extraction and calibration
magic_clean: dictionary (1 entry per MAGIC telescope)
Each entry is a MAGICClean object using the telescope camera geometry
calibrator_magic: CameraCalibrator (ctapipe.calib)
ctapipe object needed to calibrate the camera
Returns
-------
signal_pixels: Mask of the pixels selected by the cleaning
image: Array of number of p.e. in the camera pixels
peak_time: Array of the signal peak time in the camera pixels
"""

calibrator_magic._calibrate_dl0(event, tel_id)
calibrator_magic._calibrate_dl1(event, tel_id)

image = event.dl1.tel[tel_id].image.astype(np.float64)
peak_time = event.dl1.tel[tel_id].peak_time.astype(np.float64)
use_charge_correction = config_magic["charge_correction"]["use"]

if use_charge_correction:
# Scale the charges by the correction factor
image *= config_magic["charge_correction"]["factor"]

# Apply the image cleaning
signal_pixels, image, peak_time = magic_clean[tel_id].clean_image(
event_image=image, event_pulse_time=peak_time
)
return signal_pixels, image, peak_time





def Calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_geoms=None, magic_clean=None):
def calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_geoms=None, magic_clean=None):

"""
This function computes and returns some information for a single event of a telescope
Expand Down Expand Up @@ -187,7 +55,7 @@ def Calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_g

image = event.dl1.tel[tel_id].image.astype(np.float64)
peak_time = event.dl1.tel[tel_id].peak_time.astype(np.float64)
if LST_bool==False:
if (LST_bool==False) and (magic_clean!=None):
use_charge_correction = config["charge_correction"]["use"]

if use_charge_correction:
Expand All @@ -197,7 +65,7 @@ def Calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_g
signal_pixels, image, peak_time = magic_clean[tel_id].clean_image(
event_image=image, event_pulse_time=peak_time
)
else:
elif (LST_bool==True) and (obs_id!=None) and (camera_geoms!=None):
increase_nsb = config["increase_nsb"].pop("use")
increase_psf = config["increase_psf"]["use"]
use_time_delta_cleaning = config["time_delta_cleaning"].pop("use")
Expand Down Expand Up @@ -247,4 +115,7 @@ def Calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_g
max_island_label = np.argmax(n_pixels_on_island)
signal_pixels[island_labels != max_island_label] = False

else:
print("Check the provided parameters and the telescope type; calibration was not possible")
return
return signal_pixels, image, peak_time
8 changes: 4 additions & 4 deletions magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
from ctapipe.io import EventSource, HDF5TableWriter

from magicctapipe.image import MAGICClean
from magicctapipe.image.calib import Calibrate_LST, Calibrate_MAGIC
from magicctapipe.image.calib import calibrate
from magicctapipe.io import SimEventInfoContainer, format_object
from magicctapipe.utils import calculate_disp, calculate_impact
from traitlets.config import Config
Expand Down Expand Up @@ -230,12 +230,12 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):

for tel_id in tels_with_trigger:

if tel_id in LSTs_IDs: ##If the ID is in the LST list, we call Calibrate_LST()
if tel_id in LSTs_IDs: ##If the ID is in the LST list, we call calibrate on the LST()
# Calibrate the LST-1 event
signal_pixels, image, peak_time = Calibrate_LST(event, tel_id, obs_id, config_lst, camera_geoms, calibrator_lst)
signal_pixels, image, peak_time = calibrate(event=event, tel_id=tel_id, obs_id=obs_id, config=config_lst, camera_geoms=camera_geoms, calibrator=calibrator_lst, LST_bool=True)
elif tel_id in MAGICs_IDs:
# Calibrate the MAGIC event
signal_pixels, image, peak_time = Calibrate_MAGIC(event, tel_id, config_magic, magic_clean, calibrator_magic)
signal_pixels, image, peak_time = calibrate(event=event, tel_id=tel_id, config=config_magic, magic_clean=magic_clean, calibrator=calibrator_magic, LST_bool=False)
else:
logger.info(
f"--> Telescope ID {tel_id} not in LST list or MAGIC list. Please check if the IDs are OK in the configuration file"
Expand Down

0 comments on commit e0756c5

Please sign in to comment.