Merge pull request #154 from cta-observatory/Torino

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:
   - conda-forge
 dependencies:

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,11 +43,11 @@
     "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")
@@ -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,52 @@ 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
+
+
+@pytest.fixture(scope="session")
+def config_gen_4lst():
+    config_path = resource_file("test_config_general_4LST.yaml")
+    with open(config_path, "rb") as f:
+        config = yaml.safe_load(f)
+    return config
+
+
+@pytest.fixture(scope="session")
+def config_calib():
+    config_path = resource_file("test_config_calib.yaml")
+    with open(config_path, "rb") as f:
+        config = yaml.safe_load(f)
+    return config
+
+
+@pytest.fixture(scope="session")
+def config_check():
+    config_path = resource_file("test_check_list.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,140 @@
+import numpy as np
+from ctapipe.image import (
+    apply_time_delta_cleaning,   
+    number_of_islands,
+    tailcuts_clean,   
+)
+from ctapipe.instrument import CameraGeometry
+from lstchain.image.cleaning import apply_dynamic_cleaning
+from lstchain.image.modifier import (
+    add_noise_in_pixels,
+    random_psf_smearer,  
+    set_numba_seed  
+)
+from magicctapipe.image import MAGICClean
+
+
+__all__ = [
+    "calibrate"
+]
+
+
+def calibrate(event, tel_id, config, calibrator, is_lst, 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
+    is_lst: 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
+
+    """
+    if (not is_lst) and (magic_clean==None):
+        raise ValueError("Check the provided parameters and the telescope type; MAGIC calibration not possible if magic_clean not provided") 
+    if (is_lst) and (obs_id==None):
+        raise ValueError("Check the provided parameters and the telescope type; LST calibration not possible if obs_id not provided") 
+    if (is_lst) and (camera_geoms==None):
+        raise ValueError("Check the provided parameters and the telescope type; LST calibration not possible if gamera_geoms not provided") 
+    if (not is_lst) and (type(magic_clean[tel_id])!=MAGICClean):
+        raise ValueError("Check the provided magic_clean parameter; MAGIC calibration not possible if magic_clean not a dictionary of MAGICClean objects") 
+    if (is_lst) and (type(camera_geoms[tel_id])!=CameraGeometry):
+        raise ValueError("Check the provided camera_geoms parameter; LST calibration not possible if camera_geoms not a dictionary of CameraGeometry objects")
+
+    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 not is_lst:
+        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
+        )
+
+    else:
+        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
+
+        config["increase_nsb"]["use"]=increase_nsb
+        config["time_delta_cleaning"]["use"]=use_time_delta_cleaning
+        config["dynamic_cleaning"]["use"]=use_dynamic_cleaning        
+
+    return signal_pixels, image, peak_time