MCP upgrades: semi-automatic scripts and update of the current scripts (1 LST -> 4 LSTs); Torino team

environment.yml

@@ -1,5 +1,5 @@
 # A conda environment with all useful package for ctapipe developers
-name: magic-lst1
+name: magic-lst
 channels:
   - default
   - conda-forge

magicctapipe/conftest.py

@@ -7,6 +7,7 @@
 from math import trunc
 from ctapipe.utils.download import download_file_cached
 from magicctapipe.utils import resource_file
+import yaml
 
 maxjoint = 13000
 maxmonly = 500
@@ -50,6 +51,7 @@
 @pytest.fixture(scope="session")
 def temp_DL1_gamma(tmp_path_factory):
     return tmp_path_factory.mktemp("DL1_gammas")
+
 
 
 @pytest.fixture(scope="session")
@@ -384,6 +386,13 @@ def config_monly():
     config_path = resource_file("config_monly.yaml")
     return config_path
 
+@pytest.fixture(scope="session")
+def config_gen():
+    config_path = resource_file("config_general.yaml")
+    with open(config_path, "rb") as f:
+        config = yaml.safe_load(f)
+    return config
+
 
 """
 Data processing

magicctapipe/image/__init__.py

@@ -9,10 +9,17 @@
     get_leakage,
 )
 
+from .calib import (
+    Calibrate_LST, 
+    Calibrate_MAGIC
+)
+
 __all__ = [
     "MAGICClean",
     "PixelTreatment",
     "get_num_islands_MAGIC",
     "clean_image_params",
     "get_leakage",
+    "Calibrate_LST", 
+    "Calibrate_MAGIC"
 ]

magicctapipe/image/calib.py

@@ -0,0 +1,102 @@
+
+
+import numpy as np
+
+
+from ctapipe.image import (
+    apply_time_delta_cleaning,   
+    number_of_islands,
+    tailcuts_clean,   
+)
+
+from lstchain.image.cleaning import apply_dynamic_cleaning
+from lstchain.image.modifier import (
+    add_noise_in_pixels,
+    random_psf_smearer,    
+)
+
+__all__ = [
+    "Calibrate_LST", "Calibrate_MAGIC"
+]
+
+def Calibrate_LST(event, tel_id, rng, config_lst, camera_geoms, calibrator_lst, increase_nsb, use_time_delta_cleaning, use_dynamic_cleaning ):
+
+    """
+    This function computes and returns signal_pixels, image, and peak_time for LST
+    """
+
+    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_psf = config_lst["increase_psf"]["use"]
+    use_only_main_island = config_lst["use_only_main_island"]
+
+    if increase_nsb:
+        # Add extra noise in pixels
+        image = add_noise_in_pixels(rng, image, **config_lst["increase_nsb"])
+
+    if increase_psf:
+        # 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 signal_pixels, image, and peak_time for MAGIC
+    """
+
+    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

magicctapipe/io/__init__.py

@@ -6,24 +6,28 @@
     RealEventInfoContainer,
     SimEventInfoContainer,
 )
-from .gadf import (
-    create_event_hdu,
-    create_gh_cuts_hdu,
-    create_gti_hdu,
-    create_pointing_hdu,
-)
 from .io import (
+    telescope_combinations,
     format_object,
     get_dl2_mean,
     get_stereo_events,
+    get_stereo_events_old,
     load_dl2_data_file,
     load_irf_files,
     load_lst_dl1_data_file,
     load_magic_dl1_data_files,
     load_mc_dl2_data_file,
     load_train_data_files,
+    load_train_data_files_tel,
     save_pandas_data_in_table,
 )
+from .gadf import (
+    create_event_hdu,
+    create_gh_cuts_hdu,
+    create_gti_hdu,
+    create_pointing_hdu,
+)
+
 
 __all__ = [
     "BaseEventInfoContainer",
@@ -33,14 +37,17 @@
     "create_gh_cuts_hdu",
     "create_gti_hdu",
     "create_pointing_hdu",
+    "telescope_combinations",
     "format_object",
     "get_dl2_mean",
     "get_stereo_events",
+    "get_stereo_events_old",
     "load_dl2_data_file",
     "load_irf_files",
     "load_lst_dl1_data_file",
     "load_magic_dl1_data_files",
     "load_mc_dl2_data_file",
     "load_train_data_files",
+    "load_train_data_files_tel",
     "save_pandas_data_in_table",
 ]

magicctapipe/io/gadf.py

@@ -10,17 +10,17 @@
 from astropy.table import QTable
 from astropy.time import Time
 from magicctapipe import __version__
-from magicctapipe.io.io import TEL_COMBINATIONS
 from magicctapipe.utils.functions import HEIGHT_ORM, LAT_ORM, LON_ORM
 from pyirf.binning import split_bin_lo_hi
-
+#from .io import telescope_combinations
 __all__ = [
     "create_gh_cuts_hdu",
     "create_event_hdu",
     "create_gti_hdu",
     "create_pointing_hdu",
 ]
 
+
 logger = logging.getLogger(__name__)
 logger.addHandler(logging.StreamHandler())
 logger.setLevel(logging.INFO)
@@ -29,6 +29,7 @@
 MJDREF = Time(0, format="unix", scale="utc")
 
 
+
 @u.quantity_input
 def create_gh_cuts_hdu(
     gh_cuts, reco_energy_bins: u.TeV, fov_offset_bins: u.deg, **header_cards
@@ -88,7 +89,6 @@ def create_gh_cuts_hdu(
 
     return gh_cuts_hdu
 
-
 def create_event_hdu(
     event_table, on_time, deadc, source_name, source_ra=None, source_dec=None
 ):
@@ -125,6 +125,12 @@ def create_event_hdu(
         If the source name cannot be resolved and also either or both of
         source RA/Dec coordinate is set to None
     """
+    TEL_COMBINATIONS = {
+    "M1_M2": [2, 3],  # combo_type = 0
+    "LST1_M1": [1, 2],  # combo_type = 1
+    "LST1_M2": [1, 3],  # combo_type = 2
+    "LST1_M1_M2": [1, 2, 3],  # combo_type = 3
+    }  #####TO BE REMOVED WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     mjdreff, mjdrefi = np.modf(MJDREF.mjd)
 
@@ -228,6 +234,10 @@ def create_event_hdu(
     return event_hdu
 
 
+
+
+
+
 def create_gti_hdu(event_table):
     """
     Creates a fits binary table HDU for Good Time Interval (GTI).