From b26b5e3ca25188824cbd5142a35a60e7b260791f Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Tue, 28 Jul 2026 14:20:57 -0600 Subject: [PATCH 1/8] Add new ancillary dependency to hi l1b de processing --- imap_processing/cli.py | 12 ++++++++++-- imap_processing/tests/test_cli.py | 5 ++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 324a46378..de155efbe 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -941,9 +941,17 @@ def do_processing( # noqa: PLR0912 l1b_hk_file = dependencies.get_file_paths( source="hi", data_type="l1b", descriptor="hk" )[0] - esa_energies_csv = dependencies.get_file_paths(data_type="ancillary")[0] + esa_energies_csv = dependencies.get_file_paths( + data_type="ancillary", descriptor="esa-energies" + )[0] + gain_config_csv = dependencies.get_file_paths( + data_type="ancillary", descriptor="gain-configuration" + )[0] datasets = hi_l1b.annotate_direct_events( - load_cdf(l1a_de_file), load_cdf(l1b_hk_file), esa_energies_csv + load_cdf(l1a_de_file), + load_cdf(l1b_hk_file), + esa_energies_csv, + gain_config_csv, ) elif self.data_level == "l1c": if "pset" in self.descriptor: diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index 92acd8cbe..4584721d2 100644 --- a/imap_processing/tests/test_cli.py +++ b/imap_processing/tests/test_cli.py @@ -330,7 +330,10 @@ def test_post_processing_returns_empty_list_if_invoked_with_no_data( "imap_hi_l1a_90sensor-de_20241105_v001.cdf", "imap_hi_l1b_90sensor-hk_20241105_v001.cdf", ], - ["imap_hi_90sensor-esa-energies_20240101_v001.csv"], + [ + "imap_hi_90sensor-esa-energies_20240101_v001.csv", + "imap_hi_90sensor-gain-configuration_20240101_v001.csv", + ], 1, ), ("l1b", "sci", "housekeeping", ["imap_hi_l0_raw_20231212_v001.pkts"], [], 2), From 98132e5eb75aa3237889eb310abbfce5f1db8e48 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 10:05:18 -0600 Subject: [PATCH 2/8] Add gain-configuration dependency to hi L1C --- imap_processing/cli.py | 15 +++++++++------ imap_processing/tests/test_cli.py | 1 + 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/imap_processing/cli.py b/imap_processing/cli.py index de155efbe..18ee82c17 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -969,10 +969,10 @@ def do_processing( # noqa: PLR0912 anc_dependencies = dependencies.get_processing_inputs( data_type="ancillary" ) - if len(anc_dependencies) != 2: + if len(anc_dependencies) != 3: raise ValueError( - f"Expected two ancillary dependencies (cal-prod and " - f"backgrounds). Got " + f"Expected three ancillary dependencies (cal-prod, " + f"backgrounds, and gain-configuration). Got " f"{[anc_dep.descriptor for anc_dep in anc_dependencies]}" ) @@ -984,14 +984,16 @@ def do_processing( # noqa: PLR0912 for dep in anc_dependencies } - # Verify we have both required ancillary files + # Verify we have all required ancillary files if ( "cal-prod" not in anc_path_dict or "backgrounds" not in anc_path_dict + or "gain-configuration" not in anc_path_dict ): raise ValueError( - f"Missing required ancillary files. Expected 'cal-prod' and " - f"'backgrounds', got {list(anc_path_dict.keys())}" + f"Missing required ancillary files. Expected 'cal-prod', " + f"'backgrounds', and 'gain-configuration', got " + f"{list(anc_path_dict.keys())}" ) # Load goodtimes dependency @@ -1009,6 +1011,7 @@ def do_processing( # noqa: PLR0912 anc_path_dict["cal-prod"], load_cdf(goodtimes_paths[0]), anc_path_dict["backgrounds"], + anc_path_dict["gain-configuration"], ) elif self.data_level == "l2": science_paths = dependencies.get_file_paths(source="hi", data_type="l1c") diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index 4584721d2..0dbfbf987 100644 --- a/imap_processing/tests/test_cli.py +++ b/imap_processing/tests/test_cli.py @@ -348,6 +348,7 @@ def test_post_processing_returns_empty_list_if_invoked_with_no_data( [ "imap_hi_45sensor-cal-prod_20240101_v001.csv", "imap_hi_45sensor-backgrounds_20240101_v001.csv", + "imap_hi_45sensor-gain-configuration_20240101_v001.csv", ], 1, ), From aa83dcd510c47cc1a7ab00d29a220f89777a5846 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 14:56:33 -0600 Subject: [PATCH 3/8] Add check for gain tests and general gain configuration to hi l1b --- .../cdf/config/imap_hi_variable_attrs.yaml | 28 +- imap_processing/hi/hi_l1b.py | 355 ++++++++++++++++-- imap_processing/hi/utils.py | 172 +++++++++ imap_processing/quality_flags.py | 4 + imap_processing/tests/hi/conftest.py | 7 + imap_processing/tests/hi/test_hi_l1b.py | 325 +++++++++++++++- imap_processing/tests/hi/test_utils.py | 104 +++++ 7 files changed, 944 insertions(+), 51 deletions(-) diff --git a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml index 985353a55..c4e173967 100644 --- a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml @@ -446,8 +446,10 @@ hi_de_ccsds_qf: VALIDMIN: 0 VAR_NOTES: > Bitwise quality flag for CCSDS packet. Bit 0 (value 1): packet was full - (contained 664 events). Bit 2 (value 4): packet contained events from an - invalid spin. + (contained 664 events). Bit 1 (value 2): ESA (inner/outer) voltage did not + match any known esa_energy_step. Bit 2 (value 4): packet contained events + from an invalid spin. Bit 3 (value 8): detector high voltage (CEM/MCP/U-Can/ + deflector) did not match the pointing's classified gain configuration. # ======= L1C PSET Section ======= @@ -477,6 +479,28 @@ hi_pset_esa_energy_step: VAR_TYPE: support_data dtype: uint8 +hi_pset_geometric_factor: + <<: *default_float32 + CATDESC: Geometric factor for the detector gain configuration active during this pointing + DEPEND_0: epoch + DEPEND_1: esa_energy_step + DISPLAY_TYPE: no_plot + FIELDNAM: Geometric factor + FORMAT: E12.5 + LABLAXIS: Geometric Factor + LABL_PTR_1: esa_energy_step_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 pointing's classified gain configuration (see the + L1B DE product's gain_configuration_id global attribute). FILLVAL where the + gain configuration could not be classified. Not yet consumed by L2 + processing (deferred to a follow-on ticket that handles combining PSETs + from different gain configurations into a single map). + VAR_TYPE: support_data + hi_pset_calibration_prod: CATDESC: Coincidence types are combined into various calibration products FIELDNAM: Calibration product number diff --git a/imap_processing/hi/hi_l1b.py b/imap_processing/hi/hi_l1b.py index a7a648eb5..44f252175 100644 --- a/imap_processing/hi/hi_l1b.py +++ b/imap_processing/hi/hi_l1b.py @@ -16,8 +16,10 @@ HIAPID, CoincidenceBitmap, EsaEnergyStepLookupTable, + GainConfigLookupTable, HiConstants, create_dataset_variables, + load_gain_configuration, parse_sensor_number, ) from imap_processing.quality_flags import ImapHiL1bDeFlags @@ -96,7 +98,10 @@ def housekeeping(packet_file_path: str | Path) -> list[xr.Dataset]: def annotate_direct_events( - l1a_de_dataset: xr.Dataset, l1b_hk_dataset: xr.Dataset, esa_energies_anc: Path + l1a_de_dataset: xr.Dataset, + l1b_hk_dataset: xr.Dataset, + esa_energies_anc: Path, + gain_config_anc: Path, ) -> list[xr.Dataset]: """ Perform Hi L1B processing on direct event data. @@ -109,11 +114,16 @@ def annotate_direct_events( L1B housekeeping data coincident with the L1A DE data. esa_energies_anc : pathlib.Path Location of the esa-energies ancillary csv file. + gain_config_anc : pathlib.Path + Location of the gain-configuration ancillary csv file. Returns ------- l1b_datasets : list[xarray.Dataset] - List containing exactly one L1B direct event dataset. + List containing exactly one L1B direct event dataset. Its + "gain_configuration_id" global attribute records the pointing's + classified gain configuration (see `de_gain_configuration`); this is + `GainConfigLookupTable.NO_MATCH` if it could not be classified. """ logger.info( f"Running Hi L1B processing on dataset: " @@ -121,9 +131,17 @@ def annotate_direct_events( ) l1b_de_dataset = l1a_de_dataset.copy() + # Creates the baseline "ccsds_qf" (PACKET_FULL/BADSPIN bits) that + # de_esa_energy_step() and de_gain_configuration() build on. + l1b_de_dataset.update(de_ccsds_qf(l1b_de_dataset)) l1b_de_dataset.update( de_esa_energy_step(l1b_de_dataset, l1b_hk_dataset, esa_energies_anc) ) + # Modifies "esa_energy_step" and "ccsds_qf" in place, and sets the + # "gain_configuration_id" global attribute. + l1b_de_dataset = de_gain_configuration( + l1b_de_dataset, l1b_hk_dataset, gain_config_anc + ) 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)) l1b_de_dataset.update(compute_hae_coordinates(l1b_de_dataset)) @@ -135,7 +153,6 @@ def annotate_direct_events( ) ) l1b_de_dataset.update(de_esa_step_met(l1b_de_dataset)) - l1b_de_dataset.update(de_ccsds_qf(l1b_de_dataset)) l1b_de_dataset = l1b_de_dataset.drop_vars( [ "src_seq_ctr", @@ -398,15 +415,23 @@ def compute_hae_coordinates(dataset: xr.Dataset) -> dict[str, xr.DataArray]: def de_esa_energy_step( - l1b_de_ds: xr.Dataset, l1b_hk_ds: xr.Dataset, esa_energies_anc: Path + l1b_de_ds: xr.Dataset, + l1b_hk_ds: xr.Dataset, + esa_energies_anc: Path, ) -> dict[str, xr.DataArray]: """ - Compute esa_energy_step for each direct event. + Compute esa_energy_step for each direct event from ESA voltage measurements. + + Must be called after de_ccsds_qf(), which creates the "ccsds_qf" variable + this function modifies in place. Parameters ---------- l1b_de_ds : xarray.Dataset - The partial L1B dataset. + The partial L1B dataset. Must already contain "ccsds_qf" (see + de_ccsds_qf()). Modified in place: ImapHiL1bDeFlags.BAD_ESA_VOLTAGE is + set in "ccsds_qf" for packets whose measured ESA voltage didn't match + any esa_energy_step. l1b_hk_ds : xarray.Dataset L1B housekeeping data coincident with the L1A DE data. esa_energies_anc : pathlib.Path @@ -415,7 +440,11 @@ def de_esa_energy_step( Returns ------- new_vars : dict[str, xarray.DataArray] - Keys are variable names and values are `xarray.DataArray`. + Dictionary with the new "esa_energy_step" DataArray. + de_gain_configuration() must be called after this function to force + FILLVAL into "esa_energy_step" and set its own "ccsds_qf" bit for + events whose detector voltages don't match the pointing's gain + configuration. """ new_vars = create_dataset_variables( ["esa_energy_step"], @@ -426,7 +455,6 @@ def de_esa_energy_step( if not any_good_direct_events(l1b_de_ds): return new_vars - # Get the LUT object using the HK data and esa-energies ancillary csv esa_energies_lut = pd.read_csv(esa_energies_anc, comment="#") esa_to_esa_energy_step_lut = get_esa_to_esa_energy_step_lut( l1b_hk_ds, esa_energies_lut @@ -434,12 +462,282 @@ def de_esa_energy_step( new_vars["esa_energy_step"].values = esa_to_esa_energy_step_lut.query( l1b_de_ds["ccsds_met"].data, l1b_de_ds["esa_step"].data ) + # Set the ccsds_qf quality flag bit for packets whose measured ESA voltage + # didn't match any esa_energy_step. + esa_energy_step_fillval = new_vars["esa_energy_step"].attrs["FILLVAL"] + l1b_de_ds["ccsds_qf"].values[ + new_vars["esa_energy_step"].values == esa_energy_step_fillval + ] |= np.uint8(ImapHiL1bDeFlags.BAD_ESA_VOLTAGE) return new_vars +def de_gain_configuration( + l1b_de_ds: xr.Dataset, + l1b_hk_ds: xr.Dataset, + gain_config_anc: Path, +) -> xr.Dataset: + """ + Classify gain configuration and force FILLVAL for non-matching events. + + Must be called after de_esa_energy_step(), which sets the "esa_energy_step" + and "ccsds_qf" variables this function modifies in place. + + Parameters + ---------- + l1b_de_ds : xarray.Dataset + The partial L1B dataset. Must already contain "esa_energy_step" and + "ccsds_qf" (see de_esa_energy_step()). Modified in place: FILLVAL is + forced into "esa_energy_step" for events whose detector voltages + didn't match the pointing's classified gain configuration, + ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE is set in "ccsds_qf" for the + same events, and the "gain_configuration_id" global attribute is set + to the pointing's classified gain configuration id (or + GainConfigLookupTable.NO_MATCH if it could not be classified). The + geometric factor itself is not computed here -- it is constant for + the whole pointing, so downstream processing (L1C) looks up the + geometric factor per esa_energy_step from the gain-configuration + ancillary file using the recorded "gain_configuration_id". + l1b_hk_ds : xarray.Dataset + L1B housekeeping data coincident with the L1A DE data. + gain_config_anc : pathlib.Path + Location of the gain-configuration ancillary csv file. + + Returns + ------- + l1b_de_ds : xarray.Dataset + The same dataset passed in, modified in place as described above. + """ + # Check for no valid direct events. + if not any_good_direct_events(l1b_de_ds): + l1b_de_ds.attrs["gain_configuration_id"] = GainConfigLookupTable.NO_MATCH + return l1b_de_ds + + gain_config_df = load_gain_configuration(gain_config_anc) + gain_config_lut, config_id = get_gain_configuration_lut(l1b_hk_ds, gain_config_df) + + ccsds_met = l1b_de_ds["ccsds_met"].data + packet_config_ids = gain_config_lut.query(ccsds_met) + detector_voltage_bad_mask = packet_config_ids == GainConfigLookupTable.NO_MATCH + + esa_energy_step = l1b_de_ds["esa_energy_step"] + esa_energy_step.values = np.where( + detector_voltage_bad_mask, + esa_energy_step.attrs["FILLVAL"], + esa_energy_step.values, + ) + l1b_de_ds["ccsds_qf"].values[detector_voltage_bad_mask] |= np.uint8( + ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE + ) + + l1b_de_ds.attrs["gain_configuration_id"] = ( + config_id if config_id is not None else GainConfigLookupTable.NO_MATCH + ) + return l1b_de_ds + + +def _get_hvsci_segments(l1b_hk_ds: xr.Dataset) -> list[tuple[int, int]]: + """ + Find contiguous segments where op_mode == "HVSCI" in housekeeping data. + + Parameters + ---------- + l1b_hk_ds : xarray.Dataset + L1B housekeeping dataset. + + Returns + ------- + segments : list[tuple[int, int]] + List of (start_index, end_index) tuples. `end_index` is exclusive, + suitable for use with `Dataset.isel(epoch=slice(start, end))`. + """ + # Pad the boolean array `op_mode == HVSCI` with False values on each end. + # This treats starting or ending in HVSCI mode as a transition in the next + # step where np.diff is used to find op_mode transitions into and out of + # HVSCI + padded_mask = np.pad( + l1b_hk_ds["op_mode"].data == "HVSCI", (1, 1), constant_values=False + ) + mode_changes = np.diff(padded_mask.astype(int)) + starts = np.nonzero(mode_changes == 1)[0] + ends = np.nonzero(mode_changes == -1)[0] + return list(zip(starts, ends, strict=False)) + + +def _get_config_hv_row(gain_config_df: pd.DataFrame, config_id: int) -> pd.Series: + """ + Get a representative row of HV voltage/tolerance values for a config_id. + + The {field}_v/{field}_delta_v columns are constant across esa_energy_step + within a config_id (see utils.load_gain_configuration()), so any row for + that config_id can be used. + + Parameters + ---------- + gain_config_df : pandas.DataFrame + Gain configuration lookup table indexed by (config_id, esa_energy_step). + config_id : int + Configuration id to look up. + + Returns + ------- + pandas.Series + A single row of {field}_v/{field}_delta_v values for the given config_id. + """ + return gain_config_df.loc[config_id].iloc[0] + + +def _detector_voltage_matches_config( + segment_ds: xr.Dataset, gain_config_row: pd.Series +) -> bool: + """ + Check whether a segment's median detector voltages match a gain config. + + Parameters + ---------- + segment_ds : xarray.Dataset + A contiguous HVSCI segment of L1B housekeeping data. + gain_config_row : pandas.Series + A single row (config_id) of the gain-configuration ancillary table. + + Returns + ------- + bool + True if the median of every field in HiConstants.GAIN_CONFIG_HV_FIELDS + falls within that field's nominal +/- tolerance for this config. + """ + for field in HiConstants.GAIN_CONFIG_HV_FIELDS: + median_value = np.median(segment_ds[field].data) + if ( + abs(median_value - gain_config_row[f"{field}_v"]) + > gain_config_row[f"{field}_delta_v"] + ): + return False + return True + + +def classify_gain_configuration( + l1b_hk_ds: xr.Dataset, gain_config_df: pd.DataFrame +) -> int | None: + """ + Classify which gain configuration (config_id) a pointing is running. + + Uses the median detector high voltages (see HiConstants.GAIN_CONFIG_HV_ + FIELDS) of the first contiguous HVSCI segment in the pointing, matched + against the gain-configuration ancillary table. This assumes that when a + gain test is run during a pointing, the first HVSCI segment of that + pointing is run at the pointing's real (non-gain-test) configuration. + + Parameters + ---------- + l1b_hk_ds : xarray.Dataset + L1B housekeeping dataset. + gain_config_df : pandas.DataFrame + Gain configuration lookup table derived from ancillary file. See + utils.load_gain_configuration(). + + Returns + ------- + config_id : int or None + The classified gain configuration id, or None if the first HVSCI + segment's voltages matched zero or multiple configurations. + """ + segments = _get_hvsci_segments(l1b_hk_ds) + if not segments: + logger.critical("No HVSCI segments found; cannot classify gain configuration.") + return None + + i_start, i_end = segments[0] + first_segment_ds = l1b_hk_ds.isel(epoch=slice(i_start, i_end)) + config_ids = gain_config_df.index.get_level_values("config_id").unique() + matches = [ + config_id + for config_id in config_ids + if _detector_voltage_matches_config( + first_segment_ds, _get_config_hv_row(gain_config_df, config_id) + ) + ] + if len(matches) != 1: + interval = met_to_utc(first_segment_ds["shcoarse"].data[[0, -1]]) + medians = { + field: float(np.median(first_segment_ds[field].data)) + for field in HiConstants.GAIN_CONFIG_HV_FIELDS + } + if len(matches) == 0: + logger.critical( + f"No gain configuration matches found for first HVSCI segment " + f"during interval: ({interval}) with median detector " + f"voltages: {medians}." + ) + else: + logger.critical( + f"Multiple gain configuration matches found ({matches}) for " + f"first HVSCI segment during interval: ({interval}) with " + f"median detector voltages: {medians}." + ) + return None + return int(matches[0]) + + +def get_gain_configuration_lut( + l1b_hk_ds: xr.Dataset, + gain_config_df: pd.DataFrame, +) -> tuple[GainConfigLookupTable, int | None]: + """ + Generate a lookup table that associates MET ranges with a gain config_id. + + Classifies the pointing's overall gain configuration from its first + HVSCI segment (see classify_gain_configuration()), then walks every + contiguous HVSCI segment, recording the pointing's config_id for segments + whose detector voltages still match that configuration. Segments that + don't match (e.g. a mid-pointing gain test), or every segment if the + pointing's configuration could not be classified, are left out of the + LUT entirely and so return GainConfigLookupTable.NO_MATCH when queried. + + Parameters + ---------- + l1b_hk_ds : xarray.Dataset + L1B housekeeping dataset. + gain_config_df : pandas.DataFrame + Gain configuration lookup table derived from ancillary file. See + utils.load_gain_configuration(). + + Returns + ------- + gain_config_lut : GainConfigLookupTable + A lookup table object that can be used to query by MET time for the + pointing's config_id, or GainConfigLookupTable.NO_MATCH. + config_id : int or None + The gain configuration classified for this pointing, or None if it + could not be classified. + + Notes + ----- + This is unrelated to ESA energy step assignment (see + get_esa_to_esa_energy_step_lut()) other than both being evaluated per + contiguous HVSCI segment; the two are intentionally independent so that + they can be understood, tested, and evolved separately. + """ + gain_config_lut = GainConfigLookupTable() + config_id = classify_gain_configuration(l1b_hk_ds, gain_config_df) + + if config_id is not None: + gain_config_row = _get_config_hv_row(gain_config_df, config_id) + for i_start, i_end in _get_hvsci_segments(l1b_hk_ds): + segment_ds = l1b_hk_ds.isel(epoch=slice(i_start, i_end)) + if _detector_voltage_matches_config(segment_ds, gain_config_row): + gain_config_lut.add_entry( + segment_ds["shcoarse"].data[0], + segment_ds["shcoarse"].data[-1], + config_id, + ) + + return gain_config_lut, config_id + + def get_esa_to_esa_energy_step_lut( - l1b_hk_ds: xr.Dataset, esa_energies_lut: pd.DataFrame + l1b_hk_ds: xr.Dataset, + esa_energies_lut: pd.DataFrame, ) -> EsaEnergyStepLookupTable: """ Generate a lookup table that associates an esa_step to an esa_energy_step. @@ -454,8 +752,10 @@ def get_esa_to_esa_energy_step_lut( Returns ------- esa_energy_step_lut : EsaEnergyStepLookupTable - A lookup table object that can be used to query by MET time and esa_step - for the associated esa_energy_step values. + A lookup table object that can be used to query by MET time and + esa_step for the associated esa_energy_step values. Segments/esa_steps + where the measured ESA voltage did not match any esa_energy_step are + left out of the LUT entirely, so querying them returns FILLVAL. Notes ----- @@ -465,19 +765,12 @@ def get_esa_to_esa_energy_step_lut( esa_energy_step_lut = EsaEnergyStepLookupTable() # Get the set of esa_steps visited esa_steps = list(sorted(set(l1b_hk_ds["sci_esa_step"].data))) - # Break into contiguous segments where op_mode == "HVSCI" - # Pad the boolean array `op_mode == HVSCI` with False values on each end. - # This treats starting or ending in HVSCI mode as a transition in the next - # step where np.diff is used to find op_mode transitions into and out of - # HVSCI - padded_mask = np.pad( - l1b_hk_ds["op_mode"].data == "HVSCI", (1, 1), constant_values=False - ) - mode_changes = np.diff(padded_mask.astype(int)) - hsvsci_starts = np.nonzero(mode_changes == 1)[0] - hsvsci_ends = np.nonzero(mode_changes == -1)[0] - for i_start, i_end in zip(hsvsci_starts, hsvsci_ends, strict=False): + + for i_start, i_end in _get_hvsci_segments(l1b_hk_ds): contiguous_hvsci_ds = l1b_hk_ds.isel(dict(epoch=slice(i_start, i_end))) + segment_start = contiguous_hvsci_ds["shcoarse"].data[0] + segment_end = contiguous_hvsci_ds["shcoarse"].data[-1] + # Find median inner and outer ESA voltages for each ESA step for esa_step in esa_steps: single_esa_ds = contiguous_hvsci_ds.where( @@ -529,8 +822,8 @@ def get_esa_to_esa_energy_step_lut( continue # Set LUT to matching esa_energy_step for time range esa_energy_step_lut.add_entry( - contiguous_hvsci_ds["shcoarse"].data[0], - contiguous_hvsci_ds["shcoarse"].data[-1], + segment_start, + segment_end, esa_step, matching_esa_energy["esa_energy_step"].values[0], ) @@ -571,19 +864,23 @@ def de_esa_step_met(dataset: xr.Dataset) -> dict[str, xr.DataArray]: def de_ccsds_qf(dataset: xr.Dataset) -> dict[str, xr.DataArray]: """ - Compute ccsds_qf quality flag for each CCSDS packet. + Compute the baseline ccsds_qf quality flag for each CCSDS packet. - The ccsds_qf is a quality flag bitmask indicating packet characteristics. + Sets the PACKET_FULL and BADSPIN bits. Must be called first, before + de_esa_energy_step() and de_gain_configuration(), which add their own + bits (BAD_ESA_VOLTAGE, BAD_DETECTOR_VOLTAGE) to this same "ccsds_qf" + variable. Parameters ---------- dataset : xarray.Dataset - The L1A/B dataset containing ccsds_index for mapping events to packets. + The L1A/B dataset containing "ccsds_index" and "spin_invalids" for + mapping events to packets. Returns ------- new_vars : dict[str, xarray.DataArray] - Dictionary with "ccsds_qf" key and uint8 DataArray value. + Dictionary with the new "ccsds_qf" DataArray. """ max_events_per_packet = 664 diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 9d532f5c1..92146a60e 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -95,6 +95,9 @@ class HiConstants: in Filter 2. STAT_FILTER_2_BIN_PADDING : int Number of bins to add on each side of pulse angle range in Filter 2. + GAIN_CONFIG_HV_FIELDS : tuple[str, ...] + Detector high voltage monitor field names used to classify which gain + configuration a pointing is running. EXCESS_BACKGROUND_COUNT_RATE : float Constant rate offset (per second) to subtract from combined background rates. Corrects for excess counts from the outer ESA during background @@ -133,6 +136,22 @@ class HiConstants: STAT_FILTER_2_MAX_TIME_DELTA = 5000 * DE_CLOCK_TICK_S STAT_FILTER_2_BIN_PADDING = 1 + # Detector high voltage monitors used to classify which gain configuration + # a pointing is running. See the gain-configuration ancillary file + # (utils.load_gain_configuration()) and hi_l1b.classify_gain_configuration(). + # Note "tof" is the raw CCSDS mnemonic name for the U-Can voltage monitor + # (see IMAP-Hi Algorithm Document Section 8, Level 0 Packet Definitions). + GAIN_CONFIG_HV_FIELDS = ( + "pos_defl", + "neg_defl", + "tof", # aka U-Can + "mcp_f", + "mcp_b", + "cem_f", + "cem_bk_a", + "cem_bk_b", + ) + # Background rate correction # Constant offset to subtract from combined background rates to correct # for excess counts from the outer ESA during background testing. @@ -173,6 +192,73 @@ def parse_sensor_number(full_string: str) -> int: return int(match["sensor_num"]) +def load_gain_configuration(path: str | Path | IO[str]) -> pd.DataFrame: + """ + Load the detector gain configuration ancillary CSV file. + + Each row corresponds to one (config_id, esa_energy_step) combination and + gives that step's geometric factor. The detector high voltage nominal/ + tolerance ("{field}_v"/"{field}_delta_v") columns are constant across + esa_energy_step within a config_id (see HiConstants.GAIN_CONFIG_HV_FIELDS + for the monitored fields), so they only need to be specified on the first + esa_energy_step row of each config_id group -- subsequent rows are + forward-filled. + + Parameters + ---------- + path : str or pathlib.Path or file-like object + Location of the gain configuration ancillary CSV file. + + Returns + ------- + pandas.DataFrame + DataFrame indexed by (config_id, esa_energy_step), with columns + "geometric_factor" and "{field}_v"/"{field}_delta_v" for each field + in HiConstants.GAIN_CONFIG_HV_FIELDS. + + Raises + ------ + ValueError + If geometric_factor is null for any row, if {field}_v/{field}_delta_v + values are missing for a config_id's first esa_energy_step row, or if + they are inconsistent across esa_energy_step within a config_id. + """ + df = pd.read_csv(path, index_col=["config_id", "esa_energy_step"], comment="#") + + hv_columns = [ + f"{field}{suffix}" + for field in HiConstants.GAIN_CONFIG_HV_FIELDS + for suffix in ("_v", "_delta_v") + ] + df[hv_columns] = df.groupby(level="config_id")[hv_columns].ffill() + + if df["geometric_factor"].isna().any(): + null_rows = df.index[df["geometric_factor"].isna()].tolist() + raise ValueError(f"Null geometric_factor value(s) found for rows: {null_rows}.") + + # Validate {field}_v/{field}_delta_v columns are non-null and consistent + # across esa_energy_step for each config_id. This catches a missing first + # row (forward-fill leaves NaN) or an authoring mistake (inconsistent + # values). + for config_id, group in df.groupby(level="config_id"): + for col in hv_columns: + if group[col].isna().any(): + raise ValueError( + f"Missing {col} value(s) for config_id={config_id}. The " + f"first esa_energy_step row for each config_id must " + f"specify voltage and tolerance values for every field." + ) + if group[col].nunique() > 1: + raise ValueError( + f"Inconsistent {col} values across esa_energy_step for " + f"config_id={config_id}: {group[col].unique().tolist()}. " + f"Voltage and tolerance values must be identical across " + f"all esa_energy_step rows for a config_id." + ) + + return df + + def full_dataarray( name: str, attrs: dict, @@ -473,6 +559,92 @@ def query( return results.astype(self._esa_energy_step_dtype) +class GainConfigLookupTable: + """Class for holding a MET range to gain configuration (config_id) LUT.""" + + # Sentinel returned for MET values that don't fall within any segment + # matching the pointing's classified gain configuration (e.g. the + # configuration could not be classified for the pointing, or this + # particular segment doesn't match it -- such as during a mid-pointing + # gain test). This is purely an internal L1B processing detail, not tied + # to any CDF variable's FILLVAL. + NO_MATCH = -1 + + def __init__(self) -> None: + self.df = pd.DataFrame( + { + "start_met": pd.Series(dtype="float64"), + "end_met": pd.Series(dtype="float64"), + "config_id": pd.Series(dtype="int64"), + } + ) + self._indexed = False + + def add_entry(self, start_met: float, end_met: float, config_id: int) -> None: + """ + Add a single entry to the lookup table. + + Parameters + ---------- + start_met : float + Start mission elapsed time of the time range. + end_met : float + End mission elapsed time of the time range. + config_id : int + Gain configuration id that applies over this time range. + """ + new_row = pd.DataFrame( + { + "start_met": [start_met], + "end_met": [end_met], + "config_id": [config_id], + } + ) + self.df = pd.concat([self.df, new_row], ignore_index=True) + self._indexed = False + + def _ensure_indexed(self) -> None: + """Sort the internal DataFrame by start_met for faster queries.""" + if not self._indexed: + self.df = self.df.sort_values("start_met").reset_index(drop=True) + self._indexed = True + + def query(self, query_met: float | Iterable[float]) -> int | np.ndarray: + """ + Query MET(s) to retrieve the matching gain configuration id(s). + + Parameters + ---------- + query_met : float or array_like + Mission elapsed time value(s) to query. + + Returns + ------- + int or numpy.ndarray + - If input is scalar: returns int (config_id or NO_MATCH). + - If input is array-like: returns numpy array of config_ids with + the same length as the input, containing NO_MATCH for queries + with no matching entry. + + Notes + ----- + If multiple entries match a query, returns the first match found. + """ + self._ensure_indexed() + + is_scalar_met = np.isscalar(query_met) + query_mets = np.atleast_1d(query_met) + results = np.full(query_mets.shape, self.NO_MATCH, dtype=np.int64) + + for i, qm in enumerate(query_mets): + mask = (self.df["start_met"] <= qm) & (self.df["end_met"] >= qm) + matches = self.df[mask] + if not matches.empty: + results[i] = matches["config_id"].iloc[0] + + return int(results[0]) if is_scalar_met else results + + class _BaseConfigAccessor: """ Base class for configuration DataFrame accessors. diff --git a/imap_processing/quality_flags.py b/imap_processing/quality_flags.py index e0c48df2b..d59d2ce51 100644 --- a/imap_processing/quality_flags.py +++ b/imap_processing/quality_flags.py @@ -170,4 +170,8 @@ class ImapHiL1bDeFlags(FlagNameMixin): NONE = CommonFlags.NONE PACKET_FULL = 2**0 # bit 0, packet contained 664 events (max capacity) + BAD_ESA_VOLTAGE = 2**1 # bit 1, ESA voltage didn't match any esa_energy_step BADSPIN = ENAFlags.BADSPIN # bit 2, packet contained events from an invalid spin + BAD_DETECTOR_VOLTAGE = ( + 2**3 + ) # bit 3, detector HV didn't match the pointing's gain configuration diff --git a/imap_processing/tests/hi/conftest.py b/imap_processing/tests/hi/conftest.py index a63872ee0..2832736fc 100644 --- a/imap_processing/tests/hi/conftest.py +++ b/imap_processing/tests/hi/conftest.py @@ -28,6 +28,13 @@ def hi_test_background_config_path(hi_l1_test_data_path): return hi_l1_test_data_path / "imap_hi_90sensor-backgrounds_20240101_v001.csv" +@pytest.fixture(scope="session") +def hi_test_gain_configuration_path(hi_l1_test_data_path): + return ( + hi_l1_test_data_path / "imap_hi_90sensor-gain-configuration_20240101_v001.csv" + ) + + def create_metaevent(esa_step, met_subseconds, met_seconds): start_bitmask_data = 0 # META return ( diff --git a/imap_processing/tests/hi/test_hi_l1b.py b/imap_processing/tests/hi/test_hi_l1b.py index 91de0a9cb..39e9179d1 100644 --- a/imap_processing/tests/hi/test_hi_l1b.py +++ b/imap_processing/tests/hi/test_hi_l1b.py @@ -11,23 +11,84 @@ from imap_processing.hi.hi_l1b import ( annotate_direct_events, any_good_direct_events, + classify_gain_configuration, compute_coincidence_type_and_tofs, compute_hae_coordinates, de_ccsds_qf, de_esa_energy_step, de_esa_step_met, + de_gain_configuration, de_nominal_bin_and_spin_phase, get_esa_to_esa_energy_step_lut, + get_gain_configuration_lut, housekeeping, ) from imap_processing.hi.utils import ( CoincidenceBitmap, EsaEnergyStepLookupTable, + GainConfigLookupTable, HiConstants, ) from imap_processing.quality_flags import ImapHiL1bDeFlags from imap_processing.spice.geometry import SpiceFrame +# Nominal detector voltage config (config_id=0) matching the real +# gain-configuration ancillary file, for use as a permissive default in tests +# that don't care about gain-configuration classification. +NOMINAL_HV_VALUES = { + "pos_defl": 6300.0, + "neg_defl": -6300.0, + "tof": -8000.0, + "mcp_f": -3000.0, + "mcp_b": -2125.0, + "cem_f": -4500.0, + "cem_bk_a": -2350.0, + "cem_bk_b": -2350.0, +} +# Tolerances matching the real gain-configuration ancillary file +NOMINAL_HV_TOLERANCES = { + "pos_defl": 700.0, + "neg_defl": 700.0, + "tof": 50.0, + "mcp_f": 25.0, + "mcp_b": 50.0, + "cem_f": 50.0, + "cem_bk_a": 100.0, + "cem_bk_b": 100.0, +} + + +def _make_gain_config_df(hv_tolerance=None, n_esa_energy_steps=9, config_id=0): + """Build a gain configuration DataFrame for a single config_id. + + Indexed by (config_id, esa_energy_step), matching the real + gain-configuration ancillary file's long format -- one row per + esa_energy_step, with {field}_v/{field}_delta_v repeated across every row + for the config_id. + + If hv_tolerance is None, realistic per-field tolerances are used + (matching the real ancillary file). Otherwise, hv_tolerance is applied + as a flat tolerance for every field (useful for permissive matching). + """ + if hv_tolerance is None: + tolerances = NOMINAL_HV_TOLERANCES + else: + tolerances = dict.fromkeys(NOMINAL_HV_VALUES, hv_tolerance) + + hv_row = {} + for field, mean_value in NOMINAL_HV_VALUES.items(): + hv_row[f"{field}_v"] = mean_value + hv_row[f"{field}_delta_v"] = tolerances[field] + rows = [ + {"geometric_factor": 0.001 * step, **hv_row} + for step in range(1, n_esa_energy_steps + 1) + ] + index = pd.MultiIndex.from_product( + [[config_id], range(1, n_esa_energy_steps + 1)], + names=["config_id", "esa_energy_step"], + ) + return pd.DataFrame(rows, index=index) + def test_hi_l1b_hk(hi_l0_test_data_path): """Test coverage for imap_processing.hi.hi_l1b() with @@ -41,9 +102,11 @@ def test_hi_l1b_hk(hi_l0_test_data_path): @pytest.mark.external_kernel @pytest.mark.external_test_data +@mock.patch("imap_processing.hi.hi_l1b.get_gain_configuration_lut") @mock.patch("imap_processing.hi.hi_l1b.get_esa_to_esa_energy_step_lut") def test_hi_annotate_direct_events( mock_get_esa_lut, + mock_get_gain_config_lut, hi_l1_test_data_path, use_fake_spin_data_for_time, imap_ena_sim_metakernel, @@ -55,6 +118,14 @@ def test_hi_annotate_direct_events( mock_esa_lut.query.side_effect = lambda a, b: b mock_get_esa_lut.return_value = mock_esa_lut + # Mock the gain config LUT to always classify config_id=0 with no bad + # detector voltage segments. + mock_gain_lut = mock.MagicMock(spec=GainConfigLookupTable()) + mock_gain_lut.query.side_effect = lambda mets: np.zeros( + len(np.atleast_1d(mets)), dtype=np.int64 + ) + mock_get_gain_config_lut.return_value = (mock_gain_lut, 0) + # Start MET time of spin for simulated input data is 482372988 use_fake_spin_data_for_time(482372987.999) l1a_test_file_path = ( @@ -63,12 +134,18 @@ def test_hi_annotate_direct_events( esa_energies_csv = ( hi_l1_test_data_path / "imap_hi_90sensor-esa-energies_20240101_v001.csv" ) + gain_config_csv = ( + hi_l1_test_data_path / "imap_hi_90sensor-gain-configuration_20240101_v001.csv" + ) # Process using test data l1a_dataset = load_cdf(l1a_test_file_path) - l1b_datasets = annotate_direct_events(l1a_dataset, xr.Dataset(), esa_energies_csv) + l1b_datasets = annotate_direct_events( + l1a_dataset, xr.Dataset(), esa_energies_csv, gain_config_csv + ) assert len(l1b_datasets) == 1 assert l1b_datasets[0].attrs["Logical_source"] == "imap_hi_l1b_45sensor-de" + assert l1b_datasets[0].attrs["gain_configuration_id"] == 0 assert len(l1b_datasets[0].data_vars) == 18 @@ -117,6 +194,9 @@ def test_annotate_direct_events_with_hk( esa_energies_csv = ( hi_l1_test_data_path / "imap_hi_90sensor-esa-energies_20240101_v001.csv" ) + gain_config_csv = ( + hi_l1_test_data_path / "imap_hi_90sensor-gain-configuration_20240101_v001.csv" + ) # Process using test data l1a_dataset = load_cdf(l1a_de_file_path) hk_dataset = load_cdf(l1b_hk_file_path) @@ -128,7 +208,9 @@ def test_annotate_direct_events_with_hk( ) use_fake_spin_data_for_time(spin_start_met) - l1b_datasets = annotate_direct_events(l1a_dataset, hk_dataset, esa_energies_csv) + l1b_datasets = annotate_direct_events( + l1a_dataset, hk_dataset, esa_energies_csv, gain_config_csv + ) assert len(l1b_datasets) == 1 assert l1b_datasets[0].attrs["Logical_source"] == "imap_hi_l1b_90sensor-de" assert len(l1b_datasets[0].data_vars) == 18 @@ -340,8 +422,12 @@ def side_effect_func(et, inst_frame: SpiceFrame, to_frame): @mock.patch("imap_processing.hi.hi_l1b.get_esa_to_esa_energy_step_lut") def test_de_esa_energy_step(mock_get_esa_lut, mock_read_csv, mock_any_good_de): """Test coverage for de_esa_energy_step function.""" + esa_energy_step_fillval = 255 + # Packet at index 5 fails to find a matching esa_energy_step (FILLVAL). mock_esa_lut = mock.MagicMock(spec=EsaEnergyStepLookupTable()) - mock_esa_lut.query.side_effect = lambda a, b: np.arange(len(a))[::-1] % 9 + mock_esa_lut.query.side_effect = lambda a, b: np.where( + np.arange(len(a)) == 5, esa_energy_step_fillval, np.arange(len(a))[::-1] % 9 + ) mock_get_esa_lut.return_value = mock_esa_lut n_epoch = 20 @@ -352,13 +438,80 @@ def test_de_esa_energy_step(mock_get_esa_lut, mock_read_csv, mock_any_good_de): data_vars={ "ccsds_met": xr.DataArray(np.arange(n_epoch) % 9, dims=["epoch"]), "esa_step": xr.DataArray(np.arange(n_epoch), dims=["epoch"]), + # Pre-existing "ccsds_qf", as created by de_ccsds_qf(). + "ccsds_qf": xr.DataArray(np.zeros(n_epoch, dtype=np.uint8), dims=["epoch"]), + }, + ) + new_vars = de_esa_energy_step(fake_dataset, xr.Dataset(), "Fake path") + + expected_esa_energy_step = np.arange(n_epoch)[::-1] % 9 + expected_esa_energy_step[5] = esa_energy_step_fillval + np.testing.assert_array_equal( + new_vars["esa_energy_step"].values, expected_esa_energy_step + ) + # de_esa_energy_step modifies "ccsds_qf" on fake_dataset in place. + expected_qf_bits = np.zeros(n_epoch, dtype=np.uint8) + expected_qf_bits[5] = ImapHiL1bDeFlags.BAD_ESA_VOLTAGE + np.testing.assert_array_equal(fake_dataset["ccsds_qf"].values, expected_qf_bits) + + +@mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=True) +@mock.patch("imap_processing.hi.hi_l1b.load_gain_configuration") +@mock.patch("imap_processing.hi.hi_l1b.get_gain_configuration_lut") +def test_de_gain_configuration( + mock_get_gain_config_lut, mock_load_gain_config, mock_any_good_de +): + """Test coverage for de_gain_configuration function.""" + # Mock the gain config LUT to classify config_id=0, with the packet at + # ccsds_met == 5 failing to match (bad detector voltage). + mock_gain_lut = mock.MagicMock(spec=GainConfigLookupTable()) + mock_gain_lut.query.side_effect = lambda mets: np.where( + np.atleast_1d(mets) == 5, GainConfigLookupTable.NO_MATCH, 0 + ) + mock_get_gain_config_lut.return_value = (mock_gain_lut, 0) + mock_load_gain_config.return_value = _make_gain_config_df() + + n_epoch = 10 + esa_energy_step_fillval = 255 + fake_dataset = xr.Dataset( + coords={ + "epoch": xr.DataArray(np.arange(n_epoch), name="epoch", dims=["epoch"]) + }, + data_vars={ + "ccsds_met": xr.DataArray(np.arange(n_epoch), dims=["epoch"]), + "esa_energy_step": xr.DataArray( + np.full(n_epoch, 3, dtype=np.uint8), + dims=["epoch"], + attrs={"FILLVAL": esa_energy_step_fillval}, + ), + "ccsds_qf": xr.DataArray(np.zeros(n_epoch, dtype=np.uint8), dims=["epoch"]), }, ) - esa_energy_step_var = de_esa_energy_step(fake_dataset, xr.Dataset(), "Fake path") + result = de_gain_configuration(fake_dataset, xr.Dataset(), "Fake gain config path") + + # de_gain_configuration modifies fake_dataset in place and returns it. + assert result is fake_dataset + assert result.attrs["gain_configuration_id"] == 0 + expected_esa_energy_step = np.full(n_epoch, 3) + expected_esa_energy_step[5] = esa_energy_step_fillval np.testing.assert_array_equal( - esa_energy_step_var["esa_energy_step"].values, np.arange(n_epoch)[::-1] % 9 + result["esa_energy_step"].values, expected_esa_energy_step ) + expected_qf_bits = np.zeros(n_epoch, dtype=np.uint8) + expected_qf_bits[5] = ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE + np.testing.assert_array_equal(result["ccsds_qf"].values, expected_qf_bits) + + +@mock.patch("imap_processing.hi.hi_l1b.any_good_direct_events", return_value=False) +def test_de_gain_configuration_no_good_events(mock_any_good_de): + """Test that gain_configuration_id is NO_MATCH when there are no good DEs.""" + fake_dataset = xr.Dataset(attrs={}) + + result = de_gain_configuration(fake_dataset, xr.Dataset(), "Fake gain config path") + + assert result is fake_dataset + assert result.attrs["gain_configuration_id"] == GainConfigLookupTable.NO_MATCH class TestGetEsaToEsaEnergyStepLut: @@ -389,19 +542,24 @@ def create_mock_dataset( inner_esa_lo, outer_esa_values, shcoarse_values, + **hv_overrides, ): """Helper method to create mock L1B housekeeping dataset.""" - return xr.Dataset( - { - "op_mode": (["epoch"], op_modes), - "sci_esa_step": (["epoch"], esa_steps), - "inner_esa_state": (["epoch"], inner_esa_state), - "inner_esa_hi": (["epoch"], inner_esa_hi), - "inner_esa_lo": (["epoch"], inner_esa_lo), - "outer_esa": (["epoch"], outer_esa_values), - "shcoarse": (["epoch"], shcoarse_values), - } - ) + n = len(op_modes) + hv_values = dict(NOMINAL_HV_VALUES) + hv_values.update(hv_overrides) + data_vars = { + "op_mode": (["epoch"], op_modes), + "sci_esa_step": (["epoch"], esa_steps), + "inner_esa_state": (["epoch"], inner_esa_state), + "inner_esa_hi": (["epoch"], inner_esa_hi), + "inner_esa_lo": (["epoch"], inner_esa_lo), + "outer_esa": (["epoch"], outer_esa_values), + "shcoarse": (["epoch"], shcoarse_values), + } + for field, value in hv_values.items(): + data_vars[field] = (["epoch"], np.full(n, value)) + return xr.Dataset(data_vars) @mock.patch("imap_processing.hi.hi_l1b.EsaEnergyStepLookupTable") def test_basic_functionality_single_hvsci_segment(self, mock_lut_class): @@ -424,7 +582,7 @@ def test_basic_functionality_single_hvsci_segment(self, mock_lut_class): shcoarse_values=[1000, 1001, 1002, 1003], ) - result = get_esa_to_esa_energy_step_lut(l1b_hk_ds, self.esa_energies_lut) + lut = get_esa_to_esa_energy_step_lut(l1b_hk_ds, self.esa_energies_lut) # Verify LUT was instantiated mock_lut_class.assert_called_once() @@ -446,7 +604,7 @@ def test_basic_functionality_single_hvsci_segment(self, mock_lut_class): # Second call should be for esa_step 2 assert calls[1][0] == (1000, 1003, 2, 2) - assert result == self.mock_lut + assert lut == self.mock_lut @mock.patch("imap_processing.hi.hi_l1b.EsaEnergyStepLookupTable") def test_multiple_hvsci_segments(self, mock_lut_class): @@ -489,11 +647,11 @@ def test_no_hvsci_segments(self, mock_lut_class): shcoarse_values=[1000, 1001, 1002], ) - result = get_esa_to_esa_energy_step_lut(l1b_hk_ds, self.esa_energies_lut) + lut = get_esa_to_esa_energy_step_lut(l1b_hk_ds, self.esa_energies_lut) # No add_entry calls should be made self.mock_lut.add_entry.assert_not_called() - assert result == self.mock_lut + assert lut == self.mock_lut @mock.patch("imap_processing.hi.hi_l1b.EsaEnergyStepLookupTable") @mock.patch("imap_processing.hi.hi_l1b.logger") @@ -648,6 +806,133 @@ def test_cal_data(self, hi_l1_test_data_path): np.testing.assert_array_equal(lut.df["esa_step"].values, np.arange(9) + 1) +class TestGetGainConfigurationLut: + """Test suite for get_gain_configuration_lut function.""" + + def _make_hk_ds(self, op_modes, shcoarse_values=None, **hv_overrides): + n = len(op_modes) + hv_values = dict(NOMINAL_HV_VALUES) + hv_values.update(hv_overrides) + if shcoarse_values is None: + shcoarse_values = np.arange(n, dtype=float) + 1000.0 + data_vars = { + "op_mode": (["epoch"], op_modes), + "shcoarse": (["epoch"], np.asarray(shcoarse_values, dtype=float)), + } + for field, value in hv_values.items(): + data_vars[field] = (["epoch"], np.full(n, value)) + return xr.Dataset(data_vars) + + def test_matching_segment_is_added_to_lut(self): + """Test that a segment matching the classified config gets a LUT entry.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI", "HVSCI"]) + + lut, config_id = get_gain_configuration_lut(hk_ds, gain_config_df) + + assert config_id == 0 + assert lut.query(1000.0) == 0 + assert lut.query(1002.0) == 0 + assert lut.query(999.0) == GainConfigLookupTable.NO_MATCH + + def test_mismatched_segment_excluded_from_lut(self): + """Test that a segment not matching the classified config is excluded.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds( + ["HVSCI", "HVSCI", "OTHER", "HVSCI", "HVSCI"], + shcoarse_values=[1000, 1001, 1002, 1003, 1004], + ) + # Second HVSCI segment (indices 3-4) has a mismatched CEM front voltage, + # simulating a mid-pointing gain test excursion. + hk_ds["cem_f"].values[3:] = -4200.0 + + lut, config_id = get_gain_configuration_lut(hk_ds, gain_config_df) + + # First segment classifies the pointing's gain configuration (config 0) + assert config_id == 0 + # Only the first (matching) segment should have a LUT entry + assert lut.query(1000.0) == 0 + assert lut.query(1001.0) == 0 + assert lut.query(1003.0) == GainConfigLookupTable.NO_MATCH + assert lut.query(1004.0) == GainConfigLookupTable.NO_MATCH + + def test_no_hvsci_segments_returns_no_match_everywhere(self): + """Test that no HVSCI segments results in an empty LUT and config_id=None.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["OTHER", "LVSCI"]) + + lut, config_id = get_gain_configuration_lut(hk_ds, gain_config_df) + + assert config_id is None + assert lut.query(1000.0) == GainConfigLookupTable.NO_MATCH + + def test_classification_failure_returns_no_match_everywhere(self): + """Test that an unclassifiable pointing produces an empty LUT.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI"], cem_f=0.0) + + lut, config_id = get_gain_configuration_lut(hk_ds, gain_config_df) + + assert config_id is None + assert lut.query(1000.0) == GainConfigLookupTable.NO_MATCH + + +class TestClassifyGainConfiguration: + """Test suite for classify_gain_configuration function.""" + + def _make_hk_ds(self, op_modes, **hv_overrides): + n = len(op_modes) + hv_values = dict(NOMINAL_HV_VALUES) + hv_values.update(hv_overrides) + data_vars = { + "op_mode": (["epoch"], op_modes), + "shcoarse": (["epoch"], np.arange(n, dtype=float) + 1000.0), + } + for field, value in hv_values.items(): + data_vars[field] = (["epoch"], np.full(n, value)) + return xr.Dataset(data_vars) + + def test_classifies_nominal_config(self): + """Test that nominal voltages classify as config_id 0.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI", "HVSCI"]) + assert classify_gain_configuration(hk_ds, gain_config_df) == 0 + + def test_classifies_second_config(self): + """Test that a second configuration's voltages are classified correctly.""" + gain_config_df = pd.concat( + [_make_gain_config_df(config_id=0), _make_gain_config_df(config_id=1)] + ) + # Give config 1 a distinctly different CEM front voltage + gain_config_df.loc[1, "cem_f_v"] = -4200.0 + gain_config_df.loc[1, "cem_f_delta_v"] = 25.0 + + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI"], cem_f=-4200.0) + assert classify_gain_configuration(hk_ds, gain_config_df) == 1 + + def test_no_hvsci_segments_returns_none(self): + """Test that no HVSCI segments results in None.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["OTHER", "LVSCI"]) + assert classify_gain_configuration(hk_ds, gain_config_df) is None + + def test_no_match_returns_none(self): + """Test that voltages matching no configuration return None.""" + gain_config_df = _make_gain_config_df() + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI"], cem_f=0.0) + assert classify_gain_configuration(hk_ds, gain_config_df) is None + + def test_only_uses_first_hvsci_segment(self): + """Test that only the first HVSCI segment is used for classification.""" + gain_config_df = _make_gain_config_df() + # First segment matches config 0; a later segment (e.g. a mid-pointing + # gain test) does not match anything -- this should not affect the + # classification result, which is based only on the first segment. + hk_ds = self._make_hk_ds(["HVSCI", "HVSCI", "OTHER", "HVSCI", "HVSCI"]) + hk_ds["cem_f"].values[3:] = 0.0 + assert classify_gain_configuration(hk_ds, gain_config_df) == 0 + + class TestDeEsaStepMet: """Tests for de_esa_step_met function.""" diff --git a/imap_processing/tests/hi/test_utils.py b/imap_processing/tests/hi/test_utils.py index 78abfb560..f0db51aa5 100644 --- a/imap_processing/tests/hi/test_utils.py +++ b/imap_processing/tests/hi/test_utils.py @@ -15,6 +15,7 @@ CalibrationProductConfig, CoincidenceBitmap, EsaEnergyStepLookupTable, + HiConstants, compute_qualified_event_mask, create_dataset_variables, filter_events_by_coincidence, @@ -23,9 +24,52 @@ get_tof_window_mask, iter_background_events_by_config, iter_qualified_events_by_config, + load_gain_configuration, parse_sensor_number, ) +# Nominal HV mean/tolerance values matching the real gain-configuration +# ancillary file, for building minimal test CSVs. +_GAIN_CONFIG_HV_MEAN = { + "pos_defl": 6300.0, + "neg_defl": -6300.0, + "tof": -8000.0, + "mcp_f": -3000.0, + "mcp_b": -2125.0, + "cem_f": -4500.0, + "cem_bk_a": -2350.0, + "cem_bk_b": -2350.0, +} +_GAIN_CONFIG_HV_DELTA = dict.fromkeys(HiConstants.GAIN_CONFIG_HV_FIELDS, 50.0) + + +def _gain_config_csv_rows(second_row_overrides=None): + """Build (header, first_row, second_row) CSV lines for config_id=0. + + first_row (esa_energy_step=1) specifies all {field}_v/{field}_delta_v + values. second_row (esa_energy_step=2) leaves them blank (to be + forward-filled), unless overridden via second_row_overrides (column name + -> value string). + """ + hv_columns = [ + f"{field}{suffix}" + for field in HiConstants.GAIN_CONFIG_HV_FIELDS + for suffix in ("_v", "_delta_v") + ] + header = "config_id,esa_energy_step,geometric_factor," + ",".join(hv_columns) + + first_values = { + **{f"{f}_v": str(v) for f, v in _GAIN_CONFIG_HV_MEAN.items()}, + **{f"{f}_delta_v": str(v) for f, v in _GAIN_CONFIG_HV_DELTA.items()}, + } + first_row = "0,1,0.00055," + ",".join(first_values[col] for col in hv_columns) + + overrides = second_row_overrides or {} + second_values = {col: overrides.get(col, "") for col in hv_columns} + second_row = "0,2,0.00085," + ",".join(second_values[col] for col in hv_columns) + + return header, first_row, second_row + def test_hiapid(): """Test coverage for HIAPID class""" @@ -422,6 +466,66 @@ def test_calibration_product_numbers_arbitrary_values(self): assert isinstance(cal_prod_numbers, np.ndarray) +class TestLoadGainConfiguration: + """Test coverage for load_gain_configuration function.""" + + def test_from_csv(self, hi_test_gain_configuration_path): + """Test coverage for load_gain_configuration reading the ancillary CSV.""" + gain_config_df = load_gain_configuration(hi_test_gain_configuration_path) + + assert gain_config_df.index.names == ["config_id", "esa_energy_step"] + assert len(gain_config_df) == 18 # 2 configs x 9 esa_energy_steps + assert gain_config_df.loc[(0, 1), "geometric_factor"] == 0.00055 + assert gain_config_df.loc[(0, 5), "geometric_factor"] == 0.00340 + assert gain_config_df.loc[(1, 9), "geometric_factor"] == 0.01922 + # {field}_v/{field}_delta_v should be forward-filled to every + # esa_energy_step + for esa_energy_step in range(1, 10): + assert gain_config_df.loc[(0, esa_energy_step), "pos_defl_v"] == 6300.0 + assert gain_config_df.loc[(0, esa_energy_step), "pos_defl_delta_v"] == 700.0 + assert gain_config_df.loc[(0, esa_energy_step), "tof_v"] == -8000.0 + assert gain_config_df.loc[(0, esa_energy_step), "cem_f_v"] == -4500.0 + assert gain_config_df.loc[(1, esa_energy_step), "cem_f_v"] == -4200.0 + + def test_from_csv_string_io(self): + """Test coverage for load_gain_configuration with a file-like object.""" + header, first_row, second_row = _gain_config_csv_rows() + csv_content = ( + f"# comment line should be ignored\n{header}\n{first_row}\n{second_row}\n" + ) + gain_config_df = load_gain_configuration(io.StringIO(csv_content)) + assert gain_config_df.loc[(0, 1), "geometric_factor"] == 0.00055 + assert gain_config_df.loc[(0, 2), "geometric_factor"] == 0.00085 + # Forward-filled from the first row + assert gain_config_df.loc[(0, 2), "pos_defl_v"] == 6300.0 + assert gain_config_df.loc[(0, 2), "pos_defl_delta_v"] == 50.0 + + def test_missing_geometric_factor_raises(self): + """Test that a null geometric_factor raises ValueError.""" + header, first_row, _ = _gain_config_csv_rows() + second_row = "0,2,,,,,,,,,,,,,,,,," + csv_content = f"{header}\n{first_row}\n{second_row}\n" + with pytest.raises(ValueError, match="Null geometric_factor"): + load_gain_configuration(io.StringIO(csv_content)) + + def test_missing_first_row_hv_values_raises(self): + """Test that missing {field}_v/{field}_delta_v on the first row raises.""" + header, _, second_row = _gain_config_csv_rows() + first_row = "0,1,0.00055,,,,,,,,,,,,,,,," + csv_content = f"{header}\n{first_row}\n{second_row}\n" + with pytest.raises(ValueError, match="Missing pos_defl_v"): + load_gain_configuration(io.StringIO(csv_content)) + + def test_inconsistent_hv_values_raises(self): + """Test inconsistent {field}_v/{field}_delta_v within a config_id raises.""" + header, first_row, second_row = _gain_config_csv_rows( + second_row_overrides={"pos_defl_v": "6301.0"} + ) + csv_content = f"{header}\n{first_row}\n{second_row}\n" + with pytest.raises(ValueError, match="Inconsistent pos_defl_v"): + load_gain_configuration(io.StringIO(csv_content)) + + class TestBackgroundConfig: """ All test coverage for the pd.DataFrame accessor extension "background_config". From 11a5a78df1678e2ed9f6760300851fb2ba1bcc9a Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 14:57:52 -0600 Subject: [PATCH 4/8] Update hi goodtimes and hi_l1c to use L1B gain configuration and identification of gain tests --- imap_processing/hi/hi_goodtimes.py | 40 ++++++----- imap_processing/hi/hi_l1c.py | 69 ++++++++++++++++++- imap_processing/tests/hi/test_hi_goodtimes.py | 36 +++++----- imap_processing/tests/hi/test_hi_l1c.py | 62 ++++++++++++++++- 4 files changed, 170 insertions(+), 37 deletions(-) diff --git a/imap_processing/hi/hi_goodtimes.py b/imap_processing/hi/hi_goodtimes.py index 582151886..d67891540 100644 --- a/imap_processing/hi/hi_goodtimes.py +++ b/imap_processing/hi/hi_goodtimes.py @@ -53,7 +53,7 @@ class CullCode(IntEnum): STAT_FILTER_0 = 1 << 4 # 16 STAT_FILTER_1 = 1 << 5 # 32 STAT_FILTER_2 = 1 << 6 # 64 - BAD_ESA_VOLTAGE = 1 << 7 # 128 + BAD_HV_VALUE = 1 << 7 # 128 def hi_goodtimes( @@ -69,7 +69,8 @@ def hi_goodtimes( This is the top-level function that orchestrates all goodtimes culling operations for a single pointing. It applies the following filters in order: - 0. mark_bad_esa_voltage - Remove times with invalid ESA voltage configuration + 0. mark_bad_voltage - Remove times with invalid ESA or detector HV + configuration (e.g. gain test intervals) 1. mark_incomplete_spin_sets - Remove incomplete 8-spin histogram periods 2. mark_drf_times - Remove times during spacecraft drift restabilization 3. mark_bad_tdc_cal - Remove times with failed TDC calibration @@ -276,9 +277,9 @@ def _apply_goodtimes_filters( # === Apply culling filters === - # 0. Mark bad ESA voltage times - logger.info("Applying filter: mark_bad_esa_voltage") - mark_bad_esa_voltage(goodtimes_ds, current_l1b_de) + # 0. Mark bad ESA/detector HV voltage times + logger.info("Applying filter: mark_bad_voltage") + mark_bad_voltage(goodtimes_ds, current_l1b_de) # 1. Mark incomplete spin sets logger.info("Applying filter: mark_incomplete_spin_sets") @@ -953,22 +954,28 @@ def finalize_dataset(self) -> xr.Dataset: # ============================================================================== -def mark_bad_esa_voltage( +def mark_bad_voltage( goodtimes_ds: xr.Dataset, l1b_de: xr.Dataset, - cull_code: int = CullCode.BAD_ESA_VOLTAGE, + cull_code: int = CullCode.BAD_HV_VALUE, ) -> None: """ - Mark times when ESA voltages don't match expected values. + Mark times when ESA or detector high voltages don't match expected values. Filters out 8-spin periods where the ESA energy step is invalid, indicating - either calibration mode (esa_energy_step=0) or an ESA voltage mismatch + either calibration mode (esa_energy_step=0) or a voltage mismatch (esa_energy_step=FILLVAL). The voltage validation is performed during L1B - processing by matching measured inner/outer ESA voltages against the ESA - energies lookup table. + processing (see hi_l1b.get_esa_to_esa_energy_step_lut()) by matching + measured inner/outer ESA voltages against the ESA energies lookup table, + and matching measured detector high voltages (CEM/MCP/U-Can/deflector) + against the pointing's classified gain configuration. Both checks set + esa_energy_step=FILLVAL on mismatch; see the L1B DE ccsds_qf quality flag + (ImapHiL1bDeFlags.BAD_ESA_VOLTAGE / BAD_DETECTOR_VOLTAGE) to distinguish + which check failed. Algorithm Document Reference: - Section 2.3.2: Good times selection requiring valid ESA configuration + Section 2.2.4/2.3.2: Good times selection requiring valid ESA and + detector high voltage configuration. Parameters ---------- @@ -977,7 +984,7 @@ def mark_bad_esa_voltage( l1b_de : xarray.Dataset L1B Direct Event data containing esa_energy_step field. cull_code : int, optional - Cull code to use for marking bad times (default: CullCode.BAD_ESA_VOLTAGE). + Cull code to use for marking bad times (default: CullCode.BAD_HV_VALUE). Notes ----- @@ -986,10 +993,11 @@ def mark_bad_esa_voltage( Invalid ESA energy steps: - esa_energy_step = 0: Calibration mode (ESA stepping but not science data) - - esa_energy_step = FILLVAL (255): ESA voltage mismatch - measured voltages - didn't match any known energy step in the lookup table + - esa_energy_step = FILLVAL (255): ESA or detector voltage mismatch - + measured voltages didn't match any known energy step / gain + configuration """ - logger.info("Running mark_bad_esa_voltage culling") + logger.info("Running mark_bad_voltage culling") # Get FILLVAL from attributes (should be 255 for uint8) fillval = l1b_de["esa_energy_step"].attrs.get("FILLVAL", 255) diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 8884e3e74..7bf0d6527 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -15,11 +15,13 @@ from imap_processing.hi.utils import ( BackgroundConfig, CalibrationProductConfig, + GainConfigLookupTable, HiConstants, create_dataset_variables, full_dataarray, iter_background_events_by_config, iter_qualified_events_by_config, + load_gain_configuration, parse_sensor_number, ) from imap_processing.spice.geometry import ( @@ -47,6 +49,7 @@ def hi_l1c( calibration_prod_config_path: Path, goodtimes_ds: xr.Dataset, background_config_path: Path, + gain_config_anc: Path, ) -> list[xr.Dataset]: """ High level IMAP-Hi l1c processing function. @@ -61,6 +64,8 @@ def hi_l1c( Goodtimes dataset with cull_flags. background_config_path : pathlib.Path Background configuration file. + gain_config_anc : pathlib.Path + Location of the gain-configuration ancillary csv file. Returns ------- @@ -70,7 +75,11 @@ def hi_l1c( logger.info("Running Hi l1c processing") l1c_dataset = generate_pset_dataset( - de_dataset, calibration_prod_config_path, goodtimes_ds, background_config_path + de_dataset, + calibration_prod_config_path, + goodtimes_ds, + background_config_path, + gain_config_anc, ) return [l1c_dataset] @@ -81,6 +90,7 @@ def generate_pset_dataset( calibration_prod_config_path: Path, goodtimes_ds: xr.Dataset, background_config_path: Path, + gain_config_anc: Path, ) -> xr.Dataset: """ Generate IMAP-Hi l1c pset xarray dataset from l1b product. @@ -95,6 +105,8 @@ def generate_pset_dataset( Goodtimes dataset with cull_flags. background_config_path : pathlib.Path Background configuration file. + gain_config_anc : pathlib.Path + Location of the gain-configuration ancillary csv file. Returns ------- @@ -123,6 +135,11 @@ 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's + # gain configuration (classified at L1B) + pset_dataset.update( + pset_geometric_factor(pset_dataset.coords, de_dataset, gain_config_anc) + ) # Bin the counts into the spin-bins pset_dataset.update( pset_counts(pset_dataset.coords, config_df, de_dataset, goodtimes_ds) @@ -343,6 +360,56 @@ 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], + l1b_de_dataset: xr.Dataset, + gain_config_anc: Path, +) -> dict[str, xr.DataArray]: + """ + Look up the per-esa_energy_step geometric factor for this pointing. + + 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 a + "gain_configuration_id" global attribute. + gain_config_anc : pathlib.Path + Location of the gain-configuration ancillary csv file. + + Returns + ------- + dict[str, xarray.DataArray] + Dictionary containing the "geometric_factor" DataArray, dims + (epoch, esa_energy_step). + + Notes + ----- + The gain configuration is constant for a pointing (see + `hi_l1b.classify_gain_configuration`), so the L1B DE product only records + the classified gain configuration id as a global attribute rather than + duplicating the geometric factor across every direct event. This looks + up the geometric factor for each esa_energy_step directly from the + gain-configuration ancillary file. Not yet consumed by L2 processing + (deferred to a follow-on ticket that handles combining PSETs from + different gain configurations into a single map). + """ + geometric_factor_var = create_dataset_variables( + ["geometric_factor"], + coords=pset_coords, + att_manager_lookup_str="hi_pset_{0}", + ) + config_id = l1b_de_dataset.attrs["gain_configuration_id"] + if config_id != GainConfigLookupTable.NO_MATCH: + gain_config_df = load_gain_configuration(gain_config_anc) + for i, step in enumerate(pset_coords["esa_energy_step"].data): + geometric_factor_var["geometric_factor"].values[0, i] = gain_config_df.loc[ + (config_id, int(step)), "geometric_factor" + ] + return geometric_factor_var + + def pset_counts( pset_coords: dict[str, xr.DataArray], config_df: pd.DataFrame, diff --git a/imap_processing/tests/hi/test_hi_goodtimes.py b/imap_processing/tests/hi/test_hi_goodtimes.py index 67f815d3e..a24e44e1c 100644 --- a/imap_processing/tests/hi/test_hi_goodtimes.py +++ b/imap_processing/tests/hi/test_hi_goodtimes.py @@ -24,8 +24,8 @@ _sum_esa_counts, create_goodtimes_dataset, hi_goodtimes, - mark_bad_esa_voltage, mark_bad_tdc_cal, + mark_bad_voltage, mark_drf_times, mark_incomplete_spin_sets, mark_overflow_packets, @@ -88,7 +88,7 @@ def test_cull_code_values(self): assert CullCode.STAT_FILTER_0 == 16 assert CullCode.STAT_FILTER_1 == 32 assert CullCode.STAT_FILTER_2 == 64 - assert CullCode.BAD_ESA_VOLTAGE == 128 + assert CullCode.BAD_HV_VALUE == 128 def test_cull_code_is_int(self): """Test that CullCode values are integers.""" @@ -4199,7 +4199,7 @@ def test_loads_cal_config(self, tmp_path): patch( "imap_processing.hi.utils.CalibrationProductConfig.from_csv" ) as mock_cal_load, - patch("imap_processing.hi.hi_goodtimes.mark_bad_esa_voltage"), + patch("imap_processing.hi.hi_goodtimes.mark_bad_voltage"), patch("imap_processing.hi.hi_goodtimes.mark_incomplete_spin_sets"), patch("imap_processing.hi.hi_goodtimes.mark_drf_times"), patch("imap_processing.hi.hi_goodtimes.mark_overflow_packets"), @@ -4237,7 +4237,7 @@ def test_calls_all_filters(self, tmp_path): "imap_processing.hi.utils.CalibrationProductConfig.from_csv", return_value=mock_cal, ), - patch("imap_processing.hi.hi_goodtimes.mark_bad_esa_voltage") as mock_f0, + patch("imap_processing.hi.hi_goodtimes.mark_bad_voltage") as mock_f0, patch( "imap_processing.hi.hi_goodtimes.mark_incomplete_spin_sets" ) as mock_f1, @@ -4288,7 +4288,7 @@ def test_raises_statistical_filter_0_errors(self, tmp_path): "imap_processing.hi.utils.CalibrationProductConfig.from_csv", return_value=mock_cal, ), - patch("imap_processing.hi.hi_goodtimes.mark_bad_esa_voltage"), + patch("imap_processing.hi.hi_goodtimes.mark_bad_voltage"), patch("imap_processing.hi.hi_goodtimes.mark_incomplete_spin_sets"), patch("imap_processing.hi.hi_goodtimes.mark_drf_times"), patch("imap_processing.hi.hi_goodtimes.mark_bad_tdc_cal"), @@ -4324,7 +4324,7 @@ def test_raises_statistical_filter_1_errors(self, tmp_path): "imap_processing.hi.utils.CalibrationProductConfig.from_csv", return_value=mock_cal, ), - patch("imap_processing.hi.hi_goodtimes.mark_bad_esa_voltage"), + patch("imap_processing.hi.hi_goodtimes.mark_bad_voltage"), patch("imap_processing.hi.hi_goodtimes.mark_incomplete_spin_sets"), patch("imap_processing.hi.hi_goodtimes.mark_drf_times"), patch("imap_processing.hi.hi_goodtimes.mark_bad_tdc_cal"), @@ -4550,8 +4550,8 @@ def test_returns_datasets(self, tmp_path): assert result == [mock_finalized] -class TestMarkBadEsaVoltage: - """Tests for mark_bad_esa_voltage culling function.""" +class TestMarkBadVoltage: + """Tests for mark_bad_voltage culling function.""" @pytest.fixture def goodtimes_for_esa(self): @@ -4609,41 +4609,41 @@ def l1b_de_with_fillval(self): ds["esa_energy_step"].attrs["FILLVAL"] = 255 return ds - def test_mark_bad_esa_voltage_all_valid(self, goodtimes_for_esa, l1b_de_all_valid): + def test_mark_bad_voltage_all_valid(self, goodtimes_for_esa, l1b_de_all_valid): """Test that no times are marked when all ESA energy steps are valid.""" - mark_bad_esa_voltage(goodtimes_for_esa, l1b_de_all_valid) + mark_bad_voltage(goodtimes_for_esa, l1b_de_all_valid) assert np.all(goodtimes_for_esa["cull_flags"].values == CullCode.GOOD) - def test_mark_bad_esa_voltage_with_zero(self, goodtimes_for_esa, l1b_de_with_zero): + def test_mark_bad_voltage_with_zero(self, goodtimes_for_esa, l1b_de_with_zero): """Test that times are marked when esa_energy_step=0 (calibration).""" - mark_bad_esa_voltage(goodtimes_for_esa, l1b_de_with_zero) + mark_bad_voltage(goodtimes_for_esa, l1b_de_with_zero) # MET 1050 (index 1) should be culled assert np.all( - goodtimes_for_esa["cull_flags"].values[1, :] == CullCode.BAD_ESA_VOLTAGE + goodtimes_for_esa["cull_flags"].values[1, :] == CullCode.BAD_HV_VALUE ) # Other times should remain good assert np.all(goodtimes_for_esa["cull_flags"].values[0, :] == CullCode.GOOD) assert np.all(goodtimes_for_esa["cull_flags"].values[2, :] == CullCode.GOOD) - def test_mark_bad_esa_voltage_with_fillval( + def test_mark_bad_voltage_with_fillval( self, goodtimes_for_esa, l1b_de_with_fillval ): """Test that times are marked when esa_energy_step=FILLVAL.""" - mark_bad_esa_voltage(goodtimes_for_esa, l1b_de_with_fillval) + mark_bad_voltage(goodtimes_for_esa, l1b_de_with_fillval) # MET 1100 (index 2) should be culled assert np.all( - goodtimes_for_esa["cull_flags"].values[2, :] == CullCode.BAD_ESA_VOLTAGE + goodtimes_for_esa["cull_flags"].values[2, :] == CullCode.BAD_HV_VALUE ) # Other times should remain good assert np.all(goodtimes_for_esa["cull_flags"].values[0, :] == CullCode.GOOD) assert np.all(goodtimes_for_esa["cull_flags"].values[1, :] == CullCode.GOOD) - def test_mark_bad_esa_voltage_custom_cull_code( + def test_mark_bad_voltage_custom_cull_code( self, goodtimes_for_esa, l1b_de_with_zero ): """Test using a custom cull code.""" custom_code = 200 - mark_bad_esa_voltage(goodtimes_for_esa, l1b_de_with_zero, cull_code=custom_code) + mark_bad_voltage(goodtimes_for_esa, l1b_de_with_zero, cull_code=custom_code) assert np.all(goodtimes_for_esa["cull_flags"].values[1, :] == custom_code) diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index dae285987..d8b8da118 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -12,7 +12,7 @@ from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import load_cdf, write_cdf from imap_processing.hi import hi_l1c, utils -from imap_processing.hi.utils import HIAPID, HiConstants +from imap_processing.hi.utils import HIAPID, GainConfigLookupTable, HiConstants from imap_processing.spice.time import met_to_ttj2000ns, ttj2000ns_to_et @@ -37,6 +37,7 @@ def test_hi_l1c( mock_generate_pset_dataset, hi_test_cal_prod_config_path, hi_test_background_config_path, + hi_test_gain_configuration_path, ): """Test coverage for hi_l1c function""" mock_generate_pset_dataset.return_value = xr.Dataset() @@ -45,6 +46,7 @@ def test_hi_l1c( hi_test_cal_prod_config_path, xr.Dataset(), hi_test_background_config_path, + hi_test_gain_configuration_path, )[0] # Empty attributes, global values get added in post-processing assert pset.attrs == {} @@ -57,13 +59,18 @@ def test_generate_pset_dataset( hi_goodtimes_dataset, hi_test_cal_prod_config_path, hi_test_background_config_path, + hi_test_gain_configuration_path, use_fake_spin_data_for_time, use_fake_repoint_data_for_time, imap_ena_sim_metakernel, ): """Test coverage for generate_pset_dataset function""" use_fake_spin_data_for_time(482372987.999) - l1b_dataset = hi_l1b_de_dataset + l1b_dataset = hi_l1b_de_dataset.copy() + # The real fixture CDF predates the gain_configuration_id L1B global + # attribute; add a placeholder so pset_geometric_factor() has something + # to look up. + l1b_dataset.attrs["gain_configuration_id"] = 0 l1b_met = l1b_dataset["ccsds_met"].values[0] # Set repoint start and end times. seconds_per_day = 24 * 60 * 60 @@ -78,6 +85,7 @@ def test_generate_pset_dataset( hi_test_cal_prod_config_path, goodtimes, hi_test_background_config_path, + hi_test_gain_configuration_path, ) assert l1c_dataset.epoch.data[0] == l1b_dataset.epoch.data[0].astype(np.int64) @@ -111,6 +119,7 @@ def test_generate_pset_dataset_uses_midpoint_time( mock_pset_backgrounds, hi_test_cal_prod_config_path, hi_test_background_config_path, + hi_test_gain_configuration_path, ): """Test that generate_pset_dataset uses midpoint ET for pset_geometry.""" # Create a mock L1B dataset @@ -132,6 +141,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", + "gain_configuration_id": 0, }, ) @@ -157,6 +167,7 @@ def test_generate_pset_dataset_uses_midpoint_time( hi_test_cal_prod_config_path, xr.Dataset(), hi_test_background_config_path, + hi_test_gain_configuration_path, ) # Calculate expected midpoint ET @@ -179,6 +190,53 @@ def test_generate_pset_dataset_uses_midpoint_time( assert actual_sensor_arg == "45sensor" +@mock.patch("imap_processing.hi.hi_l1c.load_gain_configuration") +def test_pset_geometric_factor_looks_up_by_config_id(mock_load_gain_config): + """Test coverage for pset_geometric_factor when the pointing was classified.""" + gain_config_df = pd.DataFrame( + {"geometric_factor": [0.001, 0.002, 0.003]}, + index=pd.MultiIndex.from_tuples( + [(0, 1), (0, 2), (0, 3)], names=["config_id", "esa_energy_step"] + ), + ) + mock_load_gain_config.return_value = gain_config_df + + pset_coords = { + "epoch": xr.DataArray([0], dims=["epoch"]), + "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), + } + l1b_de_dataset = xr.Dataset(attrs={"gain_configuration_id": 0}) + + result = hi_l1c.pset_geometric_factor( + pset_coords, l1b_de_dataset, "Fake gain config path" + ) + + np.testing.assert_allclose( + result["geometric_factor"].values[0], [0.001, 0.002, 0.003], rtol=1e-6 + ) + mock_load_gain_config.assert_called_once_with("Fake gain config path") + + +def test_pset_geometric_factor_no_match_returns_fillval(): + """Test coverage for pset_geometric_factor when the pointing was not classified.""" + pset_coords = { + "epoch": xr.DataArray([0], dims=["epoch"]), + "esa_energy_step": xr.DataArray([1, 2, 3], dims=["esa_energy_step"]), + } + l1b_de_dataset = xr.Dataset( + attrs={"gain_configuration_id": GainConfigLookupTable.NO_MATCH} + ) + + result = hi_l1c.pset_geometric_factor( + pset_coords, l1b_de_dataset, "Fake gain config path" + ) + + fillval = np.float32(result["geometric_factor"].attrs["FILLVAL"]) + np.testing.assert_array_equal( + result["geometric_factor"].values[0], [fillval, fillval, fillval] + ) + + def test_empty_pset_dataset(use_fake_repoint_data_for_time): """Test coverage for empty_pset_dataset function""" n_energy_steps = 8 From 5d7819ed04ce121c3056b2397f36cf6a87de9b14 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 14:58:50 -0600 Subject: [PATCH 5/8] Add missing test gain configuration csv --- ...ensor-gain-configuration_20240101_v001.csv | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 imap_processing/tests/hi/data/l1/imap_hi_90sensor-gain-configuration_20240101_v001.csv diff --git a/imap_processing/tests/hi/data/l1/imap_hi_90sensor-gain-configuration_20240101_v001.csv b/imap_processing/tests/hi/data/l1/imap_hi_90sensor-gain-configuration_20240101_v001.csv new file mode 100644 index 000000000..8a6b7bd9b --- /dev/null +++ b/imap_processing/tests/hi/data/l1/imap_hi_90sensor-gain-configuration_20240101_v001.csv @@ -0,0 +1,40 @@ +# Gain Configuration table +# Valid start date: 2024-01-01 (from filename) +# This file will be used in processing until a new calibration file is uploaded +# to the SDC with either a higher version number or new date in the filename. +# For details on how the SDC selects ancillary files for processing, see: +# https://imap-processing.readthedocs.io/en/latest/data-access/calibration-files.html +# +# Each row describes one (config_id, esa_energy_step) combination. IMAP-Hi is +# expected to alternate whole pointings between config_id=0 (the current/ +# nominal configuration) and config_id=1 (a new configuration reflecting a +# deliberate gain change), starting August 2026. +# +# {field}_v/{field}_delta_v give the nominal value and tolerance (in volts) +# for each detector high voltage monitor used to classify which configuration +# a pointing is running (see IMAP-Hi Algorithm Document Section 2.2.4 and +# culling_v2.c drop_bad_voltages() for config_id=0's reference values). Note +# "tof" is the raw CCSDS mnemonic name for the U-Can voltage monitor (see +# Section 8, Level 0 Packet Definitions). These values are constant across +# esa_energy_step, so they only need to be specified on the first +# esa_energy_step row of each config_id group; subsequent rows are +# forward-filled. geometric_factor must be specified for every row. +config_id,esa_energy_step,geometric_factor,pos_defl_v,pos_defl_delta_v,neg_defl_v,neg_defl_delta_v,tof_v,tof_delta_v,mcp_f_v,mcp_f_delta_v,mcp_b_v,mcp_b_delta_v,cem_f_v,cem_f_delta_v,cem_bk_a_v,cem_bk_a_delta_v,cem_bk_b_v,cem_bk_b_delta_v +0,1,0.00055,6300.0,700.0,-6300.0,700.0,-8000.0,50.0,-3000.0,25.0,-2125.0,50.0,-4500.0,50.0,-2350.0,100.0,-2350.0,100.0 +0,2,0.00085,,,,,,,,,,,,,,,, +0,3,0.00126,,,,,,,,,,,,,,,, +0,4,0.00170,,,,,,,,,,,,,,,, +0,5,0.00340,,,,,,,,,,,,,,,, +0,6,0.00523,,,,,,,,,,,,,,,, +0,7,0.00659,,,,,,,,,,,,,,,, +0,8,0.01301,,,,,,,,,,,,,,,, +0,9,0.01830,,,,,,,,,,,,,,,, +1,1,0.00058,6300.0,700.0,-6300.0,700.0,-8000.0,50.0,-3000.0,25.0,-2125.0,50.0,-4200.0,50.0,-2500.0,100.0,-2500.0,100.0 +1,2,0.00089,,,,,,,,,,,,,,,, +1,3,0.00132,,,,,,,,,,,,,,,, +1,4,0.00179,,,,,,,,,,,,,,,, +1,5,0.00357,,,,,,,,,,,,,,,, +1,6,0.00549,,,,,,,,,,,,,,,, +1,7,0.00692,,,,,,,,,,,,,,,, +1,8,0.01366,,,,,,,,,,,,,,,, +1,9,0.01922,,,,,,,,,,,,,,,, From 62f6aa01b510488dc59454d1e7401ece81cafd89 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 15:31:41 -0600 Subject: [PATCH 6/8] Fix a couple of issues from PR --- imap_processing/cdf/config/imap_hi_variable_attrs.yaml | 6 +----- imap_processing/hi/utils.py | 5 +++-- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml index c4e173967..3bb52cdd6 100644 --- a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml @@ -494,11 +494,7 @@ hi_pset_geometric_factor: VALIDMIN: 0.0 VAR_NOTES: > Geometric factor per esa_energy_step, looked up from the gain-configuration - ancillary file using the pointing's classified gain configuration (see the - L1B DE product's gain_configuration_id global attribute). FILLVAL where the - gain configuration could not be classified. Not yet consumed by L2 - processing (deferred to a follow-on ticket that handles combining PSETs - from different gain configurations into a single map). + ancillary file using the l1b_de's classified gain configuration. VAR_TYPE: support_data hi_pset_calibration_prod: diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 92146a60e..4e495be81 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -206,8 +206,9 @@ def load_gain_configuration(path: str | Path | IO[str]) -> pd.DataFrame: Parameters ---------- - path : str or pathlib.Path or file-like object - Location of the gain configuration ancillary CSV file. + path : str or pathlib.Path + Location of the gain configuration ancillary CSV file. Anything + accepted by ``pandas.read_csv`` (e.g. a file-like object) also works. Returns ------- From b79071692aa0f23e262f398a69310f500913b73f Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Wed, 29 Jul 2026 15:35:41 -0600 Subject: [PATCH 7/8] Copilot feedback --- imap_processing/hi/utils.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/imap_processing/hi/utils.py b/imap_processing/hi/utils.py index 4e495be81..42a22ff07 100644 --- a/imap_processing/hi/utils.py +++ b/imap_processing/hi/utils.py @@ -637,11 +637,17 @@ def query(self, query_met: float | Iterable[float]) -> int | np.ndarray: query_mets = np.atleast_1d(query_met) results = np.full(query_mets.shape, self.NO_MATCH, dtype=np.int64) - for i, qm in enumerate(query_mets): - mask = (self.df["start_met"] <= qm) & (self.df["end_met"] >= qm) - matches = self.df[mask] - if not matches.empty: - results[i] = matches["config_id"].iloc[0] + if not self.df.empty: + starts = self.df["start_met"].to_numpy() + ends = self.df["end_met"].to_numpy() + config_ids = self.df["config_id"].to_numpy() + for start, end, cfg in zip(starts, ends, config_ids, strict=False): + mask = ( + (results == self.NO_MATCH) + & (query_mets >= start) + & (query_mets <= end) + ) + results[mask] = cfg return int(results[0]) if is_scalar_met else results From 4db1e73155e1e56d72e4ca3611cb3304ad14c7ba Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Thu, 30 Jul 2026 10:36:58 -0600 Subject: [PATCH 8/8] Fix documentation --- imap_processing/hi/hi_l1b.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/imap_processing/hi/hi_l1b.py b/imap_processing/hi/hi_l1b.py index 44f252175..e58837e1d 100644 --- a/imap_processing/hi/hi_l1b.py +++ b/imap_processing/hi/hi_l1b.py @@ -622,11 +622,12 @@ def classify_gain_configuration( """ Classify which gain configuration (config_id) a pointing is running. - Uses the median detector high voltages (see HiConstants.GAIN_CONFIG_HV_ - FIELDS) of the first contiguous HVSCI segment in the pointing, matched - against the gain-configuration ancillary table. This assumes that when a - gain test is run during a pointing, the first HVSCI segment of that - pointing is run at the pointing's real (non-gain-test) configuration. + Uses the median detector high voltages (see + HiConstants.GAIN_CONFIG_HV_FIELDS) of the first contiguous HVSCI segment + in the pointing, matched against the gain-configuration ancillary table. + This assumes that when a gain test is run during a pointing, the first + HVSCI segment of that pointing is run at the pointing's real + (non-gain-test) configuration. Parameters ----------