From 2f10bf80ffcf00ac665cb5b37eb3728b8e8bdcfa Mon Sep 17 00:00:00 2001 From: alex Date: Tue, 18 Aug 2026 10:47:36 -0600 Subject: [PATCH 01/15] 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 02/15] 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 03/15] 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 04/15] 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 05/15] 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 ---------- From 118d9075388d4f45c6d5125b9ff08e13528de99d Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Aug 2026 11:17:08 -0600 Subject: [PATCH 06/15] addressed round of comments --- .../config/imap_idex_l1a_variable_attrs.yaml | 4 +- .../config/imap_idex_l1b_variable_attrs.yaml | 4 +- imap_processing/idex/idex_constants.py | 9 ++- imap_processing/idex/idex_event_flags.py | 51 ++++++++++++++-- imap_processing/idex/idex_l2a.py | 45 +++++++++----- imap_processing/idex/idex_l2b.py | 4 +- .../tests/idex/test_idex_event_flags.py | 60 +++++++++++++++++++ imap_processing/tests/idex/test_idex_l2a.py | 21 +++++++ 8 files changed, 175 insertions(+), 23 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 a63b43784c..b57fc72957 100644 --- a/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml @@ -76,8 +76,10 @@ trigger_base: &trigger_base event_flag_base: &event_flag_base <<: *trigger_base CATDESC: "Boolean event classification flag." - FILLVAL: 0 + CDF_DATA_TYPE: CDF_UINT1 + FILLVAL: 255 FORMAT: I1 + UNITS: " " VALIDMAX: 1 VALIDMIN: 0 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 5393a2d17d..8c0eb99e85 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -50,8 +50,10 @@ trigger_base: &trigger_base event_flag_base: &event_flag_base <<: *trigger_base CATDESC: "Boolean event classification flag." - FILLVAL: 0 + CDF_DATA_TYPE: CDF_UINT1 + FILLVAL: 255 FORMAT: I1 + UNITS: " " VALIDMAX: 1 VALIDMIN: 0 diff --git a/imap_processing/idex/idex_constants.py b/imap_processing/idex/idex_constants.py index fa4583d8dd..02e953d989 100644 --- a/imap_processing/idex/idex_constants.py +++ b/imap_processing/idex/idex_constants.py @@ -73,7 +73,7 @@ class IdexConstants: class ConversionFactors(float, Enum): - """Conversion factor values (DN to picocoulombs) for each of the six waveforms.""" + """Conversion factor values (DN to picocoulombs) for each waveform.""" TOF_High = 2.89e-4 TOF_Low = 5.14e-1 @@ -83,6 +83,13 @@ class ConversionFactors(float, Enum): Ion_Grid = 7.46e-4 +# Constants for the ion-grid V(R) relation, where R is the ion-grid to target +# charge ratio and velocity is returned in km/s. +ION_GRID_VELOCITY_SCALE = 55.0 +ION_GRID_VELOCITY_EXPONENT = -3.2 +ION_GRID_VELOCITY_OFFSET = 1.5 + + SPICE_ARRAYS = [ "ephemeris_position_x", "ephemeris_position_y", diff --git a/imap_processing/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py index 6825dd4c2a..a0ad8432db 100644 --- a/imap_processing/idex/idex_event_flags.py +++ b/imap_processing/idex/idex_event_flags.py @@ -338,15 +338,35 @@ def _saturation_aware_width( float Full width at half maximum in microseconds, or NaN if unavailable. """ + if ( + high.size != mid.size + or high.size != low.size + or high.size != times.size + or high_corrected.size != high.size + or peak_index < 0 + or peak_index >= high.size + ): + return np.nan + peak_time = float(times[peak_index]) + if not np.isfinite(peak_time): + return np.nan if not _is_saturated(float(high[peak_index])): return _fwhm(high_corrected, times, peak_index) for waveform in (mid, low): - index = int(np.nanargmin(np.abs(times - peak_time))) + finite_times = np.isfinite(times) + if not np.any(finite_times): + continue + distances = np.where(finite_times, np.abs(times - peak_time), np.inf) + index = int(np.argmin(distances)) + sample = float(waveform[index]) + if not np.isfinite(sample) or _is_saturated(sample): + continue corrected, _ = _baseline_corrected(waveform, times) - if not _is_saturated(float(waveform[index])): - return _fwhm(corrected, times, index) + width = _fwhm(corrected, times, index) + if np.isfinite(width): + return width return np.nan @@ -383,6 +403,15 @@ def _fwhm(corrected: np.ndarray, times: np.ndarray, peak_index: int) -> float: float Full width at half maximum in microseconds, or NaN if unavailable. """ + if ( + corrected.ndim != 1 + or times.ndim != 1 + or corrected.size != times.size + or peak_index < 0 + or peak_index >= corrected.size + ): + return np.nan + peak_height = float(corrected[peak_index]) if not np.isfinite(peak_height) or peak_height <= 0.0: return np.nan @@ -397,7 +426,21 @@ def _fwhm(corrected: np.ndarray, times: np.ndarray, peak_index: int) -> float: and corrected[right] >= half_height ): right += 1 - if peak_index in (left, right): + left_bracketed = ( + left > 0 + and np.isfinite(corrected[left - 1]) + and corrected[left - 1] >= half_height + and np.isfinite(corrected[left]) + and corrected[left] < half_height + ) + right_bracketed = ( + right < corrected.size - 1 + and np.isfinite(corrected[right - 1]) + and corrected[right - 1] >= half_height + and np.isfinite(corrected[right]) + and corrected[right] < half_height + ) + if not left_bracketed or not right_bracketed: return np.nan left_time = _crossing_time(corrected, times, left, left + 1, half_height) right_time = _crossing_time(corrected, times, right - 1, right, half_height) diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 9749e8886e..dc5d6457da 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -326,17 +326,7 @@ def idex_l2a(l1b_dataset: xr.Dataset, ancillary_files: dict) -> xr.Dataset: ), ) - # 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) + _mask_saturated_derived_estimates(l2a_dataset) _mask_non_science_derived_estimates(l2a_dataset) @@ -397,6 +387,31 @@ def _mask_non_science_derived_estimates(l2a_dataset: xr.Dataset) -> None: l2a_dataset[estimate_name] = l2a_dataset[estimate_name].where(science_event) +def _mask_saturated_derived_estimates(l2a_dataset: xr.Dataset) -> None: + """Mask fitted charges and derived estimates for saturated waveforms. + + Fit parameters remain available for diagnostics. Impact charge, velocity, + and mass estimates are not scientifically valid when their source waveform + is saturated. + + Parameters + ---------- + l2a_dataset : xarray.Dataset + L2A dataset containing waveform saturation flags and derived estimates. + """ + for waveform_name in ("target_low", "target_high", "ion_grid"): + saturation_flag = f"{waveform_name}_saturation_flag" + if saturation_flag not in l2a_dataset: + continue + invalid = l2a_dataset[saturation_flag] == 1 + for estimate_name in ( + f"{waveform_name}_impact_charge", + f"{waveform_name}_velocity_estimate", + f"{waveform_name}_dust_mass_estimate", + ): + l2a_dataset[estimate_name] = l2a_dataset[estimate_name].where(~invalid) + + def calculate_velocity_and_mass( sig_amp: float, t_rise: float, t_rise_params: np.ndarray, yield_params: np.ndarray ) -> tuple[float, float]: @@ -481,10 +496,12 @@ def calculate_ion_grid_velocity_and_mass( 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 + velocity_estimate = ( + idex_constants.ION_GRID_VELOCITY_SCALE + * charge_ratio**idex_constants.ION_GRID_VELOCITY_EXPONENT + + idex_constants.ION_GRID_VELOCITY_OFFSET + ) mass_estimate = calculate_mass_from_velocity( target_charge, velocity_estimate, yield_params ) diff --git a/imap_processing/idex/idex_l2b.py b/imap_processing/idex/idex_l2b.py index 1c3f316485..435d86154d 100644 --- a/imap_processing/idex/idex_l2b.py +++ b/imap_processing/idex/idex_l2b.py @@ -92,8 +92,8 @@ def _select_target_values( 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. + are saturated receive NaN values. Finite-fit fallback selection is deferred + to a later L2B update. Parameters ---------- diff --git a/imap_processing/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py index d98ef0f522..ae0d6b6980 100644 --- a/imap_processing/tests/idex/test_idex_event_flags.py +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -1,13 +1,18 @@ """Test IDEX event classification and Dust Hit flags.""" import numpy as np +import pytest from imap_processing.idex.idex_event_flags import ( + ALL_FLAG_NAMES, EVENT_FLAG_NAMES, SATURATION_FLAG_NAMES, + _fwhm, + _saturation_aware_width, classify_event_flags, classify_saturation_flags, ) +from imap_processing.idex.idex_utils import get_idex_attrs def _waveforms(saturated: bool = False) -> tuple[np.ndarray, ...]: @@ -94,6 +99,20 @@ def test_saturation_flags_use_channel_bit_depth_and_95_percent_limit() -> None: assert flags["ion_grid_saturation_flag"] == 0 +def test_event_flags_have_istp_integer_flag_attributes() -> None: + """All event and saturation flags use explicit ISTP-compatible metadata.""" + for level in ("l1a", "l1b"): + attributes = get_idex_attrs(level) + for flag_name in ALL_FLAG_NAMES: + flag_attrs = attributes.get_variable_attributes(flag_name) + assert flag_attrs["CDF_DATA_TYPE"] == "CDF_UINT1" + assert flag_attrs["FILLVAL"] == 255 + assert flag_attrs["FORMAT"] == "I1" + assert flag_attrs["UNITS"] == " " + assert flag_attrs["VALIDMIN"] == 0 + assert flag_attrs["VALIDMAX"] == 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) @@ -111,3 +130,44 @@ def test_dust_hit_is_not_set_for_non_science_events() -> None: ) assert flags["pulser_flag"] == 1 assert flags["dust_hit_flag"] == 0 + + +def test_saturation_aware_width_falls_through_invalid_mid_gain() -> None: + """A non-finite Mid sample falls through to a usable Low waveform.""" + times = np.arange(9, dtype=float) + low = np.array([0.0, 0.0, 1.0, 3.0, 5.0, 3.0, 1.0, 0.0, 0.0]) + high = low.copy() + high[4] = 1023.0 + mid = low.copy() + mid[4] = np.nan + + width = _saturation_aware_width(4, high, mid, low, times, high - high[0]) + + assert width == pytest.approx(2.5) + + +def test_fwhm_rejects_truncated_boundary_peaks() -> None: + """A missing half-height crossing at either edge is not measurable.""" + times = np.arange(4, dtype=float) + assert np.isnan(_fwhm(np.array([2.0, 2.0, 1.0, 0.0]), times, 1)) + assert np.isnan(_fwhm(np.array([0.0, 1.0, 2.0, 2.0]), times, 2)) + + +def test_saturation_aware_width_rejects_invalid_peak_inputs() -> None: + """Invalid waveform lengths or times produce no measurable peak.""" + values = np.ones(4) + assert np.isnan( + _saturation_aware_width( + 1, values, values, values[:-1], np.arange(4, dtype=float), values + ) + ) + assert np.isnan( + _saturation_aware_width( + 1, + values, + values, + values, + np.full(4, np.nan), + values, + ) + ) diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 016c567aad..93ab4e920b 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -15,6 +15,7 @@ from imap_processing.idex.idex_l2a import ( BaselineNoiseTime, _mask_non_science_derived_estimates, + _mask_saturated_derived_estimates, analyze_peaks, butter_lowpass_filter, calculate_ion_grid_velocity_and_mass, @@ -53,6 +54,26 @@ def test_non_science_derived_estimates_are_nan() -> None: np.testing.assert_array_equal(dataset[estimate].values, [1.0, np.nan]) +def test_saturated_waveform_derived_values_are_nan() -> None: + """Saturation masks fitted charge, velocity, and mass estimates.""" + data = { + f"{waveform}_{estimate}": xr.DataArray([1.0, 2.0], dims="epoch") + for waveform in ("target_low", "target_high", "ion_grid") + for estimate in ("impact_charge", "velocity_estimate", "dust_mass_estimate") + } + dataset = xr.Dataset(data, coords={"epoch": [0, 1]}) + for waveform in ("target_low", "target_high", "ion_grid"): + dataset[f"{waveform}_saturation_flag"] = xr.DataArray([0, 1], dims="epoch") + + _mask_saturated_derived_estimates(dataset) + + for waveform in ("target_low", "target_high", "ion_grid"): + for estimate in ("impact_charge", "velocity_estimate", "dust_mass_estimate"): + np.testing.assert_array_equal( + dataset[f"{waveform}_{estimate}"].values, [1.0, np.nan] + ) + + def 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 From 69b9c0827b495d0fae6ed03aa3641ace9a375245 Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Aug 2026 11:17:49 -0600 Subject: [PATCH 07/15] addressed round of comments --- run_local_idex_chain.sh | 158 ---------------------------------------- 1 file changed, 158 deletions(-) delete mode 100755 run_local_idex_chain.sh diff --git a/run_local_idex_chain.sh b/run_local_idex_chain.sh deleted file mode 100755 index e8bcde2ec7..0000000000 --- a/run_local_idex_chain.sh +++ /dev/null @@ -1,158 +0,0 @@ -#!/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 1c80d140f348b8706080d728598ab84ea2ed997f Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Aug 2026 11:24:51 -0600 Subject: [PATCH 08/15] added log and mandatory state for saturation flags --- imap_processing/idex/idex_l2a.py | 7 ++++++- imap_processing/tests/idex/test_idex_l2a.py | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index dc5d6457da..4360105e45 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -376,6 +376,9 @@ def _mask_non_science_derived_estimates(l2a_dataset: xr.Dataset) -> None: L2A dataset containing the Science Event flag and derived estimates. """ if "science_event_flag" not in l2a_dataset: + logger.debug( + "Science event flag is not present; skipping non-science estimate masking." + ) return science_event = l2a_dataset["science_event_flag"] == 1 @@ -402,7 +405,9 @@ def _mask_saturated_derived_estimates(l2a_dataset: xr.Dataset) -> None: for waveform_name in ("target_low", "target_high", "ion_grid"): saturation_flag = f"{waveform_name}_saturation_flag" if saturation_flag not in l2a_dataset: - continue + message = f"Required L2A saturation flag is missing: {saturation_flag}" + logger.error(message) + raise KeyError(message) invalid = l2a_dataset[saturation_flag] == 1 for estimate_name in ( f"{waveform_name}_impact_charge", diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 93ab4e920b..8f5891f4ef 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -74,6 +74,27 @@ def test_saturated_waveform_derived_values_are_nan() -> None: ) +def test_saturated_derived_values_require_all_saturation_flags(caplog) -> None: + """Missing saturation metadata raises an informative error.""" + dataset = xr.Dataset( + { + "target_low_impact_charge": xr.DataArray([1.0], dims="epoch"), + "target_low_velocity_estimate": xr.DataArray([1.0], dims="epoch"), + "target_low_dust_mass_estimate": xr.DataArray([1.0], dims="epoch"), + "target_low_saturation_flag": xr.DataArray([0], dims="epoch"), + }, + coords={"epoch": [0]}, + ) + + with ( + caplog.at_level("ERROR"), + pytest.raises(KeyError, match="target_high_saturation_flag"), + ): + _mask_saturated_derived_estimates(dataset) + + assert "Required L2A saturation flag is missing" in caplog.text + + 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 From 8566c766e891e05bbe51631da9bacf47eec5e1ef Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Aug 2026 11:42:39 -0600 Subject: [PATCH 09/15] re-run checks From 0f16d72a7baf8ef11a7e5642ec928e83ac33e3fe Mon Sep 17 00:00:00 2001 From: Alex Date: Mon, 24 Aug 2026 11:55:04 -0600 Subject: [PATCH 10/15] fixed tests --- imap_processing/idex/idex_event_flags.py | 8 ++++---- imap_processing/tests/idex/test_idex_event_flags.py | 1 - 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/imap_processing/idex/idex_event_flags.py b/imap_processing/idex/idex_event_flags.py index a0ad8432db..600395e44f 100644 --- a/imap_processing/idex/idex_event_flags.py +++ b/imap_processing/idex/idex_event_flags.py @@ -281,7 +281,7 @@ def _baseline_corrected( 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) + baseline_mask = finite & (times < first_time + _BASELINE_WINDOW_US) samples = values[baseline_mask] if samples.size == 0: samples = values[finite] @@ -427,11 +427,11 @@ def _fwhm(corrected: np.ndarray, times: np.ndarray, peak_index: int) -> float: ): right += 1 left_bracketed = ( - left > 0 - and np.isfinite(corrected[left - 1]) - and corrected[left - 1] >= half_height + left < corrected.size - 1 and np.isfinite(corrected[left]) and corrected[left] < half_height + and np.isfinite(corrected[left + 1]) + and corrected[left + 1] >= half_height ) right_bracketed = ( right < corrected.size - 1 diff --git a/imap_processing/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py index ae0d6b6980..2db241b9cc 100644 --- a/imap_processing/tests/idex/test_idex_event_flags.py +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -105,7 +105,6 @@ def test_event_flags_have_istp_integer_flag_attributes() -> None: attributes = get_idex_attrs(level) for flag_name in ALL_FLAG_NAMES: flag_attrs = attributes.get_variable_attributes(flag_name) - assert flag_attrs["CDF_DATA_TYPE"] == "CDF_UINT1" assert flag_attrs["FILLVAL"] == 255 assert flag_attrs["FORMAT"] == "I1" assert flag_attrs["UNITS"] == " " From 74db9772a1c050db849bae673b03610231d1bcf0 Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Aug 2026 13:36:55 -0600 Subject: [PATCH 11/15] corrected fill values --- .../cdf/config/imap_idex_l2a_variable_attrs.yaml | 4 +++- imap_processing/idex/idex_l2a.py | 6 +++--- imap_processing/tests/idex/test_idex_l2a.py | 9 ++------- 3 files changed, 8 insertions(+), 11 deletions(-) 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 5d0713efdc..0a0e7344ec 100644 --- a/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l2a_variable_attrs.yaml @@ -3,9 +3,11 @@ double_fillval: &double_fillval -1.0E+31 event_flag_base: &event_flag_base CATDESC: Boolean event classification flag. + CDF_DATA_TYPE: CDF_UINT1 DEPEND_0: epoch DICT_KEY: SPASE>Support>SupportQuantity:QualityFlag - FILLVAL: 0 + DISPLAY_TYPE: time_series + FILLVAL: 255 FORMAT: I1 LABLAXIS: " " UNITS: " " diff --git a/imap_processing/idex/idex_l2a.py b/imap_processing/idex/idex_l2a.py index 4360105e45..f7da5df21d 100644 --- a/imap_processing/idex/idex_l2a.py +++ b/imap_processing/idex/idex_l2a.py @@ -405,9 +405,9 @@ def _mask_saturated_derived_estimates(l2a_dataset: xr.Dataset) -> None: for waveform_name in ("target_low", "target_high", "ion_grid"): saturation_flag = f"{waveform_name}_saturation_flag" if saturation_flag not in l2a_dataset: - message = f"Required L2A saturation flag is missing: {saturation_flag}" - logger.error(message) - raise KeyError(message) + raise KeyError( + f"Required L2A saturation flag is missing: {saturation_flag}" + ) invalid = l2a_dataset[saturation_flag] == 1 for estimate_name in ( f"{waveform_name}_impact_charge", diff --git a/imap_processing/tests/idex/test_idex_l2a.py b/imap_processing/tests/idex/test_idex_l2a.py index 8f5891f4ef..a6aaee7017 100644 --- a/imap_processing/tests/idex/test_idex_l2a.py +++ b/imap_processing/tests/idex/test_idex_l2a.py @@ -74,7 +74,7 @@ def test_saturated_waveform_derived_values_are_nan() -> None: ) -def test_saturated_derived_values_require_all_saturation_flags(caplog) -> None: +def test_saturated_derived_values_require_all_saturation_flags() -> None: """Missing saturation metadata raises an informative error.""" dataset = xr.Dataset( { @@ -86,14 +86,9 @@ def test_saturated_derived_values_require_all_saturation_flags(caplog) -> None: coords={"epoch": [0]}, ) - with ( - caplog.at_level("ERROR"), - pytest.raises(KeyError, match="target_high_saturation_flag"), - ): + with pytest.raises(KeyError, match="target_high_saturation_flag"): _mask_saturated_derived_estimates(dataset) - assert "Required L2A saturation flag is missing" in caplog.text - def mock_microphonics_noise(time: np.ndarray) -> np.ndarray: """Function to mock signal noise (linear and sine wave) due to microphonics.""" From c433242440233301f1ac11e337c5fab03b7ab4b4 Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Aug 2026 14:05:14 -0600 Subject: [PATCH 12/15] addressed data type issue --- .../cdf/config/imap_idex_l1a_variable_attrs.yaml | 6 +++++- imap_processing/idex/idex_l1a.py | 2 +- imap_processing/tests/idex/test_idex_event_flags.py | 3 +++ imap_processing/tests/idex/test_idex_l1a.py | 7 +++++++ 4 files changed, 16 insertions(+), 2 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 b57fc72957..0d1d5abb00 100644 --- a/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1a_variable_attrs.yaml @@ -74,14 +74,18 @@ trigger_base: &trigger_base VAR_TYPE: support_data event_flag_base: &event_flag_base - <<: *trigger_base CATDESC: "Boolean event classification flag." CDF_DATA_TYPE: CDF_UINT1 + DEPEND_0: epoch + DISPLAY_TYPE: no_plot + FIELDNAM: "" FILLVAL: 255 FORMAT: I1 + LABLAXIS: "" UNITS: " " VALIDMAX: 1 VALIDMIN: 0 + VAR_TYPE: support_data # <=== LABL_PTR_i Attributes ===> time_high_sample_rate_label: diff --git a/imap_processing/idex/idex_l1a.py b/imap_processing/idex/idex_l1a.py index 6bef36dc82..c5156531d8 100644 --- a/imap_processing/idex/idex_l1a.py +++ b/imap_processing/idex/idex_l1a.py @@ -1042,7 +1042,7 @@ def process(self) -> Dataset | None: ) for name in ALL_FLAG_NAMES: data_vars[name] = xr.DataArray( - data=[event_flags[name]], + data=np.asarray([event_flags[name]], dtype=np.uint8), dims=("epoch"), name=name, attrs=idex_attrs.get_variable_attributes(name), diff --git a/imap_processing/tests/idex/test_idex_event_flags.py b/imap_processing/tests/idex/test_idex_event_flags.py index 2db241b9cc..39526edec0 100644 --- a/imap_processing/tests/idex/test_idex_event_flags.py +++ b/imap_processing/tests/idex/test_idex_event_flags.py @@ -110,6 +110,9 @@ def test_event_flags_have_istp_integer_flag_attributes() -> None: assert flag_attrs["UNITS"] == " " assert flag_attrs["VALIDMIN"] == 0 assert flag_attrs["VALIDMAX"] == 1 + if level == "l1a": + assert flag_attrs["VAR_TYPE"] == "support_data" + assert flag_attrs["DISPLAY_TYPE"] == "no_plot" def test_dust_hit_requires_two_seven_sigma_peaks_and_is_saturation_aware() -> None: diff --git a/imap_processing/tests/idex/test_idex_l1a.py b/imap_processing/tests/idex/test_idex_l1a.py index 828423e1cc..4c7251315b 100644 --- a/imap_processing/tests/idex/test_idex_l1a.py +++ b/imap_processing/tests/idex/test_idex_l1a.py @@ -63,6 +63,13 @@ def test_idex_cdf_file(decom_test_data_sci: xr.Dataset): ) +def test_event_flags_use_uint8_data(decom_test_data_sci: xr.Dataset): + """L1A event and saturation flags use the declared unsigned byte type.""" + for name in decom_test_data_sci.data_vars: + if name.endswith("_flag"): + assert decom_test_data_sci[name].dtype == np.uint8 + + def test_bad_cdf_attributes(decom_test_data_sci: xr.Dataset): """Ensure an ``ISTPError`` is raised when using bad CDF attributes. From 6674c4c28e9774bbf3ebc554ba1a539f47a975f6 Mon Sep 17 00:00:00 2001 From: Alex Date: Tue, 25 Aug 2026 15:48:25 -0600 Subject: [PATCH 13/15] reverted l2b changes for cleaner merging --- imap_processing/idex/idex_l2b.py | 64 ++---------------- imap_processing/tests/idex/test_idex_l2b.py | 74 ++------------------- 2 files changed, 12 insertions(+), 126 deletions(-) diff --git a/imap_processing/idex/idex_l2b.py b/imap_processing/idex/idex_l2b.py index 435d86154d..42c3daa5a6 100644 --- a/imap_processing/idex/idex_l2b.py +++ b/imap_processing/idex/idex_l2b.py @@ -35,6 +35,7 @@ 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, @@ -85,55 +86,6 @@ 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. Finite-fit fallback selection is deferred - to a later L2B update. - - 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]: @@ -487,18 +439,16 @@ 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]) - 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) + 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] 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] - # 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) + # Convert units + mass_vals = FG_TO_KG * 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 570bd720a5..e8164f6366 100644 --- a/imap_processing/tests/idex/test_idex_l2b.py +++ b/imap_processing/tests/idex/test_idex_l2b.py @@ -8,6 +8,7 @@ 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, @@ -280,13 +281,7 @@ def test_compute_counts_by_charge_and_mass(): l2a_dataset = xr.Dataset( { "epoch": epochs, - "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_dust_mass_estimate": ((MASS_BIN_EDGES / FG_TO_KG)[:6] + 1e-5), "target_low_impact_charge": CHARGE_BIN_EDGES[:6], "spin_phase": np.full((6,), 0), "longitude": np.full(6, 5), @@ -352,15 +347,10 @@ 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-21, MASS_BIN_EDGES[-1] + 1e-21] - ), + [MASS_BIN_EDGES[0] - 1e-05, MASS_BIN_EDGES[-1] + 1e-05] + ) + / FG_TO_KG, "target_low_impact_charge": np.array( [CHARGE_BIN_EDGES[0] - 1e-05, CHARGE_BIN_EDGES[-1] + 1e-05] ), @@ -410,60 +400,6 @@ 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 eeff14b25765d3f920be73705f51f03dd0f51ad3 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 26 Aug 2026 09:50:00 -0600 Subject: [PATCH 14/15] Use upstream mA conversions for IDEX TOF channels --- .../config/imap_idex_l1b_variable_attrs.yaml | 4 ++ imap_processing/idex/idex_constants.py | 13 +++++-- imap_processing/idex/idex_l1b.py | 13 ++++--- imap_processing/tests/idex/test_idex_l1b.py | 37 ++++++++++++------- 4 files changed, 43 insertions(+), 24 deletions(-) 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 8c0eb99e85..6e5c5ebf1d 100644 --- a/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_idex_l1b_variable_attrs.yaml @@ -36,6 +36,7 @@ l1b_tof_base: &l1b_tof_base <<: *l1b_data_base DEPEND_1: time_high_sample_rate_index LABL_PTR_1: time_high_sample_rate_label + UNITS: mA l1b_target_base: &l1b_target_base <<: *l1b_data_base @@ -135,6 +136,7 @@ trigger_level_lg: <<: *trigger_base CATDESC: Low Gain Trigger Level threshold. FIELDNAM: Low Gain Trigger Level + UNITS: mA trigger_mode_mg: <<: *string_base @@ -145,6 +147,7 @@ trigger_level_mg: <<: *trigger_base CATDESC: Mid Gain Trigger level threshold. FIELDNAM: Mid Gain Trigger Level + UNITS: mA trigger_mode_hg: @@ -156,6 +159,7 @@ trigger_level_hg: <<: *trigger_base CATDESC: High Trigger Level threshold. FIELDNAM: High Trigger Level + UNITS: mA trigger_origin: <<: *string_base diff --git a/imap_processing/idex/idex_constants.py b/imap_processing/idex/idex_constants.py index c928a9ab9c..bca7251abd 100644 --- a/imap_processing/idex/idex_constants.py +++ b/imap_processing/idex/idex_constants.py @@ -73,10 +73,15 @@ class IdexConstants: class ConversionFactors(float, Enum): - """Conversion factor values (DN to picocoulombs) for each waveform.""" - TOF_High = 2.89e-4 - TOF_Low = 5.14e-1 - TOF_Mid = 1.13e-2 + """Conversion factors from DN to the engineering units for each waveform. + + TOF channels are reported in milliamperes (mA); target and ion-grid channels are + reported in picocoulombs (pC). + """ + + TOF_High = 7.50e-5 + TOF_Low = 1.34e-1 + TOF_Mid = 2.93e-3 Target_Low = 1.58e1 Target_High = 1.63e-1 Ion_Grid = 7.46e-4 diff --git a/imap_processing/idex/idex_l1b.py b/imap_processing/idex/idex_l1b.py index 59edbc1edb..c59ca74e2c 100644 --- a/imap_processing/idex/idex_l1b.py +++ b/imap_processing/idex/idex_l1b.py @@ -342,7 +342,7 @@ def convert_waveforms( l1a_dataset: xr.Dataset, idex_attrs: ImapCdfAttributes ) -> dict[str, xr.DataArray]: """ - Apply transformation from raw DN to picocoulombs (pC) for each of the six waveforms. + Apply the channel-specific transformation from raw DN to engineering units. Parameters ---------- @@ -355,17 +355,18 @@ def convert_waveforms( ------- waveforms_converted : dict A dictionary where the keys are the waveform array names and the values are - xr.DataArrays representing the waveforms transformed into picocoulombs. + xr.DataArrays representing the converted waveforms. TOF channels are in mA; + target and ion-grid channels are in pC. """ - waveforms_pc = {} + waveforms_converted = {} for var in ConversionFactors: - waveforms_pc[var.name] = l1a_dataset[var.name] * var.value - waveforms_pc[var.name].attrs = idex_attrs.get_variable_attributes( + waveforms_converted[var.name] = l1a_dataset[var.name] * var.value + waveforms_converted[var.name].attrs = idex_attrs.get_variable_attributes( var.name.lower() ) - return waveforms_pc + return waveforms_converted def get_trigger_mode_and_level( diff --git a/imap_processing/tests/idex/test_idex_l1b.py b/imap_processing/tests/idex/test_idex_l1b.py index 143d16eb82..303b07fff7 100644 --- a/imap_processing/tests/idex/test_idex_l1b.py +++ b/imap_processing/tests/idex/test_idex_l1b.py @@ -10,7 +10,7 @@ from imap_processing import imap_module_directory from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import write_cdf -from imap_processing.idex.idex_constants import DT_BLOCK +from imap_processing.idex.idex_constants import DT_BLOCK, ConversionFactors from imap_processing.idex.idex_event_flags import ALL_FLAG_NAMES, EVENT_FLAG_NAMES from imap_processing.idex.idex_l1b import ( TRIGGER_LABELS, @@ -103,18 +103,15 @@ def test_idex_waveform_units(l1b_dataset: xr.Dataset): assert l1b_dataset[var_name].attrs["UNITS"] == row["unit"] # Check waveform units - waveform_var_names = [ - "TOF_High", - "TOF_Low", - "TOF_Mid", - "Ion_Grid", - "Target_Low", - "Target_High", - ] + for var_name in ("TOF_High", "TOF_Low", "TOF_Mid"): + assert l1b_dataset[var_name].attrs["UNITS"] == "mA" - for var_name in waveform_var_names: + for var_name in ("Ion_Grid", "Target_Low", "Target_High"): assert l1b_dataset[var_name].attrs["UNITS"] == "pC" + for var_name in ("trigger_level_lg", "trigger_level_mg", "trigger_level_hg"): + assert l1b_dataset[var_name].attrs["UNITS"] == "mA" + def test_unpack_instrument_settings(): """ @@ -180,9 +177,9 @@ def test_get_trigger_settings_success(decom_test_data_sci): expected_modes_mg[0] = "MGThreshold" expected_levels_lg = np.full(n_epochs, np.nan) expected_levels_hg = expected_levels_lg.copy() - expected_levels_hg[1:] = 0.16762 + expected_levels_hg[1:] = 580.0 * 7.50e-5 expected_levels_mg = expected_levels_lg.copy() - expected_levels_mg[0] = 1023.0 * 1.13e-2 + expected_levels_mg[0] = 1023.0 * 2.93e-3 var_names = ["trigger_mode_lg", "trigger_mode_mg", "trigger_mode_hg"] expected_modes = [expected_modes_lg, expected_modes_mg, expected_modes_hg] @@ -346,6 +343,13 @@ def test_validate_l1b_idex_data_variables( event=np.arange(l1b_dataset.sizes["epoch"]) ) # Compare each corresponding variable + # The team validation file stores TOF waveforms using the legacy pC factors. + # Convert those reference arrays to the corrected L1B mA units before comparing. + legacy_tof_factors = { + "TOF L": 5.14e-1, + "TOF H": 2.89e-4, + "TOF M": 1.13e-2, + } for var in l1b_example_data.data_vars: if var not in arrays_to_skip: # Get the corresponding array name @@ -360,15 +364,20 @@ def test_validate_l1b_idex_data_variables( l1b_dataset[cdf_var] except KeyError: continue + expected_data = np.squeeze(l1b_example_data[var]) + if var in legacy_tof_factors: + expected_data = expected_data * ( + ConversionFactors[cdf_var].value / legacy_tof_factors[var] + ) if l1b_dataset[cdf_var].dtype == object: assert ( - l1b_dataset[cdf_var].data == np.squeeze(l1b_example_data[var]) + l1b_dataset[cdf_var].data == expected_data ).all(), warning else: np.testing.assert_array_almost_equal( l1b_dataset[cdf_var].data, - np.squeeze(l1b_example_data[var]), + expected_data, decimal=4, err_msg=warning, ) From 0165f7f4f555694199e346130af794764b9d7e98 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 26 Aug 2026 09:58:17 -0600 Subject: [PATCH 15/15] pre-commit fix --- imap_processing/tests/idex/test_idex_l1b.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/imap_processing/tests/idex/test_idex_l1b.py b/imap_processing/tests/idex/test_idex_l1b.py index 303b07fff7..151b6b9947 100644 --- a/imap_processing/tests/idex/test_idex_l1b.py +++ b/imap_processing/tests/idex/test_idex_l1b.py @@ -370,9 +370,7 @@ def test_validate_l1b_idex_data_variables( ConversionFactors[cdf_var].value / legacy_tof_factors[var] ) if l1b_dataset[cdf_var].dtype == object: - assert ( - l1b_dataset[cdf_var].data == expected_data - ).all(), warning + assert (l1b_dataset[cdf_var].data == expected_data).all(), warning else: np.testing.assert_array_almost_equal(