Add per-waveform baseline and RMS statistics to IDEX L1b

imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml

@@ -43,6 +43,16 @@ l1b_target_base: &l1b_target_base
   DEPEND_1: time_low_sample_rate_index
   LABL_PTR_1: time_low_sample_rate_label
 
+waveform_statistics_base: &waveform_statistics_base
+  <<: *l1b_data_base
+  DISPLAY_TYPE: time_series
+  VAR_NOTES: >-
+    Calculated from finite samples in the first 3 microseconds of each record,
+    using start <= time < start + 3 us, where start is the earliest finite sample
+    time on that channel's time axis. Baseline is the arithmetic mean; RMS is
+    sqrt(mean((signal - baseline)^2)), with ddof=0. Missing windows are fill values.
+  VAR_TYPE: support_data
+
 trigger_base: &trigger_base
   <<: *l1b_data_base
   DISPLAY_TYPE: no_plot
@@ -243,6 +253,96 @@ ion_grid:
         Sampled at 3.75 Msps with 12-bit resolution.
         Data is used to quantify dust charge.
 
+tof_high_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the TOF High waveform.
+  FIELDNAM: TOF High Baseline
+  LABLAXIS: TOF HG Baseline
+  UNITS: mA
+  VALIDMAX: 0.076725
+
+tof_high_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the TOF High waveform.
+  FIELDNAM: TOF High RMS
+  LABLAXIS: TOF HG RMS
+  UNITS: mA
+  VALIDMAX: 0.076725
+
+tof_mid_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the TOF Mid waveform.
+  FIELDNAM: TOF Mid Baseline
+  LABLAXIS: TOF MG Baseline
+  UNITS: mA
+  VALIDMAX: 2.99739
+
+tof_mid_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the TOF Mid waveform.
+  FIELDNAM: TOF Mid RMS
+  LABLAXIS: TOF MG RMS
+  UNITS: mA
+  VALIDMAX: 2.99739
+
+tof_low_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the TOF Low waveform.
+  FIELDNAM: TOF Low Baseline
+  LABLAXIS: TOF LG Baseline
+  UNITS: mA
+  VALIDMAX: 137.082
+
+tof_low_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the TOF Low waveform.
+  FIELDNAM: TOF Low RMS
+  LABLAXIS: TOF LG RMS
+  UNITS: mA
+  VALIDMAX: 137.082
+
+target_high_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the Target High waveform.
+  FIELDNAM: Target High Baseline
+  LABLAXIS: Target HG Baseline
+  VALIDMAX: 667.485
+
+target_high_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the Target High waveform.
+  FIELDNAM: Target High RMS
+  LABLAXIS: Target HG RMS
+  VALIDMAX: 667.485
+
+target_low_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the Target Low waveform.
+  FIELDNAM: Target Low Baseline
+  LABLAXIS: Target LG Baseline
+  VALIDMAX: 64701.0
+
+target_low_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the Target Low waveform.
+  FIELDNAM: Target Low RMS
+  LABLAXIS: Target LG RMS
+  VALIDMAX: 64701.0
+
+ion_grid_baseline:
+  <<: *waveform_statistics_base
+  CATDESC: Mean baseline of the Ion Grid waveform.
+  FIELDNAM: Ion Grid Baseline
+  LABLAXIS: Ion Grid Baseline
+  VALIDMAX: 3.05487
+
+ion_grid_rms:
+  <<: *waveform_statistics_base
+  CATDESC: Baseline-subtracted RMS noise of the Ion Grid waveform.
+  FIELDNAM: Ion Grid RMS
+  LABLAXIS: Ion Grid RMS
+  VALIDMAX: 3.05487
+
 detector_voltage:
   <<: *trigger_base
   CATDESC: Voltage reading of Detector on HVPS Board (ADC Channel 0)

imap_processing/idex/idex_l1b.py

@@ -32,7 +32,11 @@
     ConversionFactors,
 )
 from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES
-from imap_processing.idex.idex_utils import get_idex_attrs, setup_dataset
+from imap_processing.idex.idex_utils import (
+    calculate_baseline_and_rms,
+    get_idex_attrs,
+    setup_dataset,
+)
 from imap_processing.spice.geometry import (
     SpiceBody,
     SpiceFrame,
@@ -292,6 +296,7 @@ def idex_l1b_science(l1a_dataset: xr.Dataset) -> xr.Dataset:
         conversion_table_path=var_information_path,
         packet_name="IDEX_SCI",
     )
+    l1b_dataset.update(get_waveform_statistics(l1b_dataset, idex_attrs))
 
     logger.info("IDEX L1B science data processing completed.")
     return l1b_dataset
@@ -378,6 +383,44 @@ def convert_waveforms(
     return waveforms_converted
 
 
+def get_waveform_statistics(
+    dataset: xr.Dataset, idex_attrs: ImapCdfAttributes
+) -> dict[str, xr.DataArray]:
+    """
+    Create event-level baseline and RMS fields for all six converted waveforms.
+
+    Parameters
+    ----------
+    dataset : xarray.Dataset
+        L1b waveforms in engineering units and sample times in microseconds.
+    idex_attrs : ImapCdfAttributes
+        CDF attribute manager object.
+
+    Returns
+    -------
+    statistics : dict[str, xarray.DataArray]
+        Mean baseline and baseline-subtracted RMS noise in the first 3 microseconds
+        of each waveform. TOF fields are in mA; target and ion-grid fields are in pC.
+    """
+    statistics = {}
+    for channel in ConversionFactors:
+        waveform = dataset[channel.name]
+        time_name = waveform.dims[-1].removesuffix("_index")
+        baseline, rms = calculate_baseline_and_rms(
+            waveform.data, dataset[time_name].data
+        )
+        for suffix, values in (("baseline", baseline), ("rms", rms)):
+            name = f"{channel.name.lower()}_{suffix}"
+            statistics[name] = xr.DataArray(
+                values,
+                dims="epoch",
+                coords={"epoch": dataset.epoch},
+                name=name,
+                attrs=idex_attrs.get_variable_attributes(name),
+            )
+    return statistics
+
+
 def get_trigger_mode_and_level(
     l1a_dataset: xr.Dataset,
     idex_attrs: ImapCdfAttributes,

imap_processing/idex/idex_l2a.py

@@ -17,7 +17,6 @@
 """
 
 import logging
-from enum import IntEnum
 
 import numpy as np
 import pandas as pd
@@ -32,27 +31,15 @@
 from imap_processing.idex import idex_constants
 from imap_processing.idex.idex_constants import SPICE_ARRAYS
 from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES
-from imap_processing.idex.idex_utils import get_idex_attrs, setup_dataset
+from imap_processing.idex.idex_utils import (
+    calculate_baseline_and_rms,
+    get_idex_attrs,
+    setup_dataset,
+)
 
 logger = logging.getLogger(__name__)
 
 
-class BaselineNoiseTime(IntEnum):
-    """
-    Time range in microseconds that mark the baseline noise before a Dust impact.
-
-    Attributes
-    ----------
-    START: int
-         Beginning of the baseline noise window.
-    STOP: int
-        End of the baseline noise window.
-    """
-
-    START = -7
-    STOP = -5
-
-
 def load_calibration_files(ancillary_files: dict) -> tuple[NDArray, NDArray]:
     """
     Load calibration files for IDEX L2A processing.
@@ -94,7 +81,8 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset:
     Parameters
     ----------
     l1b_dataset : xarray.Dataset
-        IDEX L1a dataset to process.
+        IDEX L1b dataset to process, including the per-event waveform baseline
+        and RMS statistics.
     ancillary_files : dict
         Ancillary files containing calibration coefficients needed to estimate
         velocity and mass of the dust particles.
@@ -128,7 +116,12 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset:
         dims=("epoch", "time_high_sample_rate_index"),
         attrs=idex_attrs.get_variable_attributes("mass_scale"),
     )
-    snr = calculate_snr(tof_high, hs_time)
+    snr = calculate_snr(
+        tof_high,
+        hs_time,
+        l1b_dataset["tof_high_baseline"],
+        l1b_dataset["tof_high_rms"],
+    )
     # Find peaks for each event. The peaks represent a TOF of an ion.
     # Peaks_2d is a list of variable-length arrays
     peaks_2d = [find_peaks(tof, prominence=0.01)[0] for tof in tof_high]
@@ -194,14 +187,21 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset:
             l2a_dataset[flag_name].attrs = idex_attrs.get_variable_attributes(flag_name)
 
     for waveform in ["Target_Low", "Target_High", "Ion_Grid"]:
+        baseline = l1b_dataset[f"{waveform.lower()}_baseline"]
         # Get the dust mass estimates and fit results
         fit_results = xr.apply_ufunc(
             estimate_dust_mass,
             ls_time,
             l1b_dataset[waveform],
+            False,
+            waveform,
+            baseline,
             input_core_dims=[
                 ["time_low_sample_rate_index"],
                 ["time_low_sample_rate_index"],
+                [],
+                [],
+                [],
             ],
             output_core_dims=[
                 ["target_fit_parameter_index"],
@@ -213,7 +213,6 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset:
             vectorize=True,
             output_dtypes=[np.float64] * 6,
             keep_attrs=True,
-            kwargs={"waveform_name": waveform},
         )
         # Calculate mass and velocity estimates
         if waveform == "Ion_Grid":
@@ -784,7 +783,12 @@ def calculate_kappa(mass_scales: np.ndarray, peaks_2d: list) -> NDArray:
     return kappas
 
 
-def calculate_snr(tof_high: xr.DataArray, hs_time: xr.DataArray) -> NDArray:
+def calculate_snr(
+    tof_high: xr.DataArray,
+    hs_time: xr.DataArray,
+    baseline: xr.DataArray | NDArray | None = None,
+    rms: xr.DataArray | NDArray | None = None,
+) -> NDArray:
     """
     Calculate the signal-to-noise ratio.
 
@@ -794,33 +798,24 @@ def calculate_snr(tof_high: xr.DataArray, hs_time: xr.DataArray) -> NDArray:
         The time of flight array.
     hs_time : xarray.DataArray
         The high sampling time array.
+    baseline : xarray.DataArray or numpy.ndarray or None
+        Per-event L1b mean baseline. Calculated from the first 3 microseconds
+        when baseline or rms is omitted in standalone calls.
+    rms : xarray.DataArray or numpy.ndarray or None
+        Per-event L1b baseline-subtracted RMS noise (ddof=0).
 
     Returns
     -------
     numpy.ndarray
         Signal-to-noise ratio at each event.
     """
-    # Find indices where Time (High Sampling) is between -7 and -5 ns (no signal yet)
-    # To determine the baseline noise
-    baseline_noise = np.where(
-        np.logical_and(
-            hs_time >= BaselineNoiseTime.START, hs_time <= BaselineNoiseTime.STOP
-        ),
-        tof_high,
-        np.nan,
-    )
-    if np.all(np.isnan(baseline_noise)):
-        logger.warning(
-            "Unable to find baseline noise. "
-            f"There is no signal from {BaselineNoiseTime.START} to "
-            f"{BaselineNoiseTime.STOP} ns. Returning np.nan SNR values"
+    if baseline is None or rms is None:
+        baseline, rms = calculate_baseline_and_rms(
+            np.asarray(tof_high), np.asarray(hs_time)
         )
-        return np.full(len(hs_time), fill_value=np.nan)
-    # Get the max signal without baseline noise
-    tof_max = np.max(tof_high, axis=1) - np.nanmean(baseline_noise, axis=1)
-    tof_sigma = np.nanstd(baseline_noise, axis=1, ddof=1)
-    # Return snr ratio
-    return tof_max / tof_sigma
+    tof_max = np.max(tof_high, axis=1) - np.asarray(baseline)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        return np.asarray(tof_max / np.asarray(rms))
 
 
 def analyze_peaks(
@@ -1003,6 +998,7 @@ def estimate_dust_mass(
     target_signal: xr.DataArray,
     remove_noise: bool = False,
     waveform_name: str = "",
+    baseline_mean: float | None = None,
 ) -> tuple[NDArray, float, float, float, NDArray]:
     """
     Filter and fit the target or ion grid signals to get the total dust impact charge.
@@ -1018,6 +1014,9 @@ def estimate_dust_mass(
         signal.
     waveform_name : str
         Channel name used to select channel-specific fit bounds.
+    baseline_mean : float or None
+        L1b mean baseline used to initialize the fit. Calculated from the first
+        3 microseconds when omitted in standalone calls.
 
     Returns
     -------
@@ -1037,23 +1036,14 @@ def estimate_dust_mass(
     """
     signal = np.array(target_signal.data)
     time = np.array(low_sampling_time.data)
-    # window_start = float(np.min(time))
-    window_stop = float(np.min(time)) + 5.0
-    # good_mask = np.logical_and(time >= window_start, time <= window_stop)
-    good_mask = time <= window_stop
-    if not np.any(good_mask):
-        logger.warning(
-            "Unable to find baseline noise. "
-            f"There is no signal in the first 5 microseconds of the waveform "
-            f"(Beginning to {window_stop} us)."
-        )
+    if baseline_mean is None:
+        baseline, _ = calculate_baseline_and_rms(signal, time)
+        baseline_mean = float(baseline)
     if remove_noise:
         logger.debug(
             "estimate_dust_mass fits the raw low-rate waveform directly; "
             "remove_noise is ignored for this fit path."
         )
-    signal_before_impact = signal[good_mask]
-    baseline_mean = float(np.mean(signal_before_impact))
     channel_name = waveform_name or str(getattr(target_signal, "name", ""))
 
     # Initial Guess for the parameters of the ion grid signal