Add integration test for DL2-to-DL3 step

docs/contribute/beforepushing.rst

@@ -84,7 +84,7 @@ Same for *pyirf*.
 Integration tests
 ^^^^^^^^^^^^^^^^^
 
-These are neither unit-tests nor benchmarks, but rather functions that tests
+These are neither unit-tests nor benchmarks, but rather functions that test
 whole functionalities and not just single API functions.
 
 In the case of the pipeline, such functionalities are the scripts/tools

protopipe/perf/__init__.py

@@ -1,2 +1,2 @@
 from .temp import *
-from .utils import *
+from .utils import *

protopipe/perf/temp.py

@@ -1,9 +1,9 @@
-"""Temporary development code from pyirf master branch."""
+"""Temporary development code in-between pyirf releases."""
 
-from astropy.table import Table
 import numpy as np
-from scipy.stats import norm
+from astropy.table import Table
 import astropy.units as u
+from scipy.stats import norm
 
 from pyirf.binning import calculate_bin_indices
 from pyirf.benchmarks.energy_bias_resolution import inter_quantile_distance
@@ -47,6 +47,8 @@ def energy_bias_resolution(
     table["bin_index"] = calculate_bin_indices(
         table[f"{energy_type}_energy"].quantity, energy_bins
     )
+    n_bins = len(energy_bins) - 1
+    mask = (table["bin_index"] >= 0) & (table["bin_index"] < n_bins)
 
     result = Table()
     result[f"{energy_type}_energy_low"] = energy_bins[:-1]
@@ -56,8 +58,13 @@ def energy_bias_resolution(
     result["bias"] = np.nan
     result["resolution"] = np.nan
 
+    if not len(events):
+        # if we get an empty input (no selected events available)
+        # we return the table filled with NaNs
+        return result
+
     # use groupby operations to calculate the percentile in each bin
-    by_bin = table.group_by("bin_index")
+    by_bin = table[mask].group_by("bin_index")
 
     index = by_bin.groups.keys["bin_index"]
     result["bias"][index] = by_bin["rel_error"].groups.aggregate(bias_function)
@@ -89,17 +96,22 @@ def angular_resolution(
     Returns
     -------
     result : astropy.table.Table
-        Table containing the 68% containment of the angular
-        distance distribution per each reconstructed energy bin.
+        Table containing the 68% containment of the angular distance
+        distribution per each reconstructed energy bin.
     """
 
+    ONE_SIGMA_QUANTILE = norm.cdf(1) - norm.cdf(-1)
+
     # create a table to make use of groupby operations
     table = Table(events[[f"{energy_type}_energy", "theta"]])
 
     table["bin_index"] = calculate_bin_indices(
         table[f"{energy_type}_energy"].quantity, energy_bins
     )
 
+    n_bins = len(energy_bins) - 1
+    mask = (table["bin_index"] >= 0) & (table["bin_index"] < n_bins)
+
     result = Table()
     result[f"{energy_type}_energy_low"] = energy_bins[:-1]
     result[f"{energy_type}_energy_high"] = energy_bins[1:]
@@ -108,11 +120,10 @@ def angular_resolution(
     result["angular_resolution"] = np.nan * u.deg
 
     # use groupby operations to calculate the percentile in each bin
-    by_bin = table.group_by("bin_index")
+    by_bin = table[mask].group_by("bin_index")
 
     index = by_bin.groups.keys["bin_index"]
-    ONE_SIGMA_PERCENTILE = norm.cdf(1) - norm.cdf(-1)
     result["angular_resolution"][index] = by_bin["theta"].groups.aggregate(
-        lambda x: np.percentile(x, 100 * ONE_SIGMA_PERCENTILE)
+        lambda x: np.quantile(x, ONE_SIGMA_QUANTILE)
     )
     return result

protopipe/perf/utils.py

@@ -1,3 +1,4 @@
+import argparse
 import logging
 
 import numpy as np
@@ -11,6 +12,55 @@
 from pyirf.simulations import SimulatedEventsInfo
 
 
+def initialize_script_arguments():
+    """Initialize the parser of any protopipe.scripts.make_performance_XXX script.
+
+    Returns
+    -------
+    args : argparse.Namespace
+        Populated argparse namespace.
+    """
+
+    parser = argparse.ArgumentParser(description='Make performance files')
+    parser.add_argument('--config_file', type=str, required=True, help='')
+
+    mode_group = parser.add_mutually_exclusive_group()
+    mode_group.add_argument('--wave', dest="mode", action='store_const',
+                            const="wave", default="tail",
+                            help="if set, use wavelet cleaning")
+    mode_group.add_argument('--tail', dest="mode", action='store_const',
+                            const="tail",
+                            help="if set, use tail cleaning (default)")
+    parser.add_argument(
+        "-i",
+        "--indir",
+        type=str,
+        default=None,
+        help="Directory containing the required DL2 input files"
+    )
+    parser.add_argument(
+        "--template_input_file",
+        type=str,
+        default=None,
+        help="template for the filename of DL2 files (default: read from config file)",
+    )
+    parser.add_argument(
+        "--outdir_path",
+        type=str,
+        default=None,
+        help="Output directory for DL3 file"
+    )
+    parser.add_argument(
+        "--out_file_name",
+        type=str,
+        default=None,
+        help="Desired name for DL3 file"
+    )
+
+    args = parser.parse_args()
+
+    return args
+
 def percentiles(values, bin_values, bin_edges, percentile):
     # Seems complicated for vector defined as [inf, inf, .., inf]
     percentiles_binned = np.squeeze(

protopipe/scripts/make_performance_EventDisplay.py

@@ -5,12 +5,11 @@
 import astropy.units as u
 from astropy import table
 from astropy.io import fits
-import argparse
 import numpy as np
 import operator
 
 from protopipe.pipeline.utils import load_config
-from protopipe.perf.utils import read_DL2_pyirf
+from protopipe.perf.utils import initialize_script_arguments, read_DL2_pyirf
 
 from pyirf.binning import (
     add_overflow_bins,
@@ -43,44 +42,51 @@
     create_rad_max_hdu,
     create_background_2d_hdu,
 )
-# from pyirf.benchmarks import energy_bias_resolution, angular_resolution
 from protopipe.perf.temp import energy_bias_resolution, angular_resolution
 
 log = logging.getLogger("pyirf")
 
 
 def main():
 
-    # Read arguments
-    parser = argparse.ArgumentParser(description='Make performance files')
-    parser.add_argument('--config_file', type=str, required=True, help='')
+    # INITIALIZE CLI arguments
+    args = initialize_script_arguments()
 
-    mode_group = parser.add_mutually_exclusive_group()
-    mode_group.add_argument('--wave', dest="mode", action='store_const',
-                            const="wave", default="tail",
-                            help="if set, use wavelet cleaning")
-    mode_group.add_argument('--tail', dest="mode", action='store_const',
-                            const="tail",
-                            help="if set, use tail cleaning (default)")
-
-    args = parser.parse_args()
-
-    # Read configuration file
+    # LOAD CONFIGURATION FILE
     cfg = load_config(args.config_file)
 
-    # Create output directory if necessary
-    outdir = os.path.join(cfg['general']['outdir'], 'performance_protopipe_{}_CTA{}_{}_Zd{}_{}_Time{:.2f}{}'.format(
-        cfg['general']['prod'],
-        cfg['general']['site'],
-        cfg['general']['array'],
-        cfg['general']['zenith'],
-        cfg['general']['azimuth'],
-        cfg['analysis']['obs_time']['value'],
-        cfg['analysis']['obs_time']['unit']),
-    )
-
-    indir = cfg['general']['indir']
-    template_input_file = cfg['general']['template_input_file']
+    # Set final input parameters
+    # Command line as higher priority on configuration file
+    if args.indir is None:
+        indir = cfg["general"]["indir"]
+    else:
+        indir = args.indir
+
+    if args.outdir_path is None:
+        outdir_path = cfg["general"]["outdir"]
+    else:
+        outdir_path = args.outdir_path
+    if not os.path.exists(outdir_path):
+        os.makedirs(outdir_path)
+
+    if args.indir is None:
+        template_input_file = cfg['general']['template_input_file']
+    else:
+        template_input_file = args.template_input_file
+
+    if args.out_file_name is None:
+        out_file_name = 'performance_protopipe_{}_CTA{}_{}_Zd{}_{}_Time{:.2f}{}'.format(
+            cfg['general']['prod'],
+            cfg['general']['site'],
+            cfg['general']['array'],
+            cfg['general']['zenith'],
+            cfg['general']['azimuth'],
+            cfg['analysis']['obs_time']['value'],
+            cfg['analysis']['obs_time']['unit'])
+    else:
+        out_file_name = args.out_file_name
+
+    outdir = os.path.join(outdir_path, out_file_name)
 
     T_OBS = cfg['analysis']['obs_time']['value'] * u.Unit(cfg['analysis']['obs_time']['unit'])
 

protopipe/scripts/tests/test_config_analysis_north.yaml

@@ -59,7 +59,7 @@ ImageCleaning:
     - DigiCam: [0,0] # values left unset for future studies
     - CHEC: [0,0] # values left unset for future studies
     - SCTCam: [0,0] # values left unset for future studies
-   keep_isolated_pixels: False
+   keep_isolated_pixels: True
    min_number_picture_neighbors: 1
 
   wave:
@@ -92,7 +92,11 @@ ImageSelection:
 
 # Minimal number of telescopes to consider events
 Reconstruction:
- min_tel: 2
+  # for events with <2 LST images the single-LST image is removed
+  # before shower geometry
+  LST_stereo: True
+  # after this we check if the remaining images satisfy the min_tel condition
+  min_tel: 2 # any tel_type
 
 # Parameters for energy estimation
 EnergyRegressor:

protopipe/scripts/tests/test_config_analysis_south.yaml

@@ -59,7 +59,7 @@ ImageCleaning:
     - DigiCam: [0,0] # values left unset for future studies
     - CHEC: [0,0] # values left unset for future studies
     - SCTCam: [0,0] # values left unset for future studies
-   keep_isolated_pixels: False
+   keep_isolated_pixels: True
    min_number_picture_neighbors: 1
 
   wave:
@@ -92,7 +92,11 @@ ImageSelection:
 
 # Minimal number of telescopes to consider events
 Reconstruction:
- min_tel: 2
+  # for events with <2 LST images the single-LST image is removed
+  # before shower geometry
+  LST_stereo: False
+  # after this we check if the remaining images satisfy the min_tel condition
+  min_tel: 2 # any tel_type
 
 # Parameters for energy estimation
 EnergyRegressor:

protopipe/scripts/tests/test_performance_ED_prod3b.yaml

@@ -0,0 +1,52 @@
+general:
+ # Directory with input data file
+ # [...] = your analysis local full path OUTSIDE the Vagrant box
+ indir: '[...]/shared_folder/analyses/v0.4.0_dev1/data/DL2'
+ # Template name for output file
+ prod: 'Prod3b'
+ site: 'North'
+ array: 'baseline_full_array'
+ zenith: '20deg'
+ azimuth: '180deg' # 0deg -> north 180deg -> south
+ template_input_file: 'DL2_{}_{}_merged.h5' # filled with mode and particle type
+ # Directory for output files
+ outdir: '[...]/shared_folder/analyses/v0.4.0_dev1/data/DL3'
+
+analysis:
+ obs_time:
+   value: 50
+   unit: 'h'
+ cut_on_multiplicity: 4
+ # Normalisation between ON and OFF regions
+ alpha: 0.2
+
+ # Radius to use for calculating bg rate
+ max_bg_radius: 1.
+
+particle_information:
+ gamma:
+  num_use: 10
+  num_showers: 100000
+  e_min: 0.003
+  e_max: 330
+  gen_radius: 1400
+  gen_gamma: -2
+  diff_cone: 0
+
+ proton:
+  num_use: 20
+  num_showers: 200000
+  e_min: 0.004
+  e_max: 600
+  gen_radius: 1900
+  gen_gamma: -2
+  diff_cone: 10
+
+ electron:
+  num_use: 20
+  num_showers: 100000
+  e_min: 0.003
+  e_max: 330
+  gen_radius: 1900
+  gen_gamma: -2
+  diff_cone: 10