Merge branch 'master' into upgradelstchain

README.md

@@ -21,10 +21,7 @@ The following command will set up a conda virtual environment, add the necessary
 
 # Available scripts
 
-*magic-cta-pipe* scripts to perform the analysis of MAGIC only and MAGIC+LST1 data. Within the *scripts* folder we have the subdirectories:
-
-* *magic*: scripts for the analysis of MAGIC data (**currently under heavy restructuring in order to be used with ctapipe v0.12**)
-* *lst1_magic*: scripts for the MAGIC+LST1 analysis (updated to be used with ctapipe v0.12)
+*magic-cta-pipe* scripts to perform the analysis of MAGIC only and MAGIC+LST1 data. Both types of analysis can be performed using the scripts within the *lst1_magic* folder.
 
 <!--
 A brief description:

environment.yml

@@ -49,7 +49,7 @@ dependencies:
   - conda-forge::corsikaio
   - pip:
       - lstchain~=0.9.6
-      - ctapipe_io_magic~=0.4.1
+      - ctapipe_io_magic~=0.4.4
       - uproot~=4.1
       - ctaplot~=0.5.5
       - pyirf~=0.6.0

magicctapipe/__init__.py

@@ -1,10 +1,10 @@
+from .version import __version__
 from . import image
 from . import irfs
 from . import reco
 from . import scripts
 from . import utils
 
-from .version import __version__
 
 __all__ = [
     "image",

magicctapipe/image/__init__.py

@@ -3,8 +3,6 @@
     PixelTreatment,
     get_num_islands_MAGIC,
     clean_image_params,
-    eval_impact,
-    tailcuts_clean_lstchain,
 )
 
 from .leakage import (
@@ -16,8 +14,5 @@
     'PixelTreatment',
     'get_num_islands_MAGIC',
     'clean_image_params',
-    'eval_impact',
-    'tailcuts_clean_lstchain',
     'get_leakage',
 ]
-

magicctapipe/image/cleaning.py

@@ -1,12 +1,8 @@
 import copy
 import itertools
 import numpy as np
-import astropy.units as u
-from astropy.coordinates import SkyCoord, AltAz
 from scipy.sparse.csgraph import connected_components
 
-from ctapipe.coordinates import GroundFrame
-from ctapipe.core.container import Container, Field
 from ctapipe.image import (
     tailcuts_clean,
     number_of_islands,
@@ -21,10 +17,9 @@
     'PixelTreatment',
     'get_num_islands_MAGIC',
     'clean_image_params',
-    'eval_impact',
-    'tailcuts_clean_lstchain',
 ]
 
+
 class MAGICClean:
 
     def __init__(self, camera, configuration):
@@ -56,20 +51,21 @@ def __init__(self, camera, configuration):
         else:
             self.SumThresh4NNPerPixel = 1.0
 
+        # default WindowXNN values are converted from time slices to ns
         if 'Window2NN' in configuration:
             self.Window2NN = configuration['Window2NN']
         else:
-            self.Window2NN = 0.82
+            self.Window2NN = 0.82/1.639
 
         if 'Window3NN' in configuration:
             self.Window3NN = configuration['Window3NN']
         else:
-            self.Window3NN = 1.15
+            self.Window3NN = 1.15/1.639
 
         if 'Window4NN' in configuration:
             self.Window4NN = configuration['Window4NN']
         else:
-            self.Window4NN = 1.80
+            self.Window4NN = 1.80/1.639
 
         if 'clipping' in configuration:
             self.clipping = configuration['clipping']
@@ -662,7 +658,7 @@ def clean_image_params(geom, image, clean, peakpos):
     # Hillas parameters, same for LST and MAGIC. From ctapipe
     hillas_p = hillas_parameters(geom=geom[clean], image=image[clean])
     # Leakage, same for LST and MAGIC. From ctapipe
-    leakage_p = leakage(geom=geom, image=image, cleaning_mask=clean)
+    leakage_p = leakage_parameters(geom=geom, image=image, cleaning_mask=clean)
     # Timing parameters, same for LST and MAGIC. From ctapipe
     timing_p = timing_parameters(
         geom=geom[clean],
@@ -676,98 +672,3 @@ def clean_image_params(geom, image, clean, peakpos):
     #     time_grad[tel_id] = 1
 
     return hillas_p, leakage_p, timing_p
-
-
-def eval_impact(subarray, hillas_p, stereo_params):
-    """Copied from protopipe"""
-    # Impact parameter for energy estimation (/ tel)
-    ground_frame = GroundFrame()
-    impact_p = {}
-
-    class ImpactContainer(Container):
-        impact = Field(-1, "Impact")
-
-    for tel_id in hillas_p.keys():
-        pos = subarray.positions[tel_id]
-        tel_ground = SkyCoord(pos[0], pos[1], pos[2], frame=ground_frame)
-
-        core_ground = SkyCoord(
-            stereo_params.core_x, stereo_params.core_y, 0 * u.m, frame=ground_frame,
-        )
-        # Should be better handled (tilted frame)
-        impact_ = np.sqrt(
-            (core_ground.x - tel_ground.x) ** 2 + (core_ground.y - tel_ground.y) ** 2
-        )
-        impact_p[tel_id] = ImpactContainer(impact=impact_)
-    return impact_p
-
-
-def tailcuts_clean_lstchain(geom, image, peak_time, input_file, cleaning_parameters):
-    """Apply tailcuts cleaning lstchain mode
-
-    Parameters
-    ----------
-    geom: `ctapipe.instrument.CameraGeometry`
-        Camera geometry information
-    image: array
-        pixel values
-    peak_time : array
-        dl1.peak_time
-    cleaning_parameters : dict
-        dictionary composed by the cleaning parameters for tailcuts and the ones
-        for lstchain tailcuts ('delta_time' and 'use_only_main_island')
-
-    Returns
-    -------
-    tuple
-        signal_pixels, num_islands, island_labels
-    """
-
-    # pop delta_time and use_main_island, so we can cleaning_parameters to tailcuts
-    delta_time = cleaning_parameters.pop("delta_time", None)
-    use_main_island = cleaning_parameters.pop("use_only_main_island", True)
-
-    sigma = cleaning_parameters["sigma"]
-    pedestal_thresh = get_threshold_from_dl1_file(input_file, sigma)
-    picture_thresh_cfg = cleaning_parameters["cleaning_parameters"]
-    print(
-        f"Fraction of pixel cleaning thresholds above picture thr.:"
-        f"{np.sum(pedestal_thresh>picture_thresh_cfg) / len(pedestal_thresh):.3f}"
-    )
-    picture_thresh = np.clip(pedestal_thresh, picture_thresh_cfg, None)
-
-    signal_pixels = tailcuts_clean(
-        geom=geom,
-        image=image,
-        picture_thresh=picture_thresh,
-        boundary_thresh=cleaning_parameters["boundary_thresh"],
-        keep_isolated_pixels=cleaning_parameters["keep_isolated_pixels"],
-        min_number_picture_neighbors=cleaning_parameters["min_n_neighbors"],
-    )
-
-    n_pixels = np.count_nonzero(signal_pixels)
-    if n_pixels > 0:
-        num_islands, island_labels = number_of_islands(camera_geom, signal_pixels)
-        n_pixels_on_island = np.bincount(island_labels.astype(np.int64))
-        # first island is no-island and should not be considered
-        n_pixels_on_island[0] = 0
-        max_island_label = np.argmax(n_pixels_on_island)
-        if use_only_main_island:
-            signal_pixels[island_labels != max_island_label] = False
-
-        # if delta_time has been set, we require at least one
-        # neighbor within delta_time to accept a pixel in the image:
-        if delta_time is not None:
-            cleaned_pixel_times = peak_time
-            # makes sure only signal pixels are used in the time
-            # check:
-            cleaned_pixel_times[~signal_pixels] = np.nan
-            new_mask = apply_time_delta_cleaning(
-                camera_geom, signal_pixels, cleaned_pixel_times, 1, delta_time
-            )
-            signal_pixels = new_mask
-
-        # count the surviving pixels
-        n_pixels = np.count_nonzero(signal_pixels)
-
-    return signal_pixels, num_islands, island_labels