From b35c1ea7ed4f0dc9d0bc424c4083de360a42dd70 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 2 Oct 2026 13:26:00 -0600 Subject: [PATCH 1/2] Adjust statistical filter 0 to not include invalid or partial spin-sets. This removes them from corrupting the statistics. --- imap_processing/hi/hi_goodtimes.py | 80 +++++++++-- imap_processing/tests/hi/test_hi_goodtimes.py | 128 ++++++++++++++++++ 2 files changed, 195 insertions(+), 13 deletions(-) diff --git a/imap_processing/hi/hi_goodtimes.py b/imap_processing/hi/hi_goodtimes.py index 2196ec713..8c432c2f0 100644 --- a/imap_processing/hi/hi_goodtimes.py +++ b/imap_processing/hi/hi_goodtimes.py @@ -1528,12 +1528,40 @@ def _compute_normalized_counts_per_sweep( ------- xarray.Dataset Reshaped dataset with esa_sweep as a dimension containing: - - normalized_count: normalized AB coincidence counts per sweep + - normalized_count: normalized AB coincidence counts per sweep (NaN + for sweeps with no valid 8-spin sets) - All other variables from the input dataset (first value per sweep) + + Notes + ----- + Each sweep's AB coincidence count is normalized by the number of valid + 8-spin sets (ESA steps) in that sweep, so that full and incomplete sweeps + are comparable. An 8-spin set is valid if none of its packets has an + esa_energy_step of 0 (calibration) or FILLVAL (ESA or detector voltage + mismatch). Only AB coincidences from valid 8-spin sets are counted. """ if "esa_sweep" not in l1b_de.coords: raise ValueError("Dataset must have esa_sweep coordinate") + # Identify 8-spin sets (packets sharing the same esa_step_met). A set is + # invalid if any of its packets has an invalid esa_energy_step, consistent + # with mark_bad_voltage() which culls the entire 8-spin set in that case. + esa_energy_step = l1b_de["esa_energy_step"] + fillval = esa_energy_step.attrs.get("FILLVAL", 255) + is_invalid_packet = (esa_energy_step.values == 0) | ( + esa_energy_step.values == fillval + ) + # set_mets holds one esa_step_met per 8-spin set; return_inverse gives, + # for each packet, the index of its set in set_mets. + set_mets, packet_set_idx = np.unique( + l1b_de["esa_step_met"].values, return_inverse=True + ) + # OR each packet's invalid flag into its set's entry. The unbuffered .at + # form accumulates over repeated indices (multiple packets per set). + is_invalid_set: np.ndarray = np.zeros(len(set_mets), dtype=bool) + np.logical_or.at(is_invalid_set, packet_set_idx, is_invalid_packet) + is_valid_packet = ~is_invalid_set[packet_set_idx] + # Filter to valid AB coincidences tof_ab = l1b_de["tof_ab"] coincidence_type = l1b_de["coincidence_type"] @@ -1544,18 +1572,31 @@ def _compute_normalized_counts_per_sweep( np.abs(tof_ab) <= tof_ab_limit_ns ) - # Map events to sweeps via ccsds_index -> esa_sweep + # Map events to sweeps via ccsds_index -> esa_sweep, keeping only events + # from valid 8-spin sets event_epoch_idx = ccsds_index.values event_sweep_idx = l1b_de["esa_sweep"].values[event_epoch_idx] + is_counted = is_valid_ab.values & is_valid_packet[event_epoch_idx] # Count valid AB events per sweep n_sweeps = int(l1b_de["esa_sweep"].max().values) + 1 counts_per_sweep: np.ndarray = np.zeros(n_sweeps, dtype=np.int64) - np.add.at(counts_per_sweep, event_sweep_idx[is_valid_ab.values], 1) + np.add.at(counts_per_sweep, event_sweep_idx[is_counted], 1) - # Normalize by number of unique ESA energy steps - n_unique_esa_energy_steps = len(np.unique(l1b_de["esa_energy_step"].values)) - normalized_counts = counts_per_sweep / n_unique_esa_energy_steps + # Count valid 8-spin sets per sweep (each set lies within a single sweep) + set_sweep: np.ndarray = np.zeros(len(set_mets), dtype=np.int64) + set_sweep[packet_set_idx] = l1b_de["esa_sweep"].values + n_valid_sets_per_sweep: np.ndarray = np.bincount( + set_sweep[~is_invalid_set], minlength=n_sweeps + ) + + # Normalize by the number of valid 8-spin sets in each sweep. Sweeps with + # no valid 8-spin sets get NaN. + normalized_counts = np.full(n_sweeps, np.nan) + has_valid_sets = n_valid_sets_per_sweep > 0 + normalized_counts[has_valid_sets] = ( + counts_per_sweep[has_valid_sets] / n_valid_sets_per_sweep[has_valid_sets] + ) # Remove all variables that depend on event_met dimension ds = l1b_de.drop_dims("event_met", errors="ignore") @@ -1599,8 +1640,9 @@ def mark_statistical_filter_0( the penetrating background rate has changed drastically, compromising background subtraction accuracy. For each ESA sweep across all input Pointings, it computes the normalized AB coincidence count (total count - divided by number of ESA steps). It then marks ESA sweeps in the current - Pointing where the normalized count exceeds 150% of the median. + divided by number of valid ESA steps in that sweep). It then marks ESA + sweeps in the current Pointing where the normalized count exceeds 150% of + the median. Parameters ---------- @@ -1637,9 +1679,12 @@ def mark_statistical_filter_0( median` are marked as bad. Other sweeps remain unaffected. Algorithm: - 1. For each complete ESA sweep across all Pointings, count AB coincidences - where |tof_ab| <= tof_ab_limit_ns and divide by number of ESA steps - 2. Calculate median of all normalized sweep counts + 1. For each ESA sweep across all Pointings, count AB coincidences where + |tof_ab| <= tof_ab_limit_ns from valid 8-spin sets (esa_energy_step not + 0 or FILLVAL) and divide by the number of valid 8-spin sets in that + sweep + 2. Calculate median of all normalized sweep counts, excluding sweeps with + no valid 8-spin sets (these are never marked) 3. For each sweep in current Pointing, mark all METs in that sweep as bad if normalized count > threshold_factor * median """ @@ -1685,8 +1730,16 @@ def mark_statistical_filter_0( current_ds = reshaped_datasets[current_index] - # Calculate median from all sweep counts + # Calculate median from all sweep counts, excluding sweeps with no valid + # 8-spin sets (NaN normalized count) all_counts = np.concatenate(all_normalized_counts) + all_counts = all_counts[~np.isnan(all_counts)] + if len(all_counts) == 0: + logger.warning( + "Statistical Filter 0: No sweeps with valid 8-spin sets found; " + "skipping filter." + ) + return median_count = float(np.median(all_counts)) threshold = median_count * threshold_factor @@ -1695,7 +1748,8 @@ def mark_statistical_filter_0( f"threshold={threshold:.2f} ({len(all_counts)} sweeps)" ) - # Find and mark bad sweeps in current dataset + # Find and mark bad sweeps in current dataset. Sweeps with no valid + # 8-spin sets have NaN normalized count and are never marked. bad_sweep_mask = current_ds["normalized_count"] > threshold n_bad_sweeps = int(bad_sweep_mask.sum()) diff --git a/imap_processing/tests/hi/test_hi_goodtimes.py b/imap_processing/tests/hi/test_hi_goodtimes.py index 0d4f2a55a..893c407b4 100644 --- a/imap_processing/tests/hi/test_hi_goodtimes.py +++ b/imap_processing/tests/hi/test_hi_goodtimes.py @@ -2352,6 +2352,112 @@ def test_multiple_sweeps(self): assert len(result["normalized_count"]) == 5 assert result.sizes["esa_sweep"] == 5 + @staticmethod + def _create_sets_dataset( + esa_steps: list[int], + esa_energy_steps: list[int], + events_per_set: list[int], + ) -> xr.Dataset: + """Create a DE dataset with one packet per 8-spin set and AB events.""" + n_sets = len(esa_steps) + set_met = 1000.0 + 120.0 * np.arange(n_sets) + ccsds_index = np.repeat(np.arange(n_sets), events_per_set).astype(np.uint16) + n_events = len(ccsds_index) + ds = xr.Dataset( + { + "ccsds_met": (["epoch"], set_met), + "esa_step_met": (["epoch"], set_met), + "esa_step": (["epoch"], np.array(esa_steps, dtype=np.uint8)), + "esa_energy_step": ( + ["epoch"], + np.array(esa_energy_steps, dtype=np.uint8), + {"FILLVAL": 255}, + ), + "tof_ab": (["event_met"], np.zeros(n_events, dtype=np.int32)), + "coincidence_type": ( + ["event_met"], + np.full(n_events, 12, dtype=np.uint8), + ), + "ccsds_index": (["event_met"], ccsds_index), + }, + coords={"epoch": np.arange(n_sets), "event_met": np.arange(n_events)}, + ) + return _add_sweep_indices(ds) + + def test_normalized_per_sweep_valid_sets(self): + """Test incomplete sweeps are normalized by their own number of sets.""" + sweep = list(range(1, 10)) + # Full sweep, truncated 1-set sweep, full sweep; 10 events per set + esa_steps = [*sweep, 1, *sweep] + ds = self._create_sets_dataset(esa_steps, esa_steps, [10] * len(esa_steps)) + + result = _compute_normalized_counts_per_sweep(ds, tof_ab_limit_ns=15) + + np.testing.assert_array_equal(result["normalized_count"].values, [10, 10, 10]) + + def test_excludes_calibration_and_fillval_sets(self): + """Test that esa_energy_step 0 and FILLVAL sets are not counted.""" + sweep = list(range(1, 10)) + # Sweep 0: ESA 1-9 plus ESA 10 calibration (esa_energy_step=0) + # Sweep 1: ESA 1-9 with ESA 5 having a voltage mismatch (FILLVAL) + esa_steps = [*sweep, 10, *sweep] + esa_energy_steps = [*sweep, 0, 1, 2, 3, 4, 255, 6, 7, 8, 9] + # Invalid sets have many events which must not be counted + events = [10] * 9 + [100] + [10] * 4 + [100] + [10] * 4 + ds = self._create_sets_dataset(esa_steps, esa_energy_steps, events) + + result = _compute_normalized_counts_per_sweep(ds, tof_ab_limit_ns=15) + + np.testing.assert_array_equal(result["normalized_count"].values, [10, 10]) + + def test_set_with_any_invalid_packet_is_invalid(self): + """Test that one invalid packet invalidates its whole 8-spin set.""" + # Two sweeps of ESA 1-3, two packets per 8-spin set, 10 events/packet. + # In sweep 1, the ESA 2 set has one valid and one FILLVAL packet. + esa_step = np.repeat([1, 2, 3, 1, 2, 3], 2).astype(np.uint8) + esa_step_met = np.repeat(1000.0 + 120.0 * np.arange(6), 2) + esa_energy_step = esa_step.copy() + esa_energy_step[9] = 255 + n_packets = len(esa_step) + ccsds_index = np.repeat(np.arange(n_packets), 10).astype(np.uint16) + n_events = len(ccsds_index) + ds = xr.Dataset( + { + "ccsds_met": (["epoch"], 1000.0 + 60.0 * np.arange(n_packets)), + "esa_step_met": (["epoch"], esa_step_met), + "esa_step": (["epoch"], esa_step), + "esa_energy_step": (["epoch"], esa_energy_step, {"FILLVAL": 255}), + "tof_ab": (["event_met"], np.zeros(n_events, dtype=np.int32)), + "coincidence_type": ( + ["event_met"], + np.full(n_events, 12, dtype=np.uint8), + ), + "ccsds_index": (["event_met"], ccsds_index), + }, + coords={"epoch": np.arange(n_packets), "event_met": np.arange(n_events)}, + ) + ds = _add_sweep_indices(ds) + + result = _compute_normalized_counts_per_sweep(ds, tof_ab_limit_ns=15) + + # 20 events per valid set; the partially invalid set is fully excluded + np.testing.assert_array_equal(result["normalized_count"].values, [20, 20]) + + def test_sweep_without_valid_sets_is_nan(self): + """Test that a sweep with no valid 8-spin sets gets NaN.""" + sweep = list(range(1, 10)) + esa_steps = [*sweep, 1, *sweep] + esa_energy_steps = [*sweep, 255, *sweep] + ds = self._create_sets_dataset( + esa_steps, esa_energy_steps, [10] * len(esa_steps) + ) + + result = _compute_normalized_counts_per_sweep(ds, tof_ab_limit_ns=15) + + np.testing.assert_array_equal( + result["normalized_count"].values, [10, np.nan, 10] + ) + class TestStatisticalFilter0: """Test suite for mark_statistical_filter_0() integration tests.""" @@ -2475,6 +2581,28 @@ def test_fails_anomalous_sweep(self, goodtimes_for_filter): goodtimes_for_filter["cull_flags"].values == CullCode.STAT_FILTER_0 ) + def test_invalid_sets_do_not_trigger_cull(self, goodtimes_for_filter): + """Test that high counts in invalid 8-spin sets don't cull a sweep.""" + np.random.seed(42) + l1b_de_datasets = [ + self._create_l1b_de_dataset(n_sweeps=2, events_per_met=10) for _ in range(5) + ] + # Current pointing: second sweep has 5x the events, but all of its + # 8-spin sets have a voltage mismatch (esa_energy_step = FILLVAL) + current = self._create_l1b_de_dataset(n_sweeps=2, events_per_met=50) + current["esa_energy_step"].values[9:] = 255 + l1b_de_datasets[2] = current + + mark_statistical_filter_0( + goodtimes_for_filter, l1b_de_datasets, current_index=2 + ) + + # The first sweep is also 5x but valid, so it is culled; the second + # sweep has no valid sets, so it is not + cull_flags = goodtimes_for_filter["cull_flags"].values + assert np.all(cull_flags[:9] == CullCode.STAT_FILTER_0) + assert np.all(cull_flags[9:] == CullCode.GOOD) + def test_insufficient_pointings(self, goodtimes_for_filter): """Test that fewer than min_pointings raises ValueError.""" l1b_de_datasets = [ From fb19f31abdcbacfce28b9d15f84345682a9216f3 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 2 Oct 2026 14:45:37 -0600 Subject: [PATCH 2/2] Tighten the allowed delta on the mcp_b voltage --- imap_processing/hi/utils.py | 4 +++- imap_processing/tests/hi/test_hi_l1b.py | 20 ++++++++++++++------ 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 65d5bc791..8ba19804c 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -156,7 +156,9 @@ class HiConstants: "neg_defl": 1500.0, "tof": 50.0, "mcp_f": 10.0, - "mcp_b": 50.0, + # Below the smallest MCP back gain test step (+50 V, measured as + # +49.4 V on 2026-02-09), while nominal noise is ~1 V. + "mcp_b": 25.0, "cem_f": 10.0, "cem_bk_a": 25.0, "cem_bk_b": 25.0, diff --git a/imap_processing/tests/hi/test_hi_l1b.py b/imap_processing/tests/hi/test_hi_l1b.py index 645e4b515..3800ec41d 100644 --- a/imap_processing/tests/hi/test_hi_l1b.py +++ b/imap_processing/tests/hi/test_hi_l1b.py @@ -530,8 +530,16 @@ def test_nominal_pointing_multiple_matching_segments(self, mock_any_good_de): } assert actual_hv_deltas == pytest.approx(expected_hv_deltas) + @pytest.mark.parametrize( + "field, offset", + [ + ("cem_f", 300.0), + # Smallest MCP back gain test step, as measured on 2026-02-09. + ("mcp_b", 49.4), + ], + ) @mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=True) - def test_mid_pointing_gain_test_excluded(self, mock_any_good_de): + def test_mid_pointing_gain_test_excluded(self, mock_any_good_de, field, offset): """A HVSCI segment drifting outside tolerance is excluded as a gain test.""" hk_ds = self._make_hk_ds( [ @@ -549,11 +557,11 @@ def test_mid_pointing_gain_test_excluded(self, mock_any_good_de): ], [1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010], ) - # The middle HVSCI segment (indices 4-6) has cem_f drift far beyond - # HiConstants.GAIN_TEST_HV_DELTA_V["cem_f"], simulating a mid-pointing - # gain test. - assert HiConstants.GAIN_TEST_HV_DELTA_V["cem_f"] < 300.0 - hk_ds["cem_f"].values[4:7] = NOMINAL_HV_VALUES["cem_f"] + 300.0 + # The middle HVSCI segment (indices 4-6) has one field drift beyond + # its HiConstants.GAIN_TEST_HV_DELTA_V tolerance, simulating a + # mid-pointing gain test. + assert HiConstants.GAIN_TEST_HV_DELTA_V[field] < offset + hk_ds[field].values[4:7] = NOMINAL_HV_VALUES[field] + offset # One direct event per housekeeping packet's MET. de_ds = self._make_de_ds([1000, 1001, 1002, 1004, 1005, 1006, 1008, 1009, 1010])