From c0045972d1755dcd5769a73f046ef7e58bfd2c6c Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 17 Aug 2026 17:09:14 -0600 Subject: [PATCH 1/7] Hi L1C: look up geometric factor by gain state via extended cal-prod schema Replace the separate gain-configuration ancillary file and config_id classification (added in #3358, retired in #3394) with a gain-state- dependent geometric factor lookup folded directly into the existing imap_hi_{sensor}-cal-prod_*.csv file, keyed on measured detector voltage differences (front/back) and the U-Can voltage, per Paul Janzen's feedback (#3391). - hi/utils.py: extend CalibrationProductConfig with a new leading gain_config_id index level, GAIN_MATCH_FIELDS, gain_match_columns (forward-filled per gain_config_id group, placed as the final columns of the file), an extended _validate(), and match_gain_config_id(). - Ancillary CSV (test fixture): add gain_config_id and mcp_delta_v/cem_a_delta_v/cem_b_delta_v/tof_v (+ _tol) columns. Deltas are computed as back-voltage minus front-voltage (positive), matching real flight detector voltages. - hi/hi_l1c.py: implement pset_geometric_factor() to match the L1B DE product's gain_match_{field} attributes (from #3394) against CalibrationProductConfig.match_gain_config_id(), recording the geometric factor per (esa_energy_step, calibration_prod) pair -- the PSET geometric_factor variable gains a calibration_prod dimension instead of summing across calibration products. Fix pset_counts()'s config_row.Index positional lookup, which shifts because of the new leading index level. - hi/hi_l2.py: select a single gain_config_id before the CalibrationProductConfig DataFrame is converted to xarray and reindexed against the map dataset (also a consequence of the new index level). - imap_hi_variable_attrs.yaml: add DEPEND_2/LABL_PTR_2 (calibration_prod) to hi_pset_geometric_factor. Not yet consumed by L2 map-combination processing across different gain states -- that remains a follow-on. Part of #3391. Closes #3395. Co-Authored-By: Claude Sonnet 5 --- .../cdf/config/imap_hi_variable_attrs.yaml | 9 +- imap_processing/hi/hi_l1c.py | 75 ++++-- imap_processing/hi/hi_l2.py | 5 + imap_processing/hi/utils.py | 95 ++++++++ ...map_hi_90sensor-cal-prod_20240101_v001.csv | 60 +++-- imap_processing/tests/hi/test_hi_l1c.py | 126 ++++++++-- imap_processing/tests/hi/test_utils.py | 225 ++++++++++++++++-- 7 files changed, 521 insertions(+), 74 deletions(-) diff --git a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml index 1e7290a0a5..53d5b66083 100644 --- a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml @@ -482,20 +482,23 @@ hi_pset_esa_energy_step: hi_pset_geometric_factor: <<: *default_float32 - CATDESC: Geometric factor for the detector gain configuration active during this pointing + CATDESC: Geometric factor per calibration product and esa_energy_step for the detector gain configuration active during this pointing DEPEND_0: epoch DEPEND_1: esa_energy_step + DEPEND_2: calibration_prod DISPLAY_TYPE: no_plot FIELDNAM: Geometric factor FORMAT: E12.5 LABLAXIS: Geometric Factor LABL_PTR_1: esa_energy_step_label + LABL_PTR_2: calibration_prod_label UNITS: cm^2 sr VALIDMAX: 1.0 VALIDMIN: 0.0 VAR_NOTES: > - Geometric factor per esa_energy_step, looked up from the gain-configuration - ancillary file using the l1b_de's classified gain configuration. + Geometric factor per esa_energy_step and calibration_prod, looked up from + the cal-prod ancillary file's gain_config_id row matching the l1b_de's + reference detector gain state. VAR_TYPE: support_data hi_pset_calibration_prod: diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 87c5a2b23a..2656b7e6a9 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -126,8 +126,13 @@ def generate_pset_dataset( pset_dataset.epoch.data[0] + pset_dataset.epoch_delta.data[0] / 2 ) pset_dataset.update(pset_geometry(pset_midpoint_et, logical_source_parts["sensor"])) - # Look up the per-esa_energy_step geometric factor for this pointing. - pset_dataset.update(pset_geometric_factor(pset_dataset.coords)) + # Look up the per-esa_energy_step geometric factor for this pointing's + # gain state, matched from the L1B DE product's gain_match_{field} + # global attributes against the cal-prod ancillary file's gain_config_id + # rows. + pset_dataset.update( + pset_geometric_factor(pset_dataset.coords, de_dataset, config_df) + ) # Bin the counts into the spin-bins pset_dataset.update( pset_counts(pset_dataset.coords, config_df, de_dataset, goodtimes_ds) @@ -350,34 +355,66 @@ def pset_geometry(pset_et: float, sensor_str: str) -> dict[str, xr.DataArray]: def pset_geometric_factor( pset_coords: dict[str, xr.DataArray], + l1b_de_dataset: xr.Dataset, + config_df: pd.DataFrame, ) -> dict[str, xr.DataArray]: """ - Return a placeholder per-esa_energy_step geometric factor for this pointing. - - The previous gain-configuration ancillary file and config_id - classification mechanism has been retired (see #3391 / #3394) in favor - of gain-test filtering based on a pointing's own reference detector - voltages (see `hi_l1b.de_gain_test_filter`). A replacement geometric - factor lookup -- keyed on the pointing's detector gain state via an - extended cal-prod ancillary file -- is implemented in a follow-up (see - #3395). Until then, "geometric_factor" is left at FILLVAL. + Look up the geometric factor per esa_energy_step and calibration_prod. Parameters ---------- pset_coords : dict[str, xarray.DataArray] The PSET coordinates from the xarray.Dataset. + l1b_de_dataset : xarray.Dataset + The L1B dataset for the pointing being processed. Must have + "gain_match_{field}" global attributes (see + CalibrationProductConfig.GAIN_MATCH_FIELDS) recording the pointing's + reference detector voltage deltas (see hi_l1b.de_gain_test_filter()). + config_df : pandas.DataFrame + Calibration product configuration DataFrame (see + CalibrationProductConfig.from_csv()), indexed by (gain_config_id, + calibration_prod, esa_energy_step). Returns ------- dict[str, xarray.DataArray] - Dictionary containing the "geometric_factor" DataArray (all - FILLVAL), dims (epoch, esa_energy_step). + Dictionary containing the "geometric_factor" DataArray, dims + (epoch, esa_energy_step, calibration_prod). + + Notes + ----- + A pointing's gain state is constant for the whole pointing (see + `hi_l1b.de_gain_test_filter`), so the L1B DE product only records the + pointing's reference detector voltage deltas as global attributes rather + than duplicating the geometric factor across every direct event. This + matches those deltas against the cal-prod ancillary file's gain_config_id + rows, then records the geometric_factor value for each + (esa_energy_step, calibration_prod) pair directly from the matched + gain_config_id's rows. Not yet consumed by L2 processing (deferred to a + follow-on ticket that handles combining PSETs from different gain states + into a single map). """ - return create_dataset_variables( + geometric_factor_var = create_dataset_variables( ["geometric_factor"], coords=pset_coords, att_manager_lookup_str="hi_pset_{0}", ) + hv_deltas = { + field: l1b_de_dataset.attrs[f"gain_match_{field}"] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + if not any(np.isnan(value) for value in hv_deltas.values()): + gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) + if gain_config_id is not None: + gain_config_df = config_df.loc[gain_config_id] + for i, step in enumerate(pset_coords["esa_energy_step"].data): + for j, cal_prod in enumerate(pset_coords["calibration_prod"].data): + geometric_factor_var["geometric_factor"].values[0, i, j] = ( + gain_config_df.loc[ + (int(cal_prod), int(step)), "geometric_factor" + ] + ) + return geometric_factor_var def pset_counts( @@ -457,10 +494,12 @@ def pset_counts( # truncating to an integer gives the correct bin index spin_bin_indices = (filtered_de_ds["spin_phase"].data * N_SPIN_BINS).astype(int) # When iterating over rows of a dataframe, the names of the multi-index - # are not preserved. Below, `config_row.Index[0]` gets the - # calibration_prod value from the namedtuple representing the - # dataframe row. We map this to the array index using cal_prod_to_index. - i_cal_prod = cal_prod_to_index[config_row.Index[0]] + # are not preserved. Below, `config_row.Index[1]` gets the + # calibration_prod value (index level 1 of the (gain_config_id, + # calibration_prod, esa_energy_step) MultiIndex) from the namedtuple + # representing the dataframe row. We map this to the array index + # using cal_prod_to_index. + i_cal_prod = cal_prod_to_index[config_row.Index[1]] np.add.at( counts_var["counts"].data[0, i_esa, i_cal_prod], spin_bin_indices, diff --git a/imap_processing/hi/hi_l2.py b/imap_processing/hi/hi_l2.py index c0a31208ab..d08d0562a7 100644 --- a/imap_processing/hi/hi_l2.py +++ b/imap_processing/hi/hi_l2.py @@ -458,6 +458,11 @@ def calculate_ena_intensity( """ # read calibration product configuration file cal_prod_df = CalibrationProductConfig.from_csv(l2_ancillary_path_dict["cal-prod"]) + # L2 does not yet combine PSETs from different gain states into a single + # map (see hi_l1c.pset_geometric_factor()'s docstring), so use the first + # (and, today, only) gain_config_id present in the ancillary file. + gain_config_id = cal_prod_df.index.get_level_values("gain_config_id")[0] + cal_prod_df = cal_prod_df.loc[gain_config_id] # reindex_like removes esa_energy_steps and calibration products not in the # map_ds esa_energy_step and calibration_product coordinates geometric_factor = cal_prod_df.to_xarray().reindex_like(map_ds)["geometric_factor"] diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 0e1331a97d..3d17acbec0 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -660,9 +660,28 @@ class CalibrationProductConfig(_BaseConfigAccessor): """Register custom accessor for calibration product configuration DataFrames.""" index_columns = ( + "gain_config_id", "calibration_prod", "esa_energy_step", ) + # Detector voltage difference (and U-Can voltage) fields used to match a + # pointing's gain state to a gain_config_id row. See + # hi_l1b.compute_gain_match_values() for how a pointing's own values are + # derived, and match_gain_config_id() below for the matching logic. + GAIN_MATCH_FIELDS = ( + "mcp_delta_v", + "cem_a_delta_v", + "cem_b_delta_v", + "tof_v", + ) + # Columns holding the nominal value and tolerance for each gain match + # field. These are constant across (calibration_prod, esa_energy_step) + # within a gain_config_id, so the CSV only needs to specify them once per + # gain_config_id group -- placed as the final columns of the file, after + # the full calibration product definition. + gain_match_columns = tuple( + f"{field}{suffix}" for field in GAIN_MATCH_FIELDS for suffix in ("", "_tol") + ) required_columns = ( "coincidence_type_list", *[ @@ -670,8 +689,50 @@ class CalibrationProductConfig(_BaseConfigAccessor): for det_pair in _BaseConfigAccessor.tof_detector_pairs for limit in ["low", "high"] ], + *gain_match_columns, ) + def _validate(self, df: pd.DataFrame) -> None: + """ + Validate the calibration product configuration. + + Extends base validation to verify the gain match columns are + non-null and consistent across (calibration_prod, esa_energy_step) + for each gain_config_id. + + Parameters + ---------- + df : pandas.DataFrame + DataFrame to validate. + + Raises + ------ + AttributeError + If required columns or index levels are missing. + ValueError + If gain match values are missing or inconsistent within a + gain_config_id group. + """ + super()._validate(df) + + for gain_config_id, group in df.groupby(level="gain_config_id"): + for col in self.gain_match_columns: + if group[col].isna().any(): + raise ValueError( + f"Missing {col} value(s) for gain_config_id=" + f"{gain_config_id}. The first row for each " + f"gain_config_id must specify a value for every " + f"gain match field." + ) + if group[col].nunique() > 1: + raise ValueError( + f"Inconsistent {col} values across rows for " + f"gain_config_id={gain_config_id}: " + f"{group[col].unique().tolist()}. Gain match values " + f"must be identical across all calibration_prod/" + f"esa_energy_step rows for a gain_config_id." + ) + @classmethod def from_csv(cls, path: str | Path | IO[str]) -> pd.DataFrame: """ @@ -694,6 +755,11 @@ def from_csv(cls, path: str | Path | IO[str]) -> pd.DataFrame: converters={"coincidence_type_list": lambda s: tuple(s.split("|"))}, comment="#", ) + # Forward-fill gain match columns within each gain_config_id group. + # This allows the CSV to specify these values only on the group's + # first row. + gain_cols = list(cls.gain_match_columns) + df[gain_cols] = df.groupby(level="gain_config_id")[gain_cols].ffill() # Trigger the accessor to run validation and add coincidence_type_values _ = df.cal_prod_config.number_of_products return df @@ -711,6 +777,35 @@ def number_of_products(self) -> int: """ return len(self._obj.index.unique(level="calibration_prod")) + def match_gain_config_id(self, hv_deltas: dict[str, float]) -> int | None: + """ + Find the gain_config_id whose reference values match the given deltas. + + Parameters + ---------- + hv_deltas : dict[str, float] + Mapping of CalibrationProductConfig.GAIN_MATCH_FIELDS field names + to a pointing's derived values (see hi_l1b.compute_gain_match_values()). + + Returns + ------- + int or None + The matching gain_config_id, or None if zero or multiple + gain_config_id rows match. + """ + gain_config_ids = self._obj.index.get_level_values("gain_config_id").unique() + matches = [] + for gain_config_id in gain_config_ids: + row = self._obj.loc[gain_config_id].iloc[0] + if all( + abs(hv_deltas[field] - row[field]) <= row[f"{field}_tol"] + for field in self.GAIN_MATCH_FIELDS + ): + matches.append(int(gain_config_id)) + if len(matches) != 1: + return None + return matches[0] + @pd.api.extensions.register_dataframe_accessor("background_config") class BackgroundConfig(_BaseConfigAccessor): diff --git a/imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv b/imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv index 6c15fa2fb5..4a878d33e2 100644 --- a/imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv +++ b/imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv @@ -10,22 +10,44 @@ # DEs with coincidence type contained in the `coincidence_type_list` and time-of-flight # values that are in the inclusive range specified by the low and high entries for each # energy step are included in the angular pset bins. -calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high -0,1,0.00055,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,2,0.00085,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,3,0.00126,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,4,0.00170,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,5,0.00340,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,6,0.00523,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,7,0.00659,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,8,0.01301,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -0,9,0.01830,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25 -1,1,0.00055,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,2,0.00085,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,3,0.00126,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,4,0.00170,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,5,0.00340,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,6,0.00523,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,7,0.00659,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,8,0.01301,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 -1,9,0.01830,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25 \ No newline at end of file +# +# gain_config_id groups rows by the instrument's detector gain state (see +# IMAP-Hi processing issue #3391). geometric_factor is defined per +# (gain_config_id, calibration_prod, esa_energy_step), since a deliberate +# gain change (anticipated ~August 2026) changes the geometric factor at a +# given esa_energy_step. Only one gain state (gain_config_id=0) is +# flight-relevant today. +# +# The final columns (mcp_delta_v/_tol, cem_a_delta_v/_tol, cem_b_delta_v/_tol, +# tof_v/_tol) identify which gain_config_id a pointing is running: mcp_delta_v +# is the nominal (mcp_b_v - mcp_f_v) detector voltage difference, cem_a_delta_v +# is (cem_bk_a_v - cem_f_v), cem_b_delta_v is (cem_bk_b_v - cem_f_v) -- back +# minus front, since front voltages are more negative than back voltages on +# real flight hardware, giving a positive delta -- and tof_v is the nominal +# U-Can voltage monitor value; the "_tol" columns give the allowed +/- +# tolerance in volts for each. A pointing is matched to a gain_config_id if +# its own reference detector voltages (see hi_l1b.de_gain_test_filter()) +# fall within every field's nominal +/- tolerance. These values are constant +# across calibration_prod/esa_energy_step within a gain_config_id, so they +# only need to be specified on the group's first row -- subsequent rows are +# forward-filled. Values below are placeholders carried over from prior +# calibration work and should be confirmed with the Hi instrument team. +gain_config_id,calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high,mcp_delta_v,mcp_delta_v_tol,cem_a_delta_v,cem_a_delta_v_tol,cem_b_delta_v,cem_b_delta_v_tol,tof_v,tof_v_tol +0,0,1,0.00055,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,875.0,75.0,2150.0,150.0,2150.0,150.0,-8000.0,50.0 +0,0,2,0.00085,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,3,0.00126,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,4,0.00170,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,5,0.00340,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,6,0.00523,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,7,0.00659,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,8,0.01301,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,0,9,0.01830,ABC1C2|ABC1|AC1C2,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,1,0.00055,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,2,0.00085,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,3,0.00126,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,4,0.00170,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,5,0.00340,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,6,0.00523,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,7,0.00659,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,8,0.01301,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, +0,1,9,0.01830,BC1C2|AB|AC1,15,55,0,70,-50,10,5,25,,,,,,,, diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index 4888073f1d..4069a237a0 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -64,6 +64,14 @@ def test_generate_pset_dataset( """Test coverage for generate_pset_dataset function""" use_fake_spin_data_for_time(482372987.999) l1b_dataset = hi_l1b_de_dataset.copy() + # The real fixture CDF predates the gain_match_{field} L1B global + # attributes; add placeholders matching the test cal-prod config's + # gain_config_id=0 reference values so pset_geometric_factor() has + # something to look up. + l1b_dataset.attrs["gain_match_mcp_delta_v"] = 875.0 + l1b_dataset.attrs["gain_match_cem_a_delta_v"] = 2150.0 + l1b_dataset.attrs["gain_match_cem_b_delta_v"] = 2150.0 + l1b_dataset.attrs["gain_match_tof_v"] = -8000.0 l1b_met = l1b_dataset["ccsds_met"].values[0] # Set repoint start and end times. seconds_per_day = 24 * 60 * 60 @@ -132,6 +140,10 @@ def test_generate_pset_dataset_uses_midpoint_time( attrs={ "Logical_file_id": "imap_hi_l1b_45sensor-de_20250415_v999", "Logical_source": "imap_hi_l1b_45sensor-de", + "gain_match_mcp_delta_v": 875.0, + "gain_match_cem_a_delta_v": 2150.0, + "gain_match_cem_b_delta_v": 2150.0, + "gain_match_tof_v": -8000.0, }, ) @@ -179,24 +191,104 @@ def test_generate_pset_dataset_uses_midpoint_time( assert actual_sensor_arg == "45sensor" -def test_pset_geometric_factor_returns_fillval_placeholder(): - """Test that pset_geometric_factor() returns an all-FILLVAL placeholder. +def test_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path): + """Test pset_geometric_factor for a pointing whose gain-match values match + the test cal-prod config's gain_config_id=0 reference values (within + tolerance). The resulting geometric_factor should record the unique + geometric_factor per esa_energy_step and calibration_prod.""" + config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) + pset_coords = { + "epoch": xr.DataArray([0], dims=["epoch"]), + "esa_energy_step": xr.DataArray(np.arange(1, 10), dims=["esa_energy_step"]), + "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), + } + l1b_de_dataset = xr.Dataset( + attrs={ + "gain_match_mcp_delta_v": 875.0, + "gain_match_cem_a_delta_v": 2150.0, + "gain_match_cem_b_delta_v": 2150.0, + "gain_match_tof_v": -8000.0, + } + ) - The gain-configuration ancillary file and config_id classification - mechanism has been retired (see #3391 / #3394); a real gain-state-aware - lookup is implemented in a follow-up (#3395). Until then this function - is a no-op placeholder. - """ + result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) + + # geometric_factor per esa_energy_step (1-9) and calibration_prod (0, 1) + # for gain_config_id=0. calibration_prod 0 and 1 share identical + # geometric_factor values in the fixture. See + # imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv + per_step = np.array( + [ + 0.00055, + 0.00085, + 0.00126, + 0.00170, + 0.00340, + 0.00523, + 0.00659, + 0.01301, + 0.01830, + ] + ) + expected = np.stack([per_step, per_step], axis=1) + np.testing.assert_allclose( + result["geometric_factor"].values[0], expected, rtol=1e-6 + ) + + +def test_pset_geometric_factor_nan_gain_match_returns_fillval( + hi_test_cal_prod_config_path, +): + """If L1B could not determine reference detector voltages for the + pointing (nan gain_match_{field} attrs), geometric_factor stays FILLVAL.""" + config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) pset_coords = { "epoch": xr.DataArray([0], dims=["epoch"]), "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), + "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), } + l1b_de_dataset = xr.Dataset( + attrs={ + "gain_match_mcp_delta_v": float("nan"), + "gain_match_cem_a_delta_v": 2150.0, + "gain_match_cem_b_delta_v": 2150.0, + "gain_match_tof_v": -8000.0, + } + ) + + result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) + + fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) + np.testing.assert_array_equal( + result["geometric_factor"].values[0], + np.full((3, 2), fillval), + ) + - result = hi_l1c.pset_geometric_factor(pset_coords) +def test_pset_geometric_factor_no_match_returns_fillval(hi_test_cal_prod_config_path): + """If the pointing's gain-match values don't fall within tolerance of any + gain_config_id in the cal-prod config, geometric_factor stays FILLVAL.""" + config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) + pset_coords = { + "epoch": xr.DataArray([0], dims=["epoch"]), + "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), + "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), + } + l1b_de_dataset = xr.Dataset( + attrs={ + "gain_match_mcp_delta_v": 0.0, + "gain_match_cem_a_delta_v": 0.0, + "gain_match_cem_b_delta_v": 0.0, + "gain_match_tof_v": 0.0, + } + ) + + result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) np.testing.assert_array_equal( - result["geometric_factor"].values[0], [fillval, fillval, fillval] + result["geometric_factor"].values[0], + np.full((3, 2), fillval), ) @@ -419,11 +511,11 @@ def test_pset_counts_arbitrary_cal_prod_numbers( """Test pset_counts with non-sequential calibration product numbers.""" # Create a test calibration product config with non-sequential numbers csv_content = """\ -calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high -5,1,0.00055,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023 -5,2,0.00085,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023 -10,1,0.00055,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023 -10,2,0.00085,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023 +gain_config_id,calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high,mcp_delta_v,mcp_delta_v_tol,cem_a_delta_v,cem_a_delta_v_tol,cem_b_delta_v,cem_b_delta_v_tol,tof_v,tof_v_tol +0,5,1,0.00055,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,875.0,75.0,2150.0,150.0,2150.0,150.0,-8000.0,50.0 +0,5,2,0.00085,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,10,1,0.00055,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,10,2,0.00085,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, """ cal_config_df = utils.CalibrationProductConfig.from_csv(io.StringIO(csv_content)) @@ -542,7 +634,11 @@ def test_pset_counts_goodtimes_filtering( # Mock iter_qualified_events_by_config to mark all events as qualified # and return a single (esa_energy, config_row, mask) tuple mock_config_row = MagicMock() - mock_config_row.Index = (0, 1) # (calibration_prod, esa_energy_step) + mock_config_row.Index = ( + 0, + 0, + 1, + ) # (gain_config_id, calibration_prod, esa_energy_step) def mock_iter(de_ds, config_df, esa_energy_steps): n_remaining = len(de_ds["event_met"]) diff --git a/imap_processing/tests/hi/test_utils.py b/imap_processing/tests/hi/test_utils.py index 475c04f190..3d0e0f0a33 100644 --- a/imap_processing/tests/hi/test_utils.py +++ b/imap_processing/tests/hi/test_utils.py @@ -27,6 +27,81 @@ parse_sensor_number, ) +# Nominal gain-match values matching the real cal-prod-config ancillary +# file's gain_config_id=0 row, for building minimal test CSVs. +GAIN_MATCH_0 = { + "mcp_delta_v": 875.0, + "mcp_delta_v_tol": 75.0, + "cem_a_delta_v": 2150.0, + "cem_a_delta_v_tol": 150.0, + "cem_b_delta_v": 2150.0, + "cem_b_delta_v_tol": 150.0, + "tof_v": -8000.0, + "tof_v_tol": 50.0, +} +# A second, well-separated gain configuration used to test matching/ambiguity. +GAIN_MATCH_1 = { + "mcp_delta_v": 500.0, + "mcp_delta_v_tol": 50.0, + "cem_a_delta_v": 1000.0, + "cem_a_delta_v_tol": 50.0, + "cem_b_delta_v": 1000.0, + "cem_b_delta_v_tol": 50.0, + "tof_v": -5000.0, + "tof_v_tol": 50.0, +} + +_CAL_PROD_CSV_HEADER = ( + "gain_config_id,calibration_prod,esa_energy_step,geometric_factor," + "coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high," + "tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high," + "mcp_delta_v,mcp_delta_v_tol,cem_a_delta_v,cem_a_delta_v_tol," + "cem_b_delta_v,cem_b_delta_v_tol,tof_v,tof_v_tol" +) + +_GAIN_MATCH_COLUMNS = ( + "mcp_delta_v", + "mcp_delta_v_tol", + "cem_a_delta_v", + "cem_a_delta_v_tol", + "cem_b_delta_v", + "cem_b_delta_v_tol", + "tof_v", + "tof_v_tol", +) + + +def _cal_prod_csv_row( + gain_config_id, + calibration_prod, + esa_energy_step, + geometric_factor=0.00055, + coincidence_type_list="ABC1C2", + tof_windows=(15, 55, 0, 70, -50, 10, 5, 25), + gain_match_values=None, +): + """Build a single calibration product config CSV data row (as a string). + + Parameters + ---------- + gain_match_values : dict or None + Optional overrides/values for the gain match columns + (mcp_delta_v, mcp_delta_v_tol, cem_a_delta_v, cem_a_delta_v_tol, + cem_b_delta_v, cem_b_delta_v_tol, tof_v, tof_v_tol). Any column not + present in this dict is left blank in the row (useful for testing + forward-fill and missing-value validation). + """ + gain_match_values = gain_match_values or {} + gain_match_str = ",".join( + str(gain_match_values[col]) if col in gain_match_values else "" + for col in _GAIN_MATCH_COLUMNS + ) + tof_str = ",".join(str(v) for v in tof_windows) + return ( + f"{gain_config_id},{calibration_prod},{esa_energy_step}," + f"{geometric_factor},{coincidence_type_list},{tof_str},{gain_match_str}" + ) + def test_hiapid(): """Test coverage for HIAPID class""" @@ -358,8 +433,8 @@ def test_wrong_columns(self): df = pd.DataFrame( {col: [1, 2, 3] for col in include_columns}, index=pd.MultiIndex.from_tuples( - [(0, 0), (0, 1), (1, 0)], - names=["calibration_prod", "esa_energy_step"], + [(0, 0, 0), (0, 0, 1), (0, 1, 0)], + names=["gain_config_id", "calibration_prod", "esa_energy_step"], ), ) with pytest.raises(AttributeError, match="Required column.*"): @@ -370,7 +445,7 @@ def test_from_csv(self, hi_test_cal_prod_config_path): df = imap_processing.hi.utils.CalibrationProductConfig.from_csv( hi_test_cal_prod_config_path ) - assert isinstance(df["coincidence_type_list"][0, 1], tuple) + assert isinstance(df["coincidence_type_list"][0, 0, 1], tuple) def test_added_coincidence_type_values_column(self, hi_test_cal_prod_config_path): df = CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) @@ -405,15 +480,27 @@ def test_calibration_product_numbers(self, hi_test_cal_prod_config_path): def test_calibration_product_numbers_arbitrary_values(self): """Test calibration_product_numbers with arbitrary non-sequential values.""" # Create a temporary CSV with non-sequential calibration product numbers - csv_content = """\ -calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high -10,1,0.00055,BC1C2,15,55,0,70,-50,10,5,25 -10,2,0.00085,BC1C2,15,55,0,70,-50,10,5,25 -5,1,0.00055,ABC1C2,15,55,0,70,-50,10,5,25 -5,2,0.00085,ABC1C2,15,55,0,70,-50,10,5,25 -100,1,0.00055,AC1,15,55,0,70,-50,10,5,25 -100,2,0.00085,AC1,15,55,0,70,-50,10,5,25 - """ + rows = [ + _cal_prod_csv_row( + 0, 10, 1, coincidence_type_list="BC1C2", gain_match_values=GAIN_MATCH_0 + ), + _cal_prod_csv_row( + 0, + 10, + 2, + geometric_factor=0.00085, + coincidence_type_list="BC1C2", + ), + _cal_prod_csv_row(0, 5, 1, coincidence_type_list="ABC1C2"), + _cal_prod_csv_row( + 0, 5, 2, geometric_factor=0.00085, coincidence_type_list="ABC1C2" + ), + _cal_prod_csv_row(0, 100, 1, coincidence_type_list="AC1"), + _cal_prod_csv_row( + 0, 100, 2, geometric_factor=0.00085, coincidence_type_list="AC1" + ), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) cal_prod_numbers = df.cal_prod_config.calibration_product_numbers @@ -422,6 +509,104 @@ def test_calibration_product_numbers_arbitrary_values(self): np.testing.assert_array_equal(cal_prod_numbers, np.array([5, 10, 100])) assert isinstance(cal_prod_numbers, np.ndarray) + def test_from_csv_forward_fills_gain_match_columns(self): + """Test that gain match columns are forward-filled within a gain group.""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row(0, 0, 2, geometric_factor=0.00085), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + for col, val in GAIN_MATCH_0.items(): + assert df.loc[(0, 0, 2), col] == val + + def test_from_csv_missing_first_row_gain_match_raises(self): + """Test that missing gain match values on the first row raises.""" + rows = [ + _cal_prod_csv_row(0, 0, 1), + _cal_prod_csv_row(0, 0, 2, geometric_factor=0.00085), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + + with pytest.raises(ValueError, match="Missing mcp_delta_v"): + CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + def test_from_csv_inconsistent_gain_match_raises(self): + """Test inconsistent gain match values within a gain_config_id raises.""" + overridden = dict(GAIN_MATCH_0, mcp_delta_v=-900.0) + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row( + 0, 0, 2, geometric_factor=0.00085, gain_match_values=overridden + ), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + + with pytest.raises(ValueError, match="Inconsistent mcp_delta_v"): + CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + def test_match_gain_config_id_exact_match(self): + """Test that a pointing's HV deltas match the correct gain_config_id.""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row(0, 0, 2, geometric_factor=0.00085), + _cal_prod_csv_row(1, 0, 1, gain_match_values=GAIN_MATCH_1), + _cal_prod_csv_row(1, 0, 2, geometric_factor=0.00085), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + hv_deltas = { + "mcp_delta_v": 880.0, + "cem_a_delta_v": 2140.0, + "cem_b_delta_v": 2160.0, + "tof_v": -7990.0, + } + assert df.cal_prod_config.match_gain_config_id(hv_deltas) == 0 + + hv_deltas = { + "mcp_delta_v": 510.0, + "cem_a_delta_v": 990.0, + "cem_b_delta_v": 1010.0, + "tof_v": -4990.0, + } + assert df.cal_prod_config.match_gain_config_id(hv_deltas) == 1 + + def test_match_gain_config_id_no_match(self): + """Test that HV deltas matching no gain_config_id return None.""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row(1, 0, 1, gain_match_values=GAIN_MATCH_1), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + hv_deltas = { + "mcp_delta_v": 0.0, + "cem_a_delta_v": 0.0, + "cem_b_delta_v": 0.0, + "tof_v": 0.0, + } + assert df.cal_prod_config.match_gain_config_id(hv_deltas) is None + + def test_match_gain_config_id_ambiguous_returns_none(self): + """Test that HV deltas matching multiple gain_config_ids return None.""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row(1, 0, 1, gain_match_values=GAIN_MATCH_0), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + hv_deltas = { + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, + } + assert df.cal_prod_config.match_gain_config_id(hv_deltas) is None + class TestGoodMetRangeLookupTable: """Test suite for GoodMetRangeLookupTable class.""" @@ -997,10 +1182,11 @@ def mock_cal_product_config(self): "tof_bc1_high": [50, 50, 50, 50], "tof_c1c2_low": [20, 20, 20, 20], "tof_c1c2_high": [120, 120, 120, 120], + **{col: [val] * 4 for col, val in GAIN_MATCH_0.items()}, } index = pd.MultiIndex.from_tuples( - [(1, 1), (1, 2), (2, 1), (2, 2)], - names=["calibration_prod", "esa_energy_step"], + [(0, 1, 1), (0, 1, 2), (0, 2, 1), (0, 2, 2)], + names=["gain_config_id", "calibration_prod", "esa_energy_step"], ) df = pd.DataFrame(data, index=index) # Trigger the accessor to add coincidence_type_values column @@ -1196,10 +1382,11 @@ def mock_cal_product_config(self): "tof_bc1_high": [50, 50, 50, 50], "tof_c1c2_low": [20, 20, 20, 20], "tof_c1c2_high": [120, 120, 120, 120], + **{col: [val] * 4 for col, val in GAIN_MATCH_0.items()}, } index = pd.MultiIndex.from_tuples( - [(1, 1), (1, 2), (2, 1), (2, 2)], - names=["calibration_prod", "esa_energy_step"], + [(0, 1, 1), (0, 1, 2), (0, 2, 1), (0, 2, 2)], + names=["gain_config_id", "calibration_prod", "esa_energy_step"], ) df = pd.DataFrame(data, index=index) # Trigger the accessor to add coincidence_type_values column @@ -1313,7 +1500,7 @@ def test_filters_by_coincidence_and_tof( for esa_energy, config_row, mask in iter_qualified_events_by_config( mock_de_dataset, mock_cal_product_config, esa_energy_steps ): - if esa_energy == 1 and config_row.Index[0] == 1: + if esa_energy == 1 and config_row.Index[1] == 1: # Events with coincidence 15 or 14: indices 0, 1, 4, 5, 8 # But event 4 has bad TOF (200), so should fail # Events 3, 7 have wrong coincidence (8) @@ -1334,7 +1521,7 @@ def test_different_cal_products_different_masks( mock_de_dataset, mock_cal_product_config, esa_energy_steps ): if esa_energy == 1: # Only look at ESA 1 - cal_prod = config_row.Index[0] + cal_prod = config_row.Index[1] masks_by_cal_prod[cal_prod] = mask # Cal prod 1 accepts ABC1C2 and ABC1 @@ -1404,7 +1591,7 @@ def test_fill_values_pass_tof_check(self, mock_cal_product_config, mock_de_datas for esa_energy, config_row, mask in iter_qualified_events_by_config( mock_de_dataset, mock_cal_product_config, esa_energy_steps ): - if esa_energy == 1 and config_row.Index[0] == 1: + if esa_energy == 1 and config_row.Index[1] == 1: # Event 4 should now pass (has coincidence 15 and fill value TOF) assert mask[4] break From 126784ec72fcfa051f8b9436394fc4cb505f493d Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 19 Aug 2026 09:48:17 -0600 Subject: [PATCH 2/7] Move definition of gain match formulas and gain match fields to one unified location --- imap_processing/hi/hi_l1b.py | 38 +++------------------ imap_processing/hi/utils.py | 45 ++++++++++++++++++++++--- imap_processing/tests/hi/test_hi_l1b.py | 33 +++++++----------- imap_processing/tests/hi/test_utils.py | 32 ++++++++++++++++++ 4 files changed, 90 insertions(+), 58 deletions(-) diff --git a/imap_processing/hi/hi_l1b.py b/imap_processing/hi/hi_l1b.py index e78777791a..2f363eb866 100644 --- a/imap_processing/hi/hi_l1b.py +++ b/imap_processing/hi/hi_l1b.py @@ -14,6 +14,7 @@ from imap_processing.hi.hi_l1a import MILLISECOND_TO_S from imap_processing.hi.utils import ( HIAPID, + CalibrationProductConfig, CoincidenceBitmap, EsaEnergyStepLookupTable, GoodMetRangeLookupTable, @@ -500,37 +501,6 @@ def compute_reference_hv_values(hk_segment_ds: xr.Dataset) -> dict[str, float]: } -def compute_gain_match_values(raw_hv_values: dict[str, float]) -> dict[str, float]: - """ - Derive the back/front voltage differences used for geometric factor lookup. - - Computed as back minus front (rather than front minus back) so that the - resulting deltas are positive, consistent with real flight detector - voltages (front voltages are more negative than back voltages -- see - imap_processing/hi/gain_test_analysis.ipynb). - - Parameters - ---------- - raw_hv_values : dict[str, float] - Raw detector high voltage values keyed by field name, e.g. as - returned by compute_reference_hv_values() (must contain "mcp_f", - "mcp_b", "cem_f", "cem_bk_a", "cem_bk_b", and "tof"). - - Returns - ------- - dict[str, float] - Dictionary with keys "mcp_delta_v", "cem_a_delta_v", "cem_b_delta_v", - and "tof_v", matching CalibrationProductConfig.GAIN_MATCH_FIELDS, for - use with CalibrationProductConfig.match_gain_config_id(). - """ - return { - "mcp_delta_v": raw_hv_values["mcp_b"] - raw_hv_values["mcp_f"], - "cem_a_delta_v": raw_hv_values["cem_bk_a"] - raw_hv_values["cem_f"], - "cem_b_delta_v": raw_hv_values["cem_bk_b"] - raw_hv_values["cem_f"], - "tof_v": raw_hv_values["tof"], - } - - def de_gain_test_filter( l1b_de_ds: xr.Dataset, l1b_hk_ds: xr.Dataset, @@ -567,7 +537,7 @@ def de_gain_test_filter( here -- downstream processing (L1C) looks up the geometric factor per esa_energy_step from the cal-prod ancillary file's matching gain_config_id, using these recorded "gain_match_{field}" attributes - (see hi_l1c.pset_geometric_factor()). + (see hi_l1c.add_pset_geometric_factor()). l1b_hk_ds : xarray.Dataset L1B housekeeping data coincident with the L1A DE data. @@ -578,7 +548,7 @@ def de_gain_test_filter( """ nan_gain_match_attrs = { f"gain_match_{field}": value - for field, value in compute_gain_match_values( + for field, value in CalibrationProductConfig.compute_gain_match_values( {field: np.nan for field in HiConstants.GAIN_TEST_HV_DELTA_V} ).items() } @@ -670,7 +640,7 @@ def de_gain_test_filter( ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE ) - gain_match_values = compute_gain_match_values(reference_hv) + gain_match_values = CalibrationProductConfig.compute_gain_match_values(reference_hv) l1b_de_ds.attrs.update( {f"gain_match_{field}": value for field, value in gain_match_values.items()} ) diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 3d17acbec0..02de6d92cf 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -666,7 +666,7 @@ class CalibrationProductConfig(_BaseConfigAccessor): ) # Detector voltage difference (and U-Can voltage) fields used to match a # pointing's gain state to a gain_config_id row. See - # hi_l1b.compute_gain_match_values() for how a pointing's own values are + # compute_gain_match_values() for how a pointing's own values are # derived, and match_gain_config_id() below for the matching logic. GAIN_MATCH_FIELDS = ( "mcp_delta_v", @@ -777,6 +777,39 @@ def number_of_products(self) -> int: """ return len(self._obj.index.unique(level="calibration_prod")) + @classmethod + def compute_gain_match_values( + cls, raw_hv_values: dict[str, float] + ) -> dict[str, float]: + """ + Derive the back/front voltage differences used for geometric factor lookup. + + Computed as back minus front (rather than front minus back) so that + the resulting deltas are positive, consistent with real flight + detector voltages (front voltages are more negative than back + voltages -- see imap_processing/hi/gain_test_analysis.ipynb). + + Parameters + ---------- + raw_hv_values : dict[str, float] + Raw detector high voltage values keyed by field name, e.g. as + returned by hi_l1b.compute_reference_hv_values() (must contain + "mcp_f", "mcp_b", "cem_f", "cem_bk_a", "cem_bk_b", and "tof"). + + Returns + ------- + dict[str, float] + Dictionary with keys matching GAIN_MATCH_FIELDS, for use with + match_gain_config_id(). + """ + delta_formulas = { + "mcp_delta_v": raw_hv_values["mcp_b"] - raw_hv_values["mcp_f"], + "cem_a_delta_v": raw_hv_values["cem_bk_a"] - raw_hv_values["cem_f"], + "cem_b_delta_v": raw_hv_values["cem_bk_b"] - raw_hv_values["cem_f"], + "tof_v": raw_hv_values["tof"], + } + return {field: delta_formulas[field] for field in cls.GAIN_MATCH_FIELDS} + def match_gain_config_id(self, hv_deltas: dict[str, float]) -> int | None: """ Find the gain_config_id whose reference values match the given deltas. @@ -785,14 +818,18 @@ def match_gain_config_id(self, hv_deltas: dict[str, float]) -> int | None: ---------- hv_deltas : dict[str, float] Mapping of CalibrationProductConfig.GAIN_MATCH_FIELDS field names - to a pointing's derived values (see hi_l1b.compute_gain_match_values()). + to a pointing's derived values (see compute_gain_match_values()). Returns ------- int or None - The matching gain_config_id, or None if zero or multiple - gain_config_id rows match. + The matching gain_config_id, or None if any input value is NaN + (e.g. because a pointing's reference detector voltages could + not be determined) or if zero or multiple gain_config_id rows + match. """ + if any(np.isnan(value) for value in hv_deltas.values()): + return None gain_config_ids = self._obj.index.get_level_values("gain_config_id").unique() matches = [] for gain_config_id in gain_config_ids: diff --git a/imap_processing/tests/hi/test_hi_l1b.py b/imap_processing/tests/hi/test_hi_l1b.py index 8d99f9dfc5..54a5c41699 100644 --- a/imap_processing/tests/hi/test_hi_l1b.py +++ b/imap_processing/tests/hi/test_hi_l1b.py @@ -12,7 +12,6 @@ annotate_direct_events, any_good_direct_events, compute_coincidence_type_and_tofs, - compute_gain_match_values, compute_hae_coordinates, compute_reference_hv_values, de_ccsds_qf, @@ -24,6 +23,7 @@ housekeeping, ) from imap_processing.hi.utils import ( + CalibrationProductConfig, CoincidenceBitmap, EsaEnergyStepLookupTable, HiConstants, @@ -79,7 +79,9 @@ def gain_test_filter_side_effect(l1b_de_ds, l1b_hk_ds): l1b_de_ds.attrs.update( { f"gain_match_{field}": value - for field, value in compute_gain_match_values(NOMINAL_HV_VALUES).items() + for field, value in CalibrationProductConfig.compute_gain_match_values( + NOMINAL_HV_VALUES + ).items() } ) return l1b_de_ds @@ -101,7 +103,9 @@ def gain_test_filter_side_effect(l1b_de_ds, l1b_hk_ds): assert len(l1b_datasets) == 1 assert l1b_datasets[0].attrs["Logical_source"] == "imap_hi_l1b_45sensor-de" assert l1b_datasets[0].attrs["gain_match_mcp_delta_v"] == pytest.approx( - compute_gain_match_values(NOMINAL_HV_VALUES)["mcp_delta_v"] + CalibrationProductConfig.compute_gain_match_values(NOMINAL_HV_VALUES)[ + "mcp_delta_v" + ] ) assert len(l1b_datasets[0].data_vars) == 18 @@ -518,7 +522,9 @@ def test_nominal_pointing_multiple_matching_segments(self, mock_any_good_de): result["ccsds_qf"].values & np.uint8(ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE) == 0 ) - expected_gain_match = compute_gain_match_values(NOMINAL_HV_VALUES) + expected_gain_match = CalibrationProductConfig.compute_gain_match_values( + NOMINAL_HV_VALUES + ) for field in GAIN_MATCH_FIELDS: assert result.attrs[f"gain_match_{field}"] == pytest.approx( expected_gain_match[field] @@ -573,7 +579,9 @@ def test_mid_pointing_gain_test_excluded(self, mock_any_good_de): # gain_match_* attrs reflect the pointing's reference (first segment), # which is unaffected by the excluded mid-pointing gain test segment. - expected_gain_match = compute_gain_match_values(NOMINAL_HV_VALUES) + expected_gain_match = CalibrationProductConfig.compute_gain_match_values( + NOMINAL_HV_VALUES + ) for field in GAIN_MATCH_FIELDS: assert result.attrs[f"gain_match_{field}"] == pytest.approx( expected_gain_match[field] @@ -633,21 +641,6 @@ def test_compute_reference_hv_values(self): assert isinstance(result[field], float) -class TestComputeGainMatchValues: - """Test suite for compute_gain_match_values function.""" - - def test_compute_gain_match_values(self): - """Test that back/front voltage deltas are computed correctly.""" - result = compute_gain_match_values(NOMINAL_HV_VALUES) - - assert result == { - "mcp_delta_v": 875.0, - "cem_a_delta_v": 2150.0, - "cem_b_delta_v": 2150.0, - "tof_v": -8000.0, - } - - class TestGetEsaToEsaEnergyStepLut: """Test suite for get_esa_to_esa_energy_step_lut function.""" diff --git a/imap_processing/tests/hi/test_utils.py b/imap_processing/tests/hi/test_utils.py index 3d0e0f0a33..ea166d87ec 100644 --- a/imap_processing/tests/hi/test_utils.py +++ b/imap_processing/tests/hi/test_utils.py @@ -607,6 +607,38 @@ def test_match_gain_config_id_ambiguous_returns_none(self): } assert df.cal_prod_config.match_gain_config_id(hv_deltas) is None + def test_match_gain_config_id_nan_returns_none(self): + """Test that a NaN input value returns None without raising.""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + hv_deltas = dict(GAIN_MATCH_0, mcp_delta_v=np.nan) + assert df.cal_prod_config.match_gain_config_id(hv_deltas) is None + + def test_compute_gain_match_values(self): + """Test that back/front voltage deltas are computed correctly.""" + raw_hv_values = { + "mcp_f": -3000.0, + "mcp_b": -2125.0, + "cem_f": -4500.0, + "cem_bk_a": -2350.0, + "cem_bk_b": -2350.0, + "tof": -8000.0, + } + + result = CalibrationProductConfig.compute_gain_match_values(raw_hv_values) + + assert result == { + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, + } + assert tuple(result.keys()) == CalibrationProductConfig.GAIN_MATCH_FIELDS + class TestGoodMetRangeLookupTable: """Test suite for GoodMetRangeLookupTable class.""" From 1bff921ceacfe24c1c79982e1ada238bb14cf8f5 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 19 Aug 2026 09:50:01 -0600 Subject: [PATCH 3/7] Use match gain config in L1C refactor pset_geometric_factor function --- imap_processing/hi/hi_l1c.py | 52 ++++++++++--------- imap_processing/hi/hi_l2.py | 4 +- imap_processing/tests/hi/test_hi_l1c.py | 69 +++++++++++++++---------- 3 files changed, 72 insertions(+), 53 deletions(-) diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 2656b7e6a9..f22850e5bc 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -130,9 +130,7 @@ def generate_pset_dataset( # gain state, matched from the L1B DE product's gain_match_{field} # global attributes against the cal-prod ancillary file's gain_config_id # rows. - pset_dataset.update( - pset_geometric_factor(pset_dataset.coords, de_dataset, config_df) - ) + pset_dataset = add_pset_geometric_factor(pset_dataset, de_dataset, config_df) # Bin the counts into the spin-bins pset_dataset.update( pset_counts(pset_dataset.coords, config_df, de_dataset, goodtimes_ds) @@ -353,18 +351,19 @@ def pset_geometry(pset_et: float, sensor_str: str) -> dict[str, xr.DataArray]: return geometry_vars -def pset_geometric_factor( - pset_coords: dict[str, xr.DataArray], +def add_pset_geometric_factor( + pset_ds: xr.Dataset, l1b_de_dataset: xr.Dataset, config_df: pd.DataFrame, -) -> dict[str, xr.DataArray]: +) -> xr.Dataset: """ - Look up the geometric factor per esa_energy_step and calibration_prod. + Add the geometric_factor variable to a pset dataset in place. Parameters ---------- - pset_coords : dict[str, xarray.DataArray] - The PSET coordinates from the xarray.Dataset. + pset_ds : xarray.Dataset + The PSET dataset being built. Must have "esa_energy_step" and + "calibration_prod" coordinates. l1b_de_dataset : xarray.Dataset The L1B dataset for the pointing being processed. Must have "gain_match_{field}" global attributes (see @@ -377,9 +376,9 @@ def pset_geometric_factor( Returns ------- - dict[str, xarray.DataArray] - Dictionary containing the "geometric_factor" DataArray, dims - (epoch, esa_energy_step, calibration_prod). + xarray.Dataset + The input pset_ds, updated in place with a "geometric_factor" + variable, dims (epoch, esa_energy_step, calibration_prod). Notes ----- @@ -396,25 +395,28 @@ def pset_geometric_factor( """ geometric_factor_var = create_dataset_variables( ["geometric_factor"], - coords=pset_coords, + coords=pset_ds.coords, att_manager_lookup_str="hi_pset_{0}", ) hv_deltas = { field: l1b_de_dataset.attrs[f"gain_match_{field}"] for field in CalibrationProductConfig.GAIN_MATCH_FIELDS } - if not any(np.isnan(value) for value in hv_deltas.values()): - gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) - if gain_config_id is not None: - gain_config_df = config_df.loc[gain_config_id] - for i, step in enumerate(pset_coords["esa_energy_step"].data): - for j, cal_prod in enumerate(pset_coords["calibration_prod"].data): - geometric_factor_var["geometric_factor"].values[0, i, j] = ( - gain_config_df.loc[ - (int(cal_prod), int(step)), "geometric_factor" - ] - ) - return geometric_factor_var + gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) + if gain_config_id is not None: + # config_df.loc[gain_config_id] is indexed by (calibration_prod, + # esa_energy_step). Convert to xarray and reindex onto the pset's + # own coordinate values so it broadcasts directly into the output + # array (which only has dims, not coordinate labels, to reindex_like). + gain_factor_da = config_df.loc[gain_config_id, "geometric_factor"].to_xarray() + gain_factor_da = gain_factor_da.reindex( + esa_energy_step=pset_ds["esa_energy_step"].data, + calibration_prod=pset_ds["calibration_prod"].data, + ) + geometric_factor_var["geometric_factor"].values[0] = gain_factor_da.transpose( + "esa_energy_step", "calibration_prod" + ).values + return pset_ds.update(geometric_factor_var) def pset_counts( diff --git a/imap_processing/hi/hi_l2.py b/imap_processing/hi/hi_l2.py index d08d0562a7..d3c440afb0 100644 --- a/imap_processing/hi/hi_l2.py +++ b/imap_processing/hi/hi_l2.py @@ -459,8 +459,8 @@ def calculate_ena_intensity( # read calibration product configuration file cal_prod_df = CalibrationProductConfig.from_csv(l2_ancillary_path_dict["cal-prod"]) # L2 does not yet combine PSETs from different gain states into a single - # map (see hi_l1c.pset_geometric_factor()'s docstring), so use the first - # (and, today, only) gain_config_id present in the ancillary file. + # map (see hi_l1c.add_pset_geometric_factor()'s docstring), so use the + # first (and, today, only) gain_config_id present in the ancillary file. gain_config_id = cal_prod_df.index.get_level_values("gain_config_id")[0] cal_prod_df = cal_prod_df.loc[gain_config_id] # reindex_like removes esa_energy_steps and calibration products not in the diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index 4069a237a0..0e83c9e85b 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -66,7 +66,7 @@ def test_generate_pset_dataset( l1b_dataset = hi_l1b_de_dataset.copy() # The real fixture CDF predates the gain_match_{field} L1B global # attributes; add placeholders matching the test cal-prod config's - # gain_config_id=0 reference values so pset_geometric_factor() has + # gain_config_id=0 reference values so add_pset_geometric_factor() has # something to look up. l1b_dataset.attrs["gain_match_mcp_delta_v"] = 875.0 l1b_dataset.attrs["gain_match_cem_a_delta_v"] = 2150.0 @@ -191,17 +191,33 @@ def test_generate_pset_dataset_uses_midpoint_time( assert actual_sensor_arg == "45sensor" -def test_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path): - """Test pset_geometric_factor for a pointing whose gain-match values match - the test cal-prod config's gain_config_id=0 reference values (within - tolerance). The resulting geometric_factor should record the unique - geometric_factor per esa_energy_step and calibration_prod.""" +def _make_pset_ds_for_geometric_factor(esa_energy_steps, calibration_prods): + """Build a minimal pset dataset with a spin_angle_bin coordinate. + + The spin_angle_bin coordinate is included (unused by geometric_factor) + to guard against it leaking into the geometric_factor variable's shape. + """ + return xr.Dataset( + coords={ + "epoch": xr.DataArray([0], dims=["epoch"]), + "esa_energy_step": xr.DataArray(esa_energy_steps, dims=["esa_energy_step"]), + "calibration_prod": xr.DataArray( + calibration_prods, dims=["calibration_prod"] + ), + "spin_angle_bin": xr.DataArray(np.arange(5), dims=["spin_angle_bin"]), + } + ) + + +def test_add_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path): + """Test add_pset_geometric_factor for a pointing whose gain-match values + match the test cal-prod config's gain_config_id=0 reference values + (within tolerance). The resulting geometric_factor should record the + unique geometric_factor per esa_energy_step and calibration_prod, and + should not gain a spin_angle_bin dimension from the pset dataset's other + coordinates.""" config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) - pset_coords = { - "epoch": xr.DataArray([0], dims=["epoch"]), - "esa_energy_step": xr.DataArray(np.arange(1, 10), dims=["esa_energy_step"]), - "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), - } + pset_ds = _make_pset_ds_for_geometric_factor(np.arange(1, 10), [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ "gain_match_mcp_delta_v": 875.0, @@ -211,7 +227,14 @@ def test_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path) } ) - result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) + + assert result is pset_ds + assert result["geometric_factor"].dims == ( + "epoch", + "esa_energy_step", + "calibration_prod", + ) # geometric_factor per esa_energy_step (1-9) and calibration_prod (0, 1) # for gain_config_id=0. calibration_prod 0 and 1 share identical @@ -236,17 +259,13 @@ def test_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path) ) -def test_pset_geometric_factor_nan_gain_match_returns_fillval( +def test_add_pset_geometric_factor_nan_gain_match_returns_fillval( hi_test_cal_prod_config_path, ): """If L1B could not determine reference detector voltages for the pointing (nan gain_match_{field} attrs), geometric_factor stays FILLVAL.""" config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) - pset_coords = { - "epoch": xr.DataArray([0], dims=["epoch"]), - "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), - "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), - } + pset_ds = _make_pset_ds_for_geometric_factor([1, 2, 3], [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ "gain_match_mcp_delta_v": float("nan"), @@ -256,7 +275,7 @@ def test_pset_geometric_factor_nan_gain_match_returns_fillval( } ) - result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) np.testing.assert_array_equal( @@ -265,15 +284,13 @@ def test_pset_geometric_factor_nan_gain_match_returns_fillval( ) -def test_pset_geometric_factor_no_match_returns_fillval(hi_test_cal_prod_config_path): +def test_add_pset_geometric_factor_no_match_returns_fillval( + hi_test_cal_prod_config_path, +): """If the pointing's gain-match values don't fall within tolerance of any gain_config_id in the cal-prod config, geometric_factor stays FILLVAL.""" config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) - pset_coords = { - "epoch": xr.DataArray([0], dims=["epoch"]), - "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), - "calibration_prod": xr.DataArray([0, 1], dims=["calibration_prod"]), - } + pset_ds = _make_pset_ds_for_geometric_factor([1, 2, 3], [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ "gain_match_mcp_delta_v": 0.0, @@ -283,7 +300,7 @@ def test_pset_geometric_factor_no_match_returns_fillval(hi_test_cal_prod_config_ } ) - result = hi_l1c.pset_geometric_factor(pset_coords, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) np.testing.assert_array_equal( From ee58d000eb4b1be9ae2c367336f779a9b2bdfa59 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 19 Aug 2026 10:55:00 -0600 Subject: [PATCH 4/7] Unify the way that the gain voltages are computed and stored in the dataset --- .../cdf/config/imap_hi_variable_attrs.yaml | 4 +- imap_processing/hi/hi_l1b.py | 49 +++++++-------- imap_processing/hi/hi_l1c.py | 15 +++-- imap_processing/hi/utils.py | 7 ++- imap_processing/tests/hi/test_hi_l1b.py | 61 ++++++++++--------- imap_processing/tests/hi/test_hi_l1c.py | 54 ++++++++-------- 6 files changed, 97 insertions(+), 93 deletions(-) diff --git a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml index 53d5b66083..50921a251a 100644 --- a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml @@ -482,7 +482,9 @@ hi_pset_esa_energy_step: hi_pset_geometric_factor: <<: *default_float32 - CATDESC: Geometric factor per calibration product and esa_energy_step for the detector gain configuration active during this pointing + CATDESC: > + Geometric factor per calibration product and esa_energy_step for the detector + gain configuration active during this pointing DEPEND_0: epoch DEPEND_1: esa_energy_step DEPEND_2: calibration_prod diff --git a/imap_processing/hi/hi_l1b.py b/imap_processing/hi/hi_l1b.py index 2f363eb866..9e1dd1c81d 100644 --- a/imap_processing/hi/hi_l1b.py +++ b/imap_processing/hi/hi_l1b.py @@ -117,11 +117,11 @@ def annotate_direct_events( Returns ------- l1b_datasets : list[xarray.Dataset] - List containing exactly one L1B direct event dataset. Its - "gain_match_{field}" global attributes (see - `CalibrationProductConfig.GAIN_MATCH_FIELDS`) record the pointing's - reference detector voltage deltas (see `de_gain_test_filter`); these - are NaN if they could not be determined. + List containing exactly one L1B direct event dataset. Its global + attributes (one per `CalibrationProductConfig.GAIN_MATCH_FIELDS`, + named directly by field) record the pointing's reference detector + voltage deltas (see `de_gain_test_filter`); these are NaN if they + could not be determined. """ logger.info( f"Running Hi L1B processing on dataset: " @@ -142,7 +142,7 @@ def annotate_direct_events( # flight data and can spill a packet across a good/bad segment boundary. l1b_de_dataset.update(de_esa_step_met(l1b_de_dataset)) # Modifies "esa_energy_step" and "ccsds_qf" in place, and sets the - # "gain_match_{field}" global attributes. + # pointing's HV delta global attributes. l1b_de_dataset = de_gain_test_filter(l1b_de_dataset, l1b_hk_dataset) l1b_de_dataset.update(compute_coincidence_type_and_tofs(l1b_de_dataset)) l1b_de_dataset.update(de_nominal_bin_and_spin_phase(l1b_de_dataset)) @@ -530,14 +530,14 @@ def de_gain_test_filter( de_esa_step_met()). Modified in place: FILLVAL is forced into "esa_energy_step" for events falling outside a matching HVSCI segment, ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE is set in - "ccsds_qf" for the same events, and new "gain_match_{field}" global - attributes (one per CalibrationProductConfig.GAIN_MATCH_FIELDS) are - set to the pointing's reference voltage deltas (NaN if they could - not be determined). The geometric factor itself is not computed - here -- downstream processing (L1C) looks up the geometric factor - per esa_energy_step from the cal-prod ancillary file's matching - gain_config_id, using these recorded "gain_match_{field}" attributes - (see hi_l1c.add_pset_geometric_factor()). + "ccsds_qf" for the same events, and new global attributes (one per + CalibrationProductConfig.GAIN_MATCH_FIELDS, named directly by + field) are set to the pointing's reference voltage deltas (NaN if + they could not be determined). The geometric factor itself is not + computed here -- downstream processing (L1C) looks up the + geometric factor per esa_energy_step from the cal-prod ancillary + file's matching gain_config_id, using these recorded global + attributes (see hi_l1c.add_pset_geometric_factor()). l1b_hk_ds : xarray.Dataset L1B housekeeping data coincident with the L1A DE data. @@ -546,19 +546,16 @@ def de_gain_test_filter( l1b_de_ds : xarray.Dataset The same dataset passed in, modified in place as described above. """ - nan_gain_match_attrs = { - f"gain_match_{field}": value - for field, value in CalibrationProductConfig.compute_gain_match_values( - {field: np.nan for field in HiConstants.GAIN_TEST_HV_DELTA_V} - ).items() - } + nan_hv_deltas = CalibrationProductConfig.compute_gain_match_values( + {field: np.nan for field in HiConstants.GAIN_TEST_HV_DELTA_V} + ) # Check for no valid direct events. if not any_good_direct_events(l1b_de_ds): logger.critical( "No good direct events in dataset; skipping gain test filtering." ) - l1b_de_ds.attrs.update(nan_gain_match_attrs) + l1b_de_ds.attrs.update(nan_hv_deltas) return l1b_de_ds segments = _get_hvsci_segments(l1b_hk_ds) @@ -573,7 +570,7 @@ def de_gain_test_filter( l1b_de_ds["ccsds_qf"].values[:] |= np.uint8( ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE ) - l1b_de_ds.attrs.update(nan_gain_match_attrs) + l1b_de_ds.attrs.update(nan_hv_deltas) return l1b_de_ds # Use the first ~3 housekeeping packets of the first HVSCI segment as @@ -640,11 +637,9 @@ def de_gain_test_filter( ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE ) - gain_match_values = CalibrationProductConfig.compute_gain_match_values(reference_hv) - l1b_de_ds.attrs.update( - {f"gain_match_{field}": value for field, value in gain_match_values.items()} - ) - logger.info(f"Pointing reference gain match values set: {gain_match_values}.") + hv_deltas = CalibrationProductConfig.compute_gain_match_values(reference_hv) + l1b_de_ds.attrs.update(hv_deltas) + logger.info(f"Pointing reference HV deltas set: {hv_deltas}.") return l1b_de_ds diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index f22850e5bc..74e3ec1bd3 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -127,9 +127,8 @@ def generate_pset_dataset( ) pset_dataset.update(pset_geometry(pset_midpoint_et, logical_source_parts["sensor"])) # Look up the per-esa_energy_step geometric factor for this pointing's - # gain state, matched from the L1B DE product's gain_match_{field} - # global attributes against the cal-prod ancillary file's gain_config_id - # rows. + # gain state, matched from the L1B DE product's HV delta global + # attributes against the cal-prod ancillary file's gain_config_id rows. pset_dataset = add_pset_geometric_factor(pset_dataset, de_dataset, config_df) # Bin the counts into the spin-bins pset_dataset.update( @@ -365,10 +364,10 @@ def add_pset_geometric_factor( The PSET dataset being built. Must have "esa_energy_step" and "calibration_prod" coordinates. l1b_de_dataset : xarray.Dataset - The L1B dataset for the pointing being processed. Must have - "gain_match_{field}" global attributes (see - CalibrationProductConfig.GAIN_MATCH_FIELDS) recording the pointing's - reference detector voltage deltas (see hi_l1b.de_gain_test_filter()). + The L1B dataset for the pointing being processed. Must have global + attributes set (one per CalibrationProductConfig.GAIN_MATCH_FIELDS, + named directly by field) recording the pointing's reference + detector voltage deltas (see hi_l1b.de_gain_test_filter()). config_df : pandas.DataFrame Calibration product configuration DataFrame (see CalibrationProductConfig.from_csv()), indexed by (gain_config_id, @@ -399,7 +398,7 @@ def add_pset_geometric_factor( att_manager_lookup_str="hi_pset_{0}", ) hv_deltas = { - field: l1b_de_dataset.attrs[f"gain_match_{field}"] + field: l1b_de_dataset.attrs[field] for field in CalibrationProductConfig.GAIN_MATCH_FIELDS } gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 02de6d92cf..00830844f3 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -152,8 +152,8 @@ class HiConstants: # mnemonic name for the U-Can voltage monitor (see IMAP-Hi Algorithm # Document Section 8, Level 0 Packet Definitions). GAIN_TEST_HV_DELTA_V: ClassVar[dict[str, float]] = { - "pos_defl": 50.0, - "neg_defl": 50.0, + "pos_defl": 1500.0, + "neg_defl": 1500.0, "tof": 50.0, "mcp_f": 10.0, "mcp_b": 50.0, @@ -668,6 +668,9 @@ class CalibrationProductConfig(_BaseConfigAccessor): # pointing's gain state to a gain_config_id row. See # compute_gain_match_values() for how a pointing's own values are # derived, and match_gain_config_id() below for the matching logic. + # hi_l1b.de_gain_test_filter() sets these directly as L1B DE global + # attributes and hi_l1c.add_pset_geometric_factor() reads them back + # the same way. GAIN_MATCH_FIELDS = ( "mcp_delta_v", "cem_a_delta_v", diff --git a/imap_processing/tests/hi/test_hi_l1b.py b/imap_processing/tests/hi/test_hi_l1b.py index 54a5c41699..645e4b515f 100644 --- a/imap_processing/tests/hi/test_hi_l1b.py +++ b/imap_processing/tests/hi/test_hi_l1b.py @@ -74,15 +74,10 @@ def test_hi_annotate_direct_events( mock_get_esa_lut.return_value = mock_esa_lut # Mock de_gain_test_filter to pass the dataset through unmodified, with - # nominal gain_match_* attrs set (as it would for a matching pointing). + # nominal HV delta attrs set (as it would for a matching pointing). def gain_test_filter_side_effect(l1b_de_ds, l1b_hk_ds): l1b_de_ds.attrs.update( - { - f"gain_match_{field}": value - for field, value in CalibrationProductConfig.compute_gain_match_values( - NOMINAL_HV_VALUES - ).items() - } + CalibrationProductConfig.compute_gain_match_values(NOMINAL_HV_VALUES) ) return l1b_de_ds @@ -102,10 +97,11 @@ def gain_test_filter_side_effect(l1b_de_ds, l1b_hk_ds): l1b_datasets = annotate_direct_events(l1a_dataset, xr.Dataset(), esa_energies_csv) assert len(l1b_datasets) == 1 assert l1b_datasets[0].attrs["Logical_source"] == "imap_hi_l1b_45sensor-de" - assert l1b_datasets[0].attrs["gain_match_mcp_delta_v"] == pytest.approx( - CalibrationProductConfig.compute_gain_match_values(NOMINAL_HV_VALUES)[ - "mcp_delta_v" - ] + expected_hv_deltas = CalibrationProductConfig.compute_gain_match_values( + NOMINAL_HV_VALUES + ) + assert l1b_datasets[0].attrs["mcp_delta_v"] == pytest.approx( + expected_hv_deltas["mcp_delta_v"] ) assert len(l1b_datasets[0].data_vars) == 18 @@ -411,9 +407,6 @@ def test_de_esa_energy_step(mock_get_esa_lut, mock_read_csv, mock_any_good_de): np.testing.assert_array_equal(fake_dataset["ccsds_qf"].values, expected_qf_bits) -GAIN_MATCH_FIELDS = ("mcp_delta_v", "cem_a_delta_v", "cem_b_delta_v", "tof_v") - - class TestDeGainTestFilter: """Test suite for de_gain_test_filter function.""" @@ -471,15 +464,18 @@ def _make_de_ds( @mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=False) def test_no_good_direct_events(self, mock_any_good_de): - """gain_match_* attrs are all NaN and dataset is otherwise unmodified.""" + """HV delta attrs are all NaN and dataset is otherwise unmodified.""" fake_de_ds = xr.Dataset(attrs={"some_attr": "unchanged"}) result = de_gain_test_filter(fake_de_ds, xr.Dataset()) assert result is fake_de_ds assert result.attrs["some_attr"] == "unchanged" - for field in GAIN_MATCH_FIELDS: - assert np.isnan(result.attrs[f"gain_match_{field}"]) + actual_hv_deltas = { + field: result.attrs[field] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + assert all(np.isnan(value) for value in actual_hv_deltas.values()) @mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=True) def test_no_hvsci_segments(self, mock_any_good_de): @@ -497,8 +493,11 @@ def test_no_hvsci_segments(self, mock_any_good_de): assert np.all( result["ccsds_qf"].values & np.uint8(ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE) ) - for field in GAIN_MATCH_FIELDS: - assert np.isnan(result.attrs[f"gain_match_{field}"]) + actual_hv_deltas = { + field: result.attrs[field] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + assert all(np.isnan(value) for value in actual_hv_deltas.values()) @mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=True) def test_nominal_pointing_multiple_matching_segments(self, mock_any_good_de): @@ -522,13 +521,14 @@ def test_nominal_pointing_multiple_matching_segments(self, mock_any_good_de): result["ccsds_qf"].values & np.uint8(ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE) == 0 ) - expected_gain_match = CalibrationProductConfig.compute_gain_match_values( + expected_hv_deltas = CalibrationProductConfig.compute_gain_match_values( NOMINAL_HV_VALUES ) - for field in GAIN_MATCH_FIELDS: - assert result.attrs[f"gain_match_{field}"] == pytest.approx( - expected_gain_match[field] - ) + actual_hv_deltas = { + field: result.attrs[field] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + assert actual_hv_deltas == pytest.approx(expected_hv_deltas) @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): @@ -577,15 +577,16 @@ def test_mid_pointing_gain_test_excluded(self, mock_any_good_de): ) np.testing.assert_array_equal(result["ccsds_qf"].values, expected_qf_bits) - # gain_match_* attrs reflect the pointing's reference (first segment), + # HV delta attrs reflect the pointing's reference (first segment), # which is unaffected by the excluded mid-pointing gain test segment. - expected_gain_match = CalibrationProductConfig.compute_gain_match_values( + expected_hv_deltas = CalibrationProductConfig.compute_gain_match_values( NOMINAL_HV_VALUES ) - for field in GAIN_MATCH_FIELDS: - assert result.attrs[f"gain_match_{field}"] == pytest.approx( - expected_gain_match[field] - ) + actual_hv_deltas = { + field: result.attrs[field] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + assert actual_hv_deltas == pytest.approx(expected_hv_deltas) @mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=True) def test_uses_esa_step_met_not_ccsds_met(self, mock_any_good_de): diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index 0e83c9e85b..a1fe4342c4 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -64,14 +64,18 @@ def test_generate_pset_dataset( """Test coverage for generate_pset_dataset function""" use_fake_spin_data_for_time(482372987.999) l1b_dataset = hi_l1b_de_dataset.copy() - # The real fixture CDF predates the gain_match_{field} L1B global - # attributes; add placeholders matching the test cal-prod config's - # gain_config_id=0 reference values so add_pset_geometric_factor() has - # something to look up. - l1b_dataset.attrs["gain_match_mcp_delta_v"] = 875.0 - l1b_dataset.attrs["gain_match_cem_a_delta_v"] = 2150.0 - l1b_dataset.attrs["gain_match_cem_b_delta_v"] = 2150.0 - l1b_dataset.attrs["gain_match_tof_v"] = -8000.0 + # The real fixture CDF predates the HV delta L1B global attributes; add + # placeholders matching the test cal-prod config's gain_config_id=0 + # reference values so add_pset_geometric_factor() has something to + # look up. + l1b_dataset.attrs.update( + { + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, + } + ) l1b_met = l1b_dataset["ccsds_met"].values[0] # Set repoint start and end times. seconds_per_day = 24 * 60 * 60 @@ -140,10 +144,10 @@ def test_generate_pset_dataset_uses_midpoint_time( attrs={ "Logical_file_id": "imap_hi_l1b_45sensor-de_20250415_v999", "Logical_source": "imap_hi_l1b_45sensor-de", - "gain_match_mcp_delta_v": 875.0, - "gain_match_cem_a_delta_v": 2150.0, - "gain_match_cem_b_delta_v": 2150.0, - "gain_match_tof_v": -8000.0, + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, }, ) @@ -220,10 +224,10 @@ def test_add_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_p pset_ds = _make_pset_ds_for_geometric_factor(np.arange(1, 10), [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ - "gain_match_mcp_delta_v": 875.0, - "gain_match_cem_a_delta_v": 2150.0, - "gain_match_cem_b_delta_v": 2150.0, - "gain_match_tof_v": -8000.0, + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, } ) @@ -263,15 +267,15 @@ def test_add_pset_geometric_factor_nan_gain_match_returns_fillval( hi_test_cal_prod_config_path, ): """If L1B could not determine reference detector voltages for the - pointing (nan gain_match_{field} attrs), geometric_factor stays FILLVAL.""" + pointing (nan HV delta attrs), geometric_factor stays FILLVAL.""" config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) pset_ds = _make_pset_ds_for_geometric_factor([1, 2, 3], [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ - "gain_match_mcp_delta_v": float("nan"), - "gain_match_cem_a_delta_v": 2150.0, - "gain_match_cem_b_delta_v": 2150.0, - "gain_match_tof_v": -8000.0, + "mcp_delta_v": float("nan"), + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, } ) @@ -293,10 +297,10 @@ def test_add_pset_geometric_factor_no_match_returns_fillval( pset_ds = _make_pset_ds_for_geometric_factor([1, 2, 3], [0, 1]) l1b_de_dataset = xr.Dataset( attrs={ - "gain_match_mcp_delta_v": 0.0, - "gain_match_cem_a_delta_v": 0.0, - "gain_match_cem_b_delta_v": 0.0, - "gain_match_tof_v": 0.0, + "mcp_delta_v": 0.0, + "cem_a_delta_v": 0.0, + "cem_b_delta_v": 0.0, + "tof_v": 0.0, } ) From ae43cc70d0f31f418a8cad817087806a2e06add7 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 19 Aug 2026 11:25:04 -0600 Subject: [PATCH 5/7] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- imap_processing/hi/hi_l2.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/imap_processing/hi/hi_l2.py b/imap_processing/hi/hi_l2.py index d3c440afb0..910546c8e6 100644 --- a/imap_processing/hi/hi_l2.py +++ b/imap_processing/hi/hi_l2.py @@ -461,8 +461,12 @@ def calculate_ena_intensity( # L2 does not yet combine PSETs from different gain states into a single # map (see hi_l1c.add_pset_geometric_factor()'s docstring), so use the # first (and, today, only) gain_config_id present in the ancillary file. - gain_config_id = cal_prod_df.index.get_level_values("gain_config_id")[0] - cal_prod_df = cal_prod_df.loc[gain_config_id] + gain_config_ids = cal_prod_df.index.get_level_values("gain_config_id").unique() + if len(gain_config_ids) != 1: + raise NotImplementedError( + "L2 processing does not yet support multiple gain_config_id values." + ) + cal_prod_df = cal_prod_df.loc[gain_config_ids[0]] # reindex_like removes esa_energy_steps and calibration products not in the # map_ds esa_energy_step and calibration_product coordinates geometric_factor = cal_prod_df.to_xarray().reindex_like(map_ds)["geometric_factor"] From 602baf9569294c4d23ca1d3b8ed0fb519c55656e Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 19 Aug 2026 16:09:00 -0600 Subject: [PATCH 6/7] Address PR feedback regarding bug with multiple gain configurations in cal prod dataframe --- imap_processing/hi/hi_l1c.py | 83 +++++---- imap_processing/hi/utils.py | 28 +++ imap_processing/tests/hi/test_hi_l1c.py | 224 +++++++++++++++++++----- imap_processing/tests/hi/test_utils.py | 40 +++++ 4 files changed, 293 insertions(+), 82 deletions(-) diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 74e3ec1bd3..7a1b579acb 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -111,6 +111,12 @@ def generate_pset_dataset( logical_source_parts = parse_filename_like(de_dataset.attrs["Logical_source"]) # read calibration product configuration file config_df = CalibrationProductConfig.from_csv(calibration_prod_config_path) + # Select this pointing's matched gain state up front + hv_deltas = { + field: de_dataset.attrs[field] + for field in CalibrationProductConfig.GAIN_MATCH_FIELDS + } + gain_config_df = config_df.cal_prod_config.select_gain_config(hv_deltas) # read background configuration file background_df = BackgroundConfig.from_csv(background_config_path) @@ -127,12 +133,11 @@ def generate_pset_dataset( ) pset_dataset.update(pset_geometry(pset_midpoint_et, logical_source_parts["sensor"])) # Look up the per-esa_energy_step geometric factor for this pointing's - # gain state, matched from the L1B DE product's HV delta global - # attributes against the cal-prod ancillary file's gain_config_id rows. - pset_dataset = add_pset_geometric_factor(pset_dataset, de_dataset, config_df) + # gain state. + pset_dataset = add_pset_geometric_factor(pset_dataset, gain_config_df) # Bin the counts into the spin-bins pset_dataset.update( - pset_counts(pset_dataset.coords, config_df, de_dataset, goodtimes_ds) + pset_counts(pset_dataset.coords, gain_config_df, de_dataset, goodtimes_ds) ) # Calculate and add the exposure time to the pset_dataset pset_dataset.update(pset_exposure(pset_dataset.coords, de_dataset, goodtimes_ds)) @@ -352,8 +357,7 @@ def pset_geometry(pset_et: float, sensor_str: str) -> dict[str, xr.DataArray]: def add_pset_geometric_factor( pset_ds: xr.Dataset, - l1b_de_dataset: xr.Dataset, - config_df: pd.DataFrame, + gain_config_df: pd.DataFrame | None, ) -> xr.Dataset: """ Add the geometric_factor variable to a pset dataset in place. @@ -363,15 +367,11 @@ def add_pset_geometric_factor( pset_ds : xarray.Dataset The PSET dataset being built. Must have "esa_energy_step" and "calibration_prod" coordinates. - l1b_de_dataset : xarray.Dataset - The L1B dataset for the pointing being processed. Must have global - attributes set (one per CalibrationProductConfig.GAIN_MATCH_FIELDS, - named directly by field) recording the pointing's reference - detector voltage deltas (see hi_l1b.de_gain_test_filter()). - config_df : pandas.DataFrame - Calibration product configuration DataFrame (see - CalibrationProductConfig.from_csv()), indexed by (gain_config_id, - calibration_prod, esa_energy_step). + gain_config_df : pandas.DataFrame or None + This pointing's matched gain state configuration (see + CalibrationProductConfig.select_gain_config()), indexed by + (calibration_prod, esa_energy_step), or None if the pointing's HV + deltas didn't match exactly one gain_config_id. Returns ------- @@ -384,30 +384,23 @@ def add_pset_geometric_factor( A pointing's gain state is constant for the whole pointing (see `hi_l1b.de_gain_test_filter`), so the L1B DE product only records the pointing's reference detector voltage deltas as global attributes rather - than duplicating the geometric factor across every direct event. This - matches those deltas against the cal-prod ancillary file's gain_config_id - rows, then records the geometric_factor value for each - (esa_energy_step, calibration_prod) pair directly from the matched - gain_config_id's rows. Not yet consumed by L2 processing (deferred to a - follow-on ticket that handles combining PSETs from different gain states - into a single map). + than duplicating the geometric factor across every direct event. Records + the geometric_factor value for each (esa_energy_step, calibration_prod) + pair directly from gain_config_df's rows. Not yet consumed by L2 processing + (deferred to a follow-on ticket that handles combining PSETs from different + gain states into a single map). """ geometric_factor_var = create_dataset_variables( ["geometric_factor"], coords=pset_ds.coords, att_manager_lookup_str="hi_pset_{0}", ) - hv_deltas = { - field: l1b_de_dataset.attrs[field] - for field in CalibrationProductConfig.GAIN_MATCH_FIELDS - } - gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) - if gain_config_id is not None: - # config_df.loc[gain_config_id] is indexed by (calibration_prod, - # esa_energy_step). Convert to xarray and reindex onto the pset's - # own coordinate values so it broadcasts directly into the output - # array (which only has dims, not coordinate labels, to reindex_like). - gain_factor_da = config_df.loc[gain_config_id, "geometric_factor"].to_xarray() + if gain_config_df is not None: + # gain_config_df is indexed by (calibration_prod, esa_energy_step). + # Convert to xarray and reindex onto the pset's own coordinate + # values so it broadcasts directly into the output array (which + # only has dims, not coordinate labels, to reindex_like). + gain_factor_da = gain_config_df["geometric_factor"].to_xarray() gain_factor_da = gain_factor_da.reindex( esa_energy_step=pset_ds["esa_energy_step"].data, calibration_prod=pset_ds["calibration_prod"].data, @@ -420,7 +413,7 @@ def add_pset_geometric_factor( def pset_counts( pset_coords: dict[str, xr.DataArray], - config_df: pd.DataFrame, + gain_config_df: pd.DataFrame | None, l1b_de_dataset: xr.Dataset, goodtimes_ds: xr.Dataset, ) -> dict[str, xr.DataArray]: @@ -431,8 +424,10 @@ def pset_counts( ---------- pset_coords : dict[str, xarray.DataArray] The PSET coordinates from the xarray.Dataset. - config_df : pandas.DataFrame - The calibration product configuration dataframe. + gain_config_df : pandas.DataFrame or None + This pointing's matched gain state configuration indexed by + (calibration_prod, esa_energy_step), or None if the pointing's HV + deltas didn't match exactly one gain_config_id. l1b_de_dataset : xarray.Dataset The L1B dataset for the pointing being processed. goodtimes_ds : xarray.Dataset @@ -441,7 +436,8 @@ def pset_counts( Returns ------- dict[str, xarray.DataArray] - Dictionary containing counts DataArray. + Dictionary containing counts DataArray. All zero if gain_config_df + is None. """ # Generate counts variable filled with zeros counts_var = create_dataset_variables( @@ -450,6 +446,8 @@ def pset_counts( att_manager_lookup_str="hi_pset_{0}", fill_value=0, ) + if gain_config_df is None: + return counts_var # Create mapping from calibration product numbers to array indices cal_prod_to_index = { @@ -484,7 +482,7 @@ def pset_counts( # esa energy step combination. Use the shared generator to iterate over all # config combinations and get qualified event masks. for esa_energy, config_row, qualified_mask in iter_qualified_events_by_config( - de_ds, config_df, esa_energy_steps + de_ds, gain_config_df, esa_energy_steps ): # Filter events using the qualified mask filtered_de_ds = de_ds.isel(event_met=qualified_mask) @@ -495,12 +493,13 @@ def pset_counts( # truncating to an integer gives the correct bin index spin_bin_indices = (filtered_de_ds["spin_phase"].data * N_SPIN_BINS).astype(int) # When iterating over rows of a dataframe, the names of the multi-index - # are not preserved. Below, `config_row.Index[1]` gets the - # calibration_prod value (index level 1 of the (gain_config_id, - # calibration_prod, esa_energy_step) MultiIndex) from the namedtuple + # are not preserved. Below, `config_row.Index[0]` gets the + # calibration_prod value (index level 0 of gain_config_df's + # (calibration_prod, esa_energy_step) MultiIndex, already sliced to + # this pointing's single gain_config_id above) from the namedtuple # representing the dataframe row. We map this to the array index # using cal_prod_to_index. - i_cal_prod = cal_prod_to_index[config_row.Index[1]] + i_cal_prod = cal_prod_to_index[config_row.Index[0]] np.add.at( counts_var["counts"].data[0, i_esa, i_cal_prod], spin_bin_indices, diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 00830844f3..65d5bc7914 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -846,6 +846,34 @@ def match_gain_config_id(self, hv_deltas: dict[str, float]) -> int | None: return None return matches[0] + def select_gain_config(self, hv_deltas: dict[str, float]) -> pd.DataFrame | None: + """ + Select this configuration's rows for a pointing's matched gain state. + + A pointing's gain state is constant for the whole pointing (see + hi_l1b.de_gain_test_filter()), so this only needs to be done once + per pointing and the result shared by every consumer of the + calibration product configuration (geometric factor lookup, counts + binning, etc.) rather than each matching hv_deltas independently. + + Parameters + ---------- + hv_deltas : dict[str, float] + Mapping of CalibrationProductConfig.GAIN_MATCH_FIELDS field names + to a pointing's derived values (see compute_gain_match_values()). + + Returns + ------- + pandas.DataFrame or None + The subset of rows for the matched gain_config_id, indexed by + (calibration_prod, esa_energy_step), or None if hv_deltas don't + match exactly one gain_config_id (see match_gain_config_id()). + """ + gain_config_id = self.match_gain_config_id(hv_deltas) + if gain_config_id is None: + return None + return self._obj.loc[gain_config_id] + @pd.api.extensions.register_dataframe_accessor("background_config") class BackgroundConfig(_BaseConfigAccessor): diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index a1fe4342c4..7a8057e13e 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -15,6 +15,25 @@ from imap_processing.hi.utils import HIAPID, HiConstants from imap_processing.spice.time import met_to_ttj2000ns, ttj2000ns_to_et +# HV deltas matching the test cal-prod config's gain_config_id=0 reference +# values (within tolerance). See +# imap_processing/tests/hi/data/l1/imap_hi_90sensor-cal-prod_20240101_v001.csv +NOMINAL_HV_DELTAS = { + "mcp_delta_v": 875.0, + "cem_a_delta_v": 2150.0, + "cem_b_delta_v": 2150.0, + "tof_v": -8000.0, +} + + +def _select_gain_config_df(config_df, l1b_de_dataset): + """Mirror generate_pset_dataset()'s gain_config_df selection for tests.""" + hv_deltas = { + field: l1b_de_dataset.attrs[field] + for field in utils.CalibrationProductConfig.GAIN_MATCH_FIELDS + } + return config_df.cal_prod_config.select_gain_config(hv_deltas) + @pytest.fixture(scope="module") def hi_l1b_de_dataset(hi_l1_test_data_path): @@ -65,17 +84,9 @@ def test_generate_pset_dataset( use_fake_spin_data_for_time(482372987.999) l1b_dataset = hi_l1b_de_dataset.copy() # The real fixture CDF predates the HV delta L1B global attributes; add - # placeholders matching the test cal-prod config's gain_config_id=0 - # reference values so add_pset_geometric_factor() has something to - # look up. - l1b_dataset.attrs.update( - { - "mcp_delta_v": 875.0, - "cem_a_delta_v": 2150.0, - "cem_b_delta_v": 2150.0, - "tof_v": -8000.0, - } - ) + # placeholders so add_pset_geometric_factor() and pset_counts() have + # something to look up. + l1b_dataset.attrs.update(NOMINAL_HV_DELTAS) l1b_met = l1b_dataset["ccsds_met"].values[0] # Set repoint start and end times. seconds_per_day = 24 * 60 * 60 @@ -144,10 +155,7 @@ def test_generate_pset_dataset_uses_midpoint_time( attrs={ "Logical_file_id": "imap_hi_l1b_45sensor-de_20250415_v999", "Logical_source": "imap_hi_l1b_45sensor-de", - "mcp_delta_v": 875.0, - "cem_a_delta_v": 2150.0, - "cem_b_delta_v": 2150.0, - "tof_v": -8000.0, + **NOMINAL_HV_DELTAS, }, ) @@ -222,16 +230,10 @@ def test_add_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_p coordinates.""" config_df = utils.CalibrationProductConfig.from_csv(hi_test_cal_prod_config_path) pset_ds = _make_pset_ds_for_geometric_factor(np.arange(1, 10), [0, 1]) - l1b_de_dataset = xr.Dataset( - attrs={ - "mcp_delta_v": 875.0, - "cem_a_delta_v": 2150.0, - "cem_b_delta_v": 2150.0, - "tof_v": -8000.0, - } - ) + l1b_de_dataset = xr.Dataset(attrs=NOMINAL_HV_DELTAS) + gain_config_df = _select_gain_config_df(config_df, l1b_de_dataset) - result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, gain_config_df) assert result is pset_ds assert result["geometric_factor"].dims == ( @@ -278,8 +280,9 @@ def test_add_pset_geometric_factor_nan_gain_match_returns_fillval( "tof_v": -8000.0, } ) + gain_config_df = _select_gain_config_df(config_df, l1b_de_dataset) - result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, gain_config_df) fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) np.testing.assert_array_equal( @@ -303,8 +306,9 @@ def test_add_pset_geometric_factor_no_match_returns_fillval( "tof_v": 0.0, } ) + gain_config_df = _select_gain_config_df(config_df, l1b_de_dataset) - result = hi_l1c.add_pset_geometric_factor(pset_ds, l1b_de_dataset, config_df) + result = hi_l1c.add_pset_geometric_factor(pset_ds, gain_config_df) fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) np.testing.assert_array_equal( @@ -399,17 +403,23 @@ def test_pset_counts( hi_test_background_config_path, ): """Test coverage for pset_counts function.""" + # The real fixture CDF predates the HV delta L1B global attributes; add + # placeholders matching the test cal-prod config's gain_config_id=0 + # reference values so pset_counts() has a gain_config_id to match. + l1b_dataset = hi_l1b_de_dataset.copy() + l1b_dataset.attrs.update(NOMINAL_HV_DELTAS) cal_config_df = utils.CalibrationProductConfig.from_csv( hi_test_cal_prod_config_path ) empty_pset = hi_l1c.empty_pset_dataset( 100, - hi_l1b_de_dataset.esa_energy_step, + l1b_dataset.esa_energy_step, cal_config_df.cal_prod_config.calibration_product_numbers, HIAPID.H90_SCI_DE.sensor, ) + gain_config_df = _select_gain_config_df(cal_config_df, l1b_dataset) counts_var = hi_l1c.pset_counts( - empty_pset.coords, cal_config_df, hi_l1b_de_dataset, hi_goodtimes_dataset + empty_pset.coords, gain_config_df, l1b_dataset, hi_goodtimes_dataset ) assert "counts" in counts_var @@ -428,6 +438,10 @@ def test_pset_counts_empty_l1b( # remove all but one event and set its trigger_id to zero l1b_dataset = hi_l1b_de_dataset.isel(event_met=[0]).copy(deep=True) l1b_dataset["trigger_id"].data[0] = 0 + # The real fixture CDF predates the HV delta L1B global attributes; add + # placeholders matching the test cal-prod config's gain_config_id=0 + # reference values so pset_counts() has a gain_config_id to match. + l1b_dataset.attrs.update(NOMINAL_HV_DELTAS) cal_config_df = utils.CalibrationProductConfig.from_csv( hi_test_cal_prod_config_path ) @@ -437,8 +451,9 @@ def test_pset_counts_empty_l1b( cal_config_df.cal_prod_config.calibration_product_numbers, HIAPID.H90_SCI_DE.sensor, ) + gain_config_df = _select_gain_config_df(cal_config_df, l1b_dataset) counts_var = hi_l1c.pset_counts( - empty_pset.coords, cal_config_df, l1b_dataset, hi_goodtimes_dataset + empty_pset.coords, gain_config_df, l1b_dataset, hi_goodtimes_dataset ) assert counts_var["counts"].data.sum() == 0 @@ -541,6 +556,12 @@ def test_pset_counts_arbitrary_cal_prod_numbers( cal_config_df = utils.CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + # The real fixture CDF predates the HV delta L1B global attributes; add + # placeholders matching this test's gain_config_id=0 reference values + # so pset_counts() has a gain_config_id to match. + l1b_dataset = hi_l1b_de_dataset.copy() + l1b_dataset.attrs.update(NOMINAL_HV_DELTAS) + # Create PSET with non-sequential calibration product numbers l1b_met = 482373065 use_fake_repoint_data_for_time( @@ -549,7 +570,7 @@ def test_pset_counts_arbitrary_cal_prod_numbers( empty_pset = hi_l1c.empty_pset_dataset( l1b_met, - hi_l1b_de_dataset.esa_energy_step, + l1b_dataset.esa_energy_step, cal_config_df.cal_prod_config.calibration_product_numbers, HIAPID.H90_SCI_DE.sensor, ) @@ -557,12 +578,14 @@ def test_pset_counts_arbitrary_cal_prod_numbers( # Verify the calibration_prod coordinate has non-sequential values np.testing.assert_array_equal(empty_pset.calibration_prod.data, np.array([5, 10])) + gain_config_df = _select_gain_config_df(cal_config_df, l1b_dataset) + # Mock get_pointing_times to avoid SPICE kernel requirements with mock.patch( "imap_processing.hi.hi_l1c.get_pointing_times", return_value=(100, 200) ): counts_var = hi_l1c.pset_counts( - empty_pset.coords, cal_config_df, hi_l1b_de_dataset, hi_goodtimes_dataset + empty_pset.coords, gain_config_df, l1b_dataset, hi_goodtimes_dataset ) # Verify counts array has correct shape based on coordinates @@ -594,6 +617,131 @@ def test_pset_counts_arbitrary_cal_prod_numbers( ) +@pytest.mark.external_test_data +def test_pset_counts_restricted_to_matched_gain_config_id( + hi_l1b_de_dataset, hi_goodtimes_dataset, use_fake_repoint_data_for_time +): + """Test pset_counts only qualifies events against the pointing's own + matched gain_config_id. + + Regression test: with two gain_config_id groups that have identical + calibration product definitions, pset_counts must not iterate both + groups' rows into the same (esa_energy_step, calibration_prod) counts + cell -- doing so would double every count once the ancillary file has + more than one gain_config_id. + """ + # Two gain_config_id groups (0 and 1) with identical calibration + # product definitions but distinct HV delta reference values. + csv_content = """\ +gain_config_id,calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high,mcp_delta_v,mcp_delta_v_tol,cem_a_delta_v,cem_a_delta_v_tol,cem_b_delta_v,cem_b_delta_v_tol,tof_v,tof_v_tol +0,5,1,0.00055,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,875.0,75.0,2150.0,150.0,2150.0,150.0,-8000.0,50.0 +0,5,2,0.00085,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,10,1,0.00055,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,10,2,0.00085,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +1,5,1,0.00055,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,500.0,50.0,1000.0,100.0,1000.0,100.0,-4000.0,50.0 +1,5,2,0.00085,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +1,10,1,0.00055,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +1,10,2,0.00085,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, + """ + + cal_config_df = utils.CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + # Match the pointing to gain_config_id=1's reference values. + l1b_dataset = hi_l1b_de_dataset.copy() + l1b_dataset.attrs.update( + { + "mcp_delta_v": 500.0, + "cem_a_delta_v": 1000.0, + "cem_b_delta_v": 1000.0, + "tof_v": -4000.0, + } + ) + + l1b_met = 482373065 + use_fake_repoint_data_for_time( + np.asarray([l1b_met - 15 * 60, l1b_met + 24 * 60 * 60]) + ) + + empty_pset = hi_l1c.empty_pset_dataset( + l1b_met, + l1b_dataset.esa_energy_step, + cal_config_df.cal_prod_config.calibration_product_numbers, + HIAPID.H90_SCI_DE.sensor, + ) + + gain_config_df = _select_gain_config_df(cal_config_df, l1b_dataset) + + with mock.patch( + "imap_processing.hi.hi_l1c.get_pointing_times", return_value=(100, 200) + ): + counts_var = hi_l1c.pset_counts( + empty_pset.coords, gain_config_df, l1b_dataset, hi_goodtimes_dataset + ) + + # Expected totals match the single-gain_config_id case exercised by + # test_pset_counts_arbitrary_cal_prod_numbers. If pset_counts wrongly + # iterated both gain_config_id groups' identical definitions, these + # totals would be doubled. + esa_1_2_mask = ( + hi_l1b_de_dataset["esa_step"][hi_l1b_de_dataset["ccsds_index"]] < 3 + ).values + coincidence_15_mask = (hi_l1b_de_dataset["coincidence_type"] == 15).values + np.testing.assert_equal( + np.sum(counts_var["counts"].data[:, :, 0]), + np.sum(coincidence_15_mask & esa_1_2_mask), + ) + coincidence_7_mask = (hi_l1b_de_dataset["coincidence_type"] == 7).values + np.testing.assert_equal( + np.sum(counts_var["counts"].data[:, :, 1]), + np.sum(coincidence_7_mask & esa_1_2_mask), + ) + + +@pytest.mark.external_test_data +@mock.patch("imap_processing.hi.hi_l1c.get_pointing_times", return_value=(100, 200)) +def test_pset_counts_no_gain_match_returns_zero_counts( + mock_pointing_times, + hi_l1b_de_dataset, + hi_goodtimes_dataset, + hi_test_cal_prod_config_path, +): + """Test pset_counts returns all-zero counts when the pointing's HV + deltas don't match any gain_config_id. + + Without a unique gain_config_id match, the correct coincidence-type/TOF + window definitions for this pointing are unknown, so no events should + be counted (mirrors add_pset_geometric_factor() leaving + geometric_factor at FILLVAL in the same situation). + """ + cal_config_df = utils.CalibrationProductConfig.from_csv( + hi_test_cal_prod_config_path + ) + l1b_dataset = hi_l1b_de_dataset.copy() + l1b_dataset.attrs.update( + { + "mcp_delta_v": 0.0, + "cem_a_delta_v": 0.0, + "cem_b_delta_v": 0.0, + "tof_v": 0.0, + } + ) + empty_pset = hi_l1c.empty_pset_dataset( + 100, + l1b_dataset.esa_energy_step, + cal_config_df.cal_prod_config.calibration_product_numbers, + HIAPID.H90_SCI_DE.sensor, + ) + + gain_config_df = _select_gain_config_df(cal_config_df, l1b_dataset) + assert gain_config_df is None + + counts_var = hi_l1c.pset_counts( + empty_pset.coords, gain_config_df, l1b_dataset, hi_goodtimes_dataset + ) + + assert counts_var["counts"].data.sum() == 0 + + @mock.patch("imap_processing.hi.hi_l1c.get_pointing_times", return_value=(100, 200)) @mock.patch("imap_processing.hi.hi_l1c.iter_qualified_events_by_config") def test_pset_counts_goodtimes_filtering( @@ -655,23 +803,19 @@ def test_pset_counts_goodtimes_filtering( # Mock iter_qualified_events_by_config to mark all events as qualified # and return a single (esa_energy, config_row, mask) tuple mock_config_row = MagicMock() - mock_config_row.Index = ( - 0, - 0, - 1, - ) # (gain_config_id, calibration_prod, esa_energy_step) + mock_config_row.Index = (0, 1) # (calibration_prod, esa_energy_step) - def mock_iter(de_ds, config_df, esa_energy_steps): + def mock_iter(de_ds, gain_config_df, esa_energy_steps): n_remaining = len(de_ds["event_met"]) yield 1, mock_config_row, np.ones(n_remaining, dtype=bool) mock_iter_qualified.side_effect = mock_iter - # Use MagicMock for cal_config since it's not used with our mock - mock_cal_config = MagicMock() + # Use MagicMock for gain_config_df since it's not used with our mock + mock_gain_config_df = MagicMock() counts_var = hi_l1c.pset_counts( - empty_pset.coords, mock_cal_config, l1b_dataset, goodtimes_ds + empty_pset.coords, mock_gain_config_df, l1b_dataset, goodtimes_ds ) # Only 5 events (METs 100-104) should pass goodtimes filtering diff --git a/imap_processing/tests/hi/test_utils.py b/imap_processing/tests/hi/test_utils.py index ea166d87ec..bd9f2ff62e 100644 --- a/imap_processing/tests/hi/test_utils.py +++ b/imap_processing/tests/hi/test_utils.py @@ -618,6 +618,46 @@ def test_match_gain_config_id_nan_returns_none(self): hv_deltas = dict(GAIN_MATCH_0, mcp_delta_v=np.nan) assert df.cal_prod_config.match_gain_config_id(hv_deltas) is None + def test_select_gain_config_matching(self): + """Test that select_gain_config returns the matched gain_config_id's + rows, indexed by (calibration_prod, esa_energy_step).""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + _cal_prod_csv_row(0, 0, 2, geometric_factor=0.00085), + _cal_prod_csv_row(1, 0, 1, gain_match_values=GAIN_MATCH_1), + _cal_prod_csv_row(1, 0, 2, geometric_factor=0.00085), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + gain_config_df = df.cal_prod_config.select_gain_config(GAIN_MATCH_1) + + assert gain_config_df is not None + assert gain_config_df.index.names == ["calibration_prod", "esa_energy_step"] + np.testing.assert_array_equal( + gain_config_df.index.get_level_values("esa_energy_step"), [1, 2] + ) + # Confirm it's gain_config_id=1's rows, not gain_config_id=0's, by + # checking a gain-match column value only set for gain_config_id=1. + assert gain_config_df.loc[(0, 1), "mcp_delta_v"] == GAIN_MATCH_1["mcp_delta_v"] + + def test_select_gain_config_no_match_returns_none(self): + """Test that select_gain_config returns None when no gain_config_id + matches (mirrors match_gain_config_id's no-match behavior).""" + rows = [ + _cal_prod_csv_row(0, 0, 1, gain_match_values=GAIN_MATCH_0), + ] + csv_content = _CAL_PROD_CSV_HEADER + "\n" + "\n".join(rows) + "\n" + df = CalibrationProductConfig.from_csv(io.StringIO(csv_content)) + + hv_deltas = { + "mcp_delta_v": 0.0, + "cem_a_delta_v": 0.0, + "cem_b_delta_v": 0.0, + "tof_v": 0.0, + } + assert df.cal_prod_config.select_gain_config(hv_deltas) is None + def test_compute_gain_match_values(self): """Test that back/front voltage deltas are computed correctly.""" raw_hv_values = { From a266871f7b5ccfbbb2ab491b3989d70674d91957 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Tue, 25 Aug 2026 15:00:11 -0600 Subject: [PATCH 7/7] Update imap_processing/hi/hi_l1c.py Co-authored-by: Leo Werneck --- imap_processing/hi/hi_l1c.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 7a1b579acb..14be8f2dfe 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -408,7 +408,8 @@ def add_pset_geometric_factor( geometric_factor_var["geometric_factor"].values[0] = gain_factor_da.transpose( "esa_energy_step", "calibration_prod" ).values - return pset_ds.update(geometric_factor_var) + pset_ds.update(geometric_factor_var) + return pset_ds def pset_counts(