Implement IRF interpolation method

magicctapipe/scripts/lst1_magic/config.yaml

@@ -185,8 +185,8 @@ create_irf:
 
     energy_bins:
         start: 0.01   # unit: [TeV]
-        stop: 100   # unit: [TeV]
-        n_bins: 20   # number of bins in log space
+        stop: 1000   # unit: [TeV]
+        n_bins: 25   # number of bins in log space
 
     migration_bins:
         start: 0.2
@@ -217,3 +217,4 @@ dl2_to_dl3:
     source_name: 'Crab'
     source_ra: null   # used when the source name is null
     source_dec: null   # used when the source name is null
+    interpolation_method: 'nearest' # 'nearest' - no interpolation (closest point taken), 'linear', 'cubic' - interpolation

magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py

@@ -371,7 +371,7 @@ def create_irf(
 
     table_gamma, sim_info_gamma = load_dl2_data_file(input_file_gamma, quality_cuts, irf_type, dl2_weight)
 
-    if sim_info_gamma.viewcone.value != 0.0:
+    if sim_info_gamma.viewcone.value > 1.0:
         logger.info('\nHave not yet implemented functions to create diffuse IRFs. Exiting.')
         sys.exit()
 
@@ -592,7 +592,7 @@ def create_irf(
     # Save the data in an output file:
     Path(output_dir).mkdir(exist_ok=True, parents=True)
 
-    regex = r'dl2_gamma_(\S+)_run.*'
+    regex = r'dl2_(\S+)_run.*'
     file_name = Path(input_file_gamma).name
 
     if re.fullmatch(regex, file_name):

magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py

@@ -13,24 +13,33 @@
 --input-file-irf ./data/irf_40deg_90deg_off0.4deg_LST-1_MAGIC_software_gam_global0.6_theta_global0.2.fits.gz
 --output-dir ./data
 --config-file ./config.yaml
+
+if --input-dir-irf is used instead of --input-file-irf the IRFs are obtained from interpolation of the files in this directory
 """
 
 import re
 import sys
+import os
 import time
 import yaml
 import logging
 import argparse
 import operator
+import glob
 import numpy as np
 import pandas as pd
+import pyirf.interpolation as interp
 from pathlib import Path
 from astropy import units as u
 from astropy.io import fits
 from astropy.time import Time
 from astropy.table import QTable
 from astropy.coordinates import SkyCoord
 from pyirf.cuts import evaluate_binned_cut
+from pyirf.io import create_aeff2d_hdu, create_energy_dispersion_hdu
+from pyirf.binning import join_bin_lo_hi
+from scipy.interpolate import griddata
+
 from magicctapipe.utils import (
     get_dl2_mean,
     check_tel_combination,
@@ -55,6 +64,8 @@
     'create_event_list',
     'create_gti_table',
     'create_pointing_table',
+    'read_fits_bins_lo_hi',
+    'interpolate_irf',
     'dl2_to_dl3',
 ]
 
@@ -354,8 +365,152 @@ def create_pointing_table(event_table):
 
     return pnt_table, pnt_header
 
+def read_fits_bins_lo_hi(hdus_irfs, extname, tag):
+    """
+    Reads from a HDUS fits object two arrays of tag_LO and tag_HI.
+    It checks consistency of bins between different files before returning the value
 
-def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
+    Parameters
+    ----------
+    hdus_irfs: list
+        list of HDUS objects with IRFs
+    extname: string
+        name of the extension to read the data from in fits file
+    tag: string
+        name of the field in the extension to extract, _LO and _HI will be added
+
+    Returns
+    -------
+    bins: list of astropy.units.Quantity
+        list of ntuples (LO, HI) of bins (with size of extnames)
+    """
+
+    old_table = None
+    bins = list()
+    tag_lo=tag+'_LO'
+    tag_hi=tag+'_HI'
+    for hdus in hdus_irfs:
+        table = hdus[extname].data[0]
+        if old_table is not None:
+            if not old_table[tag_lo].shape == table[tag_lo].shape:
+                raise ValueError('Non matching bins in ' + extname)
+            if not ((old_table[tag_lo] == table[tag_lo]).all() and (old_table[tag_hi] == table[tag_hi]).all()):
+                raise ValueError('Non matching bins in ' + extname)
+        else:
+            bins.append(join_bin_lo_hi(table[tag_lo], table[tag_hi]))
+        old_table = table
+    bins = u.Quantity(np.array(bins), hdus[extname].columns[tag_lo].unit, copy=False)
+
+    return bins
+
+def interpolate_irf(input_file_dl2, input_irf_dir, method='linear'):
+    """
+    Interpolates a grid of IRFs read from a given directory to a specified DL2 file
+
+    Parameters
+    ----------
+    input_file_dl2: string
+        path to the  DL2 file
+    input_irf_dir: string
+        path to the directory with IRFs, files must follow irf*theta_<zenith>_az_<azimuth>_*.fits.gz format
+    method: 'linear’, ‘nearest’, ‘cubic’
+        interpolation method
+
+    Returns
+    -------
+    irfs: list
+        list of interpolated IRFs: effective area, energy dispersion, ghcuts
+    """
+    filepaths=glob.glob(input_irf_dir+"/irf*theta_*_az_*_*.fits.gz")
+    re_float="([-+]?(?:\d*\.\d+|\d+))"
+    regex="irf_\S+_theta_"+re_float+"_az_"+re_float+"_\S+.fits.gz"
+
+    aeff_ext_name="EFFECTIVE AREA"
+    edisp_ext_name="ENERGY DISPERSION"
+    points=[]
+    hdus_irfs=[]
+    aeff_all=[]
+    edisp_all=[]
+    ghcuts_low_last=None
+    ghcuts_high_last=None
+    ghcuts_center_last=None
+    ghcuts_all=[]
+    extra_headers_list=['TELESCOP', 'INSTRUME', 'FOVALIGN', 'QUAL_CUT', 'IRF_TYPE', 'DL2_WEIG', 'GH_CUT', 'GH_EFF', 'RAD_MAX', 'TH_EFF']
+    for file in filepaths:
+        name=os.path.basename(file)
+        logger.info("loading file: "+name)
+        if (re.fullmatch(regex, name)):
+            irf_theta, irf_az=re.findall(regex, name)[0]
+            coszd=np.cos(np.radians(float(irf_theta)))
+            irf_az=np.radians(float(irf_az))
+            points.append([coszd, irf_az])
+            hdus_irf=fits.open(file)
+            hdus_irfs.append(hdus_irf)
+
+            aeff=hdus_irf["EFFECTIVE AREA"]            
+            aeff_all.append(aeff.data['EFFAREA'][0])
+            edisp=hdus_irf["ENERGY DISPERSION"]
+            edisp_all.append(edisp.data['MATRIX'][0])
+            ghcuts=hdus_irf["GH_CUTS"] 
+            ghcuts_low=u.Quantity(ghcuts.data['low'], unit=ghcuts.columns['low'].unit)
+            ghcuts_high=u.Quantity(ghcuts.data['high'], unit=ghcuts.columns['high'].unit)
+            ghcuts_center=u.Quantity(ghcuts.data['center'], unit=ghcuts.columns['center'].unit)
+            if ghcuts_low_last is not None:
+                if (ghcuts_low_last!=ghcuts_low).any() or (ghcuts_high_last!=ghcuts_high).any() or (ghcuts_center_last!=ghcuts_center).any():
+                    raise ValueError('Non matching bins in GH_CUTS')                    
+
+            ghcuts_all.append(ghcuts.data['cut'])
+
+            ghcuts_low_last=ghcuts_low
+            ghcuts_high_last=ghcuts_high
+            ghcuts_center_last=ghcuts_center
+        else:
+            logger.warning("skipping "+ name)
+
+    # fix me! only the last file is checked, no consistency checks
+    extra_headers={key: aeff.header[key] for key in list(set(extra_headers_list) & set(aeff.header.keys()))}    
+
+    points=np.array(points)
+    aeff_energ_bins=read_fits_bins_lo_hi(hdus_irfs, aeff_ext_name, "ENERG")
+    aeff_theta_bins=read_fits_bins_lo_hi(hdus_irfs, aeff_ext_name, "THETA")
+
+    edisp_energ_bins=read_fits_bins_lo_hi(hdus_irfs, edisp_ext_name, "ENERG")
+    edisp_migra_bins=read_fits_bins_lo_hi(hdus_irfs, edisp_ext_name, "MIGRA")
+    edisp_theta_bins=read_fits_bins_lo_hi(hdus_irfs, edisp_ext_name, "THETA")
+
+
+    aeff_all= u.Quantity(np.array(aeff_all), hdus_irf[aeff_ext_name].columns["EFFAREA"].unit, copy=False)
+    edisp_all= u.Quantity(np.array(edisp_all), hdus_irf[edisp_ext_name].columns["MATRIX"].unit, copy=False)
+    ghcuts_all=np.array(ghcuts_all)
+
+    # now read the file and check for what to interpolate
+    data=pd.read_hdf(input_file_dl2, 'events/parameters')
+    data_coszd = np.mean(np.sin(data['pointing_alt']))
+    data_az = np.mean(data['pointing_az'])
+
+    target=np.array([data_coszd, data_az])
+    logger.info("interpolate for: "+str(target))
+    aeff_interp=interp.interpolate_effective_area_per_energy_and_fov (aeff_all, points, target, method=method)
+
+    # here we need to swap axes because the function expects shape: (n_grid_points, n_energy_bins, n_migration_bins, n_fov_offset_bins)
+    edisp_interp=interp.interpolate_energy_dispersion(np.swapaxes(edisp_all, 1, 3), points, target, method=method)
+
+    # now create the HDUS
+    # to have the same format we need to loose one dimention, this might need to be rewised for files with many bins in offset angle!
+    aeff_hdu=create_aeff2d_hdu(aeff_interp[...,0], aeff_energ_bins[0], aeff_theta_bins[0], extname=aeff_ext_name, **extra_headers)
+    edisp_hdu=create_energy_dispersion_hdu(edisp_interp[0,...], edisp_energ_bins[0], edisp_migra_bins[0], edisp_theta_bins[0], extname=edisp_ext_name)
+
+    ghcuts_interp = griddata(points, ghcuts_all, target, method=method)
+
+    ghcuts_table=QTable()
+    ghcuts_table['low']=ghcuts_low
+    ghcuts_table['high']=ghcuts_high
+    ghcuts_table['center']=ghcuts_center
+    ghcuts_table['cut']=ghcuts_interp[0]
+    ghcuts_hdu = fits.BinTableHDU(ghcuts_table, header=ghcuts.header, name="GH_CUTS")
+    return [aeff_hdu, edisp_hdu, ghcuts_hdu]
+
+def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config, hdus_irfs=None):
     """
     Creates a DL3 data file with input DL2 data and IRF files.
 
@@ -369,6 +524,8 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
         Path to a directory where to save an output DL3 data file
     config: dict
         Configuration for the LST-1 + MAGIC analysis
+    hdus_irfs: list
+        List of BinTableHDU with IRFs
     """
 
     config_dl3 = config['dl2_to_dl3']
@@ -378,12 +535,17 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
 
     hdus = fits.HDUList([fits.PrimaryHDU(), ])
 
-    # Load the input IRF file:
-    logger.info('\nLoading the input IRF file:')
-    logger.info(input_file_irf)
+    if hdus_irfs is None:
+        # Load the input IRF file:
+        logger.info('\nLoading the input IRF file:')
+        logger.info(input_file_irf)
 
-    hdus_irf = fits.open(input_file_irf)
-    header = hdus_irf[1].header
+        hdus_irf = fits.open(input_file_irf)
+        header = hdus_irf[1].header
+        hdus += hdus_irf[1:]
+    else:
+        header = hdus_irfs[0].header
+        hdus += hdus_irfs
 
     quality_cuts = header['QUAL_CUT']
     irf_type = header['IRF_TYPE']
@@ -393,7 +555,6 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
     logger.info(f'IRF type: {irf_type}')
     logger.info(f'DL2 weight: {dl2_weight}')
 
-    hdus += hdus_irf[1:]
 
     # Load the input DL2 data file:
     logger.info('\nLoading the input DL2 data file:')
@@ -411,7 +572,10 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
     else:
         logger.info('\nApplying dynamic gammaness cuts...')
 
-        cut_table_gh = QTable.read(input_file_irf, 'GH_CUTS')
+        if hdus_irfs is None:
+            cut_table_gh = QTable.read(input_file_irf, 'GH_CUTS')
+        else:            
+            cut_table_gh = QTable(hdus_irfs[-1].data, names=hdus_irfs[-1].columns.names, units=hdus_irfs[-1].columns.units)
 
         mask_gh_gamma = evaluate_binned_cut(
             values=event_table['gammaness'],
@@ -422,10 +586,13 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
 
         event_table = event_table[mask_gh_gamma]
 
-    # Create a event list HDU:
+    # Create an event list HDU:
     logger.info('\nCreating an event list HDU...')
 
-    event_list, event_header = create_event_list(event_table, deadc, **config_dl3)
+    source_name=config_dl3['source_name']
+    source_ra=config_dl3['source_ra']
+    source_dec=config_dl3['source_dec']
+    event_list, event_header = create_event_list(event_table, deadc, source_name, source_ra, source_dec)
 
     hdu_event = fits.BinTableHDU(event_list, header=event_header, name='EVENTS')
     hdus.append(hdu_event)
@@ -476,8 +643,13 @@ def main():
     )
 
     parser.add_argument(
-        '--input-file-irf', '-i', dest='input_file_irf', type=str, required=True,
-        help='Path to an input IRF file.',
+        '--input-file-irf', '-i', dest='input_file_irf', type=str, required=False, default=None,
+        help='Path to an input IRF file (single).',
+    )
+
+    parser.add_argument(
+        '--input-dir-irf', dest='input_dir_irf', type=str, required=False,
+        help='Path to an input IRF directory (interpolation will be applied).',
     )
 
     parser.add_argument(
@@ -491,12 +663,25 @@ def main():
     )
 
     args = parser.parse_args()
+    if (((args.input_file_irf is None) and (args.input_dir_irf is None)) or 
+        ((args.input_file_irf is not None) and (args.input_dir_irf is not None))):
+        logger.error ("Please provide either input-file-irf or input-dir-irf")
+        raise RuntimeError('Wrong IRF paths.')
 
     with open(args.config_file, 'rb') as f:
         config = yaml.safe_load(f)
 
+    config_dl3 = config['dl2_to_dl3']
+
+    hdus_irfs = None
+    if args.input_dir_irf is not None:
+        interpolation_method = 'linear'
+        if 'interpolation_method' in config_dl3:
+            interpolation_method = config_dl3['interpolation_method']
+            hdus_irfs=interpolate_irf(args.input_file_dl2, args.input_dir_irf, method=interpolation_method)
+
     # Process the input data:
-    dl2_to_dl3(args.input_file_dl2, args.input_file_irf, args.output_dir, config)
+    dl2_to_dl3(args.input_file_dl2, args.input_file_irf, args.output_dir, config, hdus_irfs)
 
     logger.info('\nDone.')
 
@@ -505,4 +690,4 @@ def main():
 
 
 if __name__ == '__main__':
-    main()
+    main()