diff --git a/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml b/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml index 4aa3ec40f..0d1d5abb0 100644 --- a/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml @@ -73,6 +73,20 @@ trigger_base: &trigger_base VALIDMIN: 0 VAR_TYPE: support_data +event_flag_base: &event_flag_base + CATDESC: "Boolean event classification flag." + CDF_DATA_TYPE: CDF_UINT1 + DEPEND_0: epoch + DISPLAY_TYPE: no_plot + FIELDNAM: "" + FILLVAL: 255 + FORMAT: I1 + LABLAXIS: "" + UNITS: " " + VALIDMAX: 1 + VALIDMIN: 0 + VAR_TYPE: support_data + # <=== LABL_PTR_i Attributes ===> time_high_sample_rate_label: CATDESC: High sample rate time steps for a dust event. @@ -91,6 +105,66 @@ time_low_sample_rate_label: VAR_TYPE: metadata # <=== Instrument Setting Attributes ===> +science_event_flag: + <<: *event_flag_base + CATDESC: Science event classification flag. + FIELDNAM: Science Event Flag + LABLAXIS: Science Event + +noise_capture_flag: + <<: *event_flag_base + CATDESC: Noise capture event classification flag. + FIELDNAM: Noise Capture Flag + LABLAXIS: Noise Capture + +pulser_flag: + <<: *event_flag_base + CATDESC: Pulser event classification flag. + FIELDNAM: Pulser Flag + LABLAXIS: Pulser + +dust_hit_flag: + <<: *event_flag_base + CATDESC: Dust Hit waveform classification flag. + FIELDNAM: Dust Hit Flag + LABLAXIS: Dust Hit + +tof_high_saturation_flag: + <<: *event_flag_base + CATDESC: TOF High waveform saturation flag. + FIELDNAM: TOF High Saturation Flag + LABLAXIS: TOF HG Saturation + +tof_mid_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Mid waveform saturation flag. + FIELDNAM: TOF Mid Saturation Flag + LABLAXIS: TOF MG Saturation + +tof_low_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Low waveform saturation flag. + FIELDNAM: TOF Low Saturation Flag + LABLAXIS: TOF LG Saturation + +target_high_saturation_flag: + <<: *event_flag_base + CATDESC: Target High waveform saturation flag. + FIELDNAM: Target High Saturation Flag + LABLAXIS: Target HG Saturation + +target_low_saturation_flag: + <<: *event_flag_base + CATDESC: Target Low waveform saturation flag. + FIELDNAM: Target Low Saturation Flag + LABLAXIS: Target LG Saturation + +ion_grid_saturation_flag: + <<: *event_flag_base + CATDESC: Ion Grid waveform saturation flag. + FIELDNAM: Ion Grid Saturation Flag + LABLAXIS: Ion Grid Saturation + low_sample_rate_attrs: <<: *sample_rate_base CATDESC: Low sample rate time steps for a dust event. 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 2c7d18fb7..6e5c5ebf1 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -48,6 +48,16 @@ trigger_base: &trigger_base DISPLAY_TYPE: no_plot UNITS: " " +event_flag_base: &event_flag_base + <<: *trigger_base + CATDESC: "Boolean event classification flag." + CDF_DATA_TYPE: CDF_UINT1 + FILLVAL: 255 + FORMAT: I1 + UNITS: " " + VALIDMAX: 1 + VALIDMIN: 0 + spice_base: &spice_base <<: *l1b_data_base DISPLAY_TYPE: time_series @@ -57,6 +67,66 @@ spice_base: &spice_base VAR_TYPE: data # <=== Instrument Setting Attributes ===> +science_event_flag: + <<: *event_flag_base + CATDESC: Science event classification flag. + FIELDNAM: Science Event Flag + LABLAXIS: Science Event + +noise_capture_flag: + <<: *event_flag_base + CATDESC: Noise capture event classification flag. + FIELDNAM: Noise Capture Flag + LABLAXIS: Noise Capture + +pulser_flag: + <<: *event_flag_base + CATDESC: Pulser event classification flag. + FIELDNAM: Pulser Flag + LABLAXIS: Pulser + +dust_hit_flag: + <<: *event_flag_base + CATDESC: Dust Hit waveform classification flag. + FIELDNAM: Dust Hit Flag + LABLAXIS: Dust Hit + +tof_high_saturation_flag: + <<: *event_flag_base + CATDESC: TOF High waveform saturation flag. + FIELDNAM: TOF High Saturation Flag + LABLAXIS: TOF HG Saturation + +tof_mid_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Mid waveform saturation flag. + FIELDNAM: TOF Mid Saturation Flag + LABLAXIS: TOF MG Saturation + +tof_low_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Low waveform saturation flag. + FIELDNAM: TOF Low Saturation Flag + LABLAXIS: TOF LG Saturation + +target_high_saturation_flag: + <<: *event_flag_base + CATDESC: Target High waveform saturation flag. + FIELDNAM: Target High Saturation Flag + LABLAXIS: Target HG Saturation + +target_low_saturation_flag: + <<: *event_flag_base + CATDESC: Target Low waveform saturation flag. + FIELDNAM: Target Low Saturation Flag + LABLAXIS: Target LG Saturation + +ion_grid_saturation_flag: + <<: *event_flag_base + CATDESC: Ion Grid waveform saturation flag. + FIELDNAM: Ion Grid Saturation Flag + LABLAXIS: Ion Grid Saturation + trigger_mode_lg: <<: *string_base CATDESC: Low Gain Trigger Mode. diff --git a/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml b/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml index f2f39c961..0a0e7344e 100644 --- a/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml @@ -1,5 +1,19 @@ int_fillval: &int_fillval -9223372036854775808 double_fillval: &double_fillval -1.0E+31 + +event_flag_base: &event_flag_base + CATDESC: Boolean event classification flag. + CDF_DATA_TYPE: CDF_UINT1 + DEPEND_0: epoch + DICT_KEY: SPASE>Support>SupportQuantity:QualityFlag + DISPLAY_TYPE: time_series + FILLVAL: 255 + FORMAT: I1 + LABLAXIS: " " + UNITS: " " + VALIDMAX: 1 + VALIDMIN: 0 + VAR_TYPE: data int_maxval: &int_maxval 9223372036854775807 # <=== Base Attributes ===> @@ -141,6 +155,56 @@ mass_index: # <=== Data Attributes ===> +science_event_flag: + <<: *event_flag_base + CATDESC: Science event classification flag. + FIELDNAM: Science Event Flag + +noise_capture_flag: + <<: *event_flag_base + CATDESC: Noise capture event classification flag. + FIELDNAM: Noise Capture Flag + +pulser_flag: + <<: *event_flag_base + CATDESC: Pulser event classification flag. + FIELDNAM: Pulser Flag + +dust_hit_flag: + <<: *event_flag_base + CATDESC: Dust Hit waveform classification flag. + FIELDNAM: Dust Hit Flag + +tof_high_saturation_flag: + <<: *event_flag_base + CATDESC: TOF High waveform saturation flag. + FIELDNAM: TOF High Saturation Flag + +tof_mid_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Mid waveform saturation flag. + FIELDNAM: TOF Mid Saturation Flag + +tof_low_saturation_flag: + <<: *event_flag_base + CATDESC: TOF Low waveform saturation flag. + FIELDNAM: TOF Low Saturation Flag + +target_high_saturation_flag: + <<: *event_flag_base + CATDESC: Target High waveform saturation flag. + FIELDNAM: Target High Saturation Flag + +target_low_saturation_flag: + <<: *event_flag_base + CATDESC: Target Low waveform saturation flag. + FIELDNAM: Target Low Saturation Flag + +ion_grid_saturation_flag: + <<: *event_flag_base + CATDESC: Ion Grid waveform saturation flag. + FIELDNAM: Ion Grid Saturation Flag + tof_snr: CATDESC: Signal-to-noise ratio of the TOF High waveform. DEPEND_0: epoch diff --git a/imap_processing/cdf/config/imap_idex_l2b_variable_attrs.yaml b/imap_processing/cdf/config/imap_idex_l2b_variable_attrs.yaml index 8fbe02b31..548149743 100644 --- a/imap_processing/cdf/config/imap_idex_l2b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l2b_variable_attrs.yaml @@ -12,10 +12,10 @@ mass_labels: VAR_TYPE: metadata charge_labels: - CATDESC: Labels for Impact Charge (fC) + CATDESC: Labels for Impact Charge (pC) DEPEND_1: impact_charge DICT_KEY: SPASE>Support>SupportQuantity:Other - FIELDNAM: Impact Charge (fC) + FIELDNAM: Impact Charge (pC) FORMAT: A8 VAR_TYPE: metadata @@ -66,7 +66,7 @@ impact_charge: LABLAXIS: Impact Charge LABL_PTR_1: charge_labels SCALETYP: log - UNITS: fC + UNITS: pC VALIDMAX: 1.00e04 VALIDMIN: 0.0 VAR_TYPE: support_data diff --git a/imap_processing/idex/idex_constants.py b/imap_processing/idex/idex_constants.py index 5ced74c30..bca7251ab 100644 --- a/imap_processing/idex/idex_constants.py +++ b/imap_processing/idex/idex_constants.py @@ -87,6 +87,13 @@ class ConversionFactors(float, Enum): Ion_Grid = 7.46e-4 +# Constants for the ion-grid V(R) relation, where R is the ion-grid to target +# charge ratio and velocity is returned in km/s. +ION_GRID_VELOCITY_SCALE = 55.0 +ION_GRID_VELOCITY_EXPONENT = -3.2 +ION_GRID_VELOCITY_OFFSET = 1.5 + + SPICE_ARRAYS = [ "ephemeris_position_x", "ephemeris_position_y", diff --git a/imap_processing/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py new file mode 100644 index 000000000..600395e44 --- /dev/null +++ b/imap_processing/idex/idex_event_flags.py @@ -0,0 +1,481 @@ +"""Classify IDEX science events and identify dust-like TOF waveforms.""" + +from __future__ import annotations + +from collections.abc import Mapping + +import numpy as np +from scipy.signal import find_peaks + +EVENT_FLAG_NAMES = ( + "science_event_flag", + "noise_capture_flag", + "pulser_flag", + "dust_hit_flag", +) + +SATURATION_FLAG_NAMES = ( + "tof_high_saturation_flag", + "tof_mid_saturation_flag", + "tof_low_saturation_flag", + "target_high_saturation_flag", + "target_low_saturation_flag", + "ion_grid_saturation_flag", +) + +ALL_FLAG_NAMES = EVENT_FLAG_NAMES + SATURATION_FLAG_NAMES + +_TOF_MAX_DN = 1023.0 +_LOW_RATE_MAX_DN = 4095.0 +_SATURATION_FRACTION = 0.95 +_PULSER_THRESHOLD_DN = 1000 +_BASELINE_WINDOW_US = 3.0 +_PEAK_THRESHOLD_SIGMA = 7.0 +_MIN_PEAK_WIDTH_US = 0.020 +_MIN_PEAK_COUNT = 2 +_MIN_PEAK_DISTANCE_US = 0.030 + +_TRIGGER_CHANNELS = { + 0: "TOF H", + 1: "TOF L", + 2: "TOF M", + 3: "Target H", +} + + +def classify_event_flags( + telemetry: Mapping[str, int], + tof_high: np.ndarray, + tof_mid: np.ndarray, + tof_low: np.ndarray, + time_high_sample_rate: np.ndarray, + target_high: np.ndarray | None = None, + target_low: np.ndarray | None = None, + ion_grid: np.ndarray | None = None, +) -> dict[str, int]: + """Return mutually exclusive event-type flags and the Dust Hit flag. + + The instrument state is assigned from the event trigger telemetry using + these conditions, in order: + + * ``noise_capture_flag`` is set when no trigger channels are active, or + when a software/external trigger is present and the only active channel + is TOF High. + * ``pulser_flag`` is set when TOF High is the only active channel, the TOF + High trigger mode is ``1``, and its trigger threshold is 1000 DN. + * ``science_event_flag`` is set for all remaining events. + + ``dust_hit_flag`` is set only for science events when the raw TOF waveform + contains at least two peaks that exceed seven baseline-noise standard + deviations and have a full width at half maximum of at least 20 ns. Dust + detection uses lower-gain waveforms only to measure saturated high-gain + peaks. + + Parameters + ---------- + telemetry : collections.abc.Mapping + Raw event trigger fields. + tof_high, tof_mid, tof_low : numpy.ndarray + Raw TOF waveforms in high, medium, and low gain. + time_high_sample_rate : numpy.ndarray + High-rate waveform times in microseconds. + target_high, target_low, ion_grid : numpy.ndarray or None + Raw low-rate waveforms used to calculate saturation flags. + + Returns + ------- + dict[str, int] + Event flags with values of zero or one. + """ + trigger_id = int(telemetry.get("idx__txhdrtrigid", 0)) + active_channels = { + channel for bit, channel in _TRIGGER_CHANNELS.items() if trigger_id & (1 << bit) + } + for gain, channel in (("hg", "TOF H"), ("mg", "TOF M"), ("lg", "TOF L")): + if int(telemetry.get(f"idx__txhdr{gain}trigmode", 0)) != 0: + active_channels.add(channel) + + has_software_or_external_trigger = bool(trigger_id & ((1 << 4) | (1 << 5))) + hg_mode = int(telemetry.get("idx__txhdrhgtrigmode", 0)) + hg_threshold = (int(telemetry.get("idx__txhdrhgtrigctrl1", 0)) >> 22) & 0x3FF + + if not active_channels or ( + has_software_or_external_trigger and active_channels <= {"TOF H"} + ): + event_type = "noise_capture_flag" + elif ( + active_channels == {"TOF H"} + and hg_mode == 1 + and hg_threshold == _PULSER_THRESHOLD_DN + ): + event_type = "pulser_flag" + else: + event_type = "science_event_flag" + + flags = {name: 0 for name in ALL_FLAG_NAMES} + flags[event_type] = 1 + if event_type == "science_event_flag" and _has_dust_hit( + tof_high, tof_mid, tof_low, time_high_sample_rate + ): + flags["dust_hit_flag"] = 1 + flags.update( + classify_saturation_flags( + tof_high, tof_mid, tof_low, target_high, target_low, ion_grid + ) + ) + return flags + + +def classify_saturation_flags( + tof_high: np.ndarray, + tof_mid: np.ndarray, + tof_low: np.ndarray, + target_high: np.ndarray | None, + target_low: np.ndarray | None, + ion_grid: np.ndarray | None, +) -> dict[str, int]: + """Return saturation flags for the six raw waveform channels. + + The low-rate channels are optional to keep the event-classification API + compatible with callers that only have the TOF waveforms. + + Parameters + ---------- + tof_high, tof_mid, tof_low : numpy.ndarray + Raw 10-bit TOF waveforms. + target_high, target_low, ion_grid : numpy.ndarray or None + Raw 12-bit low-rate waveforms. + + Returns + ------- + dict[str, int] + One zero-or-one saturation flag for each waveform channel. + """ + waveforms = { + "tof_high_saturation_flag": (tof_high, _TOF_MAX_DN), + "tof_mid_saturation_flag": (tof_mid, _TOF_MAX_DN), + "tof_low_saturation_flag": (tof_low, _TOF_MAX_DN), + "target_high_saturation_flag": (target_high, _LOW_RATE_MAX_DN), + "target_low_saturation_flag": (target_low, _LOW_RATE_MAX_DN), + "ion_grid_saturation_flag": (ion_grid, _LOW_RATE_MAX_DN), + } + return { + name: int(values is not None and _waveform_is_saturated(values, maximum)) + for name, (values, maximum) in waveforms.items() + } + + +def _waveform_is_saturated(values: np.ndarray, maximum_dn: float) -> bool: + """Return whether any finite waveform sample reaches the 95% limit. + + Parameters + ---------- + values : numpy.ndarray + Waveform samples in DN. + maximum_dn : float + Maximum representable DN for the channel. + + Returns + ------- + bool + Whether any finite sample reaches the saturation threshold. + """ + values_array = np.asarray(values, dtype=float) + finite_values = values_array[np.isfinite(values_array)] + return bool( + finite_values.size + and np.any(finite_values >= _SATURATION_FRACTION * maximum_dn) + ) + + +def _has_dust_hit( + tof_high: np.ndarray, + tof_mid: np.ndarray, + tof_low: np.ndarray, + time_high_sample_rate: np.ndarray, +) -> bool: + """Return whether TOF High contains two qualifying peaks. + + Parameters + ---------- + tof_high, tof_mid, tof_low : numpy.ndarray + Raw TOF waveforms in high, medium, and low gain. + time_high_sample_rate : numpy.ndarray + High-rate waveform times in microseconds. + + Returns + ------- + bool + Whether at least two peaks meet the sigma and FWHM requirements. + """ + # Candidate peaks are always located on High, then measured at lower gain + # when saturation prevents a reliable High-gain FWHM. + high = _as_1d_array(tof_high) + mid = _as_1d_array(tof_mid) + low = _as_1d_array(tof_low) + times = _as_1d_array(time_high_sample_rate) + length = min(high.size, mid.size, low.size, times.size) + if length == 0: + return False + high, mid, low, times = (array[:length] for array in (high, mid, low, times)) + + high_corrected, high_sigma = _baseline_corrected(high, times) + if not np.isfinite(high_sigma) or high_sigma <= 0.0: + return False + finite = np.isfinite(high_corrected) & np.isfinite(times) + if not np.any(finite): + return False + dt_us = _sample_spacing_us(times) + distance = max(1, round(_MIN_PEAK_DISTANCE_US / dt_us)) if dt_us > 0 else 1 + search = np.where(finite, high_corrected, -np.inf) + peaks, _ = find_peaks( + search, + height=_PEAK_THRESHOLD_SIGMA * high_sigma, + distance=distance, + ) + + qualifying_peaks = 0 + for peak_index in peaks: + width_us = _saturation_aware_width( + peak_index, high, mid, low, times, high_corrected + ) + if np.isfinite(width_us) and width_us >= _MIN_PEAK_WIDTH_US: + qualifying_peaks += 1 + return qualifying_peaks >= _MIN_PEAK_COUNT + + +def _as_1d_array(values: np.ndarray) -> np.ndarray: + """Convert an event waveform or time coordinate to one dimension. + + Parameters + ---------- + values : numpy.ndarray + Input waveform or time coordinate. + + Returns + ------- + numpy.ndarray + One-dimensional floating-point array. + """ + return np.asarray(values, dtype=float).reshape(-1) + + +def _baseline_corrected( + values: np.ndarray, times: np.ndarray +) -> tuple[np.ndarray, float]: + """Subtract the baseline and estimate its robust standard deviation. + + Parameters + ---------- + values : numpy.ndarray + Waveform samples. + times : numpy.ndarray + Sample times in microseconds. + + Returns + ------- + tuple[numpy.ndarray, float] + Baseline-corrected samples and estimated noise standard deviation. + """ + finite = np.isfinite(values) & np.isfinite(times) + if not np.any(finite): + return np.full(values.shape, np.nan), np.nan + first_time = float(times[finite][0]) + baseline_mask = finite & (times < first_time + _BASELINE_WINDOW_US) + samples = values[baseline_mask] + if samples.size == 0: + samples = values[finite] + baseline = float(np.nanmedian(samples)) + deviations = samples - baseline + sigma = 1.4826 * float(np.nanmedian(np.abs(deviations))) + if not np.isfinite(sigma) or sigma <= 0.0: + sigma = float(np.nanstd(samples)) + return values - baseline, sigma + + +def _sample_spacing_us(times: np.ndarray) -> float: + """Return the median finite sample spacing in microseconds. + + Parameters + ---------- + times : numpy.ndarray + Sample times in microseconds. + + Returns + ------- + float + Median sample spacing, or NaN when fewer than two samples are finite. + """ + finite_times = times[np.isfinite(times)] + if finite_times.size < 2: + return np.nan + return float(np.nanmedian(np.abs(np.diff(finite_times)))) + + +def _saturation_aware_width( + peak_index: int, + high: np.ndarray, + mid: np.ndarray, + low: np.ndarray, + times: np.ndarray, + high_corrected: np.ndarray, +) -> float: + """Measure a saturated peak width using the first usable gain. + + Parameters + ---------- + peak_index : int + High-gain peak index. + high, mid, low : numpy.ndarray + Raw TOF waveforms for the three gains. + times : numpy.ndarray + Sample times in microseconds. + high_corrected : numpy.ndarray + Baseline-corrected high-gain waveform. + + Returns + ------- + float + Full width at half maximum in microseconds, or NaN if unavailable. + """ + if ( + high.size != mid.size + or high.size != low.size + or high.size != times.size + or high_corrected.size != high.size + or peak_index < 0 + or peak_index >= high.size + ): + return np.nan + + peak_time = float(times[peak_index]) + if not np.isfinite(peak_time): + return np.nan + if not _is_saturated(float(high[peak_index])): + return _fwhm(high_corrected, times, peak_index) + + for waveform in (mid, low): + finite_times = np.isfinite(times) + if not np.any(finite_times): + continue + distances = np.where(finite_times, np.abs(times - peak_time), np.inf) + index = int(np.argmin(distances)) + sample = float(waveform[index]) + if not np.isfinite(sample) or _is_saturated(sample): + continue + corrected, _ = _baseline_corrected(waveform, times) + width = _fwhm(corrected, times, index) + if np.isfinite(width): + return width + return np.nan + + +def _is_saturated(value: float) -> bool: + """Return whether a TOF sample exceeds the 95 percent limit. + + Parameters + ---------- + value : float + Raw TOF sample in DN. + + Returns + ------- + bool + Whether the sample is saturated. + """ + return bool(value >= _SATURATION_FRACTION * _TOF_MAX_DN) + + +def _fwhm(corrected: np.ndarray, times: np.ndarray, peak_index: int) -> float: + """Measure a peak's full width at half maximum. + + Parameters + ---------- + corrected : numpy.ndarray + Baseline-corrected waveform. + times : numpy.ndarray + Sample times in microseconds. + peak_index : int + Index of the peak maximum. + + Returns + ------- + float + Full width at half maximum in microseconds, or NaN if unavailable. + """ + if ( + corrected.ndim != 1 + or times.ndim != 1 + or corrected.size != times.size + or peak_index < 0 + or peak_index >= corrected.size + ): + return np.nan + + peak_height = float(corrected[peak_index]) + if not np.isfinite(peak_height) or peak_height <= 0.0: + return np.nan + half_height = peak_height / 2.0 + left = peak_index + while left > 0 and np.isfinite(corrected[left]) and corrected[left] >= half_height: + left -= 1 + right = peak_index + while ( + right < corrected.size - 1 + and np.isfinite(corrected[right]) + and corrected[right] >= half_height + ): + right += 1 + left_bracketed = ( + left < corrected.size - 1 + and np.isfinite(corrected[left]) + and corrected[left] < half_height + and np.isfinite(corrected[left + 1]) + and corrected[left + 1] >= half_height + ) + right_bracketed = ( + right < corrected.size - 1 + and np.isfinite(corrected[right - 1]) + and corrected[right - 1] >= half_height + and np.isfinite(corrected[right]) + and corrected[right] < half_height + ) + if not left_bracketed or not right_bracketed: + return np.nan + left_time = _crossing_time(corrected, times, left, left + 1, half_height) + right_time = _crossing_time(corrected, times, right - 1, right, half_height) + if not np.isfinite(left_time) or not np.isfinite(right_time): + return np.nan + return abs(right_time - left_time) + + +def _crossing_time( + values: np.ndarray, times: np.ndarray, low: int, high: int, target: float +) -> float: + """Linearly interpolate a waveform crossing time. + + Parameters + ---------- + values : numpy.ndarray + Waveform values. + times : numpy.ndarray + Sample times in microseconds. + low : int + Index on the lower side of the crossing. + high : int + Index on the upper side of the crossing. + target : float + Crossing value. + + Returns + ------- + float + Interpolated crossing time, or NaN for invalid samples. + """ + y0, y1 = values[low], values[high] + t0, t1 = times[low], times[high] + if not all(np.isfinite(value) for value in (y0, y1, t0, t1)): + return np.nan + if y1 == y0: + return float(t0) + return float(t0 + (target - y0) * (t1 - t0) / (y1 - y0)) diff --git a/imap_processing/idex/idex_l1a.py b/imap_processing/idex/idex_l1a.py index 60e59a80f..c5156531d 100644 --- a/imap_processing/idex/idex_l1a.py +++ b/imap_processing/idex/idex_l1a.py @@ -33,6 +33,7 @@ from imap_processing.idex.decode import rice_decode from imap_processing.idex.evt_msg_decode_utils import render_event_template from imap_processing.idex.idex_constants import IDEXAPID +from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES, classify_event_flags from imap_processing.idex.idex_l0 import decom_packets from imap_processing.idex.idex_utils import get_10_day_window_end_date, get_idex_attrs from imap_processing.spice.time import ( @@ -1029,6 +1030,23 @@ def process(self) -> Dataset | None: attrs=idex_attrs.get_variable_attributes("ion_grid_attrs"), ), } + event_flags = classify_event_flags( + self.telemetry_items, + data_vars["TOF_High"].data[0], + data_vars["TOF_Mid"].data[0], + data_vars["TOF_Low"].data[0], + self._calc_high_sample_resolution(len(data_vars["TOF_High"].data[0])), + data_vars["Target_High"].data[0], + data_vars["Target_Low"].data[0], + data_vars["Ion_Grid"].data[0], + ) + for name in ALL_FLAG_NAMES: + data_vars[name] = xr.DataArray( + data=np.asarray([event_flags[name]], dtype=np.uint8), + dims=("epoch"), + name=name, + attrs=idex_attrs.get_variable_attributes(name), + ) # Determine coordinate variables coords = { "epoch": xr.DataArray( @@ -1062,7 +1080,8 @@ def process(self) -> Dataset | None: f"{name}_index": array.shape[1] for name, array in sampling_rates.items() } if any( - var.shape[1] != expected_shapes[var.dims[1]] for var in data_vars.values() + var.ndim > 1 and var.shape[1] != expected_shapes[var.dims[1]] + for var in data_vars.values() ): # The IDEX team requests that a warning be logged for incomplete events # (dropped packets) in the data, while still allowing the CDF to be created diff --git a/imap_processing/idex/idex_l1b.py b/imap_processing/idex/idex_l1b.py index 5f50edc40..c59ca74e2 100644 --- a/imap_processing/idex/idex_l1b.py +++ b/imap_processing/idex/idex_l1b.py @@ -30,6 +30,7 @@ IDEX_EVENT_REFERENCE_FRAME, 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.spice.geometry import ( SpiceBody, @@ -258,6 +259,11 @@ def idex_l1b_science(l1a_dataset: xr.Dataset) -> xr.Dataset: trigger_origin = get_trigger_origin( l1a_dataset["idx__txhdrtrigid"].data, idex_attrs ) + event_flags = { + name: l1a_dataset[name].copy() for name in ALL_FLAG_NAMES if name in l1a_dataset + } + for name, data_array in event_flags.items(): + data_array.attrs = idex_attrs.get_variable_attributes(name) # Create l1b Dataset prefixes = ["shcoarse", "shfine", "time_high_sample", "time_low_sample", "aid"] data_vars = ( @@ -267,6 +273,7 @@ def idex_l1b_science(l1a_dataset: xr.Dataset) -> xr.Dataset: | trigger_settings | spice_data | trigger_origin + | event_flags ) l1b_dataset = setup_dataset(l1a_dataset, prefixes, idex_attrs, data_vars) l1b_dataset.attrs = idex_attrs.get_global_attributes("imap_idex_l1b_sci") diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 220c9d15e..f7da5df21 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -31,6 +31,7 @@ from imap_processing import imap_module_directory 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 logger = logging.getLogger(__name__) @@ -185,6 +186,13 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: l1b_dataset, prefixes + SPICE_ARRAYS, idex_attrs, data_vars ) + # Preserve event classification and saturation flags for the event-level L2A + # products. Fits remain available for every event, including non-dust events. + for flag_name in ALL_FLAG_NAMES: + if flag_name in l1b_dataset: + l2a_dataset[flag_name] = l1b_dataset[flag_name].copy() + l2a_dataset[flag_name].attrs = idex_attrs.get_variable_attributes(flag_name) + for waveform in ["Target_Low", "Target_High", "Ion_Grid"]: # Get the dust mass estimates and fit results fit_results = xr.apply_ufunc( @@ -208,16 +216,32 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: kwargs={"waveform_name": waveform}, ) # Calculate mass and velocity estimates - velocity_mass_results = xr.apply_ufunc( - calculate_velocity_and_mass, - fit_results[1], # signal amplitude - fit_results[0].data[:, 3], # fit params - output_core_dims=[[], []], - vectorize=True, - output_dtypes=[np.float64, np.float64], - keep_attrs=True, - kwargs={"t_rise_params": t_rise_params, "yield_params": yield_params}, - ) + if waveform == "Ion_Grid": + velocity_mass_results = xr.apply_ufunc( + calculate_ion_grid_velocity_and_mass, + fit_results[1], + l2a_dataset["target_high_impact_charge"], + l2a_dataset["target_low_impact_charge"], + l2a_dataset["target_high_saturation_flag"], + l2a_dataset["target_low_saturation_flag"], + l2a_dataset["ion_grid_saturation_flag"], + output_core_dims=[[], []], + vectorize=True, + output_dtypes=[np.float64, np.float64], + keep_attrs=True, + kwargs={"yield_params": yield_params}, + ) + else: + velocity_mass_results = xr.apply_ufunc( + calculate_velocity_and_mass, + fit_results[1], # signal amplitude + fit_results[0].data[:, 3], # fit params + output_core_dims=[[], []], + vectorize=True, + output_dtypes=[np.float64, np.float64], + keep_attrs=True, + kwargs={"t_rise_params": t_rise_params, "yield_params": yield_params}, + ) waveform_name = waveform.lower() output_vars = { @@ -302,18 +326,9 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: ), ) - # We're inserting a NaN block here for the 2026 June release while the - # IDEX science team works through validating the fitting routines and - # derived values. - - # Ion Grid Fitting: - l2a_dataset["ion_grid_dust_mass_estimate"].data = np.full( - l2a_dataset["ion_grid_dust_mass_estimate"].shape, np.nan - ) + _mask_saturated_derived_estimates(l2a_dataset) - l2a_dataset["ion_grid_velocity_estimate"].data = np.full( - l2a_dataset["ion_grid_velocity_estimate"].shape, np.nan - ) + _mask_non_science_derived_estimates(l2a_dataset) # TOF / Mass-spec Fitting l2a_dataset["tof_peak_area_under_fit"].data = np.full( @@ -348,6 +363,60 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: return l2a_dataset +def _mask_non_science_derived_estimates(l2a_dataset: xr.Dataset) -> None: + """Mask velocity and mass estimates for non-science events. + + Fits and fitted charges remain available for diagnostics in all instrument + modes. The six derived velocity and mass estimates are only valid for + events classified as Science Events. + + Parameters + ---------- + l2a_dataset : xarray.Dataset + L2A dataset containing the Science Event flag and derived estimates. + """ + if "science_event_flag" not in l2a_dataset: + logger.debug( + "Science event flag is not present; skipping non-science estimate masking." + ) + return + + science_event = l2a_dataset["science_event_flag"] == 1 + for waveform_name in ("target_low", "target_high", "ion_grid"): + for estimate_name in ( + f"{waveform_name}_velocity_estimate", + f"{waveform_name}_dust_mass_estimate", + ): + l2a_dataset[estimate_name] = l2a_dataset[estimate_name].where(science_event) + + +def _mask_saturated_derived_estimates(l2a_dataset: xr.Dataset) -> None: + """Mask fitted charges and derived estimates for saturated waveforms. + + Fit parameters remain available for diagnostics. Impact charge, velocity, + and mass estimates are not scientifically valid when their source waveform + is saturated. + + Parameters + ---------- + l2a_dataset : xarray.Dataset + L2A dataset containing waveform saturation flags and derived estimates. + """ + for waveform_name in ("target_low", "target_high", "ion_grid"): + saturation_flag = f"{waveform_name}_saturation_flag" + if saturation_flag not in l2a_dataset: + raise KeyError( + f"Required L2A saturation flag is missing: {saturation_flag}" + ) + invalid = l2a_dataset[saturation_flag] == 1 + for estimate_name in ( + f"{waveform_name}_impact_charge", + f"{waveform_name}_velocity_estimate", + f"{waveform_name}_dust_mass_estimate", + ): + l2a_dataset[estimate_name] = l2a_dataset[estimate_name].where(~invalid) + + def calculate_velocity_and_mass( sig_amp: float, t_rise: float, t_rise_params: np.ndarray, yield_params: np.ndarray ) -> tuple[float, float]: @@ -381,12 +450,99 @@ def calculate_velocity_and_mass( if not np.isfinite(v_est): return np.nan, np.nan + return v_est, calculate_mass_from_velocity(sig_amp, v_est, yield_params) + + +def calculate_ion_grid_velocity_and_mass( + ion_grid_charge: float, + target_high_charge: float, + target_low_charge: float, + target_high_saturated: int, + target_low_saturated: int, + ion_grid_saturated: int, + yield_params: np.ndarray, +) -> tuple[float, float]: + """Estimate Ion Grid velocity and mass from the best target charge. + + Target High is preferred when it is unsaturated. Target Low is used only + when Target High is saturated or has no finite fitted charge. If both + target channels are saturated, or Ion Grid is saturated, both estimates + are invalid. + + Parameters + ---------- + ion_grid_charge : float + Fitted Ion Grid impact charge in pC. + target_high_charge, target_low_charge : float + Fitted target impact charges in pC. + target_high_saturated, target_low_saturated, ion_grid_saturated : int + Saturation flags for the corresponding channels. + yield_params : numpy.ndarray + Charge-yield calibration parameters. + + Returns + ------- + tuple[float, float] + Ion Grid velocity in km/s and mass in kg. + """ + if int(ion_grid_saturated) == 1: + return np.nan, np.nan + + target_charge = np.nan + if int(target_high_saturated) == 0 and np.isfinite(target_high_charge): + target_charge = float(target_high_charge) + elif int(target_low_saturated) == 0 and np.isfinite(target_low_charge): + target_charge = float(target_low_charge) + + if not (np.isfinite(ion_grid_charge) and np.isfinite(target_charge)): + return np.nan, np.nan + target_charge = abs(target_charge) + ion_grid_charge = abs(float(ion_grid_charge)) + if target_charge <= 0.0 or ion_grid_charge <= 0.0: + return np.nan, np.nan + + charge_ratio = ion_grid_charge / target_charge + velocity_estimate = ( + idex_constants.ION_GRID_VELOCITY_SCALE + * charge_ratio**idex_constants.ION_GRID_VELOCITY_EXPONENT + + idex_constants.ION_GRID_VELOCITY_OFFSET + ) + mass_estimate = calculate_mass_from_velocity( + target_charge, velocity_estimate, yield_params + ) + return velocity_estimate, mass_estimate + + +def calculate_mass_from_velocity( + sig_amp: float, velocity_estimate: float, yield_params: np.ndarray +) -> float: + """Calculate dust mass from fitted charge and an estimated velocity. + + Parameters + ---------- + sig_amp : float + Fitted signal amplitude in pC. + velocity_estimate : float + Estimated impact velocity in km/s. + yield_params : numpy.ndarray + Charge-yield calibration parameters. + + Returns + ------- + float + Estimated dust mass in kg. + """ + if not np.isfinite(sig_amp) or not np.isfinite(velocity_estimate): + return np.nan + log_a_y: float = float(yield_params[0]) - yield_val = 10 ** log_smooth_powerlaw(np.log10(v_est), log_a_y, yield_params[1:]) + yield_val = 10 ** log_smooth_powerlaw( + np.log10(velocity_estimate), log_a_y, yield_params[1:] + ) sig_amp_coulombs = sig_amp * idex_constants.PICOCOULOMB_TO_COULOMB mass_est = sig_amp_coulombs / yield_val - return v_est, mass_est + return mass_est def invert_rise_time_to_velocity(t_rise: float, t_rise_params: np.ndarray) -> float: diff --git a/imap_processing/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py new file mode 100644 index 000000000..39526edec --- /dev/null +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -0,0 +1,175 @@ +"""Test IDEX event classification and Dust Hit flags.""" + +import numpy as np +import pytest + +from imap_processing.idex.idex_event_flags import ( + ALL_FLAG_NAMES, + EVENT_FLAG_NAMES, + SATURATION_FLAG_NAMES, + _fwhm, + _saturation_aware_width, + classify_event_flags, + classify_saturation_flags, +) +from imap_processing.idex.idex_utils import get_idex_attrs + + +def _waveforms(saturated: bool = False) -> tuple[np.ndarray, ...]: + """Create two deterministic, two-peak TOF waveform events.""" + times = np.arange(2048, dtype=float) / 260.0 + baseline = 100.0 + 0.5 * np.sin(np.arange(times.size, dtype=float) / 3.0) + widths = 0.030 / 2.355 + peaks = sum( + 20.0 * np.exp(-0.5 * ((times - center) / widths) ** 2) for center in (5.0, 5.08) + ) + high = baseline + peaks + if saturated: + high = np.minimum(baseline + 60.0 * peaks, 1023.0) + medium = baseline + peaks + low = baseline + peaks + return high, medium, low, times + + +def _telemetry( + *, trigger_id: int = 0, hg_mode: int = 0, hg_threshold: int = 0 +) -> dict[str, int]: + """Return the raw trigger fields used by the classifier.""" + return { + "idx__txhdrtrigid": trigger_id, + "idx__txhdrhgtrigmode": hg_mode, + "idx__txhdrmgtrigmode": 0, + "idx__txhdrlgtrigmode": 0, + "idx__txhdrhgtrigctrl1": hg_threshold << 22, + } + + +def test_core_event_flags_are_mutually_exclusive() -> None: + """Exactly one of Science, Noise Capture, and Pulser is set.""" + waveforms = _waveforms() + cases = ( + _telemetry(), + _telemetry(trigger_id=1, hg_mode=1, hg_threshold=1000), + _telemetry(trigger_id=1 | 4, hg_mode=1), + ) + + for telemetry in cases: + flags = classify_event_flags(telemetry, *waveforms) + assert sum(flags[name] for name in EVENT_FLAG_NAMES[:3]) == 1 + assert set(flags.values()) <= {0, 1} + + +def test_event_type_classification() -> None: + """Classify noise, pulser, and science events from raw trigger fields.""" + waveforms = _waveforms() + + noise = classify_event_flags(_telemetry(), *waveforms) + assert noise["noise_capture_flag"] == 1 + + pulser = classify_event_flags( + _telemetry(trigger_id=1, hg_mode=1, hg_threshold=1000), *waveforms + ) + assert pulser["pulser_flag"] == 1 + assert pulser["dust_hit_flag"] == 0 + + science = classify_event_flags(_telemetry(trigger_id=1 | 4, hg_mode=1), *waveforms) + assert science["science_event_flag"] == 1 + + +def test_saturation_flags_use_channel_bit_depth_and_95_percent_limit() -> None: + """TOF uses 10-bit DN while low-rate channels use 12-bit DN.""" + tof = np.array([0.0, 1023.0 * 0.95]) + low_rate = np.array([0.0, 4095.0 * 0.95]) + flags = classify_saturation_flags(tof, tof, tof, low_rate, low_rate, low_rate) + + assert set(flags) == set(SATURATION_FLAG_NAMES) + assert all(value == 1 for value in flags.values()) + + flags = classify_saturation_flags( + np.array([1023.0 * 0.95 - 1.0]), + np.array([0.0]), + np.array([0.0]), + np.array([4095.0 * 0.95 - 1.0]), + None, + None, + ) + assert flags["tof_high_saturation_flag"] == 0 + assert flags["target_high_saturation_flag"] == 0 + assert flags["target_low_saturation_flag"] == 0 + assert flags["ion_grid_saturation_flag"] == 0 + + +def test_event_flags_have_istp_integer_flag_attributes() -> None: + """All event and saturation flags use explicit ISTP-compatible metadata.""" + for level in ("l1a", "l1b"): + attributes = get_idex_attrs(level) + for flag_name in ALL_FLAG_NAMES: + flag_attrs = attributes.get_variable_attributes(flag_name) + assert flag_attrs["FILLVAL"] == 255 + assert flag_attrs["FORMAT"] == "I1" + assert flag_attrs["UNITS"] == " " + assert flag_attrs["VALIDMIN"] == 0 + assert flag_attrs["VALIDMAX"] == 1 + if level == "l1a": + assert flag_attrs["VAR_TYPE"] == "support_data" + assert flag_attrs["DISPLAY_TYPE"] == "no_plot" + + +def test_dust_hit_requires_two_seven_sigma_peaks_and_is_saturation_aware() -> None: + """Two qualifying peaks set Dust Hit, including saturated High fallback.""" + saturated_waveforms = _waveforms(saturated=True) + flags = classify_event_flags( + _telemetry(trigger_id=1 | 4, hg_mode=1), *saturated_waveforms + ) + assert flags["science_event_flag"] == 1 + assert flags["dust_hit_flag"] == 1 + + +def test_dust_hit_is_not_set_for_non_science_events() -> None: + """Dust-shaped waveforms cannot turn a non-science event into Dust Hit.""" + flags = classify_event_flags( + _telemetry(trigger_id=1, hg_mode=1, hg_threshold=1000), *_waveforms() + ) + assert flags["pulser_flag"] == 1 + assert flags["dust_hit_flag"] == 0 + + +def test_saturation_aware_width_falls_through_invalid_mid_gain() -> None: + """A non-finite Mid sample falls through to a usable Low waveform.""" + times = np.arange(9, dtype=float) + low = np.array([0.0, 0.0, 1.0, 3.0, 5.0, 3.0, 1.0, 0.0, 0.0]) + high = low.copy() + high[4] = 1023.0 + mid = low.copy() + mid[4] = np.nan + + width = _saturation_aware_width(4, high, mid, low, times, high - high[0]) + + assert width == pytest.approx(2.5) + + +def test_fwhm_rejects_truncated_boundary_peaks() -> None: + """A missing half-height crossing at either edge is not measurable.""" + times = np.arange(4, dtype=float) + assert np.isnan(_fwhm(np.array([2.0, 2.0, 1.0, 0.0]), times, 1)) + assert np.isnan(_fwhm(np.array([0.0, 1.0, 2.0, 2.0]), times, 2)) + + +def test_saturation_aware_width_rejects_invalid_peak_inputs() -> None: + """Invalid waveform lengths or times produce no measurable peak.""" + values = np.ones(4) + assert np.isnan( + _saturation_aware_width( + 1, values, values, values[:-1], np.arange(4, dtype=float), values + ) + ) + assert np.isnan( + _saturation_aware_width( + 1, + values, + values, + values, + np.full(4, np.nan), + values, + ) + ) diff --git a/imap_processing/tests/idex/test_idex_l0.py b/imap_processing/tests/idex/test_idex_l0.py index 425168c74..abe4ce243 100644 --- a/imap_processing/tests/idex/test_idex_l0.py +++ b/imap_processing/tests/idex/test_idex_l0.py @@ -14,7 +14,7 @@ def test_idex_decom_length(decom_test_data_sci: xr.Dataset): decom_test_data_sci : xarray.Dataset The dataset to test with """ - assert len(decom_test_data_sci) == 109 + assert len(decom_test_data_sci) == 119 def test_idex_decom_event_num(decom_test_data_sci: xr.Dataset): diff --git a/imap_processing/tests/idex/test_idex_l1a.py b/imap_processing/tests/idex/test_idex_l1a.py index 828423e1c..4c7251315 100644 --- a/imap_processing/tests/idex/test_idex_l1a.py +++ b/imap_processing/tests/idex/test_idex_l1a.py @@ -63,6 +63,13 @@ def test_idex_cdf_file(decom_test_data_sci: xr.Dataset): ) +def test_event_flags_use_uint8_data(decom_test_data_sci: xr.Dataset): + """L1A event and saturation flags use the declared unsigned byte type.""" + for name in decom_test_data_sci.data_vars: + if name.endswith("_flag"): + assert decom_test_data_sci[name].dtype == np.uint8 + + def test_bad_cdf_attributes(decom_test_data_sci: xr.Dataset): """Ensure an ``ISTPError`` is raised when using bad CDF attributes. diff --git a/imap_processing/tests/idex/test_idex_l1b.py b/imap_processing/tests/idex/test_idex_l1b.py index debc82f69..151b6b994 100644 --- a/imap_processing/tests/idex/test_idex_l1b.py +++ b/imap_processing/tests/idex/test_idex_l1b.py @@ -11,6 +11,7 @@ from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import 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, @@ -59,6 +60,16 @@ def test_l1b_logical_source(l1b_dataset: xr.Dataset): assert l1b_dataset.attrs["Logical_source"] == expected_src +def test_event_flags_are_carried_to_l1b(l1b_dataset: xr.Dataset): + """Verify event flags are present and core flags remain exclusive.""" + for flag_name in ALL_FLAG_NAMES: + assert flag_name in l1b_dataset + assert set(np.unique(l1b_dataset[flag_name])) <= {0, 1} + + core_flags = sum(l1b_dataset[flag_name] for flag_name in EVENT_FLAG_NAMES[:3]) + np.testing.assert_array_equal(core_flags, np.ones_like(core_flags)) + + def test_idex_cdf_file(l1b_dataset: xr.Dataset): """Verify the CDF file can be created with no errors. diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 5699a6732..a6aaee701 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -11,11 +11,16 @@ from imap_processing.cdf.utils import load_cdf, write_cdf 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, butter_lowpass_filter, + calculate_ion_grid_velocity_and_mass, calculate_kappa, + calculate_mass_from_velocity, calculate_snr, calculate_velocity_and_mass, chi_square, @@ -30,6 +35,61 @@ ) +def test_non_science_derived_estimates_are_nan() -> None: + """Only Science Events retain derived velocity and mass estimates.""" + estimates = { + f"{waveform}_{estimate}": xr.DataArray([1.0, 2.0], dims="epoch") + for waveform in ("target_low", "target_high", "ion_grid") + for estimate in ("velocity_estimate", "dust_mass_estimate") + } + dataset = xr.Dataset( + estimates, + coords={"epoch": [0, 1]}, + ) + dataset["science_event_flag"] = xr.DataArray([1, 0], dims="epoch") + + _mask_non_science_derived_estimates(dataset) + + for estimate in estimates: + np.testing.assert_array_equal(dataset[estimate].values, [1.0, np.nan]) + + +def test_saturated_waveform_derived_values_are_nan() -> None: + """Saturation masks fitted charge, velocity, and mass estimates.""" + data = { + f"{waveform}_{estimate}": xr.DataArray([1.0, 2.0], dims="epoch") + for waveform in ("target_low", "target_high", "ion_grid") + for estimate in ("impact_charge", "velocity_estimate", "dust_mass_estimate") + } + dataset = xr.Dataset(data, coords={"epoch": [0, 1]}) + for waveform in ("target_low", "target_high", "ion_grid"): + dataset[f"{waveform}_saturation_flag"] = xr.DataArray([0, 1], dims="epoch") + + _mask_saturated_derived_estimates(dataset) + + for waveform in ("target_low", "target_high", "ion_grid"): + for estimate in ("impact_charge", "velocity_estimate", "dust_mass_estimate"): + np.testing.assert_array_equal( + dataset[f"{waveform}_{estimate}"].values, [1.0, np.nan] + ) + + +def test_saturated_derived_values_require_all_saturation_flags() -> None: + """Missing saturation metadata raises an informative error.""" + dataset = xr.Dataset( + { + "target_low_impact_charge": xr.DataArray([1.0], dims="epoch"), + "target_low_velocity_estimate": xr.DataArray([1.0], dims="epoch"), + "target_low_dust_mass_estimate": xr.DataArray([1.0], dims="epoch"), + "target_low_saturation_flag": xr.DataArray([0], dims="epoch"), + }, + coords={"epoch": [0]}, + ) + + with pytest.raises(KeyError, match="target_high_saturation_flag"): + _mask_saturated_derived_estimates(dataset) + + def mock_microphonics_noise(time: np.ndarray) -> np.ndarray: """Function to mock signal noise (linear and sine wave) due to microphonics.""" noise_frequency = idex_constants.TARGET_NOISE_FREQUENCY @@ -87,6 +147,12 @@ def test_l2a_logical_source_and_cdf(l2a_dataset: xr.Dataset, l1b_dataset: xr.Dat l1b_dataset[variable_name].values, ) + for flag_name in ALL_FLAG_NAMES: + assert flag_name in l2a_dataset + np.testing.assert_array_equal( + l2a_dataset[flag_name].values, l1b_dataset[flag_name].values + ) + with cdflib.CDF(file_name) as cdf_file: assert cdf_file.varattsget("mass_index")["CATDESC"] == ( "Index variable for mass-resolved TOF peak-fit products." @@ -181,8 +247,6 @@ def test_l2a_logical_source_and_cdf(l2a_dataset: xr.Dataset, l1b_dataset: xr.Dat # TODO: remove this NAN block when fitting logic is applied expected_nan_vars = [ - "ion_grid_dust_mass_estimate", - "ion_grid_velocity_estimate", "tof_peak_area_under_fit", "tof_peak_chi_square", "tof_peak_fit_parameters", @@ -443,6 +507,52 @@ def test_calculate_velocity_and_mass_at_10_km_s(): assert mass_estimate == pytest.approx(expected_mass_kg, rel=1e-12) +def test_ion_grid_velocity_uses_unsaturated_high_target_first(): + """Ion Grid velocity uses the highest-gain unsaturated target charge.""" + yield_params = np.array([0.06, 2.8, 5.9, 4.1, 13.0, 22.7, 8.2, 0.40]) + velocity, mass = calculate_ion_grid_velocity_and_mass( + 2.0, 1.0, 0.5, 0, 0, 0, yield_params + ) + + assert velocity == pytest.approx(55.0 * 2.0**-3.2 + 1.5) + assert mass == pytest.approx( + calculate_mass_from_velocity(1.0, velocity, yield_params) + ) + + +def test_ion_grid_velocity_falls_back_to_low_target_when_high_saturates(): + """Ion Grid velocity falls back to Target Low when Target High saturates.""" + yield_params = np.array([0.06, 2.8, 5.9, 4.1, 13.0, 22.7, 8.2, 0.40]) + velocity, _ = calculate_ion_grid_velocity_and_mass( + 1.0, 1.0, 0.5, 1, 0, 0, yield_params + ) + + assert velocity == pytest.approx(55.0 * 2.0**-3.2 + 1.5) + + +@pytest.mark.parametrize( + "target_high_saturated,target_low_saturated,ion_grid_saturated", + [(1, 1, 0), (0, 0, 1)], +) +def test_ion_grid_velocity_is_nan_when_required_channel_is_saturated( + target_high_saturated, target_low_saturated, ion_grid_saturated +): + """Saturation prevents publishing Ion Grid velocity and mass.""" + yield_params = np.array([0.06, 2.8, 5.9, 4.1, 13.0, 22.7, 8.2, 0.40]) + velocity, mass = calculate_ion_grid_velocity_and_mass( + 1.0, + 1.0, + 0.5, + target_high_saturated, + target_low_saturated, + ion_grid_saturated, + yield_params, + ) + + assert np.isnan(velocity) + assert np.isnan(mass) + + @pytest.mark.external_test_data def test_velocity_and_mass_estimate(ancillary_files): """Tests that the velocity and mass estimate function."""