diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index d5709274d1..ec72822d55 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -1003,6 +1003,10 @@ def estimate_dust_mass( """ Filter and fit the target or ion grid signals to get the total dust impact charge. + Fit only when a sample after the first 3 microseconds exceeds 4 times the + baseline RMS noise. Target excursions must be positive; Ion Grid excursions + may have either polarity. Skipped fits return NaN products. + Parameters ---------- low_sampling_time : xarray.DataArray @@ -1036,6 +1040,13 @@ def estimate_dust_mass( """ signal = np.array(target_signal.data) time = np.array(low_sampling_time.data) + null_result = ( + np.full(5, np.nan), + np.nan, + np.nan, + np.nan, + np.full(time.shape, np.nan), + ) if baseline_mean is None: baseline, _ = calculate_baseline_and_rms(signal, time) baseline_mean = float(baseline) @@ -1060,9 +1071,7 @@ def estimate_dust_mass( amplitude = float("nan") if channel_name != "Ion_Grid" and (not np.isfinite(amplitude) or amplitude <= 0.0): amplitude = float(np.max(signal) - baseline_mean) - if not np.isfinite(amplitude): - amplitude = float(np.max(signal)) - if channel_name != "Ion_Grid" and amplitude <= 0.0: + if not np.isfinite(amplitude) or (channel_name != "Ion_Grid" and amplitude <= 0.0): amplitude = float(np.max(signal)) rise_time_0 = 0.371 # How fast the signal rises (us) @@ -1084,6 +1093,9 @@ def estimate_dust_mass( [np.inf, np.inf, amplitude_upper_bound, np.inf, np.inf], ) + if not _has_low_rate_signal(time, signal, channel_name): + return null_result + try: with np.errstate(invalid="ignore", over="ignore"): param, _ = curve_fit( @@ -1101,13 +1113,7 @@ def estimate_dust_mass( f"Signal range: {min(signal):.2f} to {max(signal):.2f}\n" "Returning None." ) - return ( - np.full(len(p0), np.nan), - np.nan, - np.nan, - np.nan, - np.full_like(time, np.nan), - ) + return null_result impact_fit = fit_impact(time, *param) # Evaluate the analytic extremum of the fitted pulse instead of relying on the @@ -1173,6 +1179,41 @@ def fit_impact( ) +def _has_low_rate_signal(time: NDArray, signal: NDArray, channel_name: str) -> bool: + """ + Check for a sample above 4 times the baseline RMS noise. + + Parameters + ---------- + time : numpy.ndarray + Sample times in microseconds. + signal : numpy.ndarray + Target or Ion Grid waveform in engineering units. + channel_name : str + Waveform name. Ion Grid permits either polarity; targets require a + positive excursion from the baseline. + + Returns + ------- + bool + Whether a finite sample after the first 3 microseconds exceeds the + threshold. False when no finite baseline samples are available. + """ + noise_mean, noise_rms = calculate_baseline_and_rms(signal, time) + if not (np.isfinite(noise_mean) and np.isfinite(noise_rms)): + return False + + # Only finite samples outside the shared statistics' first 3 microseconds + # can establish the presence of a pulse. + finite_time = np.isfinite(time) + window_start = float(np.min(time[finite_time])) + signal_window = finite_time & (time >= window_start + 3.0) & np.isfinite(signal) + excursions = signal[signal_window] - noise_mean + if channel_name == "Ion_Grid": + excursions = np.abs(excursions) + return bool(np.any(excursions > 4.0 * noise_rms)) + + def remove_signal_noise( time: np.ndarray, signal: np.ndarray, good_mask: np.ndarray ) -> NDArray: diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 57931ee9a1..e5fc378473 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -804,6 +804,10 @@ def test_estimate_dust_mass_reuses_baseline(): time = xr.DataArray(np.arange(-6.0, 6.0)) signal = xr.DataArray(np.arange(12.0)) with ( + # Isolate fit initialization; the threshold gate computes its own RMS. + mock.patch( + "imap_processing.idex.idex_l2a._has_low_rate_signal", return_value=True + ), mock.patch( "imap_processing.idex.idex_l2a.calculate_baseline_and_rms", side_effect=AssertionError("L1b baseline must not be recalculated"), @@ -843,9 +847,9 @@ def test_estimate_dust_mass_remove_noise_logs_debug(caplog): ) -def test_estimate_dust_mass_nonfinite_signal_fallbacks(): +def test_estimate_dust_mass_skips_nonfinite_signal(): """ - Test fallback handling when the input signal is entirely non-finite. + Test that an entirely non-finite waveform is skipped without fitting. """ time = xr.DataArray(np.linspace(-60, 60, 16)) signal = xr.DataArray(np.full(16, np.nan)) @@ -854,14 +858,16 @@ def test_estimate_dust_mass_nonfinite_signal_fallbacks(): "imap_processing.idex.idex_l2a.curve_fit", return_value=(np.array([0.0, 0.0, 1.0, 0.371, 37.1]), None), ) as mocked_curve_fit: - estimate_dust_mass(time, signal) + result = estimate_dust_mass(time, signal) - assert np.isnan(mocked_curve_fit.call_args.kwargs["p0"][2]) + mocked_curve_fit.assert_not_called() + for value in result: + assert np.all(np.isnan(value)) -def test_estimate_dust_mass_non_ion_grid_negative_amplitude_fallback(): +def test_estimate_dust_mass_skips_constant_negative_signal(): """ - Test the non-Ion_Grid negative-amplitude fallback path. + Test that a constant negative baseline does not qualify as a signal. """ time = xr.DataArray(np.linspace(-60, 60, 16)) signal = xr.DataArray(np.full(16, -2.0)) @@ -870,9 +876,11 @@ def test_estimate_dust_mass_non_ion_grid_negative_amplitude_fallback(): "imap_processing.idex.idex_l2a.curve_fit", return_value=(np.array([0.0, -2.0, -2.0, 0.371, 37.1]), None), ) as mocked_curve_fit: - estimate_dust_mass(time, signal) + result = estimate_dust_mass(time, signal) - assert mocked_curve_fit.call_args.kwargs["p0"][2] == -2.0 + mocked_curve_fit.assert_not_called() + for value in result: + assert np.all(np.isnan(value)) def test_estimate_dust_mass_ion_grid_negative_amplitude_bounds(): diff --git a/imap_processing/tests/idex/test_idex_signal_threshold.py b/imap_processing/tests/idex/test_idex_signal_threshold.py new file mode 100644 index 0000000000..7fd316a92a --- /dev/null +++ b/imap_processing/tests/idex/test_idex_signal_threshold.py @@ -0,0 +1,276 @@ +"""Test the signal threshold for Target and Ion Grid waveform fits.""" + +from unittest import mock + +import numpy as np +import pytest +import xarray as xr + +from imap_processing.idex.idex_l2a import ( + _has_low_rate_signal, + estimate_dust_mass, + idex_l2a, + log_smooth_powerlaw, +) + + +@pytest.mark.parametrize("waveform", ["Target_High", "Target_Low", "Ion_Grid"]) +@pytest.mark.parametrize( + "time,signal,expected", + [ + # Invalid samples do not affect the baseline (10) or population RMS (1). + ([np.nan, 0, 1, 2, 3], [1000, 9, np.nan, 11, 14.01], True), + ([np.inf, 0, 1, 2, 3], [1000, 9, np.inf, 11, 14.0], False), + # The earliest time defines the window even if its signal is missing. + ([0, 1, 2, 3], [np.nan, 9, 11, 14.01], True), + # The time axis need not be sorted; the endpoint is outside the baseline. + ([3, 2, 0, 1], [14.01, 11, 9, 10], True), + # Nonfinite pulse samples cannot trigger a fit. + ([0, 1, 3, 4], [9, 11, np.inf, np.nan], False), + ([0, 1, np.nan, np.inf], [9, 11, 1000, 1000], False), + ([0, 1, 3], [np.nan, np.inf, 1000], False), + ([], [], False), + ], +) +def test_low_rate_signal_uses_finite_baseline_window(waveform, time, signal, expected): + """Shared statistics retain the gate's finite-sample and window rules.""" + assert _has_low_rate_signal(np.array(time), np.array(signal), waveform) is expected + + +@pytest.mark.parametrize("failure", [False, True]) +@pytest.mark.parametrize("time_dtype", [int, float]) +def test_low_rate_null_result_shapes_and_dtype(failure, time_dtype): + """Skipped and failed fits return float NaNs even with integer sample times.""" + time = xr.DataArray(np.arange(-6, 6, dtype=time_dtype)) + signal = np.zeros(time.size) + if failure: + signal[-1] = 10.0 + with mock.patch( + "imap_processing.idex.idex_l2a.curve_fit", + side_effect=RuntimeError("fit failed"), + ) as optimizer: + result = estimate_dust_mass( + time, xr.DataArray(signal), waveform_name="Target_High" + ) + + assert optimizer.call_count == int(failure) + assert len(result) == 5 + assert result[0].shape == (5,) + assert result[4].shape == signal.shape + for value in result: + assert np.asarray(value).dtype.kind == "f" + assert np.all(np.isnan(value)) + + +@pytest.mark.parametrize("waveform", ["Target_High", "Target_Low", "Ion_Grid"]) +@pytest.mark.parametrize( + "excursion", + [ + 0.0, + 1.5, + 3.375, + 3.5, + 3.88, + 3.99, + 4.0, + 4.01, + -1.5, + -3.375, + -3.5, + -3.88, + -3.99, + -4.0, + -4.01, + ], +) +@pytest.mark.parametrize("pulse_index", [6, -1]) +def test_low_rate_fit_requires_sample_above_noise(waveform, excursion, pulse_index): + """One sample must strictly exceed 4 sigma, with channel-specific polarity.""" + time = xr.DataArray(np.arange(-6.0, 6.0, 0.5)) + signal = np.full(time.size, 10.0) + # The first 3 microseconds have baseline 10 and population sigma 1. + signal[:6] += [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0] + signal[pulse_index] += excursion + target = xr.DataArray(signal, name=waveform) + parameters = np.array([0.0, 10.0, excursion, 0.371, 37.1]) + + with mock.patch( + "imap_processing.idex.idex_l2a.curve_fit", return_value=(parameters, None) + ) as optimizer: + result = estimate_dust_mass(time, target) + + should_fit = abs(excursion) > 4.0 if waveform == "Ion_Grid" else excursion > 4.0 + if should_fit: + optimizer.assert_called_once() + np.testing.assert_array_equal(result[0], parameters) + assert np.all(np.isfinite(result[4])) + else: + optimizer.assert_not_called() + assert result[0].shape == (5,) + assert result[4].shape == signal.shape + for value in result: + assert np.all(np.isnan(value)) + + +@pytest.mark.parametrize("waveform", ["Target_High", "Target_Low", "Ion_Grid"]) +@pytest.mark.parametrize("excursion", [0.0, 1.0, -1.0]) +def test_low_rate_fit_zero_noise(waveform, excursion): + """A flat baseline requires a nonzero excursion without dividing by zero.""" + time = xr.DataArray(np.arange(-6.0, 6.0, 0.5)) + signal = np.full(time.size, -2.0) + signal[-1] += excursion + with mock.patch( + "imap_processing.idex.idex_l2a.curve_fit", + return_value=(np.array([0.0, -2.0, excursion, 0.371, 37.1]), None), + ) as optimizer: + result = estimate_dust_mass(time, xr.DataArray(signal), waveform_name=waveform) + + should_fit = excursion != 0.0 if waveform == "Ion_Grid" else excursion > 0.0 + assert optimizer.call_count == int(should_fit) + if not should_fit: + for value in result: + assert np.all(np.isnan(value)) + + +@pytest.mark.parametrize("waveform", ["Target_High", "Target_Low", "Ion_Grid"]) +def test_low_rate_fit_does_not_count_baseline_excursions(waveform): + """An excursion confined to the noise window is not evidence of a pulse.""" + time = xr.DataArray(np.arange(96) / 8.0 - 6.0) + signal = np.full(time.size, 4.0 / 24.0) + # The baseline outlier exceeds 4 sigma but must not trigger a fit. + signal[:24] = 0.0 + signal[23] = 4.0 + with mock.patch("imap_processing.idex.idex_l2a.curve_fit") as optimizer: + result = estimate_dust_mass(time, xr.DataArray(signal), waveform_name=waveform) + + optimizer.assert_not_called() + for value in result: + assert np.all(np.isnan(value)) + + +@pytest.mark.parametrize("missing", ["signal", "time", "baseline"]) +def test_low_rate_fit_requires_finite_noise_samples(missing): + """Missing baseline information cannot establish a significant signal.""" + time = np.arange(-6.0, 6.0, 0.5) + signal = np.zeros(time.size) + signal[-1] = 10.0 + if missing == "signal": + signal[:] = np.nan + elif missing == "time": + time[:] = np.nan + else: + signal[:6] = np.nan + + with mock.patch("imap_processing.idex.idex_l2a.curve_fit") as optimizer: + result = estimate_dust_mass( + xr.DataArray(time), xr.DataArray(signal), waveform_name="Target_High" + ) + + optimizer.assert_not_called() + for value in result: + assert np.all(np.isnan(value)) + + +@pytest.mark.parametrize("waveform", ["Target_High", "Target_Low", "Ion_Grid"]) +def test_low_rate_signal_gate_is_per_event(waveform): + """Vectorized fits skip weak events and keep an independent strong event.""" + time = xr.DataArray(np.arange(-6.0, 6.0, 0.5), dims="sample") + signal = np.zeros((2, time.size)) + signal[:, :6] = [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0] + signal[0, -1] = 4.0 + signal[1, -1] = -5.0 if waveform == "Ion_Grid" else 5.0 + with mock.patch( + "imap_processing.idex.idex_l2a.curve_fit", + return_value=(np.array([0.0, 0.0, 2.0, 0.371, 37.1]), None), + ) as optimizer: + results = xr.apply_ufunc( + estimate_dust_mass, + time, + xr.DataArray(signal, dims=("epoch", "sample")), + input_core_dims=[["sample"], ["sample"]], + output_core_dims=[["parameter"], [], [], [], ["sample"]], + output_dtypes=[np.float64] * 5, + vectorize=True, + kwargs={"waveform_name": waveform}, + ) + + optimizer.assert_called_once() + for values in results: + assert np.all(np.isnan(values.isel(epoch=0))) + assert np.all(np.isfinite(values.isel(epoch=1))) + + +def test_l2a_signal_gate_preserves_independent_channels(): + """Skipped fits propagate NaNs while other channels retain their products.""" + time = np.tile(np.arange(-6.0, 6.0, 0.5), (2, 1)) + dimensions = ("epoch", "time_low_sample_rate_index") + high_dimensions = ("epoch", "time_high_sample_rate_index") + dataset = xr.Dataset( + { + "time_low_sample_rate": (dimensions, time), + "time_high_sample_rate": (high_dimensions, time - 1.0), + "TOF_High": (high_dimensions, np.tile([-1.0, 1.0], (2, 12))), + "science_event_flag": ("epoch", [1, 1]), + "tof_high_saturation_flag": ("epoch", [0, 0]), + "tof_high_baseline": ("epoch", [0.0, 0.0]), + "tof_high_rms": ("epoch", [1.0, 1.0]), + }, + coords={"epoch": [0, 1]}, + attrs={"Logical_source": "imap_idex_l1b_sci-10days"}, + ) + strong_events = {"Target_High": [0], "Target_Low": [1], "Ion_Grid": [0, 1]} + for waveform, events in strong_events.items(): + signal = np.zeros_like(time) + signal[:, :6] = [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0] + signal[events, -1] = 10.0 + if waveform == "Ion_Grid": + signal[0, -1] = -10.0 + dataset[waveform] = xr.DataArray(signal, dims=dimensions) + dataset[f"{waveform.lower()}_baseline"] = ("epoch", [0.0, 0.0]) + dataset[f"{waveform.lower()}_saturation_flag"] = ("epoch", [0, 0]) + + t_rise_params = np.array([1.27, -0.2, -2.1, -0.37, 5.3, 13.3, 13.3, 0.28]) + rise_time = 10 ** log_smooth_powerlaw(1.0, t_rise_params[0], t_rise_params[1:]) + + def fit_waveform(function, time, signal, p0, **kwargs): + parameters = np.asarray(p0) + parameters[3] = rise_time + return parameters, None + + with ( + mock.patch( + "imap_processing.idex.idex_l2a.load_calibration_files", + return_value=( + t_rise_params, + 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(time)), + ), + 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.curve_fit", side_effect=fit_waveform + ) as optimizer, + ): + result = idex_l2a(dataset, {}) + + assert optimizer.call_count == 4 + for waveform, events in strong_events.items(): + for product in ( + "fit_parameters", + "impact_charge", + "chi_squared", + "reduced_chi_squared", + "fit_results", + "velocity_estimate", + "dust_mass_estimate", + ): + values = result[f"{waveform.lower()}_{product}"] + assert np.all(np.isfinite(values.isel(epoch=events))) + weak_events = [event for event in range(2) if event not in events] + assert np.all(np.isnan(values.isel(epoch=weak_events)))