From 2f10bf80ffcf00ac665cb5b37eb3728b8e8bdcfa Mon Sep 17 00:00:00 2001 From: alex Date: Tue, 18 Aug 2026 10:47:36 -0600 Subject: [PATCH 1/5] Add IDEX event classification flags --- .../config/imap_idex_l1a_variable_attrs.yaml | 32 ++ .../config/imap_idex_l1b_variable_attrs.yaml | 32 ++ imap_processing/idex/idex_event_flags.py | 342 ++++++++++++++++++ imap_processing/idex/idex_l1a.py | 18 +- imap_processing/idex/idex_l1b.py | 9 + .../tests/idex/test_idex_event_flags.py | 88 +++++ imap_processing/tests/idex/test_idex_l0.py | 2 +- imap_processing/tests/idex/test_idex_l1b.py | 11 + run_local_idex_chain.sh | 158 ++++++++ 9 files changed, 690 insertions(+), 2 deletions(-) create mode 100644 imap_processing/idex/idex_event_flags.py create mode 100644 imap_processing/tests/idex/test_idex_event_flags.py create mode 100755 run_local_idex_chain.sh 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 4aa3ec40f1..9d6c60b903 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,30 @@ 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 + 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 276b45715e..dea6e9b8b5 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,30 @@ 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 + trigger_mode_lg: <<: *string_base CATDESC: Low Gain Trigger Mode. diff --git a/imap_processing/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py new file mode 100644 index 0000000000..05ad9ca05f --- /dev/null +++ b/imap_processing/idex/idex_event_flags.py @@ -0,0 +1,342 @@ +"""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", +) + +_TOF_MAX_DN = 1023.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, +) -> dict[str, int]: + """Return mutually exclusive event-type flags and the Dust Hit flag. + + Trigger classification follows the event-labeling logic used by the IDEX + quicklook. Dust detection is performed on the raw 10-bit TOF waveforms and + 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. + + 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 EVENT_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 + return flags + + +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 60e59a80f4..28533ab6fb 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 EVENT_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,20 @@ 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])), + ) + for name in EVENT_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 +1077,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 c3f9a3db2f..1ebe9a359f 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 EVENT_FLAG_NAMES from imap_processing.idex.idex_utils import get_idex_attrs, setup_dataset from imap_processing.spice.geometry import ( SpiceBody, @@ -258,6 +259,13 @@ 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 EVENT_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 +275,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/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py new file mode 100644 index 0000000000..1ddb0185d3 --- /dev/null +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -0,0 +1,88 @@ +"""Test IDEX event classification and Dust Hit flags.""" + +import numpy as np + +from imap_processing.idex.idex_event_flags import ( + EVENT_FLAG_NAMES, + classify_event_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_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 425168c74c..6ee9422247 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) == 113 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 2210d923e6..26187fbe37 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 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 EVENT_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/run_local_idex_chain.sh b/run_local_idex_chain.sh new file mode 100755 index 0000000000..e8bcde2ec7 --- /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" From 3a2b3afbbd1373164173c4efdafdf0e91afd1c07 Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 19 Aug 2026 11:32:00 -0600 Subject: [PATCH 2/5] Added saturation flags for all 6 waveform channels. Added flags to indicate instrument status (science, pulser, noise capture). Added flag to indicate if the event was an actual dust event. Resolved ion grid velocity and mass estimates to use the QI/Qt ratio and Qt respectively. Added logic to NaN fields based on flag states. --- .../config/imap_idex_l1a_variable_attrs.yaml | 36 ++++ .../config/imap_idex_l1b_variable_attrs.yaml | 36 ++++ .../config/imap_idex_l2a_variable_attrs.yaml | 62 ++++++ imap_processing/idex/idex_event_flags.py | 88 ++++++++- imap_processing/idex/idex_l1a.py | 7 +- imap_processing/idex/idex_l1b.py | 6 +- imap_processing/idex/idex_l2a.py | 180 +++++++++++++++--- .../tests/idex/test_idex_event_flags.py | 25 +++ imap_processing/tests/idex/test_idex_l0.py | 2 +- imap_processing/tests/idex/test_idex_l1b.py | 4 +- imap_processing/tests/idex/test_idex_l2a.py | 75 ++++++++ 11 files changed, 487 insertions(+), 34 deletions(-) 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 9d6c60b903..a63b43784c 100644 --- a/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml @@ -123,6 +123,42 @@ dust_hit_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 dea6e9b8b5..5393a2d17d 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -88,6 +88,42 @@ dust_hit_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 f2f39c9611..5d0713efdc 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/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py index 05ad9ca05f..d00290e821 100644 --- a/imap_processing/idex/idex_event_flags.py +++ b/imap_processing/idex/idex_event_flags.py @@ -14,7 +14,19 @@ "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 @@ -37,6 +49,9 @@ def classify_event_flags( 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. @@ -52,6 +67,8 @@ def classify_event_flags( 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 ------- @@ -83,15 +100,82 @@ def classify_event_flags( else: event_type = "science_event_flag" - flags = {name: 0 for name in EVENT_FLAG_NAMES} + 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, @@ -267,7 +351,7 @@ def _is_saturated(value: float) -> bool: bool Whether the sample is saturated. """ - return bool(value > _SATURATION_FRACTION * _TOF_MAX_DN) + return bool(value >= _SATURATION_FRACTION * _TOF_MAX_DN) def _fwhm(corrected: np.ndarray, times: np.ndarray, peak_index: int) -> float: diff --git a/imap_processing/idex/idex_l1a.py b/imap_processing/idex/idex_l1a.py index 28533ab6fb..6bef36dc82 100644 --- a/imap_processing/idex/idex_l1a.py +++ b/imap_processing/idex/idex_l1a.py @@ -33,7 +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 EVENT_FLAG_NAMES, classify_event_flags +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 ( @@ -1036,8 +1036,11 @@ def process(self) -> Dataset | None: 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 EVENT_FLAG_NAMES: + for name in ALL_FLAG_NAMES: data_vars[name] = xr.DataArray( data=[event_flags[name]], dims=("epoch"), diff --git a/imap_processing/idex/idex_l1b.py b/imap_processing/idex/idex_l1b.py index 1ebe9a359f..59edbc1edb 100644 --- a/imap_processing/idex/idex_l1b.py +++ b/imap_processing/idex/idex_l1b.py @@ -30,7 +30,7 @@ IDEX_EVENT_REFERENCE_FRAME, ConversionFactors, ) -from imap_processing.idex.idex_event_flags import EVENT_FLAG_NAMES +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, @@ -260,9 +260,7 @@ def idex_l1b_science(l1a_dataset: xr.Dataset) -> xr.Dataset: l1a_dataset["idx__txhdrtrigid"].data, idex_attrs ) event_flags = { - name: l1a_dataset[name].copy() - for name in EVENT_FLAG_NAMES - if name in l1a_dataset + 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) diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 220c9d15ef..9749e8886e 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/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py index 1ddb0185d3..d98ef0f522 100644 --- a/imap_processing/tests/idex/test_idex_event_flags.py +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -4,7 +4,9 @@ from imap_processing.idex.idex_event_flags import ( EVENT_FLAG_NAMES, + SATURATION_FLAG_NAMES, classify_event_flags, + classify_saturation_flags, ) @@ -69,6 +71,29 @@ def test_event_type_classification() -> None: 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) diff --git a/imap_processing/tests/idex/test_idex_l0.py b/imap_processing/tests/idex/test_idex_l0.py index 6ee9422247..abe4ce243c 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) == 113 + 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 26187fbe37..143d16eb82 100644 --- a/imap_processing/tests/idex/test_idex_l1b.py +++ b/imap_processing/tests/idex/test_idex_l1b.py @@ -11,7 +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 EVENT_FLAG_NAMES +from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES, EVENT_FLAG_NAMES from imap_processing.idex.idex_l1b import ( TRIGGER_LABELS, EventMessage, @@ -62,7 +62,7 @@ def test_l1b_logical_source(l1b_dataset: xr.Dataset): 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 EVENT_FLAG_NAMES: + for flag_name in ALL_FLAG_NAMES: assert flag_name in l1b_dataset assert set(np.unique(l1b_dataset[flag_name])) <= {0, 1} diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 5699a67323..147a52ac97 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." @@ -443,6 +472,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.""" From 204e3c914567bd3d7d5bc9792c3161034c9bc5c6 Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 19 Aug 2026 12:25:33 -0600 Subject: [PATCH 3/5] fixed expected nan condition in l2a for ion grid values --- imap_processing/tests/idex/test_idex_l2a.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 147a52ac97..016c567aad 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -210,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", From 76a42754f448aaa92bb3e94f23862c8c1e35bc4f Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 19 Aug 2026 15:25:03 -0600 Subject: [PATCH 4/5] added event category masking to l2b/c, left NaN block in place. --- .../config/imap_idex_l2b_variable_attrs.yaml | 6 +- imap_processing/idex/idex_l2b.py | 64 ++++++++++++++-- imap_processing/tests/idex/test_idex_l2b.py | 74 +++++++++++++++++-- 3 files changed, 129 insertions(+), 15 deletions(-) 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 8fbe02b31f..548149743b 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_l2b.py b/imap_processing/idex/idex_l2b.py index 42c3daa5a6..1c3f316485 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_l2b.py b/imap_processing/tests/idex/test_idex_l2b.py index e8164f6366..570bd720a5 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 From fec1ed1801c37405738fc6bd45fa0983f93da855 Mon Sep 17 00:00:00 2001 From: alex Date: Thu, 20 Aug 2026 09:57:01 -0600 Subject: [PATCH 5/5] fixed description in idex_event_flags.py --- imap_processing/idex/idex_event_flags.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/imap_processing/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py index d00290e821..6825dd4c2a 100644 --- a/imap_processing/idex/idex_event_flags.py +++ b/imap_processing/idex/idex_event_flags.py @@ -55,9 +55,21 @@ def classify_event_flags( ) -> dict[str, int]: """Return mutually exclusive event-type flags and the Dust Hit flag. - Trigger classification follows the event-labeling logic used by the IDEX - quicklook. Dust detection is performed on the raw 10-bit TOF waveforms and - uses lower-gain waveforms only to measure saturated high-gain peaks. + 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 ----------