3395 hi l1c look up geometric factor by gain state via extended cal prod schema

imap_processing/cdf/config/imap_hi_variable_attrs.yaml

@@ -482,20 +482,25 @@ 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:

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,
@@ -116,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: "
@@ -141,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))
@@ -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,
@@ -560,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.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.
 
@@ -576,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 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)
@@ -603,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
@@ -670,11 +637,9 @@ def de_gain_test_filter(
         ImapHiL1bDeFlags.BAD_DETECTOR_VOLTAGE
     )
 
-    gain_match_values = 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
 
 

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)
 
@@ -126,11 +132,12 @@ 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.
+    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))
@@ -348,41 +355,65 @@ 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],
-) -> dict[str, xr.DataArray]:
+def add_pset_geometric_factor(
+    pset_ds: xr.Dataset,
+    gain_config_df: pd.DataFrame | None,
+) -> xr.Dataset:
     """
-    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.
+    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.
+    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
     -------
-    dict[str, xarray.DataArray]
-        Dictionary containing the "geometric_factor" DataArray (all
-        FILLVAL), dims (epoch, esa_energy_step).
+    xarray.Dataset
+        The input pset_ds, updated in place with a "geometric_factor"
+        variable, 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. 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).
     """
-    return create_dataset_variables(
+    geometric_factor_var = create_dataset_variables(
         ["geometric_factor"],
-        coords=pset_coords,
+        coords=pset_ds.coords,
         att_manager_lookup_str="hi_pset_{0}",
     )
+    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,
+        )
+        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(
     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]:
@@ -393,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
@@ -403,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(
@@ -412,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 = {
@@ -446,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)
@@ -458,8 +494,11 @@ def pset_counts(
         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.
+        # 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[0]]
         np.add.at(
             counts_var["counts"].data[0, i_esa, i_cal_prod],

imap_processing/hi/hi_l2.py

@@ -458,6 +458,15 @@ 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.add_pset_geometric_factor()'s docstring), so use the
+    # first (and, today, only) gain_config_id present in the ancillary file.
+    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"]