From 9cc2764087ec0204e5a38c8e1d0e4cb78e68eda8 Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 26 Oct 2026 11:37:17 -0600 Subject: [PATCH 1/4] Add per-waveform baseline and RMS statistics to IDEX L1b --- .../config/imap_idex_l1b_variable_attrs.yaml | 100 ++++++++++++++ imap_processing/idex/idex_l1b.py | 45 +++++- imap_processing/idex/idex_l2a.py | 88 +++++++----- imap_processing/idex/idex_utils.py | 51 +++++++ imap_processing/tests/idex/test_idex_l1b.py | 62 ++++++++- imap_processing/tests/idex/test_idex_l2a.py | 128 +++++++++++++++--- imap_processing/tests/idex/test_idex_utils.py | 62 +++++++++ 7 files changed, 479 insertions(+), 57 deletions(-) create mode 100644 imap_processing/tests/idex/test_idex_utils.py diff --git a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml index d02b106787..5303420673 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/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_TYPE: support_data + 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. + 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) diff --git a/imap_processing/idex/idex_l1b.py b/imap_processing/idex/idex_l1b.py index 366512be13..782ec0de2c 100644 --- a/imap_processing/idex/idex_l1b.py +++ b/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, diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 0c3eff6b7a..f694046293 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -32,14 +32,20 @@ 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. + Legacy TOF baseline window, retained for compatibility. + + L1b statistics and L2a processing now use the first 3 microseconds instead. Attributes ---------- @@ -128,7 +134,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.get("tof_high_baseline"), + l1b_dataset.get("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 +205,27 @@ 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.get(f"{waveform.lower()}_baseline") + if baseline is None: + # Older L1b products do not contain the event-level statistics. + baseline_values, _ = calculate_baseline_and_rms( + l1b_dataset[waveform].data, ls_time.data + ) + baseline = xr.DataArray(baseline_values, dims="epoch") # 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 +237,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 +807,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 +822,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, optional + Per-event L1b mean baseline. Calculated from the first 3 microseconds + when baseline or rms is omitted, for compatibility with older products. + rms : xarray.DataArray or numpy.ndarray, optional + 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 +1022,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 +1038,9 @@ def estimate_dust_mass( signal. waveform_name : str Channel name used to select channel-specific fit bounds. + baseline_mean : float, optional + L1b mean baseline used to initialize the fit. Calculated from the first + 3 microseconds when omitted, for compatibility with older products. Returns ------- @@ -1037,23 +1060,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 diff --git a/imap_processing/idex/idex_utils.py b/imap_processing/idex/idex_utils.py index 924cbfa74f..1be649631b 100644 --- a/imap_processing/idex/idex_utils.py +++ b/imap_processing/idex/idex_utils.py @@ -1,11 +1,62 @@ """Contains helper functions to support IDEX processing.""" +import logging + +import numpy as np import pandas as pd import xarray as xr +from numpy.typing import NDArray from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.idex.idex_constants import IDEX_10_DAY_RANGES_PATH +logger = logging.getLogger(__name__) + + +def calculate_baseline_and_rms( + signal: NDArray, time: NDArray +) -> tuple[NDArray, NDArray]: + """ + Calculate the mean baseline and RMS noise in the first 3 microseconds. + + Parameters + ---------- + signal : numpy.ndarray + Waveform samples in engineering units, with samples on the last axis. + time : numpy.ndarray + Corresponding sample times in microseconds. Each record uses the half-open + interval from its earliest finite time to that time plus 3 microseconds. + + Returns + ------- + baseline : numpy.ndarray + Mean of finite samples in the window, or NaN if none are available. + rms : numpy.ndarray + Root mean square deviation from the baseline (population standard + deviation, ddof=0), in the waveform's units. NaN for an empty window. + """ + signal = np.asarray(signal, dtype=float) + time = np.asarray(time, dtype=float) + finite_time = np.isfinite(time) + start = np.min( + np.where(finite_time, time, np.inf), axis=-1, keepdims=True, initial=np.inf + ) + mask = finite_time & (time < start + 3.0) & np.isfinite(signal) + count = np.count_nonzero(mask, axis=-1) + baseline = np.full(signal.shape[:-1], np.nan) + np.divide( + np.sum(np.where(mask, signal, 0.0), axis=-1), + count, + out=baseline, + where=count > 0, + ) + residual = np.where(mask, signal - np.expand_dims(baseline, -1), 0.0) + variance = np.full_like(baseline, np.nan) + np.divide(np.sum(residual**2, axis=-1), count, out=variance, where=count > 0) + if np.any(count == 0): + logger.warning("Unable to find baseline noise in the first 3 microseconds.") + return baseline, np.sqrt(variance) + def get_idex_attrs(data_level: str) -> ImapCdfAttributes: """ diff --git a/imap_processing/tests/idex/test_idex_l1b.py b/imap_processing/tests/idex/test_idex_l1b.py index fc7972063d..ec578b5fb0 100644 --- a/imap_processing/tests/idex/test_idex_l1b.py +++ b/imap_processing/tests/idex/test_idex_l1b.py @@ -9,17 +9,19 @@ from imap_processing import imap_module_directory from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes -from imap_processing.cdf.utils import write_cdf +from imap_processing.cdf.utils import load_cdf, write_cdf from imap_processing.idex.idex_constants import DT_BLOCK, ConversionFactors from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES, EVENT_FLAG_NAMES from imap_processing.idex.idex_l1b import ( TRIGGER_LABELS, EventMessage, TriggerOrigin, + convert_waveforms, get_event_dead_time, get_spice_data, get_trigger_mode_and_level, get_trigger_origin, + get_waveform_statistics, idex_l1b, unpack_instrument_settings, ) @@ -82,6 +84,64 @@ def test_idex_cdf_file(l1b_dataset: xr.Dataset): file_name = write_cdf(l1b_dataset) assert file_name.exists() assert file_name.name == "imap_idex_l1b_sci-10days_20231218_v001.0001.cdf" + restored = load_cdf(file_name) + for channel in ConversionFactors: + for suffix in ("baseline", "rms"): + name = f"{channel.name.lower()}_{suffix}" + assert restored[name].dims == ("epoch",) + np.testing.assert_allclose(restored[name], l1b_dataset[name]) + assert restored[name].attrs["UNITS"] == l1b_dataset[name].attrs["UNITS"] + + +def test_waveform_statistics_all_channels(): + """Statistics use the correct time axis and converted units for every channel.""" + high_dim = "time_high_sample_rate_index" + low_dim = "time_low_sample_rate_index" + dataset = xr.Dataset( + { + "time_high_sample_rate": ( + ("epoch", high_dim), + [[-7.0, -6.0, -5.0, -4.0], [10.0, 11.0, 12.0, 13.0]], + ), + "time_low_sample_rate": ( + ("epoch", low_dim), + [[0.0, 0.8, 1.6, 2.4, 3.0], [-9.0, -8.2, -7.4, -6.6, -6.0]], + ), + }, + coords={"epoch": [11, 22]}, + ) + for channel in ConversionFactors: + is_tof = channel.name.startswith("TOF") + dataset[channel.name] = xr.DataArray( + [[10, 12, 14, 900], [30, 30, 30, 900]] + if is_tof + else [[20, 22, 24, 26, 900], [40, 40, 40, 40, 900]], + dims=("epoch", high_dim if is_tof else low_dim), + ) + attrs = get_idex_attrs("l1b") + dataset.update(convert_waveforms(dataset, attrs)) + original = dataset.copy(deep=True) + + statistics = get_waveform_statistics(dataset, attrs) + + assert len(statistics) == 12 + for channel in ConversionFactors: + is_tof = channel.name.startswith("TOF") + expected_baseline = np.array([12.0, 30.0] if is_tof else [23.0, 40.0]) + expected_rms = np.array([np.sqrt(8 / 3) if is_tof else np.sqrt(5), 0.0]) + expected_statistics = (("baseline", expected_baseline), ("rms", expected_rms)) + for suffix, expected in expected_statistics: + field = statistics[f"{channel.name.lower()}_{suffix}"] + np.testing.assert_allclose(field, expected * channel.value) + np.testing.assert_array_equal(field.epoch, dataset.epoch) + assert field.dims == ("epoch",) + assert field.attrs["DEPEND_0"] == "epoch" + assert "DEPEND_1" not in field.attrs + assert "LABL_PTR_1" not in field.attrs + assert field.attrs["UNITS"] == ("mA" if is_tof else "pC") + maximum_dn = 1023 if is_tof else 4095 + assert field.attrs["VALIDMAX"] == pytest.approx(maximum_dn * channel.value) + xr.testing.assert_identical(dataset, original) def test_idex_waveform_units(l1b_dataset: xr.Dataset): diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index dba158713c..0772de5c61 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -13,7 +13,6 @@ from imap_processing.idex import idex_constants from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES from imap_processing.idex.idex_l2a import ( - BaselineNoiseTime, _mask_non_science_derived_estimates, _mask_saturated_derived_estimates, analyze_peaks, @@ -26,6 +25,7 @@ chi_square, estimate_dust_mass, fit_impact, + idex_l2a, invert_rise_time_to_velocity, load_calibration_files, log_smooth_powerlaw, @@ -33,6 +33,68 @@ sine_fit, time_to_mass, ) +from imap_processing.idex.idex_utils import calculate_baseline_and_rms + + +@pytest.mark.parametrize("stored_statistics", [True, False]) +def test_l2a_uses_per_event_baselines(l1b_dataset, stored_statistics): + """L2a forwards each baseline once and computes missing legacy fields in bulk.""" + dataset = l1b_dataset.isel(epoch=slice(0, 2)).copy(deep=True) + expected_baselines = { + channel: dataset[f"{channel.lower()}_baseline"].values.copy() + for channel in ("Target_Low", "Target_High", "Ion_Grid") + } + if not stored_statistics: + dataset = dataset.drop_vars( + [name for name in dataset if name.endswith(("_baseline", "_rms"))] + ) + + def fit_waveform(time, signal, remove_noise, waveform_name, baseline_mean): + return ( + np.array([0.0, baseline_mean, 1.0, 0.371, 37.1]), + 1.0, + 0.0, + 0.0, + np.zeros_like(signal), + ) + + with ( + mock.patch( + "imap_processing.idex.idex_l2a.load_calibration_files", + return_value=( + np.array([1.27, -0.2, -2.1, -0.37, 5.3, 13.3, 13.3, 0.28]), + np.array([0.06, 2.8, 5.9, 4.1, 13.0, 22.7, 8.2, 0.40]), + ), + ), + mock.patch( + "imap_processing.idex.idex_l2a.time_to_mass", + return_value=(None, None, np.zeros_like(dataset.TOF_High)), + ), + mock.patch( + "imap_processing.idex.idex_l2a.analyze_peaks", + return_value=(np.zeros((2, 3)), np.ones(2), np.zeros(2), np.zeros(2)), + ), + mock.patch( + "imap_processing.idex.idex_l2a.estimate_dust_mass", + side_effect=fit_waveform, + ) as fit, + mock.patch( + "imap_processing.idex.idex_l2a.calculate_baseline_and_rms", + wraps=calculate_baseline_and_rms, + ) as calculate, + ): + result = idex_l2a(dataset, {}) + + assert calculate.call_count == (0 if stored_statistics else 4) + assert fit.call_count == 6 + for channel, expected in expected_baselines.items(): + supplied = [ + call.args[4] for call in fit.call_args_list if call.args[3] == channel + ] + np.testing.assert_allclose(supplied, expected) + np.testing.assert_allclose( + result[f"{channel.lower()}_fit_parameters"][:, 1], expected + ) def test_non_science_derived_estimates_are_nan() -> None: @@ -357,15 +419,15 @@ def test_calculate_snr(): """Tests the functionality of calculate_snr().""" step = 0.5 max_tof = 10 - time = np.arange(BaselineNoiseTime.START, 5, step) + time = np.arange(-7, 5, step) - # Create a baseline noise array with an std of 1 and mean of 1 - baseline_noise = np.asarray([0, 0, 1, 2, 2]) + # The first 3 microseconds have a mean of 1 and population RMS noise of 1. + baseline_noise = np.asarray([0, 2, 0, 2, 0, 2]) signal_length = len(time) - len(baseline_noise) tof_signal = np.full(int(signal_length), max_tof) tof = np.tile(np.append(baseline_noise, tof_signal), (3, 1)) - time = np.tile(time, (3, 1)) + time = np.tile(time, (3, 1)) + np.array([0, 10, 20])[:, None] snr = calculate_snr(tof, time) @@ -375,17 +437,39 @@ def test_calculate_snr(): def test_calculate_snr_warning(caplog): """Tests that calculate_snr() throws warning if no baseline noise is found.""" - time = np.tile(np.arange(10), (3, 1)) + time = np.full((3, 10), np.nan) tof = np.ones_like(time) with caplog.at_level("WARNING"): - calculate_snr(tof, time) + snr = calculate_snr(tof, time) + assert np.all(np.isnan(snr)) assert any( "Unable to find baseline noise" in message for message in caplog.text.splitlines() ) +def test_calculate_snr_reuses_statistics(): + tof = xr.DataArray([[0.0, 11.0], [0.0, 22.0]], dims=("epoch", "sample")) + time = xr.zeros_like(tof) + with mock.patch( + "imap_processing.idex.idex_l2a.calculate_baseline_and_rms", + side_effect=AssertionError("L1b statistics must not be recalculated"), + ): + snr = calculate_snr(tof, time, np.array([3.0, 6.0]), np.array([2.0, 4.0])) + + np.testing.assert_allclose(snr, [4.0, 4.0]) + + +def test_calculate_snr_zero_rms(): + tof = xr.DataArray([[1.0, 1.0, 1.0, 2.0], [1.0, 1.0, 1.0, 1.0]]) + time = xr.DataArray(np.tile(np.arange(4), (2, 1))) + + snr = calculate_snr(tof, time) + + np.testing.assert_array_equal(snr, [np.inf, np.nan]) + + def test_analyze_peaks_warning(caplog): """Tests that analyze_peaks() throws warning if the emg curve fit fails.""" # Create a 2d list of peak indices @@ -684,20 +768,10 @@ def test_estimate_dust_mass_logs_baseline_warning(caplog): """ Test that estimate_dust_mass() logs a warning if no baseline is found. """ - time = xr.DataArray(np.linspace(-60, 60, 16)) + time = xr.DataArray(np.full(16, np.nan)) signal = xr.DataArray(np.linspace(0, 1, 16)) - original_any = np.any - - def fake_any(arr): - fake_any.calls += 1 - if fake_any.calls == 1: - return False - return original_any(arr) - - fake_any.calls = 0 with ( - mock.patch("imap_processing.idex.idex_l2a.np.any", side_effect=fake_any), mock.patch( "imap_processing.idex.idex_l2a.curve_fit", return_value=(np.array([0.0, 0.0, 1.0, 0.371, 37.1]), None), @@ -712,6 +786,24 @@ def fake_any(arr): ) +def test_estimate_dust_mass_reuses_baseline(): + time = xr.DataArray(np.arange(-6.0, 6.0)) + signal = xr.DataArray(np.arange(12.0)) + with ( + mock.patch( + "imap_processing.idex.idex_l2a.calculate_baseline_and_rms", + side_effect=AssertionError("L1b baseline must not be recalculated"), + ), + mock.patch( + "imap_processing.idex.idex_l2a.curve_fit", + return_value=(np.array([0.0, 7.0, 10.0, 0.371, 37.1]), None), + ) as fit, + ): + estimate_dust_mass(time, signal, baseline_mean=7.0) + + assert fit.call_args.kwargs["p0"][1] == 7.0 + + def test_estimate_dust_mass_remove_noise_logs_debug(caplog): """ Test that estimate_dust_mass() logs that remove_noise is ignored. diff --git a/imap_processing/tests/idex/test_idex_utils.py b/imap_processing/tests/idex/test_idex_utils.py new file mode 100644 index 0000000000..4dd52eafba --- /dev/null +++ b/imap_processing/tests/idex/test_idex_utils.py @@ -0,0 +1,62 @@ +"""Tests of shared IDEX waveform statistics.""" + +import numpy as np +import pytest + +from imap_processing.idex.idex_utils import calculate_baseline_and_rms + + +def test_baseline_and_rms_use_first_three_microseconds_per_event(): + # Offsets and sample intervals vary between events. The endpoint is excluded. + time = np.array([[-8.0, -7.0, -6.0, -5.0], [4.0, 4.5, 6.9, 7.0]]) + signal = np.array([[2.0, 4.0, 6.0, 1000.0], [12.0, 14.0, 16.0, 2000.0]]) + original = signal.copy() + + baseline, rms = calculate_baseline_and_rms(signal, time) + + np.testing.assert_allclose(baseline, [4.0, 14.0]) + np.testing.assert_allclose(rms, np.sqrt([8 / 3, 8 / 3])) + np.testing.assert_array_equal(signal, original) + + +@pytest.mark.parametrize( + "signal,time,expected_baseline,expected_rms", + [ + ([5.0, 5.0], [0.0, 1.0], 5.0, 0.0), + ([7.0], [-2.0], 7.0, 0.0), + ([np.nan, 2.0, np.inf, 4.0], [-1.0, 0.0, 0.5, 1.0], 3.0, 1.0), + ([1000.0, 2.0, 4.0], [np.nan, 0.0, 1.0], 3.0, 1.0), + ([1000.0, 2.0, 4.0], [-np.inf, 0.0, 1.0], 3.0, 1.0), + # Missing initial signal samples must not move the window into the pulse. + ([np.nan, 2.0, 1000.0], [0.0, 1.0, 3.0], 2.0, 0.0), + ], +) +def test_baseline_and_rms_short_or_nonfinite_samples( + signal, time, expected_baseline, expected_rms +): + baseline, rms = calculate_baseline_and_rms(np.array(signal), np.array(time)) + + np.testing.assert_allclose(baseline, expected_baseline) + np.testing.assert_allclose(rms, expected_rms) + + +@pytest.mark.parametrize( + "signal,time", + [([np.nan, np.inf], [0.0, 1.0]), ([1.0, 2.0], [np.nan, np.inf]), ([], [])], +) +def test_baseline_and_rms_missing_window(signal, time, caplog): + baseline, rms = calculate_baseline_and_rms(np.array(signal), np.array(time)) + + assert np.isnan(baseline) + assert np.isnan(rms) + assert "Unable to find baseline noise" in caplog.text + + +def test_baseline_and_rms_missing_event_does_not_affect_other_events(): + baseline, rms = calculate_baseline_and_rms( + np.array([[np.nan, np.nan], [1.0, 3.0]]), + np.array([[0.0, 1.0], [0.0, 1.0]]), + ) + + np.testing.assert_allclose(baseline, [np.nan, 2.0]) + np.testing.assert_allclose(rms, [np.nan, 1.0]) From a6c2455dfc52c49ec23740e5359b8b0cfc5737e5 Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 26 Oct 2026 12:13:13 -0600 Subject: [PATCH 2/4] Alphabetize IDEX waveform statistics metadata attributes --- imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml index 5303420673..4c03c1a5eb 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -46,12 +46,12 @@ l1b_target_base: &l1b_target_base waveform_statistics_base: &waveform_statistics_base <<: *l1b_data_base DISPLAY_TYPE: time_series - VAR_TYPE: support_data 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 From ed1aeb01b90bd3d754624420024ed47efa27b0a3 Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 26 Oct 2026 12:35:29 -0600 Subject: [PATCH 3/4] Fix IDEX optional parameter documentation references --- imap_processing/idex/idex_l2a.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index f694046293..710c06d737 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -822,10 +822,10 @@ def calculate_snr( The time of flight array. hs_time : xarray.DataArray The high sampling time array. - baseline : xarray.DataArray or numpy.ndarray, optional + 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, for compatibility with older products. - rms : xarray.DataArray or numpy.ndarray, optional + rms : xarray.DataArray or numpy.ndarray or None Per-event L1b baseline-subtracted RMS noise (ddof=0). Returns @@ -1038,7 +1038,7 @@ def estimate_dust_mass( signal. waveform_name : str Channel name used to select channel-specific fit bounds. - baseline_mean : float, optional + baseline_mean : float or None L1b mean baseline used to initialize the fit. Calculated from the first 3 microseconds when omitted, for compatibility with older products. From 2e6e56c83fa3f2fc477cb6a7a6bfe5a4188e1046 Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 5 Nov 2026 12:53:04 -0700 Subject: [PATCH 4/4] Address review of IDEX baseline and RMS statistics --- imap_processing/idex/idex_l2a.py | 38 ++++----------------- imap_processing/idex/idex_utils.py | 20 ++++++++--- imap_processing/tests/idex/test_idex_l2a.py | 34 ++++++++++++------ 3 files changed, 46 insertions(+), 46 deletions(-) diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 710c06d737..d5709274d1 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -17,7 +17,6 @@ """ import logging -from enum import IntEnum import numpy as np import pandas as pd @@ -41,24 +40,6 @@ logger = logging.getLogger(__name__) -class BaselineNoiseTime(IntEnum): - """ - Legacy TOF baseline window, retained for compatibility. - - L1b statistics and L2a processing now use the first 3 microseconds instead. - - 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. @@ -100,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. @@ -137,8 +119,8 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: snr = calculate_snr( tof_high, hs_time, - l1b_dataset.get("tof_high_baseline"), - l1b_dataset.get("tof_high_rms"), + 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 @@ -205,13 +187,7 @@ 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.get(f"{waveform.lower()}_baseline") - if baseline is None: - # Older L1b products do not contain the event-level statistics. - baseline_values, _ = calculate_baseline_and_rms( - l1b_dataset[waveform].data, ls_time.data - ) - baseline = xr.DataArray(baseline_values, dims="epoch") + baseline = l1b_dataset[f"{waveform.lower()}_baseline"] # Get the dust mass estimates and fit results fit_results = xr.apply_ufunc( estimate_dust_mass, @@ -824,7 +800,7 @@ def calculate_snr( 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, for compatibility with older products. + 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). @@ -1040,7 +1016,7 @@ def estimate_dust_mass( 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, for compatibility with older products. + 3 microseconds when omitted in standalone calls. Returns ------- diff --git a/imap_processing/idex/idex_utils.py b/imap_processing/idex/idex_utils.py index 1be649631b..ba5a3eb3f8 100644 --- a/imap_processing/idex/idex_utils.py +++ b/imap_processing/idex/idex_utils.py @@ -22,10 +22,13 @@ def calculate_baseline_and_rms( Parameters ---------- signal : numpy.ndarray - Waveform samples in engineering units, with samples on the last axis. + Waveform samples in engineering units. Shape is (samples,) for one event + or (events, samples) for multiple events; each row is one event. time : numpy.ndarray - Corresponding sample times in microseconds. Each record uses the half-open - interval from its earliest finite time to that time plus 3 microseconds. + Corresponding sample times in microseconds, with shape (samples,) for a + shared time axis or (events, samples) for event-specific times. Each event + uses the half-open interval from its earliest finite time to that time + plus 3 microseconds. Returns ------- @@ -38,11 +41,19 @@ def calculate_baseline_and_rms( signal = np.asarray(signal, dtype=float) time = np.asarray(time, dtype=float) finite_time = np.isfinite(time) + # Ignore invalid times when locating each event's start. Keep the sample axis + # so the start broadcasts against that event's sample times. start = np.min( np.where(finite_time, time, np.inf), axis=-1, keepdims=True, initial=np.inf ) + # Use only finite samples in [start, start + 3 microseconds). An event with + # no finite times has start=inf and contributes no samples. mask = finite_time & (time < start + 3.0) & np.isfinite(signal) count = np.count_nonzero(mask, axis=-1) + if np.any(count == 0): + logger.warning("Unable to find baseline noise in the first 3 microseconds.") + # Masked sums exclude invalid samples; guarded division leaves empty events + # as NaN rather than dividing by zero. baseline = np.full(signal.shape[:-1], np.nan) np.divide( np.sum(np.where(mask, signal, 0.0), axis=-1), @@ -50,11 +61,10 @@ def calculate_baseline_and_rms( out=baseline, where=count > 0, ) + # Center each event on its own baseline and use the population variance. residual = np.where(mask, signal - np.expand_dims(baseline, -1), 0.0) variance = np.full_like(baseline, np.nan) np.divide(np.sum(residual**2, axis=-1), count, out=variance, where=count > 0) - if np.any(count == 0): - logger.warning("Unable to find baseline noise in the first 3 microseconds.") return baseline, np.sqrt(variance) diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 0772de5c61..57931ee9a1 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -36,18 +36,26 @@ from imap_processing.idex.idex_utils import calculate_baseline_and_rms -@pytest.mark.parametrize("stored_statistics", [True, False]) -def test_l2a_uses_per_event_baselines(l1b_dataset, stored_statistics): - """L2a forwards each baseline once and computes missing legacy fields in bulk.""" +@pytest.mark.parametrize( + "missing_statistic", + [ + None, + "tof_high_baseline", + "tof_high_rms", + "target_low_baseline", + "target_high_baseline", + "ion_grid_baseline", + ], +) +def test_l2a_uses_required_l1b_statistics(l1b_dataset, missing_statistic): + """L2a uses stored statistics and rejects inputs missing required fields.""" dataset = l1b_dataset.isel(epoch=slice(0, 2)).copy(deep=True) expected_baselines = { channel: dataset[f"{channel.lower()}_baseline"].values.copy() for channel in ("Target_Low", "Target_High", "Ion_Grid") } - if not stored_statistics: - dataset = dataset.drop_vars( - [name for name in dataset if name.endswith(("_baseline", "_rms"))] - ) + if missing_statistic is not None: + dataset = dataset.drop_vars(missing_statistic) def fit_waveform(time, signal, remove_noise, waveform_name, baseline_mean): return ( @@ -83,9 +91,15 @@ def fit_waveform(time, signal, remove_noise, waveform_name, baseline_mean): wraps=calculate_baseline_and_rms, ) as calculate, ): - result = idex_l2a(dataset, {}) - - assert calculate.call_count == (0 if stored_statistics else 4) + if missing_statistic is not None: + with pytest.raises(KeyError, match=missing_statistic): + idex_l2a(dataset, {}) + else: + result = idex_l2a(dataset, {}) + + calculate.assert_not_called() + if missing_statistic is not None: + return assert fit.call_count == 6 for channel, expected in expected_baselines.items(): supplied = [