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

@@ -6,24 +6,29 @@
     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,
     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 +38,18 @@
     "create_gh_cuts_hdu",
     "create_gti_hdu",
     "create_pointing_hdu",
+    "telescope_combinations",
     "format_object",
     "get_dl2_mean",
+    "get_stereo",
     "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,151 @@ def create_event_hdu(
     return event_hdu
 
 
+
+'''
+
+def create_event_hdu(
+    event_table, config, on_time, deadc, source_name, source_ra=None, source_dec=None
+):
+    """
+    Creates a fits binary table HDU for shower events.
+
+    Parameters
+    ----------
+    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
+        Dead time correction factor
+    source_name: str
+        Name of the observed source
+    source_ra: str
+        Right ascension of the observed source, whose format should be
+        acceptable by `astropy.coordinates.sky_coordinate.SkyCoord`
+        (Used only when the source name cannot be resolved)
+    source_dec: str
+        Declination of the observed source, whose format should be
+        acceptable by `astropy.coordinates.sky_coordinate.SkyCoord`
+        (Used only when the source name cannot be resolved)
+
+    Returns
+    -------
+    event_hdu: astropy.io.fits.hdu.table.BinTableHDU
+        Event HDU
+
+    Raises
+    ------
+    ValueError
+        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")
+    time_start_iso = time_start.to_value("iso", "date_hms")
+
+    time_end = Time(event_table["timestamp"][-1], format="unix", scale="utc")
+    time_end_iso = time_end.to_value("iso", "date_hms")
+
+    # Calculate the elapsed and effective time
+    elapsed_time = time_end - time_start
+    effective_time = on_time * deadc
+
+    # Get the instruments used for the observation
+    combo_types_unique = np.unique(event_table["combo_type"])
+    tel_combos = np.array(list(TEL_COMBINATIONS.keys()))[combo_types_unique]
+
+    tel_list = [tel_combo.split("_") for tel_combo in tel_combos]
+    tel_list_unique = np.unique(sum(tel_list, []))
+
+    instruments = "_".join(tel_list_unique)
+
+    # Transfer the RA/Dec directions to the galactic coordinate
+    event_coords = SkyCoord(
+        ra=event_table["reco_ra"], dec=event_table["reco_dec"], frame="icrs"
+    )
+
+    event_coords = event_coords.galactic
+
+    try:
+        # Try to get the source coordinate from the input name
+        source_coord = SkyCoord.from_name(source_name, frame="icrs")
+
+    except Exception:
+        logger.warning(
+            f"WARNING: The source name '{source_name}' could not be resolved. "
+            f"Setting the input RA/Dec coordinate ({source_ra}, {source_dec})..."
+        )
+
+        if (source_ra is None) or (source_dec is None):
+            raise ValueError("The input RA/Dec coordinate is set to `None`.")
+
+        source_coord = SkyCoord(ra=source_ra, dec=source_dec, frame="icrs")
+
+    # Create a table
+    qtable = QTable(
+        data={
+            "EVENT_ID": event_table["event_id"],
+            "TIME": event_table["timestamp"],
+            "RA": event_table["reco_ra"],
+            "DEC": event_table["reco_dec"],
+            "ENERGY": event_table["reco_energy"],
+            "GAMMANESS": event_table["gammaness"],
+            "MULTIP": event_table["multiplicity"],
+            "GLON": event_coords.l.to("deg"),
+            "GLAT": event_coords.b.to("deg"),
+            "ALT": event_table["reco_alt"].to("deg"),
+            "AZ": event_table["reco_az"].to("deg"),
+        }
+    )
+
+    # Create a header
+    header = fits.Header(
+        cards=[
+            ("CREATED", Time.now().utc.iso),
+            ("HDUCLAS1", "EVENTS"),
+            ("OBS_ID", np.unique(event_table["obs_id"])[0]),
+            ("DATE-OBS", time_start_iso[:10]),
+            ("TIME-OBS", time_start_iso[11:]),
+            ("DATE-END", time_end_iso[:10]),
+            ("TIME-END", time_end_iso[11:]),
+            ("TSTART", time_start.value),
+            ("TSTOP", time_end.value),
+            ("MJDREFI", mjdrefi),
+            ("MJDREFF", mjdreff),
+            ("TIMEUNIT", "s"),
+            ("TIMESYS", "UTC"),
+            ("TIMEREF", "TOPOCENTER"),
+            ("ONTIME", on_time.value),
+            ("TELAPSE", elapsed_time.to_value("s")),
+            ("DEADC", deadc),
+            ("LIVETIME", effective_time.value),
+            ("OBJECT", source_name),
+            ("OBS_MODE", "WOBBLE"),
+            ("N_TELS", np.max(event_table["multiplicity"])),
+            ("TELLIST", instruments),
+            ("INSTRUME", instruments),
+            ("RA_PNT", event_table["pointing_ra"][0].value, "deg"),
+            ("DEC_PNT", event_table["pointing_dec"][0].value, "deg"),
+            ("ALT_PNT", event_table["pointing_alt"][0].to_value("deg"), "deg"),
+            ("AZ_PNT", event_table["pointing_az"][0].to_value("deg"), "deg"),
+            ("RA_OBJ", source_coord.ra.to_value("deg"), "deg"),
+            ("DEC_OBJ", source_coord.dec.to_value("deg"), "deg"),
+            ("FOVALIGN", "RADEC"),
+        ]
+    )
+
+    # Create a HDU
+    event_hdu = fits.BinTableHDU(qtable, header=header, name="EVENTS")
+
+    return event_hdu
+'''
+
 def create_gti_hdu(event_table):
     """
     Creates a fits binary table HDU for Good Time Interval (GTI).