Repository: cta-observatory/magic-cta-pipe (language: python)

Issue:
fix to create HDU and OBS index files with lstchain0.10
after removing the old copied lstchain code for creating those index files and using the current (lstchain 0.10) code instead the script stoped working because of a small change of API in the function. Note that since the function is taken from lstchain, in the header file of the output index files 'lstchain' will be written as the creation program

@aleberti or  @Elisa-Visentin, please have a quick look and if you are fine with those changes we merge them into the "lstchain" branch  


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 04da3a49..3f6b04d3 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -23,7 +23,7 @@ jobs:
           fetch-depth: 0
       - uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - uses: pre-commit/action@v3.0.0
         with:
           extra_args: --files $(git diff origin/master --name-only)
@@ -37,7 +37,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v4
         with:
-          python-version: "3.8"
+          python-version: "3.9"
 
       - name: pyflakes
         run: |
@@ -54,7 +54,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v4
         with:
-          python-version: "3.8"
+          python-version: "3.9"
 
       - name: Install doc dependencies
         run: |
@@ -72,8 +72,8 @@ jobs:
 
     strategy:
       matrix:
-        python-version: ["3.8"]
-        ctapipe-version: ["v0.12.0"]
+        python-version: ["3.9", "3.10", "3.11"]
+        ctapipe-version: ["v0.19.2"]
 
     defaults:
       run:
diff --git a/docs/conf.py b/docs/conf.py
index 20851a8e..ccddb394 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -84,6 +84,13 @@
     ("py:meth", "BaseDescriptor.instance_init"),
     ("py:obj", "cls"),
     ("py:obj", "name"),
+    ("py:obj", "Field"),
+    ("py:obj", "deg"),
+    ("py:obj", "m"),
+    ("py:obj", "TeV"),
+    ("py:obj", "rad"),
+    ("py:obj", "s"),
+    ("py:obj", "ns"),
     ("py:class", "astropy.coordinates.baseframe.BaseCoordinateFrame"),
     ("py:class", "astropy.table.table.Table"),
     ("py:class", "eventio.simtel.simtelfile.SimTelFile"),
@@ -215,7 +222,7 @@
 html_theme = "pydata_sphinx_theme"
 
 intersphinx_mapping = {
-    "python": ("https://docs.python.org/3.8", None),
+    "python": ("https://docs.python.org/3.11", None),
     "numpy": ("https://numpy.org/doc/stable/", None),
     "scipy": ("https://docs.scipy.org/doc/scipy/", None),
     "astropy": ("https://docs.astropy.org/en/stable/", None),
@@ -225,7 +232,7 @@
     "cython": ("https://docs.cython.org/en/latest/", None),
     "iminuit": ("https://iminuit.readthedocs.io/en/latest/", None),
     "traitlets": ("https://traitlets.readthedocs.io/en/stable/", None),
-    "ctapipe": ("https://ctapipe.readthedocs.io/en/v0.12.0-rtd/", None),
+    "ctapipe": ("https://ctapipe.readthedocs.io/en/v0.19.2/", None),
     "pyirf": ("https://pyirf.readthedocs.io/en/stable/", None),
 }
 
diff --git a/environment.yml b/environment.yml
index 95dd946c..4cb1d870 100644
--- a/environment.yml
+++ b/environment.yml
@@ -3,37 +3,45 @@ name: magic-lst
 channels:
   - conda-forge
 dependencies:
-  - python=3.8
+  - python=3.11
   - pip
+  - numpy
+  - astropy=5
   - black
   - nbsphinx
-  - numpy=1.21
-  - ctapipe=0.12
-  - gammapy=0.19.0
-  - cython
+  - ctapipe=0.19.2
+  - gammapy=1.1
+  - ctaplot~=0.6.4
+  - h5py
+  - ipython
+  - jupyter
+  - matplotlib=3.7
+  - notebook
   - graphviz
   - iminuit>=2
-  - joblib
+  - joblib~=1.2.0
   - numpydoc
+  - jinja2~=3.0.2
   - jupyterlab=3.5.0
-  - pandas
+  - pandas>=2
   - pre-commit
-  - protozfits=2.0
+  - toml
+  - protobuf=3.20
+  - pydata-sphinx-theme
+  - pyparsing
   - psutil
   - pytest
   - pytest-cov
   - pytest-runner
   - pyyaml
-  - scikit-learn
+  - scikit-learn=1.2
   - sphinx
   - sphinx-automodapi
   - sphinx-design
   - sphinx_rtd_theme
-  - pydata-sphinx-theme
-  - uproot~=4.1
-  - corsikaio
-  - pip:
-      - lstchain~=0.9.6
-      - ctapipe_io_magic~=0.4.7
-      - ctaplot~=0.5.5
-      - pyirf~=0.6.0
+  - numba>=0.56
+  - uproot~=5.0
+  - pymongo
+  - pyirf~=0.10.0
+  - ctapipe-io-magic~=0.5.4
+  - lstchain~=0.10.5
diff --git a/magicctapipe/io/containers.py b/magicctapipe/io/containers.py
index b2840ac8..d6a68f7d 100644
--- a/magicctapipe/io/containers.py
+++ b/magicctapipe/io/containers.py
@@ -26,6 +26,9 @@ class BaseEventInfoContainer(Container):
     pointing_az = Field(DEFAULT_RAD, "Azimuth of the pointing direction", unit="rad")
     n_pixels = Field(DEFAULT_VALUE, "Number of pixels after the image cleaning")
     n_islands = Field(DEFAULT_VALUE, "Number of islands after the image cleaning")
+    tels_with_trigger = Field(
+        None, "Integer indicating the combination of telescopes that triggered"
+    )
 
 
 class RealEventInfoContainer(BaseEventInfoContainer):
diff --git a/magicctapipe/io/gadf.py b/magicctapipe/io/gadf.py
index f92b0317..1cf9bb19 100644
--- a/magicctapipe/io/gadf.py
+++ b/magicctapipe/io/gadf.py
@@ -11,10 +11,8 @@
 from astropy.time import Time
 from pyirf.binning import split_bin_lo_hi
 
-from magicctapipe import __version__
-from magicctapipe.utils.functions import HEIGHT_ORM, LAT_ORM, LON_ORM
-
-# from .io import telescope_combinations
+from .. import __version__ as MCP_VERSION
+from ..utils.functions import HEIGHT_ORM, LAT_ORM, LON_ORM
 
 __all__ = [
     "create_gh_cuts_hdu",
@@ -71,7 +69,7 @@ def create_gh_cuts_hdu(gh_cuts, reco_energy_bins, fov_offset_bins, **header_card
     # Create a header
     header = fits.Header(
         cards=[
-            ("CREATOR", f"magicctapipe v{__version__}"),
+            ("CREATOR", f"magicctapipe v{MCP_VERSION}"),
             ("HDUCLAS1", "RESPONSE"),
             ("HDUCLAS2", "GH_CUTS"),
             ("HDUCLAS3", "POINT-LIKE"),
@@ -126,12 +124,11 @@ def create_event_hdu(
         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_M1": [1, 2],  # combo_type = 0
+        "LST1_M1_M2": [1, 2, 3],  # combo_type = 1
         "LST1_M2": [1, 3],  # combo_type = 2
-        "LST1_M1_M2": [1, 2, 3],  # combo_type = 3
+        "M1_M2": [2, 3],  # combo_type = 3
     }  # TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
-
     mjdreff, mjdrefi = np.modf(MJDREF.mjd)
 
     time_start = Time(event_table["timestamp"][0], format="unix", scale="utc")
diff --git a/magicctapipe/io/io.py b/magicctapipe/io/io.py
index ab59df4f..4c7119b2 100644
--- a/magicctapipe/io/io.py
+++ b/magicctapipe/io/io.py
@@ -54,8 +54,8 @@
 # telescope pointing directions but have the same observation ID
 GROUP_INDEX_TRAIN = ["obs_id", "event_id", "true_alt", "true_az"]
 
-# The LST nominal and effective focal lengths
-NOMINAL_FOCLEN_LST = 28 * u.m
+# The LST equivalent and effective focal lengths
+EQUIVALENT_FOCLEN_LST = 28 * u.m
 EFFECTIVE_FOCLEN_LST = 29.30565 * u.m
 
 # The upper limit of the trigger time differences of consecutive events,
@@ -279,10 +279,10 @@ def get_stereo_events_old(
         Data frame of the stereo events surviving the quality cuts
     """
     TEL_COMBINATIONS = {
-        "M1_M2": [2, 3],  # combo_type = 0
-        "LST1_M1": [1, 2],  # combo_type = 1
+        "LST1_M1": [1, 2],  # combo_type = 0
+        "LST1_M1_M2": [1, 2, 3],  # combo_type = 1
         "LST1_M2": [1, 3],  # combo_type = 2
-        "LST1_M1_M2": [1, 2, 3],  # combo_type = 3
+        "M1_M2": [2, 3],  # combo_type = 3
     }  # TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
     event_data_stereo = event_data.copy()
 
@@ -617,7 +617,7 @@ def load_lst_dl1_data_file(input_file):
     # Read the subarray description
     subarray = SubarrayDescription.from_hdf(input_file)
 
-    if focal_length == NOMINAL_FOCLEN_LST:
+    if focal_length == EQUIVALENT_FOCLEN_LST:
         # Set the effective focal length to the subarray description
         subarray.tel[1].optics.equivalent_focal_length = EFFECTIVE_FOCLEN_LST
         subarray.tel[1].camera.geometry.frame = CameraFrame(
@@ -669,10 +669,19 @@ def load_magic_dl1_data_files(input_dir, config):
 
     data_list = []
 
+    # subarray of first file taken as reference for the check
+    subarray = SubarrayDescription.from_hdf(input_files[0])
+
     for input_file in input_files:
         logger.info(input_file)
 
         df_events = pd.read_hdf(input_file, key="events/parameters")
+
+        # check if subarrays are matching
+        subarray_other = SubarrayDescription.from_hdf(input_file)
+        if subarray_other != subarray:
+            raise IOError(f"Subarrays do not match for file: {input_file}")
+
         data_list.append(df_events)
 
     event_data = pd.concat(data_list)
@@ -696,10 +705,6 @@ def load_magic_dl1_data_files(input_dir, config):
     event_data.set_index(["obs_id_magic", "event_id_magic", "tel_id"], inplace=True)
     event_data.sort_index(inplace=True)
 
-    # Read the subarray description from the first input file, assuming
-    # that it is consistent with the others
-    subarray = SubarrayDescription.from_hdf(input_files[0])
-
     return event_data, subarray
 
 
@@ -736,10 +741,10 @@ def load_train_data_files(
         directory
     """
     TEL_COMBINATIONS = {
-        "M1_M2": [2, 3],  # combo_type = 0
-        "LST1_M1": [1, 2],  # combo_type = 1
+        "LST1_M1": [1, 2],  # combo_type = 0
+        "LST1_M1_M2": [1, 2, 3],  # combo_type = 1
         "LST1_M2": [1, 3],  # combo_type = 2
-        "LST1_M1_M2": [1, 2, 3],  # combo_type = 3
+        "M1_M2": [2, 3],  # combo_type = 3
     }  # TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
 
     # Find the input files
@@ -926,15 +931,13 @@ def load_mc_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2)
 
     if event_type == "software":
         # The events of the MAGIC-stereo combination are excluded
-        df_events.query(
-            "(combo_type > 0) & (magic_stereo == True)", inplace=True
-        )  # TODO: change in next PR
+        df_events.query("(combo_type < 3) & (magic_stereo == True)", inplace=True)
 
     elif event_type == "software_only_3tel":
-        df_events.query("combo_type == 3", inplace=True)  # TODO: change in next PR
+        df_events.query("combo_type == 1", inplace=True)
 
     elif event_type == "magic_only":
-        df_events.query("combo_type == 0", inplace=True)  # TODO: change in next PR
+        df_events.query("combo_type == 3", inplace=True)
 
     elif event_type != "hardware":
         raise ValueError(f"Unknown event type '{event_type}'.")
@@ -978,7 +981,7 @@ def load_mc_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2)
     sim_config = pd.read_hdf(input_file, key="simulation/config")
 
     n_total_showers = (
-        sim_config["num_showers"][0]
+        sim_config["n_showers"][0]
         * sim_config["shower_reuse"][0]
         * len(np.unique(event_table["obs_id"]))
     )
@@ -990,9 +993,7 @@ def load_mc_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2)
 
     if viewcone_diff < u.Quantity(0.001, unit="deg"):
         # Handle ring-wobble MCs as same as point-like MCs
-        viewcone = 0 * u.deg
-    else:
-        viewcone = max_viewcone_radius
+        min_viewcone_radius = max_viewcone_radius
 
     sim_info = SimulatedEventsInfo(
         n_showers=n_total_showers,
@@ -1000,7 +1001,8 @@ def load_mc_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2)
         energy_max=u.Quantity(sim_config["energy_range_max"][0], unit="TeV"),
         max_impact=u.Quantity(sim_config["max_scatter_range"][0], unit="m"),
         spectral_index=sim_config["spectral_index"][0],
-        viewcone=viewcone,
+        viewcone_min=min_viewcone_radius,
+        viewcone_max=max_viewcone_radius,
     )
 
     return event_table, pointing, sim_info
@@ -1052,13 +1054,13 @@ def load_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2):
 
     if event_type == "software":
         # The events of the MAGIC-stereo combination are excluded
-        event_data.query("combo_type > 0", inplace=True)  # TODO: change in next PR
+        event_data.query("combo_type < 3", inplace=True)
 
     elif event_type == "software_only_3tel":
-        event_data.query("combo_type == 3", inplace=True)  # TODO: change in next PR
+        event_data.query("combo_type == 1", inplace=True)
 
     elif event_type == "magic_only":
-        event_data.query("combo_type == 0", inplace=True)  # TODO: change in next PR
+        event_data.query("combo_type == 3", inplace=True)
 
     elif event_type == "hardware":
         logger.warning(
@@ -1075,6 +1077,14 @@ def load_dl2_data_file(input_file, quality_cuts, event_type, weight_type_dl2):
     # Get the mean DL2 parameters
     df_dl2_mean = get_dl2_mean(event_data, weight_type_dl2)
     df_dl2_mean.reset_index(inplace=True)
+    event_before_nan_filter = df_dl2_mean.shape[0]
+    df_dl2_mean.dropna(inplace=True)
+    event_after_nan_filter = df_dl2_mean.shape[0]
+
+    logger.info(
+        f"{event_before_nan_filter - event_after_nan_filter} "
+        f"stereo events removed because of NaN values."
+    )
 
     # Convert the pandas data frame to astropy QTable
     event_table = QTable.from_pandas(df_dl2_mean)
@@ -1215,11 +1225,10 @@ def load_irf_files(input_dir_irf):
     for input_file in input_files_irf:
         logger.info(input_file)
         irf_hdus = fits.open(input_file)
-        irf_data["file_names"].append(input_file)
+        irf_data["file_names"] = np.append(irf_data["file_names"], input_file)
 
         # Read the header
         header = irf_hdus["EFFECTIVE AREA"].header
-
         for key in extra_header.keys():
             if key in header:
                 extra_header[key].append(header[key])
@@ -1235,8 +1244,10 @@ def load_irf_files(input_dir_irf):
         aeff_data = irf_hdus["EFFECTIVE AREA"].data[0]
         edisp_data = irf_hdus["ENERGY DISPERSION"].data[0]
 
-        irf_data["effective_area"].append(aeff_data["EFFAREA"])
-        irf_data["energy_dispersion"].append(edisp_data["MATRIX"].T)
+        irf_data["effective_area"].append(aeff_data["EFFAREA"].T)  # ENERGY, THETA
+        irf_data["energy_dispersion"].append(
+            edisp_data["MATRIX"].T
+        )  # ENERGY, MIGRA, THETA
 
         energy_bins = join_bin_lo_hi(aeff_data["ENERG_LO"], aeff_data["ENERG_HI"])
         fov_offset_bins = join_bin_lo_hi(aeff_data["THETA_LO"], aeff_data["THETA_HI"])
diff --git a/magicctapipe/scripts/lst1_magic/config.yaml b/magicctapipe/scripts/lst1_magic/config.yaml
index 4cafd2ba..4d4e9c2f 100644
--- a/magicctapipe/scripts/lst1_magic/config.yaml
+++ b/magicctapipe/scripts/lst1_magic/config.yaml
@@ -1,5 +1,8 @@
 mc_tel_ids:
     LST-1: 1
+    LST-2: 0
+    LST-3: 0
+    LST-4: 0
     MAGIC-I: 2
     MAGIC-II: 3
 
diff --git a/magicctapipe/scripts/lst1_magic/create_dl3_index_files.py b/magicctapipe/scripts/lst1_magic/create_dl3_index_files.py
index 3a4badf4..65d7daae 100644
--- a/magicctapipe/scripts/lst1_magic/create_dl3_index_files.py
+++ b/magicctapipe/scripts/lst1_magic/create_dl3_index_files.py
@@ -59,9 +59,7 @@ def create_dl3_index_files(input_dir):
 
     for input_file in input_files:
         logger.info(input_file)
-
-        input_file_name = Path(input_file).name
-        file_names.append(input_file_name)
+        file_names.append(Path(input_file))
 
     # Create the DL3 index files
     logger.info("\nCreating DL3 index files...")
@@ -70,15 +68,13 @@ def create_dl3_index_files(input_dir):
     obs_index_file = f"{input_dir}/obs-index.fits.gz"
 
     create_hdu_index_hdu(
-        filename_list=file_names,
-        fits_dir=Path(input_dir),
+        file_list=file_names,
         hdu_index_file=Path(hdu_index_file),
         overwrite=True,
     )
 
     create_obs_index_hdu(
-        filename_list=file_names,
-        fits_dir=Path(input_dir),
+        file_list=file_names,
         obs_index_file=Path(obs_index_file),
         overwrite=True,
     )
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py
index 4e88ed70..25e81920 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py
@@ -126,8 +126,10 @@ def create_irf(
     event_table_gamma, pnt_gamma, sim_info_gamma = load_mc_dl2_data_file(
         input_file_gamma, quality_cuts, event_type, weight_type_dl2
     )
-
-    is_diffuse_mc = sim_info_gamma.viewcone.to_value("deg") > 0
+    viewcone = sim_info_gamma.viewcone_max.to_value(
+        "deg"
+    ) - sim_info_gamma.viewcone_min.to_value("deg")
+    is_diffuse_mc = viewcone > 0
     logger.info(f"\nIs diffuse MC: {is_diffuse_mc}")
 
     if is_diffuse_mc:
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_dl1_stereo_to_dl2.py b/magicctapipe/scripts/lst1_magic/lst1_magic_dl1_stereo_to_dl2.py
index ef8ff9f5..6483a917 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_dl1_stereo_to_dl2.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_dl1_stereo_to_dl2.py
@@ -39,10 +39,10 @@
 logger.addHandler(logging.StreamHandler())
 logger.setLevel(logging.INFO)
 TEL_COMBINATIONS = {
-    "M1_M2": [2, 3],  # combo_type = 0
-    "LST1_M1": [1, 2],  # combo_type = 1
+    "LST1_M1": [1, 2],  # combo_type = 0
+    "LST1_M1_M2": [1, 2, 3],  # combo_type = 1
     "LST1_M2": [1, 3],  # combo_type = 2
-    "LST1_M1_M2": [1, 2, 3],  # combo_type = 3
+    "M1_M2": [2, 3],  # combo_type = 3
 }  # TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
 
 
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py b/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py
index e1b3d3e6..3a316840 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py
@@ -34,7 +34,7 @@
 from astropy.io import fits
 from astropy.table import QTable
 from pyirf.cuts import evaluate_binned_cut
-from pyirf.interpolation import interpolate_effective_area_per_energy_and_fov
+from pyirf.interpolation import GridDataInterpolator
 from pyirf.io import (
     create_aeff2d_hdu,
     create_background_2d_hdu,
@@ -164,15 +164,23 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config):
 
     # Interpolate the effective area
     logger.info("\nInterpolating the effective area...")
-
-    aeff_interp = interpolate_effective_area_per_energy_and_fov(
-        effective_area=irf_data["effective_area"],
-        grid_points=irf_data["grid_points"],
-        target_point=target_point,
-        method=interpolation_method,
-    )
-
-    aeff_interp = aeff_interp[:, :, 0]  # Remove the dimension of the grid points
+    if len(irf_data["grid_points"]) > 2:
+        effective_area = irf_data["effective_area"].to_value("m2")
+        # due to large changes in Aeff we will interpolate in log space,
+        # meaning that we need to set a minimal value (1 m^2)
+        effective_area[effective_area < 1] = 1
+        interpolator = GridDataInterpolator(
+            grid_points=irf_data["grid_points"],
+            params=np.log(effective_area),
+            method=interpolation_method,
+        )
+        aeff_interp = np.exp(interpolator.interpolate(target_point)[0])
+        # setting values below the minimal value back to 0, allowing for 10% error margin
+        aeff_interp[aeff_interp < 1.1] = 0
+        aeff_interp *= u.Unit("m2")
+    else:
+        aeff_interp = irf_data["effective_area"][0]
+        print("skipping interpolation since only one point is given")
 
     aeff_hdu = create_aeff2d_hdu(
         effective_area=aeff_interp,
@@ -182,19 +190,21 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config):
         extname="EFFECTIVE AREA",
         **extra_header,
     )
-
     hdus.append(aeff_hdu)
 
-    # Interpolate the energy dispersion with a custom way, since there
-    # is a bug in the function of pyirf v0.6.0 about the renormalization
+    # Interpolate the energy dispersion with a custom way,
+    # TBD: use pyirf quantile interpolation instead
     logger.info("Interpolating the energy dispersion...")
 
-    edisp_interp = griddata(
-        points=irf_data["grid_points"],
-        values=irf_data["energy_dispersion"],
-        xi=target_point,
-        method=interpolation_method,
-    )
+    if len(irf_data["grid_points"]) > 2:
+        edisp_interp = griddata(
+            points=irf_data["grid_points"],
+            values=irf_data["energy_dispersion"],
+            xi=target_point,
+            method=interpolation_method,
+        )
+    else:
+        edisp_interp = irf_data["energy_dispersion"]
 
     edisp_interp = edisp_interp[0]  # Remove the dimension of the grid points
 
@@ -205,6 +215,10 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config):
         edisp_interp, norm, out=np.zeros_like(edisp_interp), where=mask_zeros
     )
 
+    # according to GDAF standard migration matrix is normalized to integral of bins
+    widths = np.diff(irf_data["migration_bins"])
+    edisp_interp /= widths[np.newaxis, :, np.newaxis]
+
     edisp_hdu = create_energy_dispersion_hdu(
         energy_dispersion=edisp_interp,
         true_energy_bins=irf_data["energy_bins"],
@@ -280,12 +294,15 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config):
         # Interpolate the dynamic gammaness cuts
         logger.info("Interpolating the dynamic gammaness cuts...")
 
-        gh_cuts_interp = griddata(
-            points=irf_data["grid_points"],
-            values=irf_data["gh_cuts"],
-            xi=target_point,
-            method=interpolation_method,
-        )
+        if len(irf_data["grid_points"]) > 2:
+            gh_cuts_interp = griddata(
+                points=irf_data["grid_points"],
+                values=irf_data["gh_cuts"],
+                xi=target_point,
+                method=interpolation_method,
+            )
+        else:
+            gh_cuts_interp = irf_data["gh_cuts"]
 
         # Remove the dimension of the grid points
         gh_cuts_interp = gh_cuts_interp[0]
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py
index cd230b5f..db259404 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py
@@ -38,7 +38,7 @@
     number_of_islands,
     timing_parameters,
 )
-from ctapipe.instrument import SubarrayDescription
+from ctapipe.instrument import FocalLengthKind, SubarrayDescription
 from ctapipe.io import EventSource, HDF5TableWriter
 from traitlets.config import Config
 
@@ -71,7 +71,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
     config : dict
         Configuration for the LST-1 + MAGIC analysis
     focal_length : str
-        Focal length choice, effective or nominal
+        Focal length choice, effective or equivalent
     """
 
     assigned_tel_ids = config[
@@ -86,6 +86,16 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
     # Load the input file
     logger.info(f"\nInput file: {input_file}")
 
+    if focal_length == "effective":
+        focal_length = FocalLengthKind.EFFECTIVE
+    elif focal_length == "equivalent":
+        focal_length = FocalLengthKind.EQUIVALENT
+    else:
+        raise ValueError(
+            f"Accepted choices for focal_length are 'effective' or 'equivalent'.\n"
+            f"Chosen value: {focal_length}."
+        )
+
     event_source = EventSource(
         input_file,
         allowed_tels=list(
@@ -171,7 +181,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
     # Prepare for saving data to an output file
     Path(output_dir).mkdir(exist_ok=True, parents=True)
 
-    sim_config = event_source.simulation_config
+    sim_config = event_source.simulation_config[obs_id]
     corsika_inputcard = event_source.file_.corsika_inputcards[0].decode()
 
     regex = r".*\nPRMPAR\s+(\d+)\s+.*"
@@ -227,7 +237,6 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
             )  # If both have trigger, then magic_stereo = True
 
             for tel_id in tels_with_trigger:
-
                 if (
                     tel_id in LSTs_IDs
                 ):  # If the ID is in the LST list, we call calibrate on the LST()
@@ -335,6 +344,22 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
                     lat2=event.simulation.shower.alt,
                 )
 
+                tels_with_trigger = np.intersect1d(
+                    tels_with_trigger,
+                    np.concatenate(
+                        (LSTs_IDs[LSTs_IDs != 0], MAGICs_IDs[MAGICs_IDs != 0])
+                    ),
+                    assume_unique=True,
+                ).tolist()
+
+                # encode tels_with_trigger as an int value
+                # that can be decoded later as a binary
+                # tels_with_trigger = sum_{tel_id} 2**tel_id
+                # where tel_id is only for those triggered
+                tels_with_trigger_binary_int = np.array(
+                    [2 ** (tel_id) for tel_id in tels_with_trigger]
+                ).sum()
+
                 # Set the event information
                 event_info = SimEventInfoContainer(
                     obs_id=event.index.obs_id,
@@ -352,6 +377,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, focal_length):
                     n_pixels=n_pixels,
                     n_islands=n_islands,
                     magic_stereo=magic_stereo,
+                    tels_with_trigger=tels_with_trigger_binary_int,
                 )
 
                 # Reset the telescope IDs
@@ -436,9 +462,9 @@ def main():
         "-f",
         dest="focal_length_choice",
         type=str,
-        choices=["nominal", "effective"],
+        choices=["equivalent", "effective"],
         default="effective",
-        help='Choice of focal length, either "effective" or "nominal". The default (and standard) value is "effective"',
+        help='Choice of focal length, either "effective" or "equivalent". The default (and standard) value is "effective"',
     )
 
     args = parser.parse_args()  # Here we select all 3 parameters collected above
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_muon_analysis.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_muon_analysis.py
index 42dd0308..9e4f0bef 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_muon_analysis.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_muon_analysis.py
@@ -245,7 +245,7 @@ def mc_dl0_to_muons(input_file, output_dir, config, plots_path):
                 if use_charge_correction:
                     # Scale the charges of the DL1 image by the correction factor:
                     event.dl1.tel[tel_id].image *= config_magic["charge_correction"][
-                        "correction_factor"
+                        "factor"
                     ]
 
                 # Apply the image cleaning:
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py
index 3341523b..d5a13e76 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_stereo_reco.py
@@ -41,6 +41,9 @@
     ArrayEventContainer,
     CameraHillasParametersContainer,
     ImageParametersContainer,
+    LeakageContainer,
+    MorphologyContainer,
+    TimingParametersContainer,
 )
 from ctapipe.instrument import SubarrayDescription
 from ctapipe.reco import HillasReconstructor
@@ -142,7 +145,6 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only_analysis=Fa
     logger.info(f"\nInput file: {input_file}")
 
     event_data = pd.read_hdf(input_file, key="events/parameters")
-
     # It sometimes happens that there are MAGIC events whose event and
     # telescope IDs are duplicated, so here we exclude those events
     event_data.drop_duplicates(
@@ -253,7 +255,9 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only_analysis=Fa
         df_evt = event_data.loc[(obs_id, event_id, slice(None))]
 
         # Loop over every telescope
-        tel_ids = df_evt.index.get_level_values("tel_id")
+        tel_ids = df_evt.index.get_level_values("tel_id").values
+
+        event.trigger.tels_with_trigger = tel_ids
 
         for tel_id in tel_ids:
 
@@ -269,17 +273,43 @@ def stereo_reconstruction(input_file, output_dir, config, magic_only_analysis=Fa
                 y=u.Quantity(df_tel["y"], unit="m"),
                 r=u.Quantity(df_tel["r"], unit="m"),
                 phi=Angle(df_tel["phi"], unit="deg"),
+                psi=Angle(df_tel["psi"], unit="deg"),
                 length=u.Quantity(df_tel["length"], unit="m"),
+                length_uncertainty=u.Quantity(df_tel["length_uncertainty"], unit="m"),
                 width=u.Quantity(df_tel["width"], unit="m"),
-                psi=Angle(df_tel["psi"], unit="deg"),
+                width_uncertainty=u.Quantity(df_tel["width_uncertainty"], unit="m"),
                 skewness=float(df_tel["skewness"]),
                 kurtosis=float(df_tel["kurtosis"]),
             )
+            time_par = TimingParametersContainer(
+                slope=u.Quantity(df_tel["slope"], unit="1/deg"),
+                intercept=float(df_tel["intercept"]),
+            )
+            leak_par = LeakageContainer(
+                intensity_width_1=float(df_tel["intensity_width_1"]),
+                intensity_width_2=float(df_tel["intensity_width_2"]),
+                pixels_width_1=float(df_tel["pixels_width_1"]),
+                pixels_width_2=float(df_tel["pixels_width_2"]),
+            )
 
-            event.dl1.tel[tel_id].parameters = ImageParametersContainer(
-                hillas=hillas_params
+            morph_par = MorphologyContainer(
+                n_islands=int(df_tel["n_islands"]),
+                # n_large_islands=int(df_morph_tel["n_large_islands"]),
+                # n_medium_islands=int(df_morph_tel["n_medium_islands"]),
+                # n_small_islands=int(df_morph_tel["n_small_islands"]),
+                n_large_islands=np.nan,
+                n_medium_islands=np.nan,
+                n_small_islands=np.nan,
+                n_pixels=int(df_tel["n_pixels"]),
             )
 
+            event.dl1.tel[tel_id].parameters = ImageParametersContainer(
+                hillas=hillas_params,
+                timing=time_par,
+                leakage=leak_par,
+                morphology=morph_par,
+            )
+            # event.dl1.tel[tel_id].is_valid=True
         # Reconstruct the stereo parameters
         hillas_reconstructor(event)
 
diff --git a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py
index 42416855..d2c52967 100644
--- a/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py
+++ b/magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py
@@ -70,6 +70,12 @@
 
 # The pedestal types to find bad RMS pixels
 PEDESTAL_TYPES = ["fundamental", "from_extractor", "from_extractor_rndm"]
+TEL_COMBINATIONS = {
+    "LST1_M1": [1, 2],  # combo_type = 0
+    "LST1_M1_M2": [1, 2, 3],  # combo_type = 1
+    "LST1_M2": [1, 3],  # combo_type = 2
+    "M1_M2": [2, 3],  # combo_type = 3
+}  # TODO: REMOVE WHEN SWITCHING TO THE NEW RFs IMPLEMENTTATION (1 RF PER TELESCOPE)
 
 
 def magic_calib_to_dl1(input_file, output_dir, config, max_events, process_run=False):
@@ -104,7 +110,6 @@ def magic_calib_to_dl1(input_file, output_dir, config, max_events, process_run=F
 
     obs_id = event_source.obs_ids[0]
     tel_id = event_source.telescope
-
     logger.info(f"\nObservation ID: {obs_id}")
     logger.info(f"Telescope ID: {tel_id}")
 
@@ -162,6 +167,43 @@ def magic_calib_to_dl1(input_file, output_dir, config, max_events, process_run=F
             subrun_id = event_source.metadata["subrun_number"][0]
             output_file = f"{output_dir}/dl1_M{tel_id}.Run{obs_id:08}.{subrun_id:03}.h5"
 
+    assigned_tel_ids = config[
+        "mc_tel_ids"
+    ]  # This variable becomes the dictionary {'LST-1': 1, 'MAGIC-I': 2, 'MAGIC-II': 3} or similar
+
+    # Reset the telescope IDs of the subarray description
+    if event_source.is_hast:
+        tel_positions_magic_lst = {
+            assigned_tel_ids["LST-1"]: [-8.09, 77.13, 0.78] * u.m,  # LST-1
+            assigned_tel_ids["MAGIC-I"]: [39.3, -62.55, -0.97] * u.m,  # MAGIC-I
+            assigned_tel_ids["MAGIC-II"]: [-31.21, -14.57, 0.2] * u.m,  # MAGIC-II
+        }
+
+        tel_descriptions_magic_lst = {
+            # dummy telescope description for LST-1, same as MAGIC-I
+            assigned_tel_ids["LST-1"]: subarray.tel[1],
+            assigned_tel_ids["MAGIC-I"]: subarray.tel[1],  # MAGIC-I
+            assigned_tel_ids["MAGIC-II"]: subarray.tel[2],  # MAGIC-II
+        }
+
+        subarray_magic = SubarrayDescription(
+            "MAGIC-LST-Array", tel_positions_magic_lst, tel_descriptions_magic_lst
+        )
+    else:
+        tel_positions_magic = {
+            assigned_tel_ids["MAGIC-I"]: subarray.positions[1],  # MAGIC-I
+            assigned_tel_ids["MAGIC-II"]: subarray.positions[2],  # MAGIC-II
+        }
+
+        tel_descriptions_magic = {
+            assigned_tel_ids["MAGIC-I"]: subarray.tel[1],  # MAGIC-I
+            assigned_tel_ids["MAGIC-II"]: subarray.tel[2],  # MAGIC-II
+        }
+
+        subarray_magic = SubarrayDescription(
+            "MAGIC-Array", tel_positions_magic, tel_descriptions_magic
+        )
+
     # Loop over every shower event
     logger.info("\nProcessing the events...")
 
@@ -303,10 +345,46 @@ def magic_calib_to_dl1(input_file, output_dir, config, max_events, process_run=F
 
             # Reset the telescope IDs
             if tel_id == 1:
-                event_info.tel_id = config["mc_tel_ids"]["MAGIC-I"]  # MAGIC-I
+                event_info.tel_id = assigned_tel_ids["MAGIC-I"]  # MAGIC-I
 
             elif tel_id == 2:
-                event_info.tel_id = config["mc_tel_ids"]["MAGIC-II"]  # MAGIC-II
+                event_info.tel_id = assigned_tel_ids["MAGIC-II"]  # MAGIC-II
+
+            # M1+M2 trigger
+            if event.trigger.tels_with_trigger == [1, 2]:
+                tels_with_trigger_magic_lst = [
+                    assigned_tel_ids["MAGIC-I"],
+                    assigned_tel_ids["MAGIC-II"],
+                ]
+            # M2+LST1 trigger
+            elif event.trigger.tels_with_trigger == [2, 3]:
+                tels_with_trigger_magic_lst = [
+                    assigned_tel_ids["LST-1"],
+                    assigned_tel_ids["MAGIC-II"],
+                ]
+            # M1+LST1 trigger
+            elif event.trigger.tels_with_trigger == [1, 3]:
+                tels_with_trigger_magic_lst = [
+                    assigned_tel_ids["LST-1"],
+                    assigned_tel_ids["MAGIC-I"],
+                ]
+            # M1+M2+LST1 trigger
+            elif event.trigger.tels_with_trigger == [1, 2, 3]:
+                tels_with_trigger_magic_lst = [
+                    assigned_tel_ids["LST-1"],
+                    assigned_tel_ids["MAGIC-I"],
+                    assigned_tel_ids["MAGIC-II"],
+                ]
+
+            # encode tels_with_trigger as an int value
+            # that can be decoded later as a binary
+            # tels_with_trigger = sum_{tel_id} 2**tel_id
+            # where tel_id is only for those triggered
+            tels_with_trigger_binary_int = np.array(
+                [2 ** (tel_id) for tel_id in tels_with_trigger_magic_lst]
+            ).sum()
+
+            event_info.tels_with_trigger = tels_with_trigger_binary_int
 
             # Save the parameters to the output file
             writer.write(
@@ -316,21 +394,6 @@ def magic_calib_to_dl1(input_file, output_dir, config, max_events, process_run=F
         n_events_processed = event.count + 1
         logger.info(f"\nIn total {n_events_processed} events are processed.")
 
-    # Reset the telescope IDs of the subarray description
-    tel_positions_magic = {
-        config["mc_tel_ids"]["MAGIC-I"]: subarray.positions[1],  # MAGIC-I
-        config["mc_tel_ids"]["MAGIC-II"]: subarray.positions[2],  # MAGIC-II
-    }
-
-    tel_descriptions_magic = {
-        config["mc_tel_ids"]["MAGIC-I"]: subarray.tel[1],  # MAGIC-I
-        config["mc_tel_ids"]["MAGIC-II"]: subarray.tel[2],  # MAGIC-II
-    }
-
-    subarray_magic = SubarrayDescription(
-        "MAGIC-Array", tel_positions_magic, tel_descriptions_magic
-    )
-
     # Save the subarray description
     subarray_magic.to_hdf(output_file)
 
diff --git a/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py
index 2e4a14f3..b085e043 100644
--- a/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py
+++ b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py
@@ -57,7 +57,7 @@ def magic_muons_from_cal(input_file, output_dir, config, process_run, plots_path
 
     event_source = EventSource(input_url=input_file)
     subarray = event_source.subarray
-    obs_id = event_source.obs_ids[0]
+    obs_id = event_source.obs_ids
     tel_id = event_source.telescope
 
     # Create the table which will contain the selected muon ring parameters
diff --git a/magicctapipe/utils/__init__.py b/magicctapipe/utils/__init__.py
index 90c63b6c..f84f0b77 100644
--- a/magicctapipe/utils/__init__.py
+++ b/magicctapipe/utils/__init__.py
@@ -1,5 +1,5 @@
 from .badpixels import MAGICBadPixelsCalc
-from .camera_geometry import reflected_camera_geometry, scale_camera_geometry
+from .camera_geometry import reflected_camera_geometry_mars
 from .functions import (
     HEIGHT_ORM,
     LAT_ORM,
@@ -19,8 +19,7 @@
 
 __all__ = [
     "MAGICBadPixelsCalc",
-    "scale_camera_geometry",
-    "reflected_camera_geometry",
+    "reflected_camera_geometry_mars",
     "identify_time_edges",
     "intersect_time_intervals",
     "GTIGenerator",
diff --git a/magicctapipe/utils/camera_geometry.py b/magicctapipe/utils/camera_geometry.py
index 4a109128..a9e9673c 100644
--- a/magicctapipe/utils/camera_geometry.py
+++ b/magicctapipe/utils/camera_geometry.py
@@ -1,48 +1,15 @@
 """
 Utilities for camera geometry transformations
 """
+import astropy.units as u
 from ctapipe.instrument import CameraGeometry
 
 __all__ = [
-    "scale_camera_geometry",
-    "reflected_camera_geometry",
+    "reflected_camera_geometry_mars",
 ]
 
 
-def scale_camera_geometry(camera_geom, factor):
-    """Scale given camera geometry of a given (constant) factor
-
-    Parameters
-    ----------
-    camera_geom : ctapipe.instrument.camera.geometry.CameraGeometry
-        Camera geometry
-    factor : float
-        Scale factor
-
-    Returns
-    -------
-    ctapipe.instrument.camera.geometry.CameraGeometry
-        Scaled camera geometry
-    """
-    pix_x_scaled = factor * camera_geom.pix_x
-    pix_y_scaled = factor * camera_geom.pix_y
-    pix_area_scaled = camera_geom.guess_pixel_area(
-        pix_x_scaled, pix_y_scaled, camera_geom.pix_type
-    )
-
-    return CameraGeometry(
-        camera_name="MAGICCam",
-        pix_id=camera_geom.pix_id,
-        pix_x=pix_x_scaled,
-        pix_y=pix_y_scaled,
-        pix_area=pix_area_scaled,
-        pix_type=camera_geom.pix_type,
-        pix_rotation=camera_geom.pix_rotation,
-        cam_rotation=camera_geom.cam_rotation,
-    )
-
-
-def reflected_camera_geometry(camera_geom):
+def reflected_camera_geometry_mars(camera_geom):
     """Reflect camera geometry (x->-y, y->-x)
 
     Parameters
@@ -57,14 +24,12 @@ def reflected_camera_geometry(camera_geom):
     """
 
     return CameraGeometry(
-        camera_name="MAGICCam",
+        name="MAGICCam",
         pix_id=camera_geom.pix_id,
         pix_x=-1.0 * camera_geom.pix_y,
         pix_y=-1.0 * camera_geom.pix_x,
-        pix_area=camera_geom.guess_pixel_area(
-            camera_geom.pix_x, camera_geom.pix_y, camera_geom.pix_type
-        ),
+        pix_area=camera_geom.pix_area,
         pix_type=camera_geom.pix_type,
-        pix_rotation=camera_geom.pix_rotation,
+        pix_rotation=30 * u.deg - camera_geom.pix_rotation,
         cam_rotation=camera_geom.cam_rotation,
     )
diff --git a/pyproject.toml b/pyproject.toml
index a5dbeee6..b86188b7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,10 @@
 [build-system]
-requires = ["setuptools >= 40.6.0", "wheel", "setuptools_scm[toml]>=3.4"]
+requires = ["setuptools>=64", "setuptools_scm[toml]>=8"]
 build-backend = "setuptools.build_meta"
 
+[tool.setuptools_scm]
+write_to = "magicctapipe/_version.py"
+
 [tool.black]
 line-length = 88
 target-version = ['py37', 'py38', 'py39', 'py310']
diff --git a/setup.cfg b/setup.cfg
index c2a31152..d03842e9 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -19,26 +19,89 @@ classifiers =
     Topic :: Scientific/Engineering :: Astronomy
     Topic :: Scientific/Engineering :: Physics
     Programming Language :: Python :: 3 :: Only
-    Programming Language :: Python :: 3.8
+    Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
+    Programming Language :: Python :: 3.11
 
 [options]
 packages = find:
-python_requires = >=3.8
+python_requires = >=3.9
 zip_safe = False
-install_requires=
-    lstchain~=0.9.6
-    ctapipe~=0.12.0
-    ctapipe_io_magic~=0.4.7
-    ctaplot~=0.5.5
-    gammapy~=0.19.0
-    uproot~=4.1
-    numpy<1.22.0a0
-    joblib
+install_requires =
+    ctapipe ~=0.19.2
+    lstchain ~=0.10.5
+    ctaplot ~=0.6.4
+    ctapipe_io_magic ~=0.5.4
+    gammapy ~=1.1
+    uproot ~=5.0
     pandas
-    pyirf~=0.6.0
+    pyirf ~=0.10.0
     seaborn
-    scikit-learn
     setuptools_scm
+    astropy
+    bokeh
+    eventio >=1.9.1, <2.0.0a0
+    h5py
+    iminuit >=2
+    joblib ~=1.2.0
+    matplotlib ~=3.0
+    numba
+    numpy
+    protobuf ~=3.20.0
+    scipy
+    scikit-learn ~=1.2
+    tables
+    toml
+    pymongo
+    pyparsing
+    jinja2 ~=3.0.2
+
+[options.packages.find]
+exclude =
+    magicctapipe._dev_version
+
+[options.extras_require]
+tests =
+    pytest
+docs =
+    sphinx
+    pydata_sphinx_theme
+    sphinx_automodapi
+    sphinx_argparse
+    nbsphinx
+    numpydoc
+    sphinx-design
+    jupyter
+    notebook
+    pandas
+    ipython
+dev =
+    setuptools_scm[toml]
+all =
+    %(tests)s
+    %(docs)s
+    %(dev)s
+
+[options.entry_points]
+console_scripts =
+    create_dl3_index_files = magicctapipe.scripts.lst1_magic.create_dl3_index_files:main
+    lst1_magic_create_irf = magicctapipe.scripts.lst1_magic.lst1_magic_create_irf:main
+    lst1_magic_dl1_stereo_to_dl2 = magicctapipe.scripts.lst1_magic.lst1_magic_dl1_stereo_to_dl2:main
+    lst1_magic_dl2_to_dl3 = magicctapipe.scripts.lst1_magic.lst1_magic_dl2_to_dl3:main
+    lst1_magic_event_coincidence = magicctapipe.scripts.lst1_magic.lst1_magic_event_coincidence:main
+    lst1_magic_mc_dl0_to_dl1 = magicctapipe.scripts.lst1_magic.lst1_magic_mc_dl0_to_dl1:main
+    lst1_magic_stereo_reco = magicctapipe.scripts.lst1_magic.lst1_magic_stereo_reco:main
+    lst1_magic_train_rfs = magicctapipe.scripts.lst1_magic.lst1_magic_train_rfs:main
+    magic_calib_to_dl1 = magicctapipe.scripts.lst1_magic.magic_calib_to_dl1:main
+    merge_hdf_files = magicctapipe.scripts.lst1_magic.merge_hdf_files:main
+
+[tool:pytest]
+minversion=3.0
+norecursedirs=build docs/_build
+addopts = -v
+
+[aliases]
+test = pytest
 
 [flake8]
 exclude = build,docs,.eggs
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 3ad62a26..00000000
--- a/setup.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#!/usr/bin/env python
-# Licensed under a 3-clause BSD style license - see LICENSE.rst
-
-import os
-
-from setuptools import find_packages, setup
-
-entry_points = {}
-entry_points["console_scripts"] = [
-    "create_dl3_index_files = magicctapipe.scripts.lst1_magic.create_dl3_index_files:main",
-    "lst1_magic_create_irf = magicctapipe.scripts.lst1_magic.lst1_magic_create_irf:main",
-    "lst1_magic_dl1_stereo_to_dl2 = magicctapipe.scripts.lst1_magic.lst1_magic_dl1_stereo_to_dl2:main",
-    "lst1_magic_dl2_to_dl3 = magicctapipe.scripts.lst1_magic.lst1_magic_dl2_to_dl3:main",
-    "lst1_magic_event_coincidence = magicctapipe.scripts.lst1_magic.lst1_magic_event_coincidence:main",
-    "lst1_magic_mc_dl0_to_dl1 = magicctapipe.scripts.lst1_magic.lst1_magic_mc_dl0_to_dl1:main",
-    "lst1_magic_stereo_reco = magicctapipe.scripts.lst1_magic.lst1_magic_stereo_reco:main",
-    "lst1_magic_train_rfs = magicctapipe.scripts.lst1_magic.lst1_magic_train_rfs:main",
-    "magic_calib_to_dl1 = magicctapipe.scripts.lst1_magic.magic_calib_to_dl1:main",
-    "merge_hdf_files = magicctapipe.scripts.lst1_magic.merge_hdf_files:main",
-]
-
-tests_require = ["pytest", "pandas>=0.24.0", "importlib_resources;python_version<'3.9'"]
-
-docs_require = [
-    "sphinx",
-    "sphinx-automodapi",
-    "sphinx_argparse",
-    "sphinx_rtd_theme",
-    "pydata_sphinx_theme",
-    "numpydoc",
-    "nbsphinx",
-    "sphinx-design",
-]
-
-setup(
-    use_scm_version={"write_to": os.path.join("magicctapipe", "_version.py")},
-    packages=find_packages(),
-    install_requires=[
-        "lstchain~=0.9.6",
-        "ctapipe~=0.12.0",
-        "ctapipe_io_magic~=0.4.7",
-        "ctaplot~=0.5.5",
-        "gammapy~=0.19.0",
-        "uproot~=4.1",
-        "numba<=0.56",
-        "numpy<1.22.0a0",
-        "joblib",
-        "pandas",
-        "pyirf~=0.6.0",
-        "seaborn",
-        "scikit-learn",
-        "setuptools_scm",
-    ],
-    entry_points=entry_points,
-    extras_require={
-        "all": tests_require + docs_require,
-        "tests": tests_require,
-        "docs": docs_require,
-    },
-)


Test patch:
diff --git a/magicctapipe/conftest.py b/magicctapipe/conftest.py
index cbf40fb9..2c82ce82 100644
--- a/magicctapipe/conftest.py
+++ b/magicctapipe/conftest.py
@@ -11,7 +11,7 @@
 
 from magicctapipe.io import resource_file
 
-maxjoint = 13000
+maxjoint = 4000
 maxmonly = 500
 
 DL0_gamma_data = [
@@ -35,16 +35,21 @@
     # "20201216_M1_05093711.002_Y_CrabNebula-W0.40+035.root",
     # "20201216_M1_05093711.003_Y_CrabNebula-W0.40+035.root",
     # "20201216_M1_05093711.004_Y_CrabNebula-W0.40+035.root",
-    "20201216_M1_05093711.014_Y_CrabNebula-W0.40+035.root",
+    # "20201216_M1_05093711.014_Y_CrabNebula-W0.40+035.root",
+    "20231109_M1_05110669.001_Y_CrabNebula-W0.40+215.root"
 ]
 DL0_M2_data = [
     # "20201216_M2_05093711.001_Y_CrabNebula-W0.40+035.root",
     # "20201216_M2_05093711.002_Y_CrabNebula-W0.40+035.root",
     # "20201216_M2_05093711.003_Y_CrabNebula-W0.40+035.root",
     # "20201216_M2_05093711.004_Y_CrabNebula-W0.40+035.root",
-    "20201216_M2_05093711.014_Y_CrabNebula-W0.40+035.root",
+    # "20201216_M2_05093711.014_Y_CrabNebula-W0.40+035.root",
+    "20231109_M2_05110669.001_Y_CrabNebula-W0.40+215.root"
+]
+DL1_LST_data = [
+    # "dl1_LST-1.Run15337.0001.h5",
+    "dl1_LST-1.Run15337.0002.h5"
 ]
-DL1_LST_data = ["dl1_LST-1.Run03265.0094.h5"]
 
 """
 Temporary paths
diff --git a/magicctapipe/image/tests/test_calib.py b/magicctapipe/image/tests/test_calib.py
index 6426e0df..530be0ab 100644
--- a/magicctapipe/image/tests/test_calib.py
+++ b/magicctapipe/image/tests/test_calib.py
@@ -1,5 +1,6 @@
 import pytest
 from ctapipe.calib import CameraCalibrator
+from ctapipe.instrument import FocalLengthKind
 from ctapipe.io import EventSource
 from traitlets.config import Config
 
@@ -18,11 +19,12 @@ def tel_id_MAGIC():
 
 
 def test_calibrate_LST(dl0_gamma, config_calib, tel_id_LST):
-
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
 
         obs_id = event_source.obs_ids[0]
@@ -67,11 +69,12 @@ def test_calibrate_LST(dl0_gamma, config_calib, tel_id_LST):
 
 
 def test_calibrate_MAGIC(dl0_gamma, config_calib, tel_id_MAGIC):
-
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
 
         subarray = event_source.subarray
@@ -89,7 +92,6 @@ def test_calibrate_MAGIC(dl0_gamma, config_calib, tel_id_MAGIC):
         config_extractor_magic = {extractor_type_magic: config_magic["image_extractor"]}
         magic_clean = {}
         for k in [1, 2]:
-
             magic_clean[k] = MAGICClean(camera_geoms[k], config_magic["magic_clean"])
         calibrator_magic = CameraCalibrator(
             image_extractor_type=extractor_type_magic,
@@ -122,7 +124,9 @@ def test_calibrate_exc_1(dl0_gamma, config_calib, tel_id_MAGIC):
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
         subarray = event_source.subarray
         config_magic = config_calib["MAGIC"]
@@ -157,7 +161,9 @@ def test_calibrate_exc_2(dl0_gamma, config_calib, tel_id_LST):
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
 
         subarray = event_source.subarray
@@ -200,7 +206,9 @@ def test_calibrate_exc_3(dl0_gamma, config_calib, tel_id_LST):
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
 
         obs_id = event_source.obs_ids[0]
@@ -239,7 +247,9 @@ def test_calibrate_exc_4(dl0_gamma, config_calib, tel_id_MAGIC):
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
         subarray = event_source.subarray
         tel_descriptions = subarray.tel
@@ -280,7 +290,9 @@ def test_calibrate_exc_5(dl0_gamma, config_calib, tel_id_LST):
     assigned_tel_ids = [1, 2, 3]
     for input_file in dl0_gamma:
         event_source = EventSource(
-            input_file, allowed_tels=assigned_tel_ids, focal_length_choice="effective"
+            input_file,
+            allowed_tels=assigned_tel_ids,
+            focal_length_choice=FocalLengthKind.EFFECTIVE,
         )
 
         obs_id = event_source.obs_ids[0]
diff --git a/magicctapipe/io/tests/test_gadf.py b/magicctapipe/io/tests/test_gadf.py
index 5853e08a..c5b31edc 100644
--- a/magicctapipe/io/tests/test_gadf.py
+++ b/magicctapipe/io/tests/test_gadf.py
@@ -62,10 +62,10 @@ def event():
 def test_create_gh_cuts_hdu(gammaness_cut, energy_bins, fov_bins, header):
     g_cuts_fits = create_gh_cuts_hdu(gammaness_cut, energy_bins, fov_bins, **header)
     g_cuts = g_cuts_fits.data
-    assert np.array_equal(g_cuts["ENERG_LO"][0], np.array([0.1, 1, 10]))
-    assert np.array_equal(g_cuts["ENERG_HI"][0], np.array([1, 10, 100]))
+    assert np.allclose(g_cuts["ENERG_LO"][0], np.array([0.1, 1.0, 10.0]))
+    assert np.allclose(g_cuts["ENERG_HI"][0], np.array([1.0, 10.0, 100.0]))
     assert g_cuts["THETA_LO"] == 0
-    assert np.array_equal(g_cuts["GH_CUTS"][0], np.array([0.7, 0.8, 0.9]))
+    assert np.allclose(g_cuts["GH_CUTS"][0], np.array([0.7, 0.8, 0.9]))
     g_head = g_cuts_fits.header
     assert g_head["TELESCOP"] == "CTA-N"
     assert g_head["INSTRUME"] == "LST-1_MAGIC"
diff --git a/magicctapipe/io/tests/test_io.py b/magicctapipe/io/tests/test_io.py
index 84a0905a..e7100ba9 100644
--- a/magicctapipe/io/tests/test_io.py
+++ b/magicctapipe/io/tests/test_io.py
@@ -255,14 +255,9 @@ def test_load_train_data_files_p(p_stereo):
     """
 
     events = load_train_data_files(str(p_stereo[0]))
-    assert list(events.keys()) == [
-        "M1_M2",
-        "LST1_M1",
-        "LST1_M2",
-        "LST1_M1_M2",
-    ]  # TODO: change in next PR
+    assert list(events.keys()) == ["LST1_M1", "LST1_M1_M2", "LST1_M2", "M1_M2"]
     data = events["LST1_M1"]
-    assert np.all(data["combo_type"] == 1)  # TODO: change in next PR
+    assert np.all(data["combo_type"] == 0)
     assert "off_axis" in data.columns
     assert "true_event_class" not in data.columns
 
@@ -273,14 +268,9 @@ def test_load_train_data_files_g(gamma_stereo):
     """
 
     events = load_train_data_files(str(gamma_stereo[0]))
-    assert list(events.keys()) == [
-        "M1_M2",
-        "LST1_M1",
-        "LST1_M2",
-        "LST1_M1_M2",
-    ]  # TODO: change in next PR
+    assert list(events.keys()) == ["LST1_M1", "LST1_M1_M2", "LST1_M2", "M1_M2"]
     data = events["LST1_M1"]
-    assert np.all(data["combo_type"] == 1)  # TODO: change in next PR
+    assert np.all(data["combo_type"] == 0)
     assert "off_axis" in data.columns
     assert "true_event_class" not in data.columns
 
@@ -394,7 +384,7 @@ def test_load_mc_dl2_data_file_opt(p_dl2, gamma_dl2):
     dl2_mc = [p for p in gamma_dl2.glob("*")] + [p for p in p_dl2.glob("*")]
     for file in dl2_mc:
         data_s, _, _ = load_mc_dl2_data_file(str(file), "width>0", "software", "simple")
-        assert np.all(data_s["combo_type"] > 0)  # TODO: change in next PR
+        assert np.all(data_s["combo_type"] < 3)
 
 
 def test_load_mc_dl2_data_file_exc(p_dl2, gamma_dl2):
@@ -612,7 +602,7 @@ def test_load_dl2_data_file_opt(real_dl2):
     """
     for file in real_dl2.glob("*"):
         data_s, _, _ = load_dl2_data_file(str(file), "width>0", "software", "simple")
-        assert np.all(data_s["combo_type"] > 0)  # TODO: change in next PR
+        assert np.all(data_s["combo_type"] < 3)
 
 
 def test_load_dl2_data_file_exc(real_dl2):
diff --git a/magicctapipe/io/tests/test_io_monly.py b/magicctapipe/io/tests/test_io_monly.py
index b26696b6..29e3f710 100644
--- a/magicctapipe/io/tests/test_io_monly.py
+++ b/magicctapipe/io/tests/test_io_monly.py
@@ -84,7 +84,7 @@ def test_load_train_data_files_p(p_stereo_monly):
     events = load_train_data_files(str(p_stereo_monly[0]))
     assert list(events.keys()) == ["M1_M2"]
     data = events["M1_M2"]
-    assert np.all(data["combo_type"] == 0)  # TODO: change in next PR
+    assert np.all(data["combo_type"] == 3)
     assert "off_axis" in data.columns
     assert "true_event_class" not in data.columns
 
@@ -97,7 +97,7 @@ def test_load_train_data_files_g(gamma_stereo_monly):
     events = load_train_data_files(str(gamma_stereo_monly[0]))
     assert list(events.keys()) == ["M1_M2"]
     data = events["M1_M2"]
-    assert np.all(data["combo_type"] == 0)  # TODO: change in next PR
+    assert np.all(data["combo_type"] == 3)
     assert "off_axis" in data.columns
     assert "true_event_class" not in data.columns
 
@@ -215,7 +215,7 @@ def test_load_mc_dl2_data_file_opt(p_dl2_monly, gamma_dl2_monly):
             str(file), "width>0", "magic_only", "simple"
         )
 
-        assert np.all(data_m["combo_type"] == 0)  # TODO: change in next PR
+        assert np.all(data_m["combo_type"] == 3)
 
 
 def test_load_mc_dl2_data_file_exc(p_dl2_monly, gamma_dl2_monly):
@@ -437,7 +437,7 @@ def test_load_dl2_data_file_opt(real_dl2_monly):
 
         data_m, _, _ = load_dl2_data_file(str(file), "width>0", "magic_only", "simple")
 
-        assert np.all(data_m["combo_type"] == 0)  # TODO: change in next PR
+        assert np.all(data_m["combo_type"] == 3)
 
 
 def test_load_dl2_data_file_exc(real_dl2_monly):
diff --git a/magicctapipe/resources/test_config.yaml b/magicctapipe/resources/test_config.yaml
index 8394a031..04b04a73 100644
--- a/magicctapipe/resources/test_config.yaml
+++ b/magicctapipe/resources/test_config.yaml
@@ -71,7 +71,7 @@ MAGIC:
 event_coincidence:
     timestamp_type_lst: "dragon_time"  # select "dragon_time", "tib_time" or "ucts_time"
     window_half_width: "300 ns"
-    pre_offset_search: false
+    pre_offset_search: true
     n_pre_offset_search_events: 100
     time_offset:
         start: "-10 us"
diff --git a/magicctapipe/resources/test_config_general.yaml b/magicctapipe/resources/test_config_general.yaml
index 0a15083c..2a5dcc4e 100644
--- a/magicctapipe/resources/test_config_general.yaml
+++ b/magicctapipe/resources/test_config_general.yaml
@@ -5,7 +5,7 @@ mc_tel_ids:
     LST-4: 0
     MAGIC-I: 2
     MAGIC-II: 3
-    
-general: 
-    focal_length           : "effective" #effective #nominal
-    
+
+general:
+    focal_length           : "effective" #effective #equivalent
+
diff --git a/magicctapipe/resources/test_config_general_4LST.yaml b/magicctapipe/resources/test_config_general_4LST.yaml
index bdca6480..47f1c259 100644
--- a/magicctapipe/resources/test_config_general_4LST.yaml
+++ b/magicctapipe/resources/test_config_general_4LST.yaml
@@ -5,7 +5,7 @@ mc_tel_ids:
     LST-4: 5
     MAGIC-I: 0
     MAGIC-II: 0
-    
+
 general:
-    focal_length           : "effective" #effective #nominal
-    
+    focal_length           : "effective" #effective #equivalent
+
diff --git a/magicctapipe/resources/test_config_monly.yaml b/magicctapipe/resources/test_config_monly.yaml
index 5cc7a085..8fba5ac2 100644
--- a/magicctapipe/resources/test_config_monly.yaml
+++ b/magicctapipe/resources/test_config_monly.yaml
@@ -72,7 +72,7 @@ MAGIC:
 event_coincidence:
     timestamp_type_lst: "dragon_time"  # select "dragon_time", "tib_time" or "ucts_time"
     window_half_width: "300 ns"
-    pre_offset_search: false
+    pre_offset_search: true
     n_pre_offset_search_events: 100
     time_offset:
         start: "-10 us"