Add reading of calibration files in fits format

src/ctapipe_io_lst/calibration.py

@@ -4,6 +4,7 @@
 import astropy.units as u
 from numba import njit
 import tables
+from astropy.io import fits
 
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
@@ -14,6 +15,7 @@
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
 from ctapipe.containers import FlatFieldContainer, MonitoringCameraContainer, MonitoringContainer, PedestalContainer, PixelStatusContainer, WaveformCalibrationContainer
 from ctapipe.io import HDF5TableReader, read_table
+from astropy.table import QTable
 from traitlets import Enum
 
 from .compat import CTAPIPE_GE_0_21
@@ -34,6 +36,12 @@
     'LSTR0Corrections',
 ]
 
+def to_native(data):
+    """Convert a numpy array to native byteorder."""
+    if not data.dtype.isnative:
+        data = data.byteswap()
+        data = data.view(data.dtype.newbyteorder("="))
+    return data
 
 def get_first_capacitors_for_pixels(first_capacitor_id, expected_pixel_id=None):
     '''
@@ -340,6 +348,7 @@
             # time shift from flat fielding
             if self.mon_data is not None and self.add_calibration_timeshift:
                 time_corr = self.mon_data.tel[tel_id].calibration.time_correction
+
                 # time_shift is subtracted in ctapipe,
                 # but time_correction should be added
                 if CTAPIPE_GE_0_21 or r1.selected_gain_channel is None:
@@ -356,35 +365,74 @@
         Read the correction from hdf5 calibration file
         """
 
-        with tables.open_file(path) as f:
-            tel_groups = [
-                key for key in f.root._v_children.keys()
-                if key.startswith('tel_')
-            ]
-
         mon = MonitoringContainer()
 
-        for base in tel_groups:
-            with HDF5TableReader(path) as h5_table:
-                # read the calibration data
-                tel_id = int(base[4:])
-                mon.tel[tel_id] = MonitoringCameraContainer(
-                    calibration=next(h5_table.read(f'/{base}/calibration', WaveformCalibrationContainer)),
-                    pedestal=next(h5_table.read(f'/{base}/pedestal', PedestalContainer)),
-                    flatfield=next(h5_table.read(f'/{base}/flatfield', FlatFieldContainer)),
-                    pixel_status=next(h5_table.read(f"/{base}/pixel_status", PixelStatusContainer)),
-                )
+        if path.name.endswith(".h5"):
+            # get keys in file
+            with tables.open_file(path) as f:
+                tel_groups = [
+                    key for key in f.root._v_children.keys()
+                    if key.startswith('tel_')
+                ] 
+
+            for base in tel_groups:
+                with HDF5TableReader(path) as h5_table:
+                    # read the calibration data
+                    tel_id = int(base[4:])
+                    mon.tel[tel_id] = MonitoringCameraContainer(
+                        calibration=next(h5_table.read(f'/{base}/calibration', WaveformCalibrationContainer)),
+                        pedestal=next(h5_table.read(f'/{base}/pedestal', PedestalContainer)),
+                        flatfield=next(h5_table.read(f'/{base}/flatfield', FlatFieldContainer)),
+                        pixel_status=next(h5_table.read(f"/{base}/pixel_status", PixelStatusContainer)),
+                    )
+
+        elif path.name.endswith(".fits") or path.name.endswith(".fits.gz"):
+
+            CALIB_CONTAINERS = {
+                "calibration": WaveformCalibrationContainer,
+                "flatfield": FlatFieldContainer,
+                "pedestal": PedestalContainer,
+                "pixel_status": PixelStatusContainer,
+            } 
+
+            mon_data = MonitoringCameraContainer()
+
+            with fits.open(path) as f:
+                tel_id = f['PRIMARY'].header['TEL_ID']
+                for key in CALIB_CONTAINERS.keys():
+                    table = QTable.read(f, hdu=key)
+                    row = table[0]
+
+                    for col in row.keys():  
+                        mon_data[key][col] = row[col]
+
+            mon.tel[tel_id] = mon_data
+
+        else:
+            raise ValueError("Not supported file format : %s", path)
+
+
         return mon
+
+
 
     @staticmethod
     def load_drs4_time_calibration_file(path):
         """
         Function to load calibration file.
         """
-        with tables.open_file(path, 'r') as f:
-            fan = f.root.fan[:]
-            fbn = f.root.fbn[:]
-
+        if path.name.endswith(".h5"):
+            with tables.open_file(path, 'r') as f:
+                fan = f.root.fan[:]
+                fbn = f.root.fbn[:]
+
+        elif path.name.endswith(".fits") or path.name.endswith(".fits.gz"):
+            with fits.open(path) as f:
+                fan = to_native(f["fan"].data)
+                fbn = to_native(f["fbn"].data)
+        else:
+            raise ValueError("Not supported file format : %s", path)
+
         return fan, fbn
 
     def load_drs4_time_calibration_file_for_tel(self, tel_id):
@@ -396,7 +444,7 @@
         """
         Return pulse time after time correction.
         """
-
+ 
         if self.drs4_time_calibration_path.tel[tel_id] is None:
             if CTAPIPE_GE_0_21 or selected_gain_channel is None:
                 return np.zeros((N_GAINS, N_PIXELS))
@@ -407,7 +455,7 @@
         if tel_id not in self.fan:
             self.load_drs4_time_calibration_file_for_tel(tel_id)
 
-        if CTAPIPE_GE_0_21 or selected_gain_channel is None:
+        if CTAPIPE_GE_0_21 or selected_gain_channel is None: 
             return calc_drs4_time_correction_both_gains(
                 first_capacitors,
                 self.fan[tel_id],
@@ -439,14 +487,24 @@
                 " but no file provided for telescope"
             )
 
-        table = read_table(path, f'/r1/monitoring/drs4_baseline/tel_{tel_id:03d}')
-
         pedestal_data = np.empty(
-            (N_GAINS, N_PIXELS_MODULE * N_MODULES, N_CAPACITORS_PIXEL + N_SAMPLES),
-            dtype=np.float32
-        )
-        pedestal_data[:, :, :N_CAPACITORS_PIXEL] = table[0]['baseline_mean']
-        pedestal_data[:, :, N_CAPACITORS_PIXEL:] = pedestal_data[:, :, :N_SAMPLES]
+                (N_GAINS, N_PIXELS_MODULE * N_MODULES, N_CAPACITORS_PIXEL + N_SAMPLES),
+                dtype=np.float32
+            )
+
+        if path.name.endswith(".h5"):
+            table = read_table(path, f'/r1/monitoring/drs4_baseline/tel_{tel_id:03d}')
+
+            pedestal_data[:, :, :N_CAPACITORS_PIXEL] = table[0]['baseline_mean']
+            pedestal_data[:, :, N_CAPACITORS_PIXEL:] = pedestal_data[:, :, :N_SAMPLES]
+
+        elif path.name.endswith(".fits") or path.name.endswith(".fits.gz"):   
+            with fits.open(path) as f:
+                pedestal_data[:, :, :N_CAPACITORS_PIXEL] = to_native(f["baseline_mean"].data)
+                pedestal_data[:, :, N_CAPACITORS_PIXEL:] = pedestal_data[:, :, :N_SAMPLES]
+
+        else:
+            raise ValueError("Not supported file format : %s", path)
 
         return pedestal_data
 
@@ -457,9 +515,14 @@
                 "DRS4 spike correction requested"
                 " but no pedestal file provided for telescope"
             )
+        if path.name.endswith(".h5"):
+            table = read_table(path, f'/r1/monitoring/drs4_baseline/tel_{tel_id:03d}')
+            spike_height = np.array(table[0]['spike_height'])
+
+        elif path.name.endswith(".fits") or path.name.endswith(".fits.gz"):  
+            with fits.open(path) as f:
+                spike_height = to_native(f["spike_height"].data)
 
-        table = read_table(path, f'/r1/monitoring/drs4_baseline/tel_{tel_id:03d}')
-        spike_height = np.array(table[0]['spike_height'])
         return spike_height
 
     def subtract_pedestal(self, event, tel_id):
@@ -471,9 +534,10 @@
         event : `ctapipe` event-container
         tel_id : id of the telescope
         """
+
         pedestal = self._get_drs4_pedestal_data(
             self.drs4_pedestal_path.tel[tel_id],
-            tel_id,
+            tel_id
         )
 
         if event.r1.tel[tel_id].selected_gain_channel is None:
@@ -580,7 +644,7 @@
         """
         run_id = event.lst.tel[tel_id].svc.configuration_id
         spike_height = self._get_spike_heights(self.drs4_pedestal_path.tel[tel_id], tel_id)
-
+        
         r1 = event.r1.tel[tel_id]
         if r1.selected_gain_channel is None:
             subtract_spikes(
@@ -1123,7 +1187,6 @@
     first_capacitors, fan, fbn
 ):
     time = np.zeros((N_GAINS, N_PIXELS))
-
     for gain in range(N_GAINS):
         for pixel in range(N_PIXELS):
             first_capacitor = first_capacitors[gain, pixel]

src/ctapipe_io_lst/tests/test_calib.py

@@ -22,10 +22,28 @@
 calib_version = "ctapipe-v0.17"
 calib_path = test_data / 'real/monitoring/PixelCalibration/Cat-A/'
 test_calib_path = calib_path / f'calibration/20200218/{calib_version}/calibration_filters_52.Run02006.0000.h5'
+test_calib_path_fits = calib_path / f'calibration/20200218/{calib_version}/calibration_filters_52.Run02006.0000.fits.gz'
 test_drs4_pedestal_path = calib_path / f'drs4_baseline/20200218/{calib_version}/drs4_pedestal.Run02005.0000.h5'
+test_drs4_pedestal_path_fits = calib_path / f'drs4_baseline/20200218/{calib_version}/drs4_pedestal.Run02005.0000.fits.gz'
 test_time_calib_path = calib_path / f'drs4_time_sampling_from_FF/20191124/{calib_version}/time_calibration.Run01625.0000.h5'
+test_time_calib_path_fits = calib_path / f'drs4_time_sampling_from_FF/20191124/{calib_version}/time_calibration.Run01625.0000.fits.gz'
 
 
+@pytest.mark.parametrize(
+    "array",
+    [
+        np.array([1, 2], dtype=">i4"),
+        np.array([1, 2], dtype="<i4"),
+        np.array([1, 2], dtype="i4"),
+    ]
+)
+def test_to_native(array):
+    from ctapipe_io_lst.calibration import to_native
+
+    converted = to_native(array)
+    assert converted.dtype.isnative
+    np.testing.assert_array_equal(converted, array)
+
 def test_get_first_capacitor():
     from ctapipe_io_lst import LSTEventSource
     from ctapipe_io_lst.calibration import (
@@ -56,28 +74,50 @@ def test_get_first_capacitor():
     assert np.all(first_caps == expected)
 
 
-def test_read_calib_file():
+
+@pytest.mark.parametrize(
+    "path",
+    [
+        pytest.param(test_calib_path, id="hdf5"),
+        pytest.param(test_calib_path_fits, id="fits"),
+    ]
+)
+def test_read_calib_file(path):
     from ctapipe_io_lst.calibration import LSTR0Corrections
 
-    mon = LSTR0Corrections._read_calibration_file(test_calib_path)
+    mon = LSTR0Corrections._read_calibration_file(path)
     # only one telescope in that file
     assert mon.tel.keys() == {1, }
 
 
-def test_read_drs4_pedestal_file():
+@pytest.mark.parametrize(
+    "path",
+    [
+        pytest.param(test_drs4_pedestal_path, id="hdf5"),
+        pytest.param(test_drs4_pedestal_path_fits, id="fits"),
+    ]
+)
+def test_read_drs4_pedestal_file(path):
     from ctapipe_io_lst.calibration import LSTR0Corrections, N_CAPACITORS_PIXEL, N_SAMPLES
 
-    pedestal = LSTR0Corrections._get_drs4_pedestal_data(test_drs4_pedestal_path, tel_id=1)
+    pedestal = LSTR0Corrections._get_drs4_pedestal_data(path, tel_id=1)
 
     assert pedestal.shape[-1] == N_CAPACITORS_PIXEL + N_SAMPLES
     # check circular boundary
     assert np.all(pedestal[..., :N_SAMPLES] == pedestal[..., N_CAPACITORS_PIXEL:])
 
+@pytest.mark.parametrize(
+    "path",
+    [
+        pytest.param(test_time_calib_path, id="hdf5"),
+        pytest.param(test_time_calib_path_fits, id="fits"),
+    ]
+)
 
-def test_read_drs_time_calibration_file():
+def test_read_drs_time_calibration_file(path):
     from ctapipe_io_lst.calibration import LSTR0Corrections, N_GAINS, N_PIXELS
 
-    fan, fbn = LSTR0Corrections.load_drs4_time_calibration_file(test_time_calib_path)
+    fan, fbn = LSTR0Corrections.load_drs4_time_calibration_file(path)
 
     assert fan.shape == fbn.shape
     assert fan.shape[0] == N_GAINS
@@ -139,20 +179,33 @@ def test_source_with_calibration():
         for event in source:
             assert event.r1.tel[1].waveform is not None
 
-
+fits_calibration = {
+    "drs4_pedestal_path": test_drs4_pedestal_path,
+    "drs4_time_calibration_path": test_time_calib_path,
+    "calibration_path": test_calib_path,
+}
+hdf5_calibration = {
+    "drs4_pedestal_path": test_drs4_pedestal_path_fits,
+    "drs4_time_calibration_path": test_time_calib_path_fits,
+    "calibration_path": test_calib_path_fits,
+}
 @pytest.mark.parametrize("trigger_information", [True, False])
-def test_source_with_all(trigger_information):
+@pytest.mark.parametrize(
+    "calib_config",
+    (
+        pytest.param(hdf5_calibration, id="hdf5"),
+        pytest.param(fits_calibration, id="fits"),
+    )
+)
+
+def test_calibration(trigger_information,calib_config):
     from ctapipe_io_lst import LSTEventSource
 
     config = Config({
         'LSTEventSource': {
             'pointing_information': False,
-            'trigger_information': trigger_information,
-            'LSTR0Corrections': {
-                'drs4_pedestal_path': test_drs4_pedestal_path,
-                'drs4_time_calibration_path': test_time_calib_path,
-                'calibration_path': test_calib_path,
-            },
+            'trigger_information': trigger_information, 
+            'LSTR0Corrections': calib_config
         },
     })
 
@@ -167,7 +220,6 @@ def test_source_with_all(trigger_information):
             assert event.r1.tel[1].waveform is not None
             assert np.any(event.calibration.tel[1].dl1.time_shift != 0)
 
-
 def test_missing_module():
     from ctapipe_io_lst import LSTEventSource
     from ctapipe_io_lst.constants import N_PIXELS_MODULE, N_SAMPLES