From ecaf02ca27ca0218194b5e7d7290961bef5ff051 Mon Sep 17 00:00:00 2001 From: Sean Hoyt Date: Mon, 5 Jan 2026 15:07:01 -0700 Subject: [PATCH 1/5] resweep_histogram_data working for all hist fields --- imap_processing/lo/l1b/lo_l1b.py | 248 ++++++++-------- ...eep-table-small_20250101_20260301_v001.csv | 1 + imap_processing/tests/lo/test_lo_l1b.py | 264 ++++++++++++++++-- 3 files changed, 376 insertions(+), 137 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index ce995eb9e7..25dd2bdb0a 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -134,21 +134,21 @@ def lo_l1b(sci_dependencies: dict, anc_dependencies: list) -> list[Path]: spin_data = sci_dependencies["imap_lo_l1a_spin"] # initialize the L1B Histogram Rates dataset from the L1A Histogram Rates # This carries over the epoch and count fields from L1A - l1b_histrates = initialize_l1b_histrates(l1a_hist, attr_mgr_l1b, logical_source) + l1b_all_rates = initialize_all_rates(l1a_hist, attr_mgr_l1b) # set spin cycle and remove invalid spin ASCs l1b_histrates = set_spin_cycle_from_spin_data( - l1a_hist, l1b_histrates, spin_data + l1a_hist, l1b_all_rates, spin_data ) pointing_start_met, pointing_end_met = get_pointing_times( ttj2000ns_to_met(l1a_hist["epoch"].values[0].item()) ) - l1b_histrates = set_esa_mode( - pointing_start_met, pointing_end_met, anc_dependencies, l1b_histrates + l1b_all_rates = set_esa_mode( + pointing_start_met, pointing_end_met, anc_dependencies, l1b_all_rates ) # resweep the histogram data - l1b_histrates, exposure_factor = resweep_histogram_data( - l1b_histrates, anc_dependencies + l1b_all_rates, exposure_factor = resweep_histogram_data( + l1b_all_rates, anc_dependencies ) # Get the start and end times for each spin epoch acq_start, acq_end = convert_start_end_acq_times(spin_data) @@ -156,13 +156,13 @@ def lo_l1b(sci_dependencies: dict, anc_dependencies: list) -> list[Path]: avg_spin_durations_per_cycle = get_avg_spin_durations_per_cycle( acq_start, acq_end ) - l1b_histrates = calculate_histogram_rates( - l1b_histrates, - acq_start, - acq_end, - avg_spin_durations_per_cycle, - exposure_factor, - ) + # l1b_histrates = calculate_histogram_rates( + # l1b_histrates, + # acq_start, + # acq_end, + # avg_spin_durations_per_cycle, + # exposure_factor, + # ) datasets_to_return.append(l1b_histrates) return datasets_to_return @@ -1370,8 +1370,8 @@ def create_badtimes_dataset() -> xr.Dataset: return thruster_ds -def initialize_l1b_histrates( - l1a_hist: xr.Dataset, attr_mgr_l1b: ImapCdfAttributes, logical_source: str +def initialize_all_rates( + l1a_hist: xr.Dataset, attr_mgr_l1b: ImapCdfAttributes ) -> xr.Dataset: """ Initialize the L1B histogram rates dataset. @@ -1382,15 +1382,13 @@ def initialize_l1b_histrates( The L1A histogram rates dataset. attr_mgr_l1b : ImapCdfAttributes Attribute manager used to get the L1B histogram rates dataset attributes. - logical_source : str - The logical source of the data product that's being created. Returns ------- - l1b_histrates : xr.Dataset - The initialized L1B histogram rates dataset. + l1b_all_rates : xr.Dataset + The initialized L1B histogram and monitor rates dataset. """ - l1b_histrates = xr.Dataset( + l1b_all_rates = xr.Dataset( coords={ "epoch": xr.DataArray(l1a_hist["epoch"].values, dims=["epoch"]), "esa_step": l1a_hist["esa_step"], @@ -1398,40 +1396,60 @@ def initialize_l1b_histrates( l1a_hist["azimuth_6"].values, dims=["spin_bin_6"], ), + "spin_bin_60": xr.DataArray( + l1a_hist["azimuth_60"].values, + dims=["spin_bin_60"], + ), }, - attrs=attr_mgr_l1b.get_global_attributes(logical_source), ) + # Define field mappings from L1A to the L1B dataset + # according to the algorithm document, only the histogram rate product should have + # counts (hydrogen and oxygen). The rest of the fields will be dropped after + # processing when the monitor rates product is created. These fields are needed for + # processing though and are added to the all_rates dataset temporarily. + # {"l1a_field_name": "l1b_all_rates_field_name"} + spin_bin_6_fields = { + "hydrogen": "h_counts", + "oxygen": "o_counts", + "tof0_tof1": "tof0_tof1_counts", + "tof0_tof2": "tof0_tof2_counts", + "tof1_tof3": "tof1_tof3_counts", + "silver_triple": "silver_triple_counts", + } + spin_bin_60_fields = { + "start_a": "start_a_counts", + "start_c": "start_c_counts", + "stop_b0": "stop_b0_counts", + "stop_b3": "stop_b3_counts", + "tof0": "tof0_counts", + "tof1": "tof1_counts", + "tof2": "tof2_counts", + "tof3": "tof3_counts", + } + + for l1a_field, l1b_field in spin_bin_6_fields.items(): + l1b_all_rates[l1b_field] = xr.DataArray( + l1a_hist[l1a_field].values, + dims=["epoch", "esa_step", "spin_bin_6"], + ) - # l1b_histrates["epoch"] = xr.DataArray( - # l1a_hist["epoch"].values, - # dims=["epoch"], - # attrs=attr_mgr_l1b.get_variable_attributes("epoch"), - # ) - # Copy over fields from L1A DE that will not change in L1B processing - l1b_histrates["h_counts"] = xr.DataArray( - l1a_hist["hydrogen"].values, - dims=["epoch", "esa_step", "spin_bin_6"], - # TODO: Add hydrogen to YAML file - # attrs=attr_mgr.get_variable_attributes("hydrogen"), - ) - l1b_histrates["o_counts"] = xr.DataArray( - l1a_hist["oxygen"].values, - dims=["epoch", "esa_step", "spin_bin_6"], - # TODO: Add oxygen to YAML file - # attrs=attr_mgr.get_variable_attributes("oxygen"), - ) + for l1a_field, l1b_field in spin_bin_60_fields.items(): + l1b_all_rates[l1b_field] = xr.DataArray( + l1a_hist[l1a_field].values, + dims=["epoch", "esa_step", "spin_bin_60"], + ) - return l1b_histrates + return l1b_all_rates def resweep_histogram_data( l1b_histrates: xr.Dataset, anc_dependencies: list, -) -> tuple[xr.Dataset, np.ndarray]: +) -> tuple[xr.Dataset, dict[str, np.ndarray]]: """ Correct energy steps in histogram data based on sweep and LUT tables. - Returns the updated dataset and a 3D array of reswept counts + Returns the updated dataset and a dictionary of 3D arrays of reswept counts (epoch, azimuth, esa_step) indicating how many original steps were reswept into each final step. @@ -1446,12 +1464,10 @@ def resweep_histogram_data( ------- l1b_histrates : xr.Dataset The updated L1B histogram rates dataset with reswept counts. - exposure_factor : np.ndarray - 3D array of exposure factors (epoch, azimuth, esa_step) indicating how many - ESA steps were reswept during resweeping. + exposure_factors : dict[str, np.ndarray] + Dictionary mapping field names to their 3D exposure factor arrays + (epoch, azimuth, esa_step). """ - # The sweep table contains the mapping of dates to the LUT table which shows how - # the ESA steps should be reswept. sweep_df = lo_ancillary.read_ancillary_file( next(str(s) for s in anc_dependencies if "sweep-table" in str(s)) ) @@ -1459,110 +1475,114 @@ def resweep_histogram_data( next(str(s) for s in anc_dependencies if "esa-mode-lut" in str(s)) ) - # Get the time information to compare the epochs to the sweep table dates + # Define field groups with their spin bin dimensions + spin_bin_6_fields = [ + "h_counts", + "o_counts", + "tof0_tof1_counts", + "tof0_tof2_counts", + "tof1_tof3_counts", + "silver_triple_counts", + ] + spin_bin_60_fields = [ + "start_a_counts", + "start_c_counts", + "stop_b0_counts", + "stop_b3_counts", + "tof0_counts", + "tof1_counts", + "tof2_counts", + "tof3_counts", + ] + sweep_dates = sweep_df["Date"].astype(str) epochs = l1b_histrates["epoch"].values epoch_utc = et_to_utc(ttj2000ns_to_et(epochs)) - # initialize the reswept counts arrays - h_counts_reswept = np.zeros_like(l1b_histrates["h_counts"].values) - o_counts_reswept = np.zeros_like(l1b_histrates["o_counts"].values) + # Initialize reswept arrays for all fields + reswept_data = {} + for field in spin_bin_6_fields + spin_bin_60_fields: + reswept_data[field] = np.zeros_like(l1b_histrates[field].values) - # Get the number of azimuth bins from the l1b_histrates dataset - num_azimuth = l1b_histrates.sizes["spin_bin_6"] - # initialize exposure factor to 1 as this will be used to scale (multiply) - # the exposure time later - exposure_factor = np.full( - (len(epochs), l1b_histrates.sizes["esa_step"], num_azimuth), 1, dtype=int - ) + # Initialize exposure factors for each field type + num_azimuth_6 = l1b_histrates.sizes["spin_bin_6"] + num_azimuth_60 = l1b_histrates.sizes["spin_bin_60"] + exposure_factors = {} + + for field in spin_bin_6_fields: + exposure_factors[field] = np.full( + (len(epochs), l1b_histrates.sizes["esa_step"], num_azimuth_6), 1, dtype=int + ) + + for field in spin_bin_60_fields: + exposure_factors[field] = np.full( + (len(epochs), l1b_histrates.sizes["esa_step"], num_azimuth_60), 1, dtype=int + ) for epoch_idx, epoch in enumerate(epoch_utc): - # Get only the date portion of the epoch string for comparison with the - # sweep table epoch_date_only = epoch.split("T")[0] - # if the epoch dat is not in the sweep table, raise an error if epoch_date_only not in sweep_dates.values: - raise ValueError( - f"No sweep table entry found for date {epoch} at epoch idx {epoch_idx}" + logger.warning( + f"Epoch {epoch} at index {epoch_idx} not found in sweep table" ) + continue - # Get the matching sweep table entry for the epoch date and its LUT table index matching_sweep = sweep_df[sweep_dates == epoch_date_only] unique_lut_tables = matching_sweep["LUT_table"].unique() - # There should only be one unique LUT table for each date if len(unique_lut_tables) != 1: - raise ValueError( - f"Expected exactly 1 unique LUT_table value for date {epoch_date_only}," - f" but found {len(unique_lut_tables)}: {unique_lut_tables}" + logger.warning( + f"Multiple LUT tables found for epoch {epoch} at index {epoch_idx}, " + f"but found tables {unique_lut_tables}." ) + continue - # Get the LUT entries for the identified LUT table index lut_table_idx = unique_lut_tables[0] lut_entries = lut_df[lut_df["Tbl_Idx"] == lut_table_idx].copy() - # If there are no LUT entries for the identified LUT table, log a warning - # and skip resweeping for this epoch if len(lut_entries) == 0: - logger.warning(f"No LUT entries found for table index {lut_table_idx}") - h_counts_reswept[epoch_idx] = l1b_histrates["h_counts"].values[epoch_idx] - o_counts_reswept[epoch_idx] = l1b_histrates["o_counts"].values[epoch_idx] + logger.warning( + f"No LUT entries for epoch {epoch} at index {epoch_idx}. Looking" + f"for table index {lut_table_idx}." + ) continue - # Sort the LUT entries by E-Step_Idx to ensure correct mapping order lut_entries = lut_entries.sort_values("E-Step_Idx") - # Create a mapping of original ESA step index to true ESA step energy_step_mapping = {} - # Loop through the LUT entries and populate the mapping for _, row in lut_entries.iterrows(): - # Original ESA step index is 1-based, convert to 0-based esa_idx = int(row["E-Step_Idx"]) - 1 - # True ESA step is 1-based true_esa_step = int(row["E-Step_lvl"]) - # Populate the mapping energy_step_mapping[esa_idx] = true_esa_step - # TODO: Change all instances of azimuth to spin bin - # Resweep the counts for each spin bin using the energy step mapping - for az_idx in range(num_azimuth): - h_original = l1b_histrates["h_counts"].values[epoch_idx, :, az_idx] - o_original = l1b_histrates["o_counts"].values[epoch_idx, :, az_idx] - - # Loop through the original ESA step indices and map to the true ESA steps - for orig_idx, true_esa_step in energy_step_mapping.items(): - # Check that the original index and true ESA step are within bounds - if orig_idx < len(h_original) and 1 <= true_esa_step <= 7: - # Resweep the counts into the true ESA step - # (convert to 0-based index) - reswept_idx = true_esa_step - 1 - h_counts_reswept[epoch_idx, reswept_idx, az_idx] += h_original[ + # Process spin_bin_6 fields + for field in spin_bin_6_fields: + for az_idx in range(num_azimuth_6): + original = l1b_histrates[field].values[epoch_idx, :, az_idx] + for orig_idx, true_esa_step in energy_step_mapping.items(): + target_idx = true_esa_step - 1 + reswept_data[field][epoch_idx, target_idx, az_idx] += original[ orig_idx ] - o_counts_reswept[epoch_idx, reswept_idx, az_idx] += o_original[ + exposure_factors[field][epoch_idx, target_idx, az_idx] += 1 + + # Process spin_bin_60 fields + for field in spin_bin_60_fields: + for az_idx in range(num_azimuth_60): + original = l1b_histrates[field].values[epoch_idx, :, az_idx] + for orig_idx, true_esa_step in energy_step_mapping.items(): + target_idx = true_esa_step - 1 + reswept_data[field][epoch_idx, target_idx, az_idx] += original[ orig_idx ] - # If a reswept was needed for this index, increment the exposure - # factor to so the exposure time can be scaled accordingly - if orig_idx != reswept_idx: - exposure_factor[epoch_idx, reswept_idx, az_idx] += 1 - - else: - logger.warning( - f"Original ESA index {orig_idx} or " - f"true ESA step {true_esa_step}" - f" out of bounds at epoch idx {epoch_idx}, " - f"spin bin idx {az_idx}" - ) - - l1b_histrates["h_counts"].values = h_counts_reswept - l1b_histrates["o_counts"].values = o_counts_reswept - l1b_histrates.attrs["energy_step_correction"] = ( - "Applied LUT table energy step mapping" - ) + exposure_factors[field][epoch_idx, target_idx, az_idx] += 1 + + # Update dataset with reswept data + for field in spin_bin_6_fields + spin_bin_60_fields: + l1b_histrates[field].values = reswept_data[field] - return l1b_histrates, exposure_factor + return l1b_histrates, exposure_factors def calculate_histogram_rates( diff --git a/imap_processing/tests/lo/test_anc/imap_lo_sweep-table-small_20250101_20260301_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_sweep-table-small_20250101_20260301_v001.csv index 4ceadfe573..6b4d507a81 100644 --- a/imap_processing/tests/lo/test_anc/imap_lo_sweep-table-small_20250101_20260301_v001.csv +++ b/imap_processing/tests/lo/test_anc/imap_lo_sweep-table-small_20250101_20260301_v001.csv @@ -1,4 +1,5 @@ YYYYDDD,ISN/ENA,Operation Mode,GoodTime_start,GoodTime_end,Instrument,LUT_table,ESA_Mode,E-Step1,E-Step2,E-Step3,E-Step4,E-Step5,E-Step6,E-Step7 +2024001,ENA,Nominal Mode,469024800,469084200,Lo,20,HiRes,1,2,3,4,5,6,7 2025001,ISN,HiRes Mode,473389200,473407618,Lo,1,HiRes,1,2,3,4,5,6,7 2025001,ISN,HiRes Mode,473407618,473420896,Lo,1,HiRes,1,2,3,4,5,6,7 2025001,ISN,HiRes Mode,473420896,473472000,Lo,1,HiRes,1,2,3,4,5,6,7 diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index b2baeb00ab..4dc48e136b 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1,3 +1,4 @@ +import logging from collections import namedtuple from unittest.mock import patch @@ -118,11 +119,24 @@ def l1a_hist(): { "hydrogen": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), "oxygen": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "tof0_tof1": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "tof0_tof2": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "tof1_tof3": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "silver_triple": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "start_a": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "start_c": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "stop_b0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "stop_b3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof1": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof2": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), }, coords={ "epoch": epoch_date, "esa_step": np.arange(1, 8), "azimuth_6": np.arange(60), + "azimuth_60": np.arange(6), }, attrs={"Logical_source": "imap_lo_l1a_histogram"}, ) @@ -802,17 +816,75 @@ def test_resweep_histogram_success(anc_dependencies): ) l1b_histrate = xr.Dataset( { - "h_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), + "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "tof0_tof1_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof0_tof2_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof1_tof3_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "silver_triple_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "start_a_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "start_c_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b0_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b3_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), }, coords={ "epoch": epoch_date, "esa_step": np.arange(1, 8), "spin_bin_6": np.arange(60), + "spin_bin_60": np.arange(6), }, ) - exposure_factor_expected = np.full((2, 7, 60), 1) - exposure_factor_expected[:, 0, :] = 2 + spin_bin_6_fields = [ + "h_counts", + "o_counts", + "tof0_tof1_counts", + "tof0_tof2_counts", + "tof1_tof3_counts", + "silver_triple_counts", + ] + spin_bin_60_fields = [ + "start_a_counts", + "start_c_counts", + "stop_b0_counts", + "stop_b3_counts", + "tof0_counts", + "tof1_counts", + "tof2_counts", + "tof3_counts", + ] + exposure_factor_init = np.full((2, 7, 60), 1) + exposure_factor_init[:, 0, :] = 2 + exposure_factor_expected = dict() + for field in spin_bin_6_fields + spin_bin_60_fields: + exposure_factor_expected[field] = exposure_factor_init l1b_histrate.h_counts[0, 0, 0] = 5 l1b_histrate.h_counts[0, 1, 0] = 10 @@ -834,23 +906,61 @@ def test_resweep_histogram_success(anc_dependencies): assert l1b_histrates.o_counts[1, 1, 0] == 0 assert l1b_histrates.o_counts[1, 2, 0] == 4 - assert np.array_equal(exposure_factor, exposure_factor_expected) + for field in spin_bin_6_fields + spin_bin_60_fields: + assert np.array_equal(l1b_histrates[field], l1b_histrate[field]) -def test_resweep_histogram_no_date(anc_dependencies): +def test_resweep_histogram_no_date_in_sweep(anc_dependencies, caplog): # Arrange epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-25T02:00:00", "2025-04-25T03:00:00"]) ) l1b_histrate = xr.Dataset( { - "h_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), + "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "tof0_tof1_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof0_tof2_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof1_tof3_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "silver_triple_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "start_a_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "start_c_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b0_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b3_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), }, coords={ "epoch": epoch_date, "esa_step": np.arange(1, 8), "spin_bin_6": np.arange(60), + "spin_bin_60": np.arange(6), }, ) @@ -858,36 +968,144 @@ def test_resweep_histogram_no_date(anc_dependencies): l1b_histrate.h_counts[0, 1, 0] = 10 l1b_histrate.h_counts[0, 2, 0] = 2 - with pytest.raises( - ValueError, - match="No sweep table entry found for date " - "2025-04-25T02:00:00.000 at epoch idx 0", - ): - resweep_histogram_data(l1b_histrate, anc_dependencies) + with caplog.at_level(logging.WARNING): + result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + + # Check that warning was logged + assert any( + "not found in sweep table" in record.message for record in caplog.records + ) -def test_resweep_histogram_multiple_lut(anc_dependencies): +def test_resweep_histogram_no_table_in_lut(anc_dependencies, caplog): + # Arrange + epoch_date = et_to_ttj2000ns( + str_to_et(["2024-01-01T02:00:00", "2024-01-01T03:00:00"]) + ) + l1b_histrate = xr.Dataset( + { + "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "tof0_tof1_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof0_tof2_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof1_tof3_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "silver_triple_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "start_a_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "start_c_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b0_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b3_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + }, + coords={ + "epoch": epoch_date, + "esa_step": np.arange(1, 8), + "spin_bin_6": np.arange(60), + "spin_bin_60": np.arange(6), + }, + ) + + l1b_histrate.h_counts[0, 0, 0] = 5 + l1b_histrate.h_counts[0, 1, 0] = 10 + l1b_histrate.h_counts[0, 2, 0] = 2 + + with caplog.at_level(logging.WARNING): + result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + + # Check that warning was logged + assert any( + "No LUT entries for epoch" in record.message for record in caplog.records + ) + + +def test_resweep_histogram_multiple_lut(anc_dependencies, caplog): epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-16T02:00:00", "2025-04-16T03:00:00"]) ) l1b_histrate = xr.Dataset( { - "h_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), + "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), + "tof0_tof1_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof0_tof2_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "tof1_tof3_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "silver_triple_counts": ( + ("epoch", "esa_step", "spin_bin_6"), + np.zeros((2, 7, 60)), + ), + "start_a_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "start_c_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b0_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "stop_b3_counts": ( + ("epoch", "esa_step", "spin_bin_60"), + np.zeros((2, 7, 6)), + ), + "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), + "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), }, coords={ "epoch": epoch_date, "esa_step": np.arange(1, 8), "spin_bin_6": np.arange(60), + "spin_bin_60": np.arange(6), }, ) - with pytest.raises( - ValueError, - match=f"Expected exactly 1 unique LUT_table " - f"value for date 2025-04-16, but found 2:{[1, 2]}", - ): - resweep_histogram_data(l1b_histrate, anc_dependencies) + with caplog.at_level(logging.WARNING): + result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + + # Check that warning was logged + assert any( + "Multiple LUT tables found for epoch" in record.message + for record in caplog.records + ) + assert any("but found tables" in record.message for record in caplog.records) def test_calculate_histogram_rates(l1b_histrates): From 0f075ccb34db5ea8c3e5c6e8212809d9f630bf6a Mon Sep 17 00:00:00 2001 From: Sean Hoyt Date: Wed, 7 Jan 2026 10:17:40 -0700 Subject: [PATCH 2/5] temp --- imap_processing/lo/l1b/lo_l1b.py | 125 +++++++++++++++++-------------- 1 file changed, 67 insertions(+), 58 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 25dd2bdb0a..6287835c24 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -1586,87 +1586,96 @@ def resweep_histogram_data( def calculate_histogram_rates( - l1b_histrates: xr.Dataset, - acq_start: xr.DataArray, - acq_end: xr.DataArray, - avg_spin_durations_per_cycle: xr.DataArray, - exposure_factor: np.ndarray, + l1b_histrates: xr.Dataset, + acq_start: xr.DataArray, + acq_end: xr.DataArray, + avg_spin_durations_per_cycle: xr.DataArray, + exposure_factors: dict[str, np.ndarray], ) -> xr.Dataset: """ Calculate histogram rates by dividing reswept counts by exposure time. For each epoch in l1b_histrates, this function finds the corresponding - spin interval, calculates the exposure time for 6-degree bins, + spin interval, calculates the exposure time for each bin type, and divides the counts by the exposure time. The exposure time is scaled by the number of ESA steps that were reswept during resweeping. Parameters ---------- l1b_histrates : xr.Dataset - The L1B histogram rates dataset containing reswept h_counts and o_counts. + The L1B histogram rates dataset containing reswept counts. acq_start : xr.DataArray Start times for each spin cycle in MET seconds. acq_end : xr.DataArray End times for each spin cycle in MET seconds. avg_spin_durations_per_cycle : xr.DataArray Average spin duration for each cycle in seconds. - exposure_factor : np.ndarray - 3D array of exposure factors (epoch, azimuth, esa_step) indicating how many - ESA steps were reswept during resweeping. + exposure_factors : dict[str, np.ndarray] + Dictionary mapping field names to their 3D exposure factor arrays + (epoch, esa_step, azimuth) indicating how many ESA steps were + reswept during resweeping. Returns ------- l1b_histrates : xr.Dataset - Updated dataset with h_rates and o_rates added. + The L1B histogram rates dataset with rates calculated. """ - epochs = l1b_histrates["epoch"].values - h_counts = l1b_histrates["h_counts"].values - o_counts = l1b_histrates["o_counts"].values - - h_rates = np.zeros_like(h_counts, dtype=float) - o_rates = np.zeros_like(o_counts, dtype=float) - num_azimuth = h_counts.shape[2] - exposure_times = np.zeros((len(epochs), 7, num_azimuth), dtype=float) - - # Calculate rates for each epoch - for epoch_idx, epoch in enumerate(epochs): - # Find the spin cycle that contains the current epoch - spin_cycle_mask = (epoch >= met_to_ttj2000ns(acq_start.values)) & ( - epoch <= met_to_ttj2000ns(acq_end.values) - ) - spin_cycle_indices = np.nonzero(spin_cycle_mask)[0] + epochs_ttj2000 = l1b_histrates["epoch"].values + epochs_met = ttj2000ns_to_met(epochs_ttj2000) - # If no matching spin cycle is found, log a warning and set rates to NaN - if len(spin_cycle_indices) == 0: - logger.warning(f"Epoch {epoch_idx} not found in any spin_cycle interval") - h_rates[epoch_idx] = np.nan - o_rates[epoch_idx] = np.nan - continue + # Match each histogram epoch to its corresponding spin cycle + closest_spin_idx = np.abs(epochs_met[:, None] - acq_start.values).argmin(axis=1) - spin_cycle_idx = spin_cycle_indices[0] - # Calculate the base exposure time for the spin cycle in minutes - base_exposure_time = ( - 4 * avg_spin_durations_per_cycle.values[spin_cycle_idx] / 60 - ) - # Scale the exposure time by the exposure factor from resweeping - scaled_exposure = base_exposure_time * exposure_factor[epoch_idx, ...] - # Avoid division by zero by setting zero exposure times to NaN - exposure_times[epoch_idx, ...] = scaled_exposure - with np.errstate(divide="ignore"): - h_rates[epoch_idx, ...] = h_counts[epoch_idx, ...] / scaled_exposure - o_rates[epoch_idx, ...] = o_counts[epoch_idx, ...] / scaled_exposure - - l1b_histrates["exposure_time"] = xr.DataArray( - exposure_times, - dims=["epoch", "esa_step", "spin_bin_6"], - ) - l1b_histrates["h_rates"] = xr.DataArray( - h_rates, - dims=l1b_histrates["h_counts"].dims, - ) - l1b_histrates["o_rates"] = xr.DataArray( - o_rates, - dims=l1b_histrates["o_counts"].dims, - ) + # Get spin durations for each epoch + spin_durations = avg_spin_durations_per_cycle.values[closest_spin_idx] + + # Calculate exposure time for 6-degree bins (60 bins per spin) + exposure_time_6deg = spin_durations / 60 + + # Calculate exposure time for 0.6-degree bins (600 bins per spin) + exposure_time_06deg = spin_durations / 600 + + # Process all fields + spin_bin_6_fields = ["h_counts", "o_counts", "tof0_tof1_counts", + "tof0_tof2_counts", "tof1_tof3_counts", + "silver_triple_counts"] + spin_bin_60_fields = ["start_a_counts", "start_c_counts", "stop_b0_counts", + "stop_b3_counts", "tof0_counts", "tof1_counts", + "tof2_counts", "tof3_counts"] + + # Process 6-degree bin fields + for field in spin_bin_6_fields: + counts = l1b_histrates[field].values # (epoch, esa_step, spin_bin_6) + exp_factor = exposure_factors[field] # (epoch, esa_step, spin_bin_6) + + # Calculate effective exposure time with broadcasting + # Shape: (epoch, 1, 1) * (epoch, esa_step, spin_bin_6) + effective_exposure = exposure_time_6deg[:, None, None] * exp_factor + + # Avoid division by zero + with np.errstate(divide='ignore', invalid='ignore'): + rates = np.where(effective_exposure > 0, + counts / effective_exposure, + 0) + + l1b_histrates[field].values = rates + + # Process 0.6-degree bin fields + for field in spin_bin_60_fields: + counts = l1b_histrates[field].values # (epoch, esa_step, spin_bin_60) + exp_factor = exposure_factors[field] # (epoch, esa_step, spin_bin_60) + + # Calculate effective exposure time with broadcasting + # Shape: (epoch, 1, 1) * (epoch, esa_step, spin_bin_60) + effective_exposure = exposure_time_06deg[:, None, None] * exp_factor + + # Avoid division by zero + with np.errstate(divide='ignore', invalid='ignore'): + rates = np.where(effective_exposure > 0, + counts / effective_exposure, + 0) + + l1b_histrates[field].values = rates return l1b_histrates + From 69ad6e51d714c1f52adb8b27e6f6b18fe40cafba Mon Sep 17 00:00:00 2001 From: Sean Hoyt Date: Thu, 15 Jan 2026 15:43:26 -0700 Subject: [PATCH 3/5] working monitor rate code --- imap_processing/lo/l1b/lo_l1b.py | 327 ++++++++++++++------ imap_processing/tests/lo/test_lo_l1b.py | 384 ++++++++---------------- 2 files changed, 356 insertions(+), 355 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 6287835c24..1b5a639610 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -33,6 +33,134 @@ logger = logging.getLogger(__name__) +# ------------------------------------------------------------------- +# Centralized field definitions to avoid repetition across functions +# ------------------------------------------------------------------- +# spin-bin fields (count fields used in multiple places) +SPIN_BIN_6_FIELDS = [ + "h_counts", + "o_counts", + "tof0_tof1_counts", + "tof0_tof2_counts", + "tof1_tof2_counts", + "silver_triple_counts", +] + +SPIN_BIN_60_FIELDS = [ + "start_a_counts", + "start_c_counts", + "stop_b0_counts", + "stop_b3_counts", + "tof0_counts", + "tof1_counts", + "tof2_counts", + "tof3_counts", + "disc_tof0_counts", + "disc_tof1_counts", + "disc_tof2_counts", + "disc_tof3_counts", + "pos0_counts", + "pos1_counts", + "pos2_counts", + "pos3_counts", +] + +# Mapping from L1A field names to L1B count field names used in initialize_all_rates +SPIN_BIN_6_L1A_TO_L1B = { + "hydrogen": "h_counts", + "oxygen": "o_counts", + "tof0_tof1": "tof0_tof1_counts", + "tof0_tof2": "tof0_tof2_counts", + "tof1_tof2": "tof1_tof2_counts", + "silver": "silver_triple_counts", +} + +SPIN_BIN_60_L1A_TO_L1B = { + "start_a": "start_a_counts", + "start_c": "start_c_counts", + "stop_b0": "stop_b0_counts", + "stop_b3": "stop_b3_counts", + "tof0_count": "tof0_counts", + "tof1_count": "tof1_counts", + "tof2_count": "tof2_counts", + "tof3_count": "tof3_counts", + "disc_tof0": "disc_tof0_counts", + "disc_tof1": "disc_tof1_counts", + "disc_tof2": "disc_tof2_counts", + "disc_tof3": "disc_tof3_counts", + "pos0": "pos0_counts", + "pos1": "pos1_counts", + "pos2": "pos2_counts", + "pos3": "pos3_counts", +} + +# Count-field -> rate-field mappings used by calculate_histogram_rates +SPIN_BIN_6_COUNT_TO_RATE = { + "h_counts": "h_rates", + "o_counts": "o_rates", + "tof0_tof1_counts": "tof0_tof1_rates", + "tof0_tof2_counts": "tof0_tof2_rates", + "tof1_tof2_counts": "tof1_tof2_rates", + "silver_triple_counts": "silver_triple_rates", +} + +SPIN_BIN_60_COUNT_TO_RATE = { + "start_a_counts": "start_a_rates", + "start_c_counts": "start_c_rates", + "stop_b0_counts": "stop_b0_rates", + "stop_b3_counts": "stop_b3_rates", + "tof0_counts": "tof0_rates", + "tof1_counts": "tof1_rates", + "tof2_counts": "tof2_rates", + "tof3_counts": "tof3_rates", + "disc_tof0_counts": "disc_tof0_rates", + "disc_tof1_counts": "disc_tof1_rates", + "disc_tof2_counts": "disc_tof2_rates", + "disc_tof3_counts": "disc_tof3_rates", + "pos0_counts": "pos0_rates", + "pos1_counts": "pos1_rates", + "pos2_counts": "pos2_rates", + "pos3_counts": "pos3_rates", +} + +# Fields to include in the split hist/monitor rate datasets +HIST_RATE_FIELDS = [ + "h_rates", + "o_rates", + "h_counts", + "o_counts", + "esa_mode", + "exposure_time_6deg", + "spin_cycle", +] +MONITOR_RATE_FIELDS = [ + "tof0_tof1_rates", + "tof0_tof2_rates", + "tof1_tof2_rates", + "silver_triple_rates", + "start_a_rates", + "start_c_rates", + "stop_b0_rates", + "stop_b3_rates", + "tof0_rates", + "tof1_rates", + "tof2_rates", + "tof3_rates", + "disc_tof0_rates", + "disc_tof1_rates", + "disc_tof2_rates", + "disc_tof3_rates", + "pos0_rates", + "pos1_rates", + "pos2_rates", + "pos3_rates", + "esa_mode", + "exposure_time_60deg", + "exposure_time_6deg", + "spin_cycle", +] +# ------------------------------------------------------------------- + def lo_l1b(sci_dependencies: dict, anc_dependencies: list) -> list[Path]: """ @@ -136,7 +264,7 @@ def lo_l1b(sci_dependencies: dict, anc_dependencies: list) -> list[Path]: # This carries over the epoch and count fields from L1A l1b_all_rates = initialize_all_rates(l1a_hist, attr_mgr_l1b) # set spin cycle and remove invalid spin ASCs - l1b_histrates = set_spin_cycle_from_spin_data( + l1b_all_rates = set_spin_cycle_from_spin_data( l1a_hist, l1b_all_rates, spin_data ) @@ -156,14 +284,18 @@ def lo_l1b(sci_dependencies: dict, anc_dependencies: list) -> list[Path]: avg_spin_durations_per_cycle = get_avg_spin_durations_per_cycle( acq_start, acq_end ) - # l1b_histrates = calculate_histogram_rates( - # l1b_histrates, - # acq_start, - # acq_end, - # avg_spin_durations_per_cycle, - # exposure_factor, - # ) - datasets_to_return.append(l1b_histrates) + l1b_all_rates = calculate_histogram_rates( + l1b_all_rates, + acq_start, + acq_end, + avg_spin_durations_per_cycle, + exposure_factor, + ) + + l1b_hist_rates, l1b_monitor_rates = split_rate_dataset( + l1b_all_rates, attr_mgr_l1b + ) + datasets_to_return.extend([l1b_hist_rates, l1b_monitor_rates]) return datasets_to_return @@ -1402,38 +1534,14 @@ def initialize_all_rates( ), }, ) - # Define field mappings from L1A to the L1B dataset - # according to the algorithm document, only the histogram rate product should have - # counts (hydrogen and oxygen). The rest of the fields will be dropped after - # processing when the monitor rates product is created. These fields are needed for - # processing though and are added to the all_rates dataset temporarily. - # {"l1a_field_name": "l1b_all_rates_field_name"} - spin_bin_6_fields = { - "hydrogen": "h_counts", - "oxygen": "o_counts", - "tof0_tof1": "tof0_tof1_counts", - "tof0_tof2": "tof0_tof2_counts", - "tof1_tof3": "tof1_tof3_counts", - "silver_triple": "silver_triple_counts", - } - spin_bin_60_fields = { - "start_a": "start_a_counts", - "start_c": "start_c_counts", - "stop_b0": "stop_b0_counts", - "stop_b3": "stop_b3_counts", - "tof0": "tof0_counts", - "tof1": "tof1_counts", - "tof2": "tof2_counts", - "tof3": "tof3_counts", - } - - for l1a_field, l1b_field in spin_bin_6_fields.items(): + # Use centralized mappings for field definitions + for l1a_field, l1b_field in SPIN_BIN_6_L1A_TO_L1B.items(): l1b_all_rates[l1b_field] = xr.DataArray( l1a_hist[l1a_field].values, dims=["epoch", "esa_step", "spin_bin_6"], ) - for l1a_field, l1b_field in spin_bin_60_fields.items(): + for l1a_field, l1b_field in SPIN_BIN_60_L1A_TO_L1B.items(): l1b_all_rates[l1b_field] = xr.DataArray( l1a_hist[l1a_field].values, dims=["epoch", "esa_step", "spin_bin_60"], @@ -1475,33 +1583,13 @@ def resweep_histogram_data( next(str(s) for s in anc_dependencies if "esa-mode-lut" in str(s)) ) - # Define field groups with their spin bin dimensions - spin_bin_6_fields = [ - "h_counts", - "o_counts", - "tof0_tof1_counts", - "tof0_tof2_counts", - "tof1_tof3_counts", - "silver_triple_counts", - ] - spin_bin_60_fields = [ - "start_a_counts", - "start_c_counts", - "stop_b0_counts", - "stop_b3_counts", - "tof0_counts", - "tof1_counts", - "tof2_counts", - "tof3_counts", - ] - sweep_dates = sweep_df["Date"].astype(str) epochs = l1b_histrates["epoch"].values epoch_utc = et_to_utc(ttj2000ns_to_et(epochs)) # Initialize reswept arrays for all fields reswept_data = {} - for field in spin_bin_6_fields + spin_bin_60_fields: + for field in SPIN_BIN_6_FIELDS + SPIN_BIN_60_FIELDS: reswept_data[field] = np.zeros_like(l1b_histrates[field].values) # Initialize exposure factors for each field type @@ -1509,12 +1597,12 @@ def resweep_histogram_data( num_azimuth_60 = l1b_histrates.sizes["spin_bin_60"] exposure_factors = {} - for field in spin_bin_6_fields: + for field in SPIN_BIN_6_FIELDS: exposure_factors[field] = np.full( (len(epochs), l1b_histrates.sizes["esa_step"], num_azimuth_6), 1, dtype=int ) - for field in spin_bin_60_fields: + for field in SPIN_BIN_60_FIELDS: exposure_factors[field] = np.full( (len(epochs), l1b_histrates.sizes["esa_step"], num_azimuth_60), 1, dtype=int ) @@ -1557,7 +1645,7 @@ def resweep_histogram_data( energy_step_mapping[esa_idx] = true_esa_step # Process spin_bin_6 fields - for field in spin_bin_6_fields: + for field in SPIN_BIN_6_FIELDS: for az_idx in range(num_azimuth_6): original = l1b_histrates[field].values[epoch_idx, :, az_idx] for orig_idx, true_esa_step in energy_step_mapping.items(): @@ -1568,7 +1656,7 @@ def resweep_histogram_data( exposure_factors[field][epoch_idx, target_idx, az_idx] += 1 # Process spin_bin_60 fields - for field in spin_bin_60_fields: + for field in SPIN_BIN_60_FIELDS: for az_idx in range(num_azimuth_60): original = l1b_histrates[field].values[epoch_idx, :, az_idx] for orig_idx, true_esa_step in energy_step_mapping.items(): @@ -1579,18 +1667,18 @@ def resweep_histogram_data( exposure_factors[field][epoch_idx, target_idx, az_idx] += 1 # Update dataset with reswept data - for field in spin_bin_6_fields + spin_bin_60_fields: + for field in SPIN_BIN_6_FIELDS + SPIN_BIN_60_FIELDS: l1b_histrates[field].values = reswept_data[field] return l1b_histrates, exposure_factors def calculate_histogram_rates( - l1b_histrates: xr.Dataset, - acq_start: xr.DataArray, - acq_end: xr.DataArray, - avg_spin_durations_per_cycle: xr.DataArray, - exposure_factors: dict[str, np.ndarray], + l1b_histrates: xr.Dataset, + acq_start: xr.DataArray, + acq_end: xr.DataArray, + avg_spin_durations_per_cycle: xr.DataArray, + exposure_factors: dict[str, np.ndarray], ) -> xr.Dataset: """ Calculate histogram rates by dividing reswept counts by exposure time. @@ -1633,49 +1721,98 @@ def calculate_histogram_rates( exposure_time_6deg = spin_durations / 60 # Calculate exposure time for 0.6-degree bins (600 bins per spin) - exposure_time_06deg = spin_durations / 600 + exposure_time_60deg = spin_durations / 6 # Process all fields - spin_bin_6_fields = ["h_counts", "o_counts", "tof0_tof1_counts", - "tof0_tof2_counts", "tof1_tof3_counts", - "silver_triple_counts"] - spin_bin_60_fields = ["start_a_counts", "start_c_counts", "stop_b0_counts", - "stop_b3_counts", "tof0_counts", "tof1_counts", - "tof2_counts", "tof3_counts"] - # Process 6-degree bin fields - for field in spin_bin_6_fields: - counts = l1b_histrates[field].values # (epoch, esa_step, spin_bin_6) - exp_factor = exposure_factors[field] # (epoch, esa_step, spin_bin_6) + for count_field, rate_field in SPIN_BIN_6_COUNT_TO_RATE.items(): + counts = l1b_histrates[count_field].values # (epoch, esa_step, spin_bin_6) + exp_factor = exposure_factors[count_field] # (epoch, esa_step, spin_bin_6) # Calculate effective exposure time with broadcasting - # Shape: (epoch, 1, 1) * (epoch, esa_step, spin_bin_6) effective_exposure = exposure_time_6deg[:, None, None] * exp_factor # Avoid division by zero - with np.errstate(divide='ignore', invalid='ignore'): - rates = np.where(effective_exposure > 0, - counts / effective_exposure, - 0) + with np.errstate(divide="ignore", invalid="ignore"): + rates = np.where(effective_exposure > 0, counts / effective_exposure, 0) - l1b_histrates[field].values = rates + l1b_histrates[rate_field] = xr.DataArray( + rates, + dims=l1b_histrates[count_field].dims, + ) - # Process 0.6-degree bin fields - for field in spin_bin_60_fields: - counts = l1b_histrates[field].values # (epoch, esa_step, spin_bin_60) - exp_factor = exposure_factors[field] # (epoch, esa_step, spin_bin_60) + l1b_histrates["exposure_time_6deg"] = xr.DataArray( + effective_exposure, + dims=["epoch", "esa_step", "spin_bin_6"], + ) + + # Process 60-degree bin fields + for count_field, rate_field in SPIN_BIN_60_COUNT_TO_RATE.items(): + counts = l1b_histrates[count_field].values # (epoch, esa_step, spin_bin_60) + exp_factor = exposure_factors[count_field] # (epoch, esa_step, spin_bin_60) # Calculate effective exposure time with broadcasting - # Shape: (epoch, 1, 1) * (epoch, esa_step, spin_bin_60) - effective_exposure = exposure_time_06deg[:, None, None] * exp_factor + effective_exposure = exposure_time_60deg[:, None, None] * exp_factor # Avoid division by zero - with np.errstate(divide='ignore', invalid='ignore'): - rates = np.where(effective_exposure > 0, - counts / effective_exposure, - 0) + with np.errstate(divide="ignore", invalid="ignore"): + rates = np.where(effective_exposure > 0, counts / effective_exposure, 0) - l1b_histrates[field].values = rates + l1b_histrates[rate_field] = xr.DataArray( + rates, + dims=l1b_histrates[count_field].dims, + ) + + l1b_histrates["exposure_time_60deg"] = xr.DataArray( + effective_exposure, + dims=["epoch", "esa_step", "spin_bin_60"], + ) return l1b_histrates + +def split_rate_dataset( + l1b_all_rates: xr.Dataset, attr_mgr_l1b: ImapCdfAttributes +) -> tuple[xr.Dataset, xr.Dataset]: + """ + Split the L1B all rates dataset into histogram rates and monitor rates datasets. + + Parameters + ---------- + l1b_all_rates : xr.Dataset + The L1B all rates dataset containing both histogram and monitor rates. + attr_mgr_l1b : ImapCdfAttributes + Attribute manager used to get the L1B histogram and monitor rates dataset + attributes. + + Returns + ------- + l1b_hist_rates : xr.Dataset + The L1B histogram rates dataset. + l1b_monitor_rates : xr.Dataset + The L1B monitor rates dataset. + """ + # Use centralized lists for fields to include in split datasets + l1b_hist_rates = xr.Dataset( + attrs=attr_mgr_l1b.get_global_attributes("imap_lo_l1b_histrates"), + coords=l1b_all_rates.coords, + ) + l1b_monitor_rates = xr.Dataset( + attrs=attr_mgr_l1b.get_global_attributes("imap_lo_l1b_monitorrates"), + coords=l1b_all_rates.coords, + ) + + for field in HIST_RATE_FIELDS: + l1b_hist_rates[field] = xr.DataArray( + l1b_all_rates[field].values, + dims=l1b_all_rates[field].dims, + # attrs=attr_mgr_l1b.get_variable_attributes(field) + ) + for field in MONITOR_RATE_FIELDS: + l1b_monitor_rates[field] = xr.DataArray( + l1b_all_rates[field].values, + dims=l1b_all_rates[field].dims, + # attrs=attr_mgr_l1b.get_variable_attributes(field) + ) + + return l1b_hist_rates, l1b_monitor_rates diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 4dc48e136b..e298b648e9 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -44,6 +44,34 @@ ttj2000ns_to_met, ) +SPIN_BIN_6_FIELDS = [ + "h_counts", + "o_counts", + "tof0_tof1_counts", + "tof0_tof2_counts", + "tof1_tof2_counts", + "silver_triple_counts", +] + +SPIN_BIN_60_FIELDS = [ + "start_a_counts", + "start_c_counts", + "stop_b0_counts", + "stop_b3_counts", + "tof0_counts", + "tof1_counts", + "tof2_counts", + "tof3_counts", + "disc_tof0_counts", + "disc_tof1_counts", + "disc_tof2_counts", + "disc_tof3_counts", + "pos0_counts", + "pos1_counts", + "pos2_counts", + "pos3_counts", +] + @pytest.fixture def dependencies(): @@ -97,15 +125,21 @@ def l1b_histrates(): epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-15T02:00:00", "2025-04-15T03:00:00"]) ) + + # Build dataset with all expected fields + data_vars = {} + for f in SPIN_BIN_6_FIELDS: + data_vars[f] = (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))) + for f in SPIN_BIN_60_FIELDS: + data_vars[f] = (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))) + l1b_histrates = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "azimuth_6"), np.zeros((2, 7, 60))), - }, + data_vars, coords={ "epoch": epoch_date, "esa_step": np.arange(1, 8), - "azimuth_6": np.arange(60), + "spin_bin_6": np.arange(60), + "spin_bin_60": np.arange(6), }, ) @@ -121,16 +155,24 @@ def l1a_hist(): "oxygen": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), "tof0_tof1": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), "tof0_tof2": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), - "tof1_tof3": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), - "silver_triple": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "tof1_tof2": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), + "silver": (("epoch", "esa_step", "azimuth_6"), np.zeros((1, 7, 60))), "start_a": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), "start_c": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), "stop_b0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), "stop_b3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), - "tof0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), - "tof1": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), - "tof2": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), - "tof3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof0_count": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof1_count": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof2_count": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "tof3_count": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "disc_tof0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "disc_tof1": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "disc_tof2": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "disc_tof3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "pos0": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "pos1": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "pos2": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), + "pos3": (("epoch", "esa_step", "azimuth_60"), np.zeros((1, 7, 6))), }, coords={ "epoch": epoch_date, @@ -222,11 +264,11 @@ def test_lo_l1b_histogram_rates( l1b_datasets = lo_l1b(sci_dependencies, anc_dependencies) # Assert - assert "h_rates" in l1b_datasets[-1].data_vars - assert "o_rates" in l1b_datasets[-1].data_vars - assert "exposure_time" in l1b_datasets[-1].data_vars - assert "h_counts" in l1b_datasets[-1].data_vars - assert "o_counts" in l1b_datasets[-1].data_vars + assert "h_rates" in l1b_datasets[-2].data_vars + assert "o_rates" in l1b_datasets[-2].data_vars + assert "exposure_time_6deg" in l1b_datasets[-2].data_vars + assert "h_counts" in l1b_datasets[-2].data_vars + assert "o_counts" in l1b_datasets[-2].data_vars # @pytest.mark.external_kernel @@ -809,93 +851,28 @@ def test_badtimes_with_spin(spice_test_data_path, use_test_spin_data_csv): np.testing.assert_array_equal(badtimes_ds["badtime_flag"], 1) -def test_resweep_histogram_success(anc_dependencies): +def test_resweep_histogram_success(l1b_histrates, anc_dependencies): # Arrange epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-15T02:00:00", "2025-04-15T03:00:00"]) ) - l1b_histrate = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "tof0_tof1_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof0_tof2_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof1_tof3_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "silver_triple_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "start_a_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "start_c_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b0_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b3_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - }, - coords={ - "epoch": epoch_date, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - "spin_bin_60": np.arange(6), - }, - ) - spin_bin_6_fields = [ - "h_counts", - "o_counts", - "tof0_tof1_counts", - "tof0_tof2_counts", - "tof1_tof3_counts", - "silver_triple_counts", - ] - spin_bin_60_fields = [ - "start_a_counts", - "start_c_counts", - "stop_b0_counts", - "stop_b3_counts", - "tof0_counts", - "tof1_counts", - "tof2_counts", - "tof3_counts", - ] + l1b_histrates["epoch"] = epoch_date exposure_factor_init = np.full((2, 7, 60), 1) exposure_factor_init[:, 0, :] = 2 exposure_factor_expected = dict() - for field in spin_bin_6_fields + spin_bin_60_fields: + for field in SPIN_BIN_6_FIELDS + SPIN_BIN_60_FIELDS: exposure_factor_expected[field] = exposure_factor_init - l1b_histrate.h_counts[0, 0, 0] = 5 - l1b_histrate.h_counts[0, 1, 0] = 10 - l1b_histrate.h_counts[0, 2, 0] = 2 + l1b_histrates.h_counts[0, 0, 0] = 5 + l1b_histrates.h_counts[0, 1, 0] = 10 + l1b_histrates.h_counts[0, 2, 0] = 2 - l1b_histrate.o_counts[1, 0, 0] = 2 - l1b_histrate.o_counts[1, 1, 0] = 3 - l1b_histrate.o_counts[1, 2, 0] = 4 + l1b_histrates.o_counts[1, 0, 0] = 2 + l1b_histrates.o_counts[1, 1, 0] = 3 + l1b_histrates.o_counts[1, 2, 0] = 4 l1b_histrates, exposure_factor = resweep_histogram_data( - l1b_histrate, anc_dependencies + l1b_histrates, anc_dependencies ) assert l1b_histrates.h_counts[0, 0, 0] == 15 @@ -906,70 +883,23 @@ def test_resweep_histogram_success(anc_dependencies): assert l1b_histrates.o_counts[1, 1, 0] == 0 assert l1b_histrates.o_counts[1, 2, 0] == 4 - for field in spin_bin_6_fields + spin_bin_60_fields: - assert np.array_equal(l1b_histrates[field], l1b_histrate[field]) + for field in SPIN_BIN_6_FIELDS + SPIN_BIN_60_FIELDS: + assert np.array_equal(l1b_histrates[field], l1b_histrates[field]) -def test_resweep_histogram_no_date_in_sweep(anc_dependencies, caplog): +def test_resweep_histogram_no_date_in_sweep(l1b_histrates, anc_dependencies, caplog): # Arrange epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-25T02:00:00", "2025-04-25T03:00:00"]) ) - l1b_histrate = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "tof0_tof1_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof0_tof2_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof1_tof3_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "silver_triple_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "start_a_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "start_c_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b0_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b3_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - }, - coords={ - "epoch": epoch_date, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - "spin_bin_60": np.arange(6), - }, - ) + l1b_histrates["epoch"] = epoch_date - l1b_histrate.h_counts[0, 0, 0] = 5 - l1b_histrate.h_counts[0, 1, 0] = 10 - l1b_histrate.h_counts[0, 2, 0] = 2 + l1b_histrates.h_counts[0, 0, 0] = 5 + l1b_histrates.h_counts[0, 1, 0] = 10 + l1b_histrates.h_counts[0, 2, 0] = 2 with caplog.at_level(logging.WARNING): - result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + result, _ = resweep_histogram_data(l1b_histrates, anc_dependencies) # Check that warning was logged assert any( @@ -977,66 +907,19 @@ def test_resweep_histogram_no_date_in_sweep(anc_dependencies, caplog): ) -def test_resweep_histogram_no_table_in_lut(anc_dependencies, caplog): +def test_resweep_histogram_no_table_in_lut(l1b_histrates, anc_dependencies, caplog): # Arrange epoch_date = et_to_ttj2000ns( str_to_et(["2024-01-01T02:00:00", "2024-01-01T03:00:00"]) ) - l1b_histrate = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "tof0_tof1_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof0_tof2_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof1_tof3_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "silver_triple_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "start_a_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "start_c_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b0_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b3_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - }, - coords={ - "epoch": epoch_date, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - "spin_bin_60": np.arange(6), - }, - ) + l1b_histrates["epoch"] = epoch_date - l1b_histrate.h_counts[0, 0, 0] = 5 - l1b_histrate.h_counts[0, 1, 0] = 10 - l1b_histrate.h_counts[0, 2, 0] = 2 + l1b_histrates.h_counts[0, 0, 0] = 5 + l1b_histrates.h_counts[0, 1, 0] = 10 + l1b_histrates.h_counts[0, 2, 0] = 2 with caplog.at_level(logging.WARNING): - result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + result, _ = resweep_histogram_data(l1b_histrates, anc_dependencies) # Check that warning was logged assert any( @@ -1044,61 +927,15 @@ def test_resweep_histogram_no_table_in_lut(anc_dependencies, caplog): ) -def test_resweep_histogram_multiple_lut(anc_dependencies, caplog): +def test_resweep_histogram_multiple_lut(l1b_histrates, anc_dependencies, caplog): epoch_date = et_to_ttj2000ns( str_to_et(["2025-04-16T02:00:00", "2025-04-16T03:00:00"]) ) - l1b_histrate = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), np.zeros((2, 7, 60))), - "tof0_tof1_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof0_tof2_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "tof1_tof3_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "silver_triple_counts": ( - ("epoch", "esa_step", "spin_bin_6"), - np.zeros((2, 7, 60)), - ), - "start_a_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "start_c_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b0_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "stop_b3_counts": ( - ("epoch", "esa_step", "spin_bin_60"), - np.zeros((2, 7, 6)), - ), - "tof0_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof1_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof2_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - "tof3_counts": (("epoch", "esa_step", "spin_bin_60"), np.zeros((2, 7, 6))), - }, - coords={ - "epoch": epoch_date, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - "spin_bin_60": np.arange(6), - }, - ) + + l1b_histrates["epoch"] = epoch_date with caplog.at_level(logging.WARNING): - result, _ = resweep_histogram_data(l1b_histrate, anc_dependencies) + result, _ = resweep_histogram_data(l1b_histrates, anc_dependencies) # Check that warning was logged assert any( @@ -1122,8 +959,17 @@ def test_calculate_histogram_rates(l1b_histrates): ] ) avg_spin_durations_per_cycle = xr.DataArray([30, 15]) - exposure_factor = np.zeros((2, 7, 60)) - exposure_factor[0, 0, 0] = 1 + + exposure_factors = {} + # default zeros then set a sample exposure as in original test intent + for field in SPIN_BIN_6_FIELDS: + arr = np.zeros((2, 7, 60)) + arr[0, 0, 0] = 1 + exposure_factors[field] = arr + for field in SPIN_BIN_60_FIELDS: + exposure_factors[field] = np.zeros((2, 7, 6)) + + # Populate counts used by assertions l1b_histrates.h_counts[0, 0, 0] = 30 l1b_histrates.h_counts[0, 1, 0] = 10 l1b_histrates.h_counts[0, 2, 0] = 2 @@ -1139,7 +985,11 @@ def test_calculate_histogram_rates(l1b_histrates): l1b_histrates.o_counts[1, 2, 0] = 4 l1b_histrate = calculate_histogram_rates( - l1b_histrates, acq_start, acq_end, avg_spin_durations_per_cycle, exposure_factor + l1b_histrates, + acq_start, + acq_end, + avg_spin_durations_per_cycle, + exposure_factors, ) hist_rates_h_epoch_0 = l1b_histrate["h_rates"] @@ -1177,13 +1027,20 @@ def test_calculate_histogram_rates_no_interval_found(l1b_histrates): ] ) avg_spin_durations_per_cycle = xr.DataArray([30, 15]) - exposure_factor = np.zeros((2, 7, 60)) + + exposure_factors = {f: np.zeros((2, 7, 60)) for f in SPIN_BIN_6_FIELDS} + exposure_factors.update({f: np.zeros((2, 7, 6)) for f in SPIN_BIN_60_FIELDS}) + l1b_histrate = calculate_histogram_rates( - l1b_histrates, acq_start, acq_end, avg_spin_durations_per_cycle, exposure_factor + l1b_histrates, + acq_start, + acq_end, + avg_spin_durations_per_cycle, + exposure_factors, ) - np.testing.assert_array_equal(l1b_histrate["h_rates"], np.full((2, 7, 60), np.nan)) - np.testing.assert_array_equal(l1b_histrate["o_rates"], np.full((2, 7, 60), np.nan)) + np.testing.assert_array_equal(l1b_histrate["h_rates"], np.zeros((2, 7, 60))) + np.testing.assert_array_equal(l1b_histrate["o_rates"], np.zeros((2, 7, 60))) def test_calculate_histogram_rates_zero_exposure_time(l1b_histrates): @@ -1200,13 +1057,20 @@ def test_calculate_histogram_rates_zero_exposure_time(l1b_histrates): ] ) avg_spin_durations_per_cycle = xr.DataArray([0, 15]) - exposure_factor = np.zeros((2, 7, 60)) + + exposure_factors = {f: np.zeros((2, 7, 60)) for f in SPIN_BIN_6_FIELDS} + exposure_factors.update({f: np.zeros((2, 7, 6)) for f in SPIN_BIN_60_FIELDS}) + l1b_histrate = calculate_histogram_rates( - l1b_histrates, acq_start, acq_end, avg_spin_durations_per_cycle, exposure_factor + l1b_histrates, + acq_start, + acq_end, + avg_spin_durations_per_cycle, + exposure_factors, ) - np.testing.assert_array_equal(l1b_histrate["h_rates"], np.full((2, 7, 60), np.nan)) - np.testing.assert_array_equal(l1b_histrate["o_rates"], np.full((2, 7, 60), np.nan)) + np.testing.assert_array_equal(l1b_histrate["h_rates"], np.zeros((2, 7, 60))) + np.testing.assert_array_equal(l1b_histrate["o_rates"], np.zeros((2, 7, 60))) def test_set_spin_cycle_from_spin_data_histogram(): From 0531eeeb060829a0dba8f1e03a5e7bd0564f1000 Mon Sep 17 00:00:00 2001 From: Sean Hoyt Date: Tue, 27 Jan 2026 17:26:15 -0700 Subject: [PATCH 4/5] fixed docstrings --- imap_processing/lo/l1b/lo_l1b.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 0d3bb6c028..06c9446613 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -1572,9 +1572,10 @@ def resweep_histogram_data( ------- l1b_histrates : xr.Dataset The updated L1B histogram rates dataset with reswept counts. - exposure_factor : np.ndarray - 3D array of exposure factors (epoch, azimuth, esa_step) indicating how many - ESA steps were reswept during resweeping. + exposure_factor : dict[str, np.ndarray] + Dictionary mapping bin types to their 3D exposure factor arrays + (epoch, esa_step, azimuth) indicating how many ESA steps were + reswept during resweeping. """ epochs = l1b_histrates["epoch"].values energy_mapping = _get_esa_level_indices(epochs, anc_dependencies=anc_dependencies) @@ -1590,12 +1591,15 @@ def resweep_histogram_data( ) l1b_histrates[field].values = reswept + # Calculate exposure factors for each bin type exposure_factor_6deg = np.zeros_like(l1b_histrates["h_counts"].values, dtype=int) exposure_factor_60deg = np.zeros_like( l1b_histrates["start_a_counts"].values, dtype=int ) np.add.at(exposure_factor_6deg, (slice(None), energy_mapping, slice(None)), 1) np.add.at(exposure_factor_60deg, (slice(None), energy_mapping, slice(None)), 1) + + # Create a dictionary to hold exposure factors for both bin types exposure_factors = {} exposure_factors["6deg"] = exposure_factor_6deg exposure_factors["60deg"] = exposure_factor_60deg @@ -1629,7 +1633,7 @@ def calculate_histogram_rates( avg_spin_durations_per_cycle : xr.DataArray Average spin duration for each cycle in seconds. exposure_factors : dict[str, np.ndarray] - Dictionary mapping field names to their 3D exposure factor arrays + Dictionary mapping bin types to their 3D exposure factor arrays (epoch, esa_step, azimuth) indicating how many ESA steps were reswept during resweeping. From 5f1e8b5b8d6fa147832487daa976211d040938c4 Mon Sep 17 00:00:00 2001 From: Greg Lucas Date: Thu, 29 Jan 2026 09:29:45 -0700 Subject: [PATCH 5/5] Lo L1b change to all-rates to include histogram and monitor --- imap_processing/lo/l1b/lo_l1b.py | 39 ++++++++----------------- imap_processing/tests/lo/test_lo_l1b.py | 2 +- 2 files changed, 13 insertions(+), 28 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index dee6afde8f..020c11b0f2 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -214,9 +214,9 @@ def lo_l1b( datasets_to_return.append(ds) # If dependencies are used to create Histogram Rates - if descriptor == "histrates": - logger.info("\nProcessing IMAP-Lo L1B Histogram Rates...") - ds = l1b_histrates(sci_dependencies, anc_dependencies, attr_mgr_l1b) + if descriptor == "all-rates": + logger.info("\nProcessing IMAP-Lo L1B Hist and Monitor Rates...") + ds = l1b_allrates(sci_dependencies, anc_dependencies, attr_mgr_l1b) datasets_to_return.extend(ds) if descriptor == "derates": @@ -313,7 +313,7 @@ def l1b_de( return l1b_de -def l1b_histrates( +def l1b_allrates( sci_dependencies: dict, anc_dependencies: list, attr_mgr_l1b: ImapCdfAttributes ) -> xr.Dataset: """ @@ -330,8 +330,8 @@ def l1b_histrates( Returns ------- - l1b_histrates : xr.Dataset - The IMAP-Lo L1B Histogram Rates dataset. + [xr.Dataset, xr.Dataset] + The IMAP-Lo L1B Histogram and Monitor Rates datasets. """ datasets_to_return = [] # get the dependency dataset for l1b histogram rates @@ -1738,7 +1738,7 @@ def calculate_histogram_rates( exposure_time_6deg[:, None, None] * exposure_factors["6deg"] ) - # Calculate exposure time for 0.6-degree bins (600 bins per spin) + # Calculate exposure time for 60-degree bins (6 bins per spin) exposure_time_60deg = spin_durations / 6 effective_exposure_60deg = ( exposure_time_60deg[:, None, None] * exposure_factors["60deg"] @@ -2080,28 +2080,13 @@ def split_rate_dataset( The L1B monitor rates dataset. """ # Use centralized lists for fields to include in split datasets - l1b_hist_rates = xr.Dataset( - attrs=attr_mgr_l1b.get_global_attributes("imap_lo_l1b_histrates"), - coords=l1b_all_rates.coords, - ) - l1b_monitor_rates = xr.Dataset( - attrs=attr_mgr_l1b.get_global_attributes("imap_lo_l1b_monitorrates"), - coords=l1b_all_rates.coords, + l1b_hist_rates = l1b_all_rates[HIST_RATE_FIELDS] + l1b_hist_rates.attrs = attr_mgr_l1b.get_global_attributes("imap_lo_l1b_histrates") + l1b_monitor_rates = l1b_all_rates[MONITOR_RATE_FIELDS] + l1b_monitor_rates.attrs = attr_mgr_l1b.get_global_attributes( + "imap_lo_l1b_monitorrates" ) - for field in HIST_RATE_FIELDS: - l1b_hist_rates[field] = xr.DataArray( - l1b_all_rates[field].values, - dims=l1b_all_rates[field].dims, - # attrs=attr_mgr_l1b.get_variable_attributes(field) - ) - for field in MONITOR_RATE_FIELDS: - l1b_monitor_rates[field] = xr.DataArray( - l1b_all_rates[field].values, - dims=l1b_all_rates[field].dims, - # attrs=attr_mgr_l1b.get_variable_attributes(field) - ) - return l1b_hist_rates, l1b_monitor_rates diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 4561e1a618..9f69f3cc03 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -279,7 +279,7 @@ def test_lo_l1b_histogram_rates( } # Act - l1b_datasets = lo_l1b(sci_dependencies, anc_dependencies, descriptor="histrates") + l1b_datasets = lo_l1b(sci_dependencies, anc_dependencies, descriptor="all-rates") # Assert assert "h_rates" in l1b_datasets[-2].data_vars