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
@@ -42,14 +43,14 @@
     "20201216_M2_05093711.014_Y_CrabNebula-W0.40+035.root",
 ]
 DL1_LST_data = ["dl1_LST-1.Run03265.0094.h5"]
+
 """
 Temporary paths
 """
 
-
 @pytest.fixture(scope="session")
 def temp_DL1_gamma(tmp_path_factory):
-    return tmp_path_factory.mktemp("DL1_gammas")
+    return tmp_path_factory.mktemp("DL1_gammas")       
 
 
 @pytest.fixture(scope="session")
@@ -231,12 +232,10 @@ def temp_DL2_real_monly(tmp_path_factory):
 def temp_DL3_monly(tmp_path_factory):
     return tmp_path_factory.mktemp("DL3_monly")
 
-
 """
 Custom data
 """
 
-
 @pytest.fixture(scope="session")
 def dl2_test(temp_DL2_test):
     """
@@ -270,12 +269,10 @@ def pd_test():
     df = pd.DataFrame(np.array([[1, 2], [3, 4], [5, 6]]), columns=["a", "b"])
     return df
 
-
 """
 Remote paths (to download test files)
 """
 
-
 @pytest.fixture(scope="session")
 def base_url():
     return "http://www.magic.iac.es/mcp-testdata"
@@ -286,12 +283,10 @@ def env_prefix():
     # ENVIRONMENT VARIABLES TO BE CREATED
     return "MAGIC_CTA_DATA_"
 
-
 """
 Downloads: files
 """
 
-
 @pytest.fixture(scope="session")
 def dl0_gamma(base_url, env_prefix):
     gamma_dl0 = []
@@ -375,21 +370,26 @@ def dl1_lst(base_url, env_prefix):
 
 @pytest.fixture(scope="session")
 def config():
-    config_path = resource_file("config.yaml")
+    config_path = resource_file("test_config.yaml")
     return config_path
 
 
 @pytest.fixture(scope="session")
 def config_monly():
-    config_path = resource_file("config_monly.yaml")
-    return config_path
+    config_path = resource_file("test_config_monly.yaml")
+    return config_path 
 
+@pytest.fixture(scope="session")
+def config_gen():
+    config_path = resource_file("test_config_general.yaml")
+    with open(config_path, "rb") as f:
+        config = yaml.safe_load(f)
+    return config
 
 """
 Data processing
 """
 
-
 @pytest.fixture(scope="session")
 def gamma_l1(temp_DL1_gamma, dl0_gamma, config):
     """

magicctapipe/image/__init__.py

@@ -4,15 +4,15 @@
     get_num_islands_MAGIC,
     clean_image_params,
 )
+from .leakage import get_leakage
+from .calib import calibrate
 
-from .leakage import (
-    get_leakage,
-)
 
 __all__ = [
     "MAGICClean",
     "PixelTreatment",
     "get_num_islands_MAGIC",
+    "calibrate",
     "clean_image_params",
     "get_leakage",
 ]

magicctapipe/image/calib.py

@@ -0,0 +1,123 @@
+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,  
+    set_numba_seed  
+)
+
+
+__all__ = [
+    "calibrate"
+]
+
+
+def calibrate(event, tel_id, config, calibrator, LST_bool, obs_id=None, camera_geoms=None, magic_clean=None):
+
+    """
+    This function calibrates the camera image for a single event of a telescope 
+
+    Parameters
+    ----------
+    event: event 
+        From an EventSource
+    tel_id: int
+        Telescope ID     
+    config: dictionary
+        Parameters for image extraction and calibration    
+    calibrator: CameraCalibrator (ctapipe.calib)
+        ctapipe object needed to calibrate the camera
+    LST_bool: bool
+        Whether the telescope is a LST
+    obs_id: int
+        Observation ID. Unsed in case of LSt telescope
+    camera_geoms: telescope.camera.geometry
+        Camera geometry. Used in case of LST telescope
+    magic_clean: dictionary (1 entry per MAGIC telescope)
+        Each entry is a MAGICClean object using the telescope camera geometry. Used in case of MAGIC telescope
+
+
+    Returns
+    -------
+    signal_pixels: boolean mask
+        Mask of the pixels selected by the cleaning
+    image: numpy array
+        Array of number of p.e. in the camera pixels
+    peak_time: numpy array
+        Array of the signal peak time in the camera pixels
+
+    """
+
+    calibrator._calibrate_dl0(event, tel_id)
+    calibrator._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)
+    if (LST_bool==False) and (magic_clean!=None): 
+        use_charge_correction = config["charge_correction"]["use"]
+
+        if use_charge_correction:
+            # Scale the charges by the correction factor
+            image *= config["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
+        )
+    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")
+        use_dynamic_cleaning = config["dynamic_cleaning"].pop("use")
+        use_only_main_island = config["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["increase_nsb"])
+
+        if increase_psf:
+            set_numba_seed(obs_id)
+            # Smear the image
+            image = random_psf_smearer(
+                image=image,
+                fraction=config["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["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["time_delta_cleaning"],
+            )
+
+        if use_dynamic_cleaning:
+            signal_pixels = apply_dynamic_cleaning(
+                image, signal_pixels, **config["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
+
+    else: 
+        raise ValueError("Check the provided parameters and the telescope type; calibration was not possible") 
+    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 (
     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,
+    telescope_combinations,
 )
+from .gadf import (
+    create_event_hdu,
+    create_gh_cuts_hdu,
+    create_gti_hdu,
+    create_pointing_hdu,
+)
+
 
 __all__ = [
     "BaseEventInfoContainer",
@@ -36,11 +40,14 @@
     "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",
+    "telescope_combinations",
 ]

magicctapipe/io/gadf.py

@@ -2,17 +2,16 @@
 # coding: utf-8
 
 import logging
-
 import numpy as np
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 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",
@@ -125,6 +124,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
+    }  #TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
 
     mjdreff, mjdrefi = np.modf(MJDREF.mjd)