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/io/__init__.py

@@ -13,6 +13,7 @@
     create_pointing_hdu,
 )
 from .io import (
+    telescope_combinations,
     format_object,
     get_dl2_mean,
     get_stereo_events,
@@ -22,6 +23,7 @@
     load_magic_dl1_data_files,
     load_mc_dl2_data_file,
     load_train_data_files,
+    load_train_data_files_tel,
     save_pandas_data_in_table,
 )
 
@@ -33,6 +35,7 @@
     "create_gh_cuts_hdu",
     "create_gti_hdu",
     "create_pointing_hdu",
+    "telescope_combinations",
     "format_object",
     "get_dl2_mean",
     "get_stereo_events",
@@ -42,5 +45,6 @@
     "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,7 +10,6 @@
 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
 
@@ -21,6 +20,7 @@
     "create_pointing_hdu",
 ]
 
+
 logger = logging.getLogger(__name__)
 logger.addHandler(logging.StreamHandler())
 logger.setLevel(logging.INFO)
@@ -29,6 +29,55 @@
 MJDREF = Time(0, format="unix", scale="utc")
 
 
+def telescope_combinations(config):
+    """
+    Generates all possible telescope combinations without repetition. E.g.: "LST1_M1", "LST2_LST4_M2", "LST1_LST2_LST3_M1" and so on.
+
+    Parameters
+    ----------
+    config: dict
+        yaml file with information about the telescope IDs. Typically evoked from "config_general.yaml" in the main scripts.
+
+    Returns
+    -------
+    TEL_NAMES: dict
+        Dictionary with telescope IDs and names.
+    TEL_COMBINATIONS: dict
+        Dictionary with all telescope combinations with no repetions.
+    """
+
+
+    TEL_NAMES = {}
+    for k, v in config["mc_tel_ids"].items():    #Here we swap the dictionary keys and values just for convenience.
+        if v > 0:
+            TEL_NAMES[v] =  k
+
+    TEL_COMBINATIONS = {}
+    keys = list(TEL_NAMES.keys())
+
+    def recursive_solution(current_tel, current_comb):
+
+        if current_tel == len(keys):     #The function stops once we reach the last telescope
+            return 
+
+        current_comb_name = current_comb[0] + '_' + TEL_NAMES[keys[current_tel]]  #Name of the combo (at this point it can even be a single telescope)
+        current_comb_list = current_comb[1] + [keys[current_tel]]                 #List of telescopes (including individual telescopes)
+
+        if len(current_comb_list) > 1:                          #We save them in the new dictionary excluding the single-telescope values
+            TEL_COMBINATIONS[current_comb_name[1:]] = current_comb_list;
+
+        current_comb = [current_comb_name, current_comb_list]   #We save the current results in this varible to recal the function recursively ("for" loop below)
+
+        for i in range(1, len(keys)-current_tel):               
+            recursive_solution(current_tel+i, current_comb)
+
+
+    for key in range(len(keys)):
+        recursive_solution(key, ['',[]])
+
+
+    return TEL_NAMES, TEL_COMBINATIONS
+
 @u.quantity_input
 def create_gh_cuts_hdu(
     gh_cuts, reco_energy_bins: u.TeV, fov_offset_bins: u.deg, **header_cards
@@ -90,7 +139,7 @@ def create_gh_cuts_hdu(
 
 
 def create_event_hdu(
-    event_table, on_time, deadc, source_name, source_ra=None, source_dec=None
+    event_table, config, on_time, deadc, source_name, source_ra=None, source_dec=None
 ):
     """
     Creates a fits binary table HDU for shower events.
@@ -99,6 +148,8 @@ def create_event_hdu(
     ----------
     event_table: astropy.table.table.QTable
         Table of the DL2 events surviving gammaness cuts
+    config: dict
+        yaml file with information about the telescope IDs. Typically called evoked from "config_DL3.yaml" in the main scripts.
     on_time: astropy.table.table.QTable
         ON time of the input data
     deadc: float
@@ -125,7 +176,8 @@ 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 = telescope_combinations(config)
+
     mjdreff, mjdrefi = np.modf(MJDREF.mjd)
 
     time_start = Time(event_table["timestamp"][0], format="unix", scale="utc")