Merge pull request #121 from cta-observatory/spike_subtraction

Implement and use spike subtraction by default, fixes #122

ctapipe_io_lst/calibration.py

@@ -1,7 +1,6 @@
 from functools import lru_cache
 
 import numpy as np
-from astropy.io import fits
 import astropy.units as u
 from numba import njit
 import tables
@@ -14,8 +13,9 @@
 
 from ctapipe.calib.camera.gainselection import ThresholdGainSelector
 from ctapipe.containers import MonitoringContainer
-from ctapipe.io import HDF5TableReader
+from ctapipe.io import HDF5TableReader, read_table
 from .containers import LSTArrayEventContainer
+from traitlets import Enum
 
 
 from .constants import (
@@ -114,6 +114,7 @@ class LSTR0Corrections(TelescopeComponent):
         help=(
             'Path to the LST pedestal file'
             ', required when `apply_drs4_pedestal_correction=True`'
+            ' or when using spike subtraction'
         ),
     ).tag(config=True)
 
@@ -175,6 +176,12 @@ class LSTR0Corrections(TelescopeComponent):
         help='Threshold for the ThresholdGainSelector.'
     ).tag(config=True)
 
+    spike_correction_method = Enum(
+        values=['subtraction', 'interpolation'],
+        default_value='subtraction',
+        help='Wheter to use spike subtraction (default) or interpolation',
+    ).tag(config=True)
+
     def __init__(self, subarray, config=None, parent=None, **kwargs):
         """
         The R0 calibrator for LST data.
@@ -234,7 +241,11 @@ def apply_drs4_corrections(self, event: LSTArrayEventContainer):
                 self.time_lapse_corr(event, tel_id)
 
             if self.apply_spike_correction:
-                self.interpolate_spikes(event, tel_id)
+                if self.spike_correction_method == 'subtraction':
+                    self.subtract_spikes(event, tel_id)
+                else:
+                    self.interpolate_spikes(event, tel_id)
+
 
             # remove samples at beginning / end of waveform
             start = self.r1_sample_start.tel[tel_id]
@@ -404,7 +415,7 @@ def get_drs4_time_correction(self, tel_id, first_capacitors, selected_gain_chann
 
     @staticmethod
     @lru_cache(maxsize=4)
-    def _get_drs4_pedestal_data(path):
+    def _get_drs4_pedestal_data(path, tel_id):
         """
         Function to load pedestal file.
 
@@ -420,17 +431,29 @@ def _get_drs4_pedestal_data(path):
                 " 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.int16
+            dtype=np.float32
         )
-        with fits.open(path) as f:
-            pedestal_data[:, :, :N_CAPACITORS_PIXEL] = f[1].data
-
+        pedestal_data[:, :, :N_CAPACITORS_PIXEL] = table[0]['baseline_mean']
         pedestal_data[:, :, N_CAPACITORS_PIXEL:] = pedestal_data[:, :, :N_SAMPLES]
 
         return pedestal_data
 
+    @lru_cache(maxsize=4)
+    def _get_spike_heights(self, path, tel_id):
+        if path is None:
+            raise ValueError(
+                "DRS4 spike correction requested"
+                " but no pedestal file provided for telescope"
+            )
+
+        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):
         """
         Subtract cell offset using pedestal file.
@@ -441,7 +464,8 @@ def subtract_pedestal(self, event, tel_id):
         tel_id : id of the telescope
         """
         pedestal = self._get_drs4_pedestal_data(
-            self.drs4_pedestal_path.tel[tel_id]
+            self.drs4_pedestal_path.tel[tel_id],
+            tel_id,
         )
 
         if event.r1.tel[tel_id].selected_gain_channel is None:
@@ -509,12 +533,9 @@ def time_lapse_corr(self, event, tel_id):
 
     def interpolate_spikes(self, event, tel_id):
         """
-        Interpolates spike A & B.
-        Fill the R1 container.
-        Parameters
-        ----------
-        event : `ctapipe` event-container
-        tel_id : id of the telescope
+        Interpolate spikes at known positions from their neighboring values
+
+        Mutates the R1 waveform.
         """
         run_id = event.lst.tel[tel_id].svc.configuration_id
 
@@ -535,6 +556,34 @@ def interpolate_spikes(self, event, tel_id):
                 run_id=run_id,
             )
 
+    def subtract_spikes(self, event, tel_id):
+        """
+        Subtract mean spike height from known spike positions
+
+        Mutates the R1 waveform.
+        """
+        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(
+                waveform=r1.waveform,
+                first_capacitors=self.first_cap[tel_id],
+                previous_first_capacitors=self.first_cap_old[tel_id],
+                run_id=run_id,
+                spike_height=spike_height
+            )
+        else:
+            subtract_spikes_gain_selected(
+                waveform=r1.waveform,
+                first_capacitors=self.first_cap[tel_id],
+                previous_first_capacitors=self.first_cap_old[tel_id],
+                selected_gain_channel=r1.selected_gain_channel,
+                run_id=run_id,
+                spike_height=spike_height
+            )
+
 
 def convert_to_pe(waveform, calibration, selected_gain_channel):
     if selected_gain_channel is None:
@@ -544,6 +593,7 @@ def convert_to_pe(waveform, calibration, selected_gain_channel):
         waveform -= calibration.pedestal_per_sample[selected_gain_channel, PIXEL_INDEX, np.newaxis]
         waveform *= calibration.dc_to_pe[selected_gain_channel, PIXEL_INDEX, np.newaxis]
 
+
 @njit(cache=True)
 def interpolate_spike_A(waveform, position):
     """
@@ -720,6 +770,12 @@ def interpolate_spikes_gain_selected(waveform, first_capacitors, previous_first_
         Value of first capacitor from previous event
         stored in a numpy array of shape
         (N_GAINS, N_PIXELS).
+    selected_gain_channel: ndarray
+        ndarray of shape (N_PIXELS, ) containing the selected gain channel
+        for each pixel
+    run_id: int
+        Run id of the run, used to determine if code for new firmware
+        or old firmware has to be used
     """
 
     for pixel in range(N_PIXELS):
@@ -738,6 +794,117 @@ def interpolate_spikes_gain_selected(waveform, first_capacitors, previous_first_
         )
 
 
+@njit(cache=True)
+def subtract_spikes_at_positions(waveform, positions, spike_height):
+    '''Subtract the spikes at given positions in waveform'''
+    for spike_position in positions:
+        for i in range(3):
+            sample = spike_position + i
+            if 0 <= sample < N_SAMPLES:
+                waveform[sample] -= spike_height[i]
+
+
+@njit(cache=True)
+def subtract_spikes(
+    waveform,
+    first_capacitors,
+    previous_first_capacitors,
+    run_id,
+    spike_height,
+):
+    """
+    Subtract mean spike heights for spike type A.
+
+    Modifies waveform in place.
+
+    Parameters
+    ----------
+    waveform : ndarray
+        Waveform stored in a numpy array of shape
+        (N_GAINS, N_PIXELS, N_SAMPLES).
+    first_capacitors : ndarray
+        Value of first capacitor stored in a numpy array of shape
+        (N_GAINS, N_PIXELS).
+    previous_first_capacitors : ndarray
+        Value of first capacitor from previous event
+        stored in a numpy array of shape
+        (N_GAINS, N_PIXELS).
+    run_id: int
+        Run id of the run, used to determine if code for new firmware
+        or old firmware has to be used
+    spike_height: ndarray
+        ndarry of shape (N_GAINS, N_PIXELS, 3) of the three spike_heights
+    """
+    for gain in range(N_GAINS):
+        for pixel in range(N_PIXELS):
+            current_fc = first_capacitors[gain, pixel]
+            last_fc = previous_first_capacitors[gain, pixel]
+
+            if run_id > LAST_RUN_WITH_OLD_FIRMWARE:
+                positions = get_spike_A_positions(current_fc, last_fc)
+            else:
+                positions = get_spike_A_positions_old_firmware(current_fc, last_fc)
+
+            subtract_spikes_at_positions(
+                waveform=waveform[gain, pixel],
+                positions=positions,
+                spike_height=spike_height[gain, pixel],
+            )
+
+
+@njit(cache=True)
+def subtract_spikes_gain_selected(
+    waveform,
+    first_capacitors,
+    previous_first_capacitors,
+    selected_gain_channel,
+    run_id,
+    spike_height,
+):
+    """
+    Subtract mean spike heights for spike type A for gain selected input data
+
+    Modifies waveform in place.
+
+    Parameters
+    ----------
+    waveform : ndarray
+        Waveform stored in a numpy array of shape
+        (N_GAINS, N_PIXELS, N_SAMPLES).
+    first_capacitors : ndarray
+        Value of first capacitor stored in a numpy array of shape
+        (N_GAINS, N_PIXELS).
+    previous_first_capacitors : ndarray
+        Value of first capacitor from previous event
+        stored in a numpy array of shape
+        (N_GAINS, N_PIXELS).
+    selected_gain_channel: ndarray
+        ndarray of shape (N_PIXELS, ) containing the selected gain channel
+        for each pixel
+    run_id: int
+        Run id of the run, used to determine if code for new firmware
+        or old firmware has to be used
+    spike_height: ndarray
+        ndarry of shape (N_GAINS, N_PIXELS, 3) of the three spike_heights
+    """
+
+    for pixel in range(N_PIXELS):
+        gain = selected_gain_channel[pixel]
+        current_fc = first_capacitors[gain, pixel]
+        last_fc = previous_first_capacitors[gain, pixel]
+
+        if run_id > LAST_RUN_WITH_OLD_FIRMWARE:
+            positions = get_spike_A_positions(current_fc, last_fc)
+        else:
+            positions = get_spike_A_positions_old_firmware(current_fc, last_fc)
+
+        subtract_spikes_at_positions(
+            waveform=waveform[pixel],
+            positions=positions,
+            spike_height=spike_height[gain, pixel],
+        )
+
+
 @njit(cache=True)
 def subtract_pedestal(
     waveform,

ctapipe_io_lst/tests/test_calib.py

@@ -14,9 +14,9 @@
 test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data'))
 test_r0_path = test_data / 'real/R0/20200218/LST-1.1.Run02008.0000_first50.fits.fz'
 test_r0_calib_path = test_data / 'real/R0/20200218/LST-1.1.Run02006.0004.fits.fz'
-test_calib_path = test_data / 'real/calibration/20200218/v05/calibration.Run2006.0000.hdf5'
-test_drs4_pedestal_path = test_data / 'real/calibration/20200218/v05/drs4_pedestal.Run2005.0000.fits'
-test_time_calib_path = test_data / 'real/calibration/20200218/v05/time_calibration.Run2006.0000.hdf5'
+test_calib_path = test_data / 'real/monitoring/PixelCalibration/LevelA/calibration/20200218/v0.7.6.dev568+gdecfb58c/calibration_filters_52.Run02006.0000.h5'
+test_drs4_pedestal_path = test_data / 'real/monitoring/PixelCalibration/LevelA/drs4_baseline/20200218/v0.7.6.dev568+gdecfb58c/drs4_pedestal.Run02005.0000.h5'
+test_time_calib_path = test_data / 'real/monitoring/PixelCalibration/LevelA/drs4_time_sampling_from_FF/20191124/v0.7.6.dev568+gdecfb58c/time_calibration.Run01625.0000.h5'
 test_missing_module_path = test_data / 'real/R0/20210215/LST-1.1.Run03669.0000_first50.fits.fz'
 test_r0_gainselected_path = test_data / 'real/R0/20200218/LST-1.1.Run02008.0000_first50_gainselected.fits.fz'
 
@@ -63,7 +63,7 @@ def test_read_calib_file():
 def test_read_drs4_pedestal_file():
     from ctapipe_io_lst.calibration import LSTR0Corrections, N_CAPACITORS_PIXEL, N_SAMPLES
 
-    pedestal = LSTR0Corrections._get_drs4_pedestal_data(test_drs4_pedestal_path)
+    pedestal = LSTR0Corrections._get_drs4_pedestal_data(test_drs4_pedestal_path, tel_id=1)
 
     assert pedestal.shape[-1] == N_CAPACITORS_PIXEL + N_SAMPLES
     # check circular boundary

ctapipe_io_lst/tests/test_stage1.py

@@ -11,9 +11,9 @@
 
 test_data = Path(os.getenv('LSTCHAIN_TEST_DATA', 'test_data'))
 test_r0_path = test_data / 'real/R0/20200218/LST-1.1.Run02008.0000_first50.fits.fz'
-test_calib_path = test_data / 'real/calibration/20200218/v05/calibration.Run2006.0000.hdf5'
-test_drs4_pedestal_path = test_data / 'real/calibration/20200218/v05/drs4_pedestal.Run2005.0000.fits'
-test_time_calib_path = test_data / 'real/calibration/20200218/v05/time_calibration.Run2006.0000.hdf5'
+test_calib_path = test_data / 'real/monitoring/PixelCalibration/LevelA/calibration/20200218/v0.7.6.dev568+gdecfb58c/calibration_filters_52.Run02006.0000.h5'
+test_drs4_pedestal_path = test_data / 'real/monitoring/PixelCalibration/LevelA/drs4_baseline/20200218/v0.7.6.dev568+gdecfb58c/drs4_pedestal.Run02005.0000.h5'
+test_time_calib_path = test_data / 'real/monitoring/PixelCalibration/LevelA/drs4_time_sampling_from_FF/20191124/v0.7.6.dev568+gdecfb58c/time_calibration.Run01625.0000.h5'
 test_drive_report = test_data / 'real/monitoring/DrivePositioning/drive_log_20200218.txt'
 test_run_summary = test_data / 'real/monitoring/RunSummary/RunSummary_20200218.ecsv'
 

download_test_data.sh

@@ -21,6 +21,7 @@ wget \
 	-R "*.html*,*.gif" \
 	--no-host-directories --cut-dirs=1 \
 	--no-parent \
+	--level=inf \
 	--user="$TEST_DATA_USER" \
 	--password="$TEST_DATA_PASSWORD" \
 	--no-verbose \