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..a63b43784 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,14 @@ trigger_base: &trigger_base VALIDMIN: 0 VAR_TYPE: support_data +event_flag_base: &event_flag_base + <<: *trigger_base + CATDESC: "Boolean event classification flag." + FILLVAL: 0 + FORMAT: I1 + VALIDMAX: 1 + VALIDMIN: 0 + # <=== LABL_PTR_i Attributes ===> time_high_sample_rate_label: CATDESC: High sample rate time steps for a dust event. @@ -91,6 +99,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 276b45715..5393a2d17 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -47,6 +47,14 @@ trigger_base: &trigger_base DISPLAY_TYPE: no_plot UNITS: " " +event_flag_base: &event_flag_base + <<: *trigger_base + CATDESC: "Boolean event classification flag." + FILLVAL: 0 + FORMAT: I1 + VALIDMAX: 1 + VALIDMIN: 0 + spice_base: &spice_base <<: *l1b_data_base DISPLAY_TYPE: time_series @@ -56,6 +64,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..5d0713efd 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,17 @@ int_fillval: &int_fillval -9223372036854775808 double_fillval: &double_fillval -1.0E+31 + +event_flag_base: &event_flag_base + CATDESC: Boolean event classification flag. + DEPEND_0: epoch + DICT_KEY: SPASE>Support>SupportQuantity:QualityFlag + FILLVAL: 0 + FORMAT: I1 + LABLAXIS: " " + UNITS: " " + VALIDMAX: 1 + VALIDMIN: 0 + VAR_TYPE: data int_maxval: &int_maxval 9223372036854775807 # <=== Base Attributes ===> @@ -141,6 +153,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_event_flags.py b/imap_processing/idex/idex_event_flags.py new file mode 100644 index 000000000..6825dd4c2 --- /dev/null +++ b/imap_processing/idex/idex_event_flags.py @@ -0,0 +1,438 @@ +"""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. + """ + peak_time = float(times[peak_index]) + if not _is_saturated(float(high[peak_index])): + return _fwhm(high_corrected, times, peak_index) + + for waveform in (mid, low): + index = int(np.nanargmin(np.abs(times - peak_time))) + corrected, _ = _baseline_corrected(waveform, times) + if not _is_saturated(float(waveform[index])): + return _fwhm(corrected, times, index) + 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. + """ + 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 + if peak_index in (left, right): + 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..6bef36dc8 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=[event_flags[name]], + 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 c3f9a3db2..59edbc1ed 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..9749e8886 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,19 @@ 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. + # Saturated channels retain their fits for diagnostics, but their derived + # velocity and mass products are not scientifically valid. + for waveform_name in ("target_low", "target_high", "ion_grid"): + saturation_flag = f"{waveform_name}_saturation_flag" + if saturation_flag in l2a_dataset: + invalid = l2a_dataset[saturation_flag] == 1 + for estimate_name in ( + f"{waveform_name}_velocity_estimate", + f"{waveform_name}_dust_mass_estimate", + ): + l2a_dataset[estimate_name] = l2a_dataset[estimate_name].where(~invalid) - # Ion Grid Fitting: - l2a_dataset["ion_grid_dust_mass_estimate"].data = np.full( - l2a_dataset["ion_grid_dust_mass_estimate"].shape, np.nan - ) - - 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 +373,30 @@ 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: + 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 calculate_velocity_and_mass( sig_amp: float, t_rise: float, t_rise_params: np.ndarray, yield_params: np.ndarray ) -> tuple[float, float]: @@ -381,12 +430,97 @@ 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 + + # SPECIAL NOTE: The constants for the V(R) equation may need to be moved + # to idex_constants.py once the calibration is finalized. + charge_ratio = ion_grid_charge / target_charge + velocity_estimate = 55.0 * charge_ratio**-3.2 + 1.5 + 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/idex/idex_l2b.py b/imap_processing/idex/idex_l2b.py index 42c3daa5a..1c3f31648 100644 --- a/imap_processing/idex/idex_l2b.py +++ b/imap_processing/idex/idex_l2b.py @@ -35,7 +35,6 @@ from imap_processing.ena_maps.ena_maps import SkyTilingType from imap_processing.ena_maps.utils.spatial_utils import AzElSkyGrid from imap_processing.idex.idex_constants import ( - FG_TO_KG, IDEX_EVENT_REFERENCE_FRAME, IDEX_SPACING_DEG, SECONDS_IN_DAY, @@ -86,6 +85,55 @@ IDEX_INT_FILLVAL = np.iinfo(np.int64).min +def _select_target_values( + l2a_dataset: xr.Dataset, event_indices: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """Select target charge and mass from the highest-gain valid channel. + + Target high is preferred whenever it is not saturated. Target low is used + only when target high is saturated. Events for which both target channels + are saturated receive NaN values; no additional validity filtering is + performed here. + + Parameters + ---------- + l2a_dataset : xarray.Dataset + IDEX L2A dataset containing target estimates and saturation flags. + event_indices : np.ndarray + Indices of events for which values should be selected. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + Selected target mass and charge values. + """ + target_high_saturated = ( + l2a_dataset["target_high_saturation_flag"].data[event_indices] == 1 + ) + target_low_saturated = ( + l2a_dataset["target_low_saturation_flag"].data[event_indices] == 1 + ) + + use_target_high = ~target_high_saturated + use_target_low = target_high_saturated & ~target_low_saturated + + mass_values = np.full(event_indices.size, np.nan) + charge_values = np.full(event_indices.size, np.nan) + mass_values[use_target_high] = l2a_dataset["target_high_dust_mass_estimate"].data[ + event_indices[use_target_high] + ] + charge_values[use_target_high] = l2a_dataset["target_high_impact_charge"].data[ + event_indices[use_target_high] + ] + mass_values[use_target_low] = l2a_dataset["target_low_dust_mass_estimate"].data[ + event_indices[use_target_low] + ] + charge_values[use_target_low] = l2a_dataset["target_low_impact_charge"].data[ + event_indices[use_target_low] + ] + return mass_values, charge_values + + def idex_l2b( l2a_datasets: list[xr.Dataset], msg_data_l1b: list[xr.Dataset] ) -> list[xr.Dataset]: @@ -439,16 +487,18 @@ def compute_counts_by_charge_and_mass( ] # Set the epoch for the current day to be the mean epoch of the day. daily_epoch[i] = np.mean(l2a_dataset["epoch"].data[current_day_indices]) - mass_vals = l2a_dataset["target_low_dust_mass_estimate"].data[ - current_day_indices - ] - charge_vals = l2a_dataset["target_low_impact_charge"].data[current_day_indices] + science_and_dust = ( + l2a_dataset["science_event_flag"].data[current_day_indices] == 1 + ) & (l2a_dataset["dust_hit_flag"].data[current_day_indices] == 1) + current_day_indices = current_day_indices[science_and_dust] + mass_vals, charge_vals = _select_target_values(l2a_dataset, current_day_indices) spin_phase_angles = l2a_dataset["spin_phase"].data[current_day_indices] # Make sure longitude values are in the range [0, 360) longitude = np.mod(l2a_dataset["longitude"].data[current_day_indices], 360) latitude = l2a_dataset["latitude"].data[current_day_indices] - # Convert units - mass_vals = FG_TO_KG * np.atleast_1d(mass_vals) + # L2A masses and the L2B bin edges are both in kg. L2A charges and the + # L2B bin edges are both in pC, so neither quantity needs conversion. + mass_vals = np.atleast_1d(mass_vals) # Bin spin phases binned_spin_phase = bin_spin_phases(spin_phase_angles) # Clip arrays to ensure that the values are within the valid range of bins. 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..d98ef0f52 --- /dev/null +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -0,0 +1,113 @@ +"""Test IDEX event classification and Dust Hit flags.""" + +import numpy as np + +from imap_processing.idex.idex_event_flags import ( + EVENT_FLAG_NAMES, + SATURATION_FLAG_NAMES, + classify_event_flags, + classify_saturation_flags, +) + + +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_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 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_l1b.py b/imap_processing/tests/idex/test_idex_l1b.py index 2210d923e..143d16eb8 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 +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..016c567aa 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -11,11 +11,15 @@ 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, 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 +34,25 @@ ) +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 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 +110,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 +210,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 +470,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.""" diff --git a/imap_processing/tests/idex/test_idex_l2b.py b/imap_processing/tests/idex/test_idex_l2b.py index e8164f636..570bd720a 100644 --- a/imap_processing/tests/idex/test_idex_l2b.py +++ b/imap_processing/tests/idex/test_idex_l2b.py @@ -8,7 +8,6 @@ from imap_processing.cdf.utils import write_cdf from imap_processing.idex.idex_constants import ( - FG_TO_KG, IDEX_SPACING_DEG, NANOSECONDS_IN_DAY, SECONDS_IN_DAY, @@ -281,7 +280,13 @@ def test_compute_counts_by_charge_and_mass(): l2a_dataset = xr.Dataset( { "epoch": epochs, - "target_low_dust_mass_estimate": ((MASS_BIN_EDGES / FG_TO_KG)[:6] + 1e-5), + "science_event_flag": np.ones(6, dtype=np.uint8), + "dust_hit_flag": np.ones(6, dtype=np.uint8), + "target_high_saturation_flag": np.ones(6, dtype=np.uint8), + "target_low_saturation_flag": np.zeros(6, dtype=np.uint8), + "target_high_dust_mass_estimate": np.full(6, np.nan), + "target_high_impact_charge": np.full(6, np.nan), + "target_low_dust_mass_estimate": MASS_BIN_EDGES[:6] + 1e-21, "target_low_impact_charge": CHARGE_BIN_EDGES[:6], "spin_phase": np.full((6,), 0), "longitude": np.full(6, 5), @@ -347,10 +352,15 @@ def test_compute_counts_by_charge_and_mass_out_of_bounds(): l2a_dataset = xr.Dataset( { "epoch": epochs, + "science_event_flag": np.ones(2, dtype=np.uint8), + "dust_hit_flag": np.ones(2, dtype=np.uint8), + "target_high_saturation_flag": np.ones(2, dtype=np.uint8), + "target_low_saturation_flag": np.zeros(2, dtype=np.uint8), + "target_high_dust_mass_estimate": np.full(2, np.nan), + "target_high_impact_charge": np.full(2, np.nan), "target_low_dust_mass_estimate": np.array( - [MASS_BIN_EDGES[0] - 1e-05, MASS_BIN_EDGES[-1] + 1e-05] - ) - / FG_TO_KG, + [MASS_BIN_EDGES[0] - 1e-21, MASS_BIN_EDGES[-1] + 1e-21] + ), "target_low_impact_charge": np.array( [CHARGE_BIN_EDGES[0] - 1e-05, CHARGE_BIN_EDGES[-1] + 1e-05] ), @@ -400,6 +410,60 @@ def test_compute_counts_by_charge_and_mass_out_of_bounds(): np.testing.assert_array_equal(mass_map, expected_map_array) +def test_compute_counts_masks_non_science_and_non_dust_events(): + """Only science dust-hit events contribute to L2B/L2C counts.""" + epochs = np.array([1, 1, 1]) * NANOSECONDS_IN_DAY + l2a_dataset = xr.Dataset( + { + "epoch": epochs, + "science_event_flag": [1, 0, 1], + "dust_hit_flag": [1, 1, 0], + "target_high_saturation_flag": [0, 0, 0], + "target_low_saturation_flag": [0, 0, 0], + "target_high_dust_mass_estimate": [MASS_BIN_EDGES[0], 0, 0], + "target_high_impact_charge": [CHARGE_BIN_EDGES[0], 0, 0], + "target_low_dust_mass_estimate": [0, 0, 0], + "target_low_impact_charge": [0, 0, 0], + "spin_phase": [0, 0, 0], + "longitude": [5, 5, 5], + "latitude": [0, 0, 0], + } + ) + results = compute_counts_by_charge_and_mass(l2a_dataset, np.array([2])) + + assert results[0].sum() == 1 + assert results[1].sum() == 1 + assert results[2].sum() == 1 + assert results[3].sum() == 1 + + +def test_compute_counts_prefers_unsaturated_target_high(): + """Target-high values are selected when target-low is also available.""" + epochs = np.array([1, 1]) * NANOSECONDS_IN_DAY + l2a_dataset = xr.Dataset( + { + "epoch": epochs, + "science_event_flag": [1, 1], + "dust_hit_flag": [1, 1], + "target_high_saturation_flag": [0, 1], + "target_low_saturation_flag": [0, 0], + "target_high_dust_mass_estimate": [MASS_BIN_EDGES[0], 0], + "target_high_impact_charge": [CHARGE_BIN_EDGES[0], 0], + "target_low_dust_mass_estimate": [MASS_BIN_EDGES[-2], MASS_BIN_EDGES[1]], + "target_low_impact_charge": [CHARGE_BIN_EDGES[-2], CHARGE_BIN_EDGES[1]], + "spin_phase": [0, 0], + "longitude": [5, 5], + "latitude": [0, 0], + } + ) + results = compute_counts_by_charge_and_mass(l2a_dataset, np.array([2])) + + assert results[0].sum() == 2 + assert results[1].sum() == 2 + assert results[0][0, 0, 0] == 1 + assert results[0][0, 1, 0] == 1 + + def test_compute_rates_by_charge_and_mass(): """Test the compute_rates_by_charge_and_mass function.""" # Mock example inputs diff --git a/run_local_idex_chain.sh b/run_local_idex_chain.sh new file mode 100755 index 000000000..e8bcde2ec --- /dev/null +++ b/run_local_idex_chain.sh @@ -0,0 +1,158 @@ +#!/usr/bin/env bash +# Run a local IDEX L1A -> L1B -> L2A processing chain. +# +# The script deliberately fails before processing if an input is not already +# present under DATA_ROOT. This prevents imap_cli from downloading a missing +# dependency while testing local products. + +set -Eeuo pipefail + +SCRIPT_DIR=$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")" && pwd) +DATA_ROOT=${IMAP_DATA_DIR:-"$SCRIPT_DIR/data/imap"} +START_DATE=${IDEX_START_DATE:-20260719} +LOCAL_MAJOR_VERSION=${IDEX_LOCAL_MAJOR_VERSION:-2} +LOCAL_MINOR_VERSION=${IDEX_LOCAL_MINOR_VERSION:-0} +WORK_DIR=${IDEX_LOCAL_WORK_DIR:-"$SCRIPT_DIR/.local_dependencies/$START_DATE"} + +L1A_TEMPLATE=${L1A_TEMPLATE:-"$DATA_ROOT/dependency/idex/l1a/2026/07/imap_idex_l1a_all-70eae46c-eb87af55_20260719_v001.0004.json"} +L1B_TEMPLATE=${L1B_TEMPLATE:-"$DATA_ROOT/dependency/idex/l1b/2026/07/imap_idex_l1b_sci-10days-3d1c847d-5eded07e_20260719_v001.0009.json"} +L2A_TEMPLATE=${L2A_TEMPLATE:-"$DATA_ROOT/dependency/idex/l2a/2026/07/imap_idex_l2a_sci-10days-11542dff-5eded07e_20260719_v001.0009.json"} + +L1A_DEPENDENCY="$WORK_DIR/l1a.json" +L1B_DEPENDENCY="$WORK_DIR/l1b.json" +L2A_DEPENDENCY="$WORK_DIR/l2a.json" +L1A_PRODUCT="$DATA_ROOT/idex/l1a/2026/07/imap_idex_l1a_sci-10days_${START_DATE}_v${LOCAL_MAJOR_VERSION}.$(printf '%04d' "$LOCAL_MINOR_VERSION").cdf" +L1B_PRODUCT="$DATA_ROOT/idex/l1b/2026/07/imap_idex_l1b_sci-10days_${START_DATE}_v${LOCAL_MAJOR_VERSION}.$(printf '%04d' "$LOCAL_MINOR_VERSION").cdf" + +die() { + echo "ERROR: $*" >&2 + exit 1 +} + +command -v imap_cli >/dev/null 2>&1 || die "imap_cli is not on PATH. Activate the IDEX environment first." +command -v python >/dev/null 2>&1 || die "python is not on PATH. Activate the IDEX environment first." + +[[ -d "$DATA_ROOT" ]] || die "IMAP data directory does not exist: $DATA_ROOT" +[[ -f "$L1A_TEMPLATE" ]] || die "Missing L1A dependency template: $L1A_TEMPLATE" +[[ -f "$L1B_TEMPLATE" ]] || die "Missing L1B dependency template: $L1B_TEMPLATE" +[[ -f "$L2A_TEMPLATE" ]] || die "Missing L2A dependency template: $L2A_TEMPLATE" + +mapfile -t TEMPLATE_FILES < <( + python - "$L1A_TEMPLATE" "$L1B_TEMPLATE" "$L2A_TEMPLATE" <<'PY' +import json +import sys + +for index, filename in enumerate(sys.argv[1:]): + with open(filename) as stream: + document = json.load(stream) + for dependency in document.get("dependency", []): + # L1B and L2A science inputs are replaced by the local predecessor. + if index > 0 and dependency.get("type") == "science": + continue + for item in dependency.get("files", []): + print(item) +PY +) + +for filename in "${TEMPLATE_FILES[@]}"; do + [[ -n "$filename" ]] || continue + if [[ -z "$(find "$DATA_ROOT" -type f -name "$filename" -print -quit)" ]]; then + die "Required local dependency is missing: $filename" + fi +done + +rm -rf -- "$WORK_DIR" +mkdir -p -- "$WORK_DIR" +cp -- "$L1A_TEMPLATE" "$WORK_DIR/template_l1a.json" +cp -- "$L1B_TEMPLATE" "$WORK_DIR/template_l1b.json" +cp -- "$L2A_TEMPLATE" "$WORK_DIR/template_l2a.json" +L1A_TEMPLATE="$WORK_DIR/template_l1a.json" +L1B_TEMPLATE="$WORK_DIR/template_l1b.json" +L2A_TEMPLATE="$WORK_DIR/template_l2a.json" + +echo "Removing generated IDEX L1A/L1B/L2A products for $START_DATE" +for level in l1a l1b l2a; do + level_dir="$DATA_ROOT/idex/$level" + [[ -d "$level_dir" ]] && find "$level_dir" -type f -name "*${START_DATE}*" -delete +done + +rewrite_dependency() { + local template=$1 + local output=$2 + local replacement=${3:-} + python - "$template" "$output" "$replacement" "$LOCAL_MAJOR_VERSION" "$LOCAL_MINOR_VERSION" <<'PY' +import json +import sys + +template, output, replacement, major, minor = sys.argv[1:] +with open(template) as stream: + document = json.load(stream) + +if replacement: + science_entries = [ + entry for entry in document.get("dependency", []) if entry.get("type") == "science" + ] + if len(science_entries) != 1: + raise SystemExit( + f"Expected exactly one science dependency in {template}; found {len(science_entries)}" + ) + science_entries[0]["files"] = [replacement] + +for versions in document.get("version", {}).values(): + versions["major_version"] = int(major) + versions["minor_version"] = int(minor) + +with open(output, "w") as stream: + json.dump(document, stream, indent=2) + stream.write("\n") +PY +} + +rewrite_dependency "$L1A_TEMPLATE" "$L1A_DEPENDENCY" + +echo "Processing L1A locally" +imap_cli --instrument idex --data-level l1a --descriptor all \ + --start-date "$START_DATE" --version v000 --dependency "$L1A_DEPENDENCY" + +[[ -f "$L1A_PRODUCT" ]] || die "Expected L1A product was not created: $L1A_PRODUCT" + +rewrite_dependency "$L1B_TEMPLATE" "$L1B_DEPENDENCY" "$(basename "$L1A_PRODUCT")" + +echo "Processing L1B locally" +imap_cli --instrument idex --data-level l1b --descriptor sci-10days \ + --start-date "$START_DATE" --version v000 --dependency "$L1B_DEPENDENCY" + +[[ -f "$L1B_PRODUCT" ]] || die "Expected L1B product was not created: $L1B_PRODUCT" + +python - "$L1A_PRODUCT" "$L1B_PRODUCT" <<'PY' +import sys + +import numpy as np +from cdflib import CDF + +flag_names = ( + "science_event_flag", + "noise_capture_flag", + "pulser_flag", + "dust_hit_flag", +) +l1a = CDF(sys.argv[1]) +l1b = CDF(sys.argv[2]) +for name in flag_names: + if name not in l1a.cdf_info().zVariables or name not in l1b.cdf_info().zVariables: + raise SystemExit(f"Missing {name} in L1A or L1B product") + if not np.array_equal(l1a.varget(name), l1b.varget(name)): + raise SystemExit(f"L1A/L1B values differ for {name}") +print("Verified all four event flags are preserved from L1A to L1B") +PY + +rewrite_dependency "$L2A_TEMPLATE" "$L2A_DEPENDENCY" "$(basename "$L1B_PRODUCT")" + +echo "Processing L2A locally" +imap_cli --instrument idex --data-level l2a --descriptor sci-10days \ + --start-date "$START_DATE" --version v000 --dependency "$L2A_DEPENDENCY" + +echo "Local IDEX processing chain completed for $START_DATE" +echo "L1A: $L1A_PRODUCT" +echo "L1B: $L1B_PRODUCT" +echo "L2A: $DATA_ROOT/idex/l2a/2026/07/imap_idex_l2a_sci-10days_${START_DATE}_v${LOCAL_MAJOR_VERSION}.$(printf '%04d' "$LOCAL_MINOR_VERSION").cdf"