diff --git a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml index 10b405ad51..ddead25a9a 100644 --- a/imap_processing/cdf/config/imap_hi_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_hi_variable_attrs.yaml @@ -509,7 +509,7 @@ hi_pset_calibration_prod: CATDESC: Coincidence types are combined into various calibration products FIELDNAM: Calibration product number FILLVAL: 255 - FORMAT: I1 + FORMAT: I3 LABLAXIS: Calibration Prod MONOTON: INCREASE SCALETYP: linear diff --git a/imap_processing/hi/hi_l1c.py b/imap_processing/hi/hi_l1c.py index 14be8f2dfe..8d8dbb8dd0 100644 --- a/imap_processing/hi/hi_l1c.py +++ b/imap_processing/hi/hi_l1c.py @@ -116,14 +116,30 @@ def generate_pset_dataset( 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) + gain_config_id = config_df.cal_prod_config.match_gain_config_id(hv_deltas) + gain_config_df = None if gain_config_id is None else config_df.loc[gain_config_id] + if gain_config_id is None: + logger.warning( + "Pointing's detector gain state matches zero or multiple " + "gain_config_ids in the cal-prod config. PSET will have a single " + "fill calibration_prod, zero counts, fill geometric_factor and " + "backgrounds, and a NaN gain_config_id global attribute." + ) + cal_prod_numbers = None + else: + logger.info(f"Pointing gain state matches gain_config_id={gain_config_id}") + # The calibration_prod coordinate holds the matched gain_config_id's + # calibration products only. + cal_prod_numbers = np.sort( + gain_config_df.index.unique(level="calibration_prod").values + ) # read background configuration file background_df = BackgroundConfig.from_csv(background_config_path) pset_dataset = empty_pset_dataset( de_dataset.ccsds_met.data.mean(), de_dataset.esa_energy_step, - config_df.cal_prod_config.calibration_product_numbers, + cal_prod_numbers, logical_source_parts["sensor"], ) # Calculate and add despun_z, hae_latitude, and hae_longitude variables to @@ -132,6 +148,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"])) + # Record which gain configuration this pointing ran so that L2 can + # combine PSETs from different gain configurations. NaN when unmatched, + # matching the L1B DE convention for the gain-state global attributes. + pset_dataset.attrs["gain_config_id"] = ( + np.nan if gain_config_id is None else gain_config_id + ) # 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) @@ -159,7 +181,7 @@ def generate_pset_dataset( def empty_pset_dataset( l1b_met: float, l1b_energy_steps: xr.DataArray, - cal_prod_numbers: npt.NDArray[np.int_], + cal_prod_numbers: npt.NDArray[np.int_] | None, sensor_str: str, ) -> xr.Dataset: """ @@ -172,9 +194,11 @@ def empty_pset_dataset( repoint-table data to get the start and end times of the pointing. l1b_energy_steps : xarray.DataArray The array of esa_energy_step data from the L1B DE product. - cal_prod_numbers : numpy.ndarray + cal_prod_numbers : numpy.ndarray or None Array of calibration product numbers from the configuration file. - These can be arbitrary integers, not necessarily starting at 0. + These can be arbitrary integers, not necessarily starting at 0. If + None (pointing matched no gain configuration), the calibration_prod + coordinate is a single FILLVAL entry. sensor_str : str '45sensor' or '90sensor'. @@ -227,6 +251,8 @@ def empty_pset_dataset( "hi_pset_calibration_prod", check_schema=False ).copy() dtype = attrs.pop("dtype") + if cal_prod_numbers is None: + cal_prod_numbers = np.array([attrs["FILLVAL"]]) coords["calibration_prod"] = xr.DataArray( cal_prod_numbers.astype(dtype), name="calibration_prod", @@ -651,13 +677,19 @@ def pset_backgrounds( excess counts from the outer ESA during background testing. The result is clipped to zero to prevent negative rates. + If the calibration_prod coordinate is the single FILLVAL entry used for a + pointing with no matched gain configuration, the background variables are + returned filled with FILLVAL. + Parameters ---------- pset_coords : dict[str, xarray.DataArray] The PSET coordinates from the xarray.Dataset. background_config_df : pandas.DataFrame Background configuration DataFrame with MultiIndex - (calibration_prod, background_index, esa_energy_step). + (calibration_prod, background_index, esa_energy_step). Must contain + every calibration_prod in pset_coords and may contain others, which + are ignored. l1b_de_dataset : xarray.Dataset The L1B dataset for the pointing being processed. goodtimes_ds : xarray.Dataset @@ -685,6 +717,13 @@ def pset_backgrounds( for var_name in ["background_rates", "background_rates_uncertainty"] } + # A pointing with no matched gain configuration has a single FILLVAL + # calibration_prod, so there are no backgrounds to compute and they are + # left at FILLVAL. + cal_prod = pset_coords["calibration_prod"] + if np.all(cal_prod.values == cal_prod.attrs["FILLVAL"]): + return output_vars + # Get total exposure time total_exposure_time = float(exposure_times.sum()) @@ -705,15 +744,19 @@ def pset_backgrounds( # Config now has dims: (calibration_prod, background_index, esa_energy_step) config_ds = background_config_df.to_xarray() - # Validate calibration products match (compare values, not DataArray metadata) - if not np.array_equal( - config_ds["calibration_prod"].values, pset_coords["calibration_prod"].values - ): + # Backgrounds are independent of gain configuration, so the background + # config may define more calibration products than this PSET (which holds + # only its matched gain configuration's products). Require every PSET + # calibration product to be present, then select just those. + pset_cal_prods = pset_coords["calibration_prod"].values + missing_cal_prods = np.setdiff1d(pset_cal_prods, config_ds["calibration_prod"]) + if missing_cal_prods.size > 0: raise ValueError( - f"Calibration products in pset_coords and background_config_df " - f"do not match. pset_coords: {pset_coords['calibration_prod'].values}, " + f"Calibration products in pset_coords are missing from " + f"background_config_df. pset_coords: {pset_cal_prods}, " f"background_config_df: {config_ds['calibration_prod'].values}" ) + config_ds = config_ds.sel(calibration_prod=pset_cal_prods) # Validate ESA energy steps match (compare values, not DataArray metadata) if not np.array_equal( diff --git a/imap_processing/tests/hi/test_hi_l1c.py b/imap_processing/tests/hi/test_hi_l1c.py index 7a8057e13e..21a2447f5b 100644 --- a/imap_processing/tests/hi/test_hi_l1c.py +++ b/imap_processing/tests/hi/test_hi_l1c.py @@ -71,7 +71,9 @@ def test_hi_l1c( @pytest.mark.external_kernel @pytest.mark.external_test_data +@pytest.mark.parametrize("gain_matched", [True, False]) def test_generate_pset_dataset( + gain_matched, hi_l1b_de_dataset, hi_goodtimes_dataset, hi_test_cal_prod_config_path, @@ -87,6 +89,9 @@ def test_generate_pset_dataset( # placeholders so add_pset_geometric_factor() and pset_counts() have # something to look up. l1b_dataset.attrs.update(NOMINAL_HV_DELTAS) + if not gain_matched: + # HV deltas matching no gain_config_id in the test cal-prod config + l1b_dataset.attrs.update({field: 0.0 for field in NOMINAL_HV_DELTAS}) l1b_met = l1b_dataset["ccsds_met"].values[0] # Set repoint start and end times. seconds_per_day = 24 * 60 * 60 @@ -110,15 +115,30 @@ def test_generate_pset_dataset( np.testing.assert_array_equal(l1c_dataset.hae_latitude.data.shape, (1, 3600)) np.testing.assert_array_equal(l1c_dataset.hae_longitude.data.shape, (1, 3600)) np.testing.assert_array_equal(l1c_dataset.exposure_times.data.shape, (1, 9, 3600)) + n_cal_prods = 2 if gain_matched else 1 for var in [ "counts", "background_rates", "background_rates_uncertainty", ]: - np.testing.assert_array_equal(l1c_dataset[var].data.shape, (1, 9, 2, 3600)) - - # Test ISTP compliance by writing CDF - write_cdf(l1c_dataset) + np.testing.assert_array_equal( + l1c_dataset[var].data.shape, (1, 9, n_cal_prods, 3600) + ) + if gain_matched: + assert l1c_dataset.attrs["gain_config_id"] == 0 + else: + assert np.isnan(l1c_dataset.attrs["gain_config_id"]) + np.testing.assert_array_equal(l1c_dataset["calibration_prod"].values, [255]) + assert np.all(l1c_dataset["counts"].values == 0) + + # Test ISTP compliance by writing CDF, and that the gain_config_id global + # attribute survives the round trip + cdf_path = write_cdf(l1c_dataset) + round_trip_id = load_cdf(cdf_path).attrs["gain_config_id"] + if gain_matched: + assert round_trip_id == 0 + else: + assert np.isnan(round_trip_id) @mock.patch("imap_processing.hi.hi_l1c.pset_backgrounds") @@ -221,6 +241,127 @@ def _make_pset_ds_for_geometric_factor(esa_energy_steps, calibration_prods): ) +# Two gain configurations whose calibration products differ, to verify that +# a PSET's calibration_prod coordinate holds only its matched gain +# configuration's products. +_TWO_GAIN_CONFIG_CSV = """\ +gain_config_id,calibration_prod,esa_energy_step,geometric_factor,coincidence_type_list,tof_ab_low,tof_ab_high,tof_ac1_low,tof_ac1_high,tof_bc1_low,tof_bc1_high,tof_c1c2_low,tof_c1c2_high,mcp_delta_v,mcp_delta_v_tol,cem_a_delta_v,cem_a_delta_v_tol,cem_b_delta_v,cem_b_delta_v_tol,tof_v,tof_v_tol +0,0,1,0.00055,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,875.0,75.0,2150.0,150.0,2150.0,150.0,-8000.0,50.0 +0,0,2,0.00085,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,1,1,0.00055,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +0,1,2,0.00085,BC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +1,0,1,0.00066,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,500.0,50.0,1000.0,100.0,1000.0,100.0,-4000.0,50.0 +1,0,2,0.00099,ABC1C2,0,1023,-1023,1023,-1023,1023,0,1023,,,,,,,, +""" + + +@pytest.mark.parametrize( + "hv_deltas, expected_gain_config_id, expected_cal_prods", + [ + (NOMINAL_HV_DELTAS, 0, [0, 1]), + ( + { + "mcp_delta_v": 500.0, + "cem_a_delta_v": 1000.0, + "cem_b_delta_v": 1000.0, + "tof_v": -4000.0, + }, + 1, + [0], + ), + ( + { + "mcp_delta_v": 0.0, + "cem_a_delta_v": 0.0, + "cem_b_delta_v": 0.0, + "tof_v": 0.0, + }, + None, + [255], + ), + ], +) +@mock.patch("imap_processing.hi.hi_l1c.pset_backgrounds", return_value={}) +@mock.patch("imap_processing.hi.hi_l1c.pset_exposure") +@mock.patch("imap_processing.hi.hi_l1c.pset_counts", return_value={}) +@mock.patch("imap_processing.hi.hi_l1c.pset_geometry", return_value={}) +@mock.patch( + "imap_processing.hi.hi_l1c.get_pointing_times", + return_value=(482372065.0, 482374065.0), +) +def test_generate_pset_dataset_gain_config( + mock_get_pointing_times, + mock_pset_geometry, + mock_pset_counts, + mock_pset_exposure, + mock_pset_backgrounds, + hv_deltas, + expected_gain_config_id, + expected_cal_prods, + hi_test_background_config_path, + tmp_path, +): + """Test generate_pset_dataset records the matched gain_config_id and + builds the calibration_prod coordinate from that gain configuration's + products only (a single FILLVAL entry if the pointing matched no gain + config).""" + cal_prod_path = tmp_path / "imap_hi_90sensor-cal-prod_20240101_v001.csv" + cal_prod_path.write_text(_TWO_GAIN_CONFIG_CSV) + n_energy_steps = 2 + mock_l1b_dataset = xr.Dataset( + coords={"epoch": xr.DataArray(np.arange(6), dims=["epoch"])}, + data_vars={ + "ccsds_met": xr.DataArray(np.full(6, 482373065.0), dims=["epoch"]), + "esa_energy_step": xr.DataArray( + np.arange(n_energy_steps + 1).repeat(2), attrs={"FILLVAL": 255} + ), + }, + attrs={ + "Logical_file_id": "imap_hi_l1b_90sensor-de_20250415_v999", + "Logical_source": "imap_hi_l1b_90sensor-de", + **hv_deltas, + }, + ) + mock_pset_exposure.return_value = { + "exposure_times": xr.DataArray( + np.ones((1, n_energy_steps, 3600), dtype=np.float32), + dims=["epoch", "esa_energy_step", "spin_angle_bin"], + ) + } + + pset = hi_l1c.generate_pset_dataset( + mock_l1b_dataset, cal_prod_path, xr.Dataset(), hi_test_background_config_path + ) + + np.testing.assert_array_equal(pset["calibration_prod"].values, expected_cal_prods) + np.testing.assert_array_equal( + pset["calibration_prod_label"].values, np.array(expected_cal_prods).astype(str) + ) + if expected_gain_config_id is None: + assert np.isnan(pset.attrs["gain_config_id"]) + else: + assert pset.attrs["gain_config_id"] == expected_gain_config_id + mock_pset_backgrounds.assert_called_once() + # The gain_config_df passed to pset_counts is the matched config's rows + gain_config_df = mock_pset_counts.call_args[0][1] + if expected_gain_config_id is None: + assert gain_config_df is None + assert ( + pset["calibration_prod"].values[0] + == (pset["calibration_prod"].attrs["FILLVAL"]) + ) + fillval = np.float32(pset["geometric_factor"].attrs["FILLVAL"]) + assert np.all(pset["geometric_factor"].values == fillval) + else: + assert gain_config_df.index.names == ["calibration_prod", "esa_energy_step"] + expected_gf = 0.00066 if expected_gain_config_id == 1 else 0.00055 + np.testing.assert_allclose( + pset["geometric_factor"].sel(esa_energy_step=1, calibration_prod=0), + [expected_gf], + rtol=1e-6, + ) + + def test_add_pset_geometric_factor_matching_gain_state(hi_test_cal_prod_config_path): """Test add_pset_geometric_factor for a pointing whose gain-match values match the test cal-prod config's gain_config_id=0 reference values @@ -993,6 +1134,98 @@ def test_compute_background_counts_missing_cal_prod_raises_error( ) +@mock.patch("imap_processing.hi.hi_l1c._compute_background_counts") +def test_pset_backgrounds_unmatched_gain_config_returns_fillval( + mock_compute_background_counts, +): + """Test pset_backgrounds returns FILLVAL for a fill calibration_prod.""" + n_energy = 2 + n_spin_bins = 3600 + pset_coords = { + "epoch": xr.DataArray(np.array([0], dtype=np.int64), dims=["epoch"]), + "esa_energy_step": xr.DataArray( + np.arange(n_energy) + 1, dims=["esa_energy_step"] + ), + "calibration_prod": xr.DataArray( + np.array([255], dtype=np.uint8), + dims=["calibration_prod"], + attrs={"FILLVAL": 255}, + ), + "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), + } + # Nonzero exposure so the early return is not the zero-exposure one + exposure_times = xr.DataArray( + np.ones((1, n_energy, n_spin_bins), dtype=np.float32), + dims=["epoch", "esa_energy_step", "spin_angle_bin"], + ) + + result = hi_l1c.pset_backgrounds( + pset_coords, pd.DataFrame(), xr.Dataset(), xr.Dataset(), exposure_times + ) + + mock_compute_background_counts.assert_not_called() + for var in ["background_rates", "background_rates_uncertainty"]: + assert result[var].shape == (1, n_energy, 1, n_spin_bins) + fillval = np.float32(result[var].attrs["FILLVAL"]) + assert np.all(result[var].values == fillval) + + +@mock.patch("imap_processing.hi.hi_l1c._compute_background_counts") +def test_pset_backgrounds_cal_prod_subset_of_background_config( + mock_compute_background_counts, hi_test_background_config_path +): + """Test pset_backgrounds when the PSET has fewer cal prods than the config. + + Backgrounds are gain-independent, so the background config defines every + calibration product while a PSET holds only its matched gain config's. + """ + background_df = utils.BackgroundConfig.from_csv(hi_test_background_config_path) + assert background_df.index.unique(level="calibration_prod").tolist() == [0, 1] + esa_steps = background_df.index.unique(level="esa_energy_step").values + background_indices = background_df.index.unique(level="background_index").values + n_spin_bins = 3600 + pset_coords = { + "epoch": xr.DataArray(np.array([0], dtype=np.int64), dims=["epoch"]), + "esa_energy_step": xr.DataArray(esa_steps, dims=["esa_energy_step"]), + "calibration_prod": xr.DataArray( + np.array([0], dtype=np.uint8), + dims=["calibration_prod"], + attrs={"FILLVAL": 255}, + ), + "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), + } + exposure_times = xr.DataArray( + np.ones((1, len(esa_steps), n_spin_bins), dtype=np.float32), + dims=["epoch", "esa_energy_step", "spin_angle_bin"], + ) + mock_compute_background_counts.return_value = xr.DataArray( + np.full((1, 1, len(background_indices)), 1000.0), + dims=["epoch", "calibration_prod", "background_index"], + coords={ + "epoch": pset_coords["epoch"], + "calibration_prod": [0], + "background_index": background_indices, + }, + ) + + result = hi_l1c.pset_backgrounds( + pset_coords, background_df, xr.Dataset(), xr.Dataset(), exposure_times + ) + + # Same result as passing a config pre-filtered to the PSET's cal prods + expected = hi_l1c.pset_backgrounds( + pset_coords, + background_df.loc[[0]], + xr.Dataset(), + xr.Dataset(), + exposure_times, + ) + for var in ["background_rates", "background_rates_uncertainty"]: + assert result[var].shape == (1, len(esa_steps), 1, n_spin_bins) + assert np.all(np.isfinite(result[var].values)) + np.testing.assert_array_equal(result[var].values, expected[var].values) + + @mock.patch("imap_processing.hi.hi_l1c._compute_background_counts") def test_pset_backgrounds_cal_prod_mismatch_raises_error( mock_compute_background_counts, @@ -1014,6 +1247,7 @@ def test_pset_backgrounds_cal_prod_mismatch_raises_error( "calibration_prod": xr.DataArray( np.array([0, 1], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), } @@ -1072,8 +1306,8 @@ def test_pset_backgrounds_cal_prod_mismatch_raises_error( # Verify that pset_backgrounds raises ValueError with expected message with pytest.raises( ValueError, - match="Calibration products in pset_coords and " - "background_config_df do not match", + match="Calibration products in pset_coords are missing from " + "background_config_df", ): hi_l1c.pset_backgrounds( pset_coords, @@ -1103,6 +1337,7 @@ def test_pset_backgrounds_esa_energy_step_mismatch_raises_error( "calibration_prod": xr.DataArray( np.array([0], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), } @@ -1185,6 +1420,7 @@ def test_pset_backgrounds_applies_offset_correction(mock_compute_background_coun "calibration_prod": xr.DataArray( np.array([0], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), } @@ -1274,6 +1510,7 @@ def test_pset_backgrounds_offset_does_not_go_negative(mock_compute_background_co "calibration_prod": xr.DataArray( np.array([0], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), } @@ -1368,6 +1605,7 @@ def test_pset_backgrounds_uncertainty_includes_constant_terms( "calibration_prod": xr.DataArray( np.array([0], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), } @@ -1464,6 +1702,7 @@ def test_pset_backgrounds_esa_7_8_9_extra_uncertainty(mock_compute_background_co "calibration_prod": xr.DataArray( np.array([0], dtype=np.int64), dims=["calibration_prod"], + attrs={"FILLVAL": 255}, ), "spin_angle_bin": xr.DataArray(np.arange(n_spin_bins), dims=["spin_angle_bin"]), }