3529 bug hi goodtimes statistical filter 0 normalizes sweep counts by per pointing esa step count

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())
 

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,

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 = [

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])