From ded2f4414e004749fbdab6d0e0e4c1c92b4677a6 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 13 May 2026 17:16:03 -0400 Subject: [PATCH 1/5] psets now take goodtimes/background information from science dependencies instead of ancillary files --- .../cdf/config/imap_lo_global_cdf_attrs.yaml | 2 +- imap_processing/cli.py | 6 + imap_processing/lo/l1b/lo_l1b.py | 105 ------- imap_processing/lo/l1c/lo_l1c.py | 197 ++++++------- imap_processing/tests/lo/test_lo_l1b.py | 42 --- imap_processing/tests/lo/test_lo_l1c.py | 266 ++++++++++-------- 6 files changed, 232 insertions(+), 386 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml index 5f0b47af38..b6e2b8448d 100644 --- a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml @@ -110,7 +110,7 @@ imap_lo_l1b_bgrates: imap_lo_l1b_goodtimes: <<: *instrument_base Data_type: L1B_goodtimes>Level-1B Goodtimes List - Logical_source: imap_lo_l1b_good-times + Logical_source: imap_lo_l1b_goodtimes Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1B Data imap_lo_l1c_goodtimes: diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 9f992e3391..a3ffccd987 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -1230,6 +1230,12 @@ def do_processing( source="lo", data_type="ancillary" ) science_files = dependencies.get_file_paths(source="lo", descriptor="de") + science_files += dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="goodtimes" + ) + science_files += dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="bgrates" + ) for file in science_files: dataset = load_cdf(file) data_dict[dataset.attrs["Logical_source"]] = dataset diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 467b93924f..02b88968ed 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -6,7 +6,6 @@ from pathlib import Path import numpy as np -import pandas as pd import spiceypy import xarray as xr @@ -323,8 +322,6 @@ def l1b_de( # calculate and set the pointing bin based on the spin phase # pointing bin is 3600 x 40 bins l1b_de = set_pointing_bin(l1b_de) - # set the badtimes - l1b_de = set_bad_times(l1b_de, anc_dependencies) return l1b_de @@ -1125,108 +1122,6 @@ def identify_species(l1b_de: xr.Dataset) -> xr.Dataset: return l1b_de -def set_bad_times(l1b_de: xr.Dataset, anc_dependencies: list) -> xr.Dataset: - """ - Set the bad times for each direct event. - - Parameters - ---------- - l1b_de : xarray.Dataset - The L1B DE dataset. - anc_dependencies : list - List of ancillary file paths. - - Returns - ------- - l1b_de : xarray.Dataset - The L1B DE dataset with the bad times added. - """ - badtimes_df = lo_ancillary.read_ancillary_file( - next(str(s) for s in anc_dependencies if "bad-times" in str(s)) - ) - - esa_steps = l1b_de["esa_step"].values - epochs = l1b_de["epoch"].values - spin_bins = l1b_de["spin_bin"].values - - badtimes = set_bad_or_goodtimes(badtimes_df, epochs, esa_steps, spin_bins) - - # 1 = badtime, 0 = not badtime - l1b_de["badtimes"] = xr.DataArray( - badtimes, - dims=["epoch"], - # TODO: Add to yaml - # attrs=attr_mgr.get_variable_attributes("bad_times"), - ) - - return l1b_de - - -def set_bad_or_goodtimes( - times_df: pd.DataFrame, - epochs: np.ndarray, - esa_steps: np.ndarray, - spin_bins: np.ndarray, -) -> np.ndarray: - """ - Find the good/bad time flags for each epoch based on the provided times DataFrame. - - Parameters - ---------- - times_df : pd.DataFrame - Good or Bad times dataframe containing time ranges and corresponding flags. - epochs : np.ndarray - Array of epochs in TTJ2000ns format. - esa_steps : np.ndarray - Array of ESA steps corresponding to each epoch. - spin_bins : np.ndarray - Array of spin bins corresponding to each epoch. - - Returns - ------- - time_flags : np.ndarray - Array of time good or bad time flags for each epoch. - """ - if "BadTime_start" in times_df.columns and "BadTime_end" in times_df.columns: - times_start = met_to_ttj2000ns(times_df["BadTime_start"]) - times_end = met_to_ttj2000ns(times_df["BadTime_end"]) - elif "GoodTime_start" in times_df.columns and "GoodTime_end" in times_df.columns: - times_start = met_to_ttj2000ns(times_df["GoodTime_start"]) - times_end = met_to_ttj2000ns(times_df["GoodTime_end"]) - else: - raise ValueError("DataFrame must contain either BadTime or GoodTime columns.") - - # Create masks for time and bin ranges using broadcasting - # the bin_start and bin_end are 6 degree bins and need to be converted to - # 0.1 degree bins to align with the spin_bins, so multiply by 60 - time_mask = (epochs[:, None] >= times_start) & (epochs[:, None] <= times_end) - # The ancillary file binning uses 0-59 for the 6 degree bins, so add 1 to bin_end - # so the upper bound is inclusive of the full bin range. - bin_mask = (spin_bins[:, None] >= times_df["bin_start"].values * 60) & ( - spin_bins[:, None] < (times_df["bin_end"].values + 1) * 60 - ) - - # Combined mask for epochs that fall within the time and bin ranges - combined_mask = time_mask & bin_mask - - # TODO: Handle the case where no matching rows are found, because - # otherwise, the bacgkround rates will be set to 0 for those epochs, - # which is not correct. - - # Get the time flags for each epoch's esa_step from matching rows - time_flags: np.ndarray = np.zeros(len(epochs), dtype=int) - for epoch_idx in range(len(epochs)): - matching_rows = np.where(combined_mask[epoch_idx])[0] - if len(matching_rows) > 0: - # Use the first matching row - row_idx = matching_rows[0] - esa_step = esa_steps[epoch_idx] - if f"E-Step{esa_step}" in times_df.columns: - time_flags[epoch_idx] = times_df[f"E-Step{esa_step}"].iloc[row_idx] - - return time_flags - - def set_pointing_direction(l1b_de: xr.Dataset) -> xr.Dataset: """ Set the pointing direction for each direct event. diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 70ca1ed379..544b8b373a 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -12,8 +12,6 @@ from imap_processing.ena_maps.utils.corrections import ( add_spacecraft_position_and_velocity_to_pset, ) -from imap_processing.lo import lo_ancillary -from imap_processing.lo.l1b.lo_l1b import set_bad_or_goodtimes from imap_processing.spice.geometry import ( SpiceFrame, frame_transform_az_el, @@ -90,7 +88,9 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: if "imap_lo_l1b_de" in sci_dependencies: logical_source = "imap_lo_l1c_pset" l1b_de = sci_dependencies["imap_lo_l1b_de"] - l1b_goodtimes_only = filter_goodtimes(l1b_de, anc_dependencies) + l1b_goodtimes_only = filter_goodtimes( + l1b_de, sci_dependencies["imap_lo_l1b_goodtimes"] + ) # Handle case where no good times are found after filtering, # which would lead to an empty dataset @@ -174,15 +174,12 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: pset["h_counts"] = create_pset_counts(l1b_goodtimes_only, FilterType.HYDROGEN) pset["o_counts"] = create_pset_counts(l1b_goodtimes_only, FilterType.OXYGEN) - # Read good-times for exposure time calculation - goodtimes_df = lo_ancillary.read_ancillary_file( - next(str(s) for s in anc_dependencies if "good-times" in str(s)) - ) - # Set the exposure time using statistical off-pointing sampling # with good-times filtering applied pset["exposure_time"] = calculate_exposure_times( - pointing_start_met, pointing_end_met, goodtimes_df + pointing_start_met, + pointing_end_met, + sci_dependencies["imap_lo_l1b_goodtimes"], ) # Set backgrounds @@ -191,10 +188,8 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: pset["h_background_rates_stat_uncert"], pset["h_background_rates_sys_err"], ) = set_background_rates( - pset["pointing_start_met"].item(), - pset["pointing_end_met"].item(), FilterType.HYDROGEN, - anc_dependencies, + sci_dependencies, attr_mgr, ) @@ -203,10 +198,8 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: pset["o_background_rates_stat_uncert"], pset["o_background_rates_sys_err"], ) = set_background_rates( - pset["pointing_start_met"].item(), - pset["pointing_end_met"].item(), FilterType.OXYGEN, - anc_dependencies, + sci_dependencies, attr_mgr, ) @@ -240,43 +233,38 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: return [pset] -def filter_goodtimes(l1b_de: xr.Dataset, anc_dependencies: list) -> xr.Dataset: +def filter_goodtimes(l1b_de: xr.Dataset, goodtimes_ds: xr.Dataset) -> xr.Dataset: """ Filter the L1B Direct Event dataset to only include good times. - The good times are read from the sweep table ancillary file. + The good times are read from the L1B goodtimes dataset produced by + l1b_bgrates_and_goodtimes. Parameters ---------- l1b_de : xarray.Dataset L1B Direct Event dataset. - - anc_dependencies : list - Ancillary files needed for L1C data product creation. + goodtimes_ds : xarray.Dataset + L1B goodtimes dataset containing gt_start_met and gt_end_met variables + that define good time windows in MET seconds. Returns ------- - l1b_de : xarray.Dataset - Filtered L1B Direct Event dataset. + xarray.Dataset + Filtered L1B Direct Event dataset containing only events within good + time windows. """ - # The goodtimes are one of several ancillary files needed for L1C processing - goodtimes_table_df = lo_ancillary.read_ancillary_file( - next(str(s) for s in anc_dependencies if "good-times" in str(s)) - ) - - esa_steps = l1b_de["esa_step"].values epochs = l1b_de["epoch"].values - spin_bins = l1b_de["spin_bin"].values + gt_starts = met_to_ttj2000ns(goodtimes_ds["gt_start_met"].values) + gt_ends = met_to_ttj2000ns(goodtimes_ds["gt_end_met"].values) - # Get array of bools for each epoch 1 = good time, 0 not good time - goodtimes_mask = set_bad_or_goodtimes( - goodtimes_table_df, epochs, esa_steps, spin_bins + # Keep events that fall within any goodtime window + in_goodtime = np.any( + (epochs[:, None] >= gt_starts) & (epochs[:, None] <= gt_ends), + axis=1, ) - # Filter the dataset using the mask - filtered_epochs = l1b_de.sel(epoch=goodtimes_mask.astype(bool)) - - return filtered_epochs + return l1b_de.sel(epoch=in_goodtime) def get_triple_coincidences(de: xr.Dataset) -> xr.Dataset: @@ -671,7 +659,7 @@ def calculate_bin_weights(off_angles: np.ndarray) -> np.ndarray: def create_goodtimes_fraction( - goodtimes_df: pd.DataFrame, + goodtimes_ds: xr.Dataset, pointing_start_met: float, pointing_end_met: float, ) -> np.ndarray: @@ -685,9 +673,9 @@ def create_goodtimes_fraction( Parameters ---------- - goodtimes_df : pandas.DataFrame - DataFrame containing the good-times ancillary data with columns: - GoodTime_start, GoodTime_end, bin_start, bin_end, E-Step1 through E-Step7. + goodtimes_ds : xarray.Dataset + Dataset containing the good-times data with variables: + gt_start_met, gt_end_met, bin_start, bin_end, E-Step1 - E-Step7. pointing_start_met : float The start MET time of the pointing. pointing_end_met : float @@ -712,12 +700,12 @@ def create_goodtimes_fraction( return goodtimes_fraction # Filter good-times to only those overlapping with the pointing period - pointing_goodtimes = goodtimes_df[ - (goodtimes_df["GoodTime_start"] < pointing_end_met) - & (goodtimes_df["GoodTime_end"] > pointing_start_met) - ] + pointing_goodtimes_mask = ( + (goodtimes_ds["gt_start_met"] < pointing_end_met) + & (goodtimes_ds["gt_end_met"] > pointing_start_met) + ).values - if len(pointing_goodtimes) == 0: + if not pointing_goodtimes_mask.any(): logging.warning( f"No good-times found for pointing period " f"[{pointing_start_met}, {pointing_end_met}]. " @@ -725,11 +713,15 @@ def create_goodtimes_fraction( ) return goodtimes_fraction + pointing_goodtimes = goodtimes_ds.isel(epoch=pointing_goodtimes_mask) + # Process each good-time row and accumulate fractional coverage - for _, row in pointing_goodtimes.iterrows(): + for i in range(len(pointing_goodtimes["epoch"])): + row = pointing_goodtimes.isel(epoch=i) + # Calculate the overlap between this good-time period and the pointing - goodtime_start = max(row["GoodTime_start"], pointing_start_met) - goodtime_end = min(row["GoodTime_end"], pointing_end_met) + goodtime_start = max(float(row["gt_start_met"]), pointing_start_met) + goodtime_end = min(float(row["gt_end_met"]), pointing_end_met) overlap_duration = goodtime_end - goodtime_start if overlap_duration <= 0: @@ -738,25 +730,15 @@ def create_goodtimes_fraction( # Calculate fraction of pointing duration covered by this good-time time_fraction = overlap_duration / total_pointing_duration - # Convert bin_start/bin_end from 6-degree to 0.1-degree resolution - # bin_start and bin_end are in units of 6-degree bins (0-59), inclusive - # We need to convert to 0.1-degree bins (0-3599) - bin_start_6deg = int(row["bin_start"]) - bin_end_6deg = int(row["bin_end"]) - + # spin_bin_start and spin_bin_end cover the entire 360 degree range # Convert to 0.1-degree resolution (multiply by 60) - # bin_end is inclusive, so add 1 after scaling for Python slice indexing - spin_bin_start = bin_start_6deg * 60 - spin_bin_end = (bin_end_6deg + 1) * 60 # +1 because bin_end is inclusive + spin_bin_start = 0 + spin_bin_end = 3600 # For each ESA step, accumulate the fractional coverage for esa_idx in range(N_ESA_ENERGY_STEPS): - esa_step_col = f"E-Step{esa_idx + 1}" - if row[esa_step_col] == 1: - # Add this time fraction to the affected bins - goodtimes_fraction[esa_idx, spin_bin_start:spin_bin_end] += ( - time_fraction - ) + # Add this time fraction to the affected bins + goodtimes_fraction[esa_idx, spin_bin_start:spin_bin_end] += time_fraction # Clip to [0, 1] in case of overlapping good-time periods goodtimes_fraction = np.clip(goodtimes_fraction, 0.0, 1.0) @@ -775,7 +757,7 @@ def create_goodtimes_fraction( def calculate_exposure_times( pointing_start_met: float, pointing_end_met: float, - goodtimes_df: pd.DataFrame | None = None, + goodtimes_ds: xr.Dataset | None = None, n_representative_spins: int = DEFAULT_N_REPRESENTATIVE_SPINS, ) -> xr.DataArray: """ @@ -793,8 +775,8 @@ def calculate_exposure_times( The start MET time of the pointing. pointing_end_met : float The end MET time of the pointing. - goodtimes_df : pandas.DataFrame, optional - DataFrame containing the good-times ancillary data. If provided, + goodtimes_ds : xarray.Dataset, optional + Dataset containing the good-times data. If provided, exposure times will be zeroed for invalid spin_angle bins and ESA steps. n_representative_spins : int, optional Number of representative spins to sample. Default is 5. @@ -876,9 +858,9 @@ def calculate_exposure_times( ).copy() # Apply good-times fraction if provided - if goodtimes_df is not None: + if goodtimes_ds is not None: goodtimes_fraction = create_goodtimes_fraction( - goodtimes_df, pointing_start_met, pointing_end_met + goodtimes_ds, pointing_start_met, pointing_end_met ) # Expand fraction to include off_angle dimension # (fraction is same for all off_angles) @@ -1046,27 +1028,27 @@ def create_datasets( def set_background_rates( - pointing_start_met: float, - pointing_end_met: float, species: FilterType, - anc_dependencies: list, + sci_dependencies: dict, attr_mgr: ImapCdfAttributes, ) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray]: """ Set the background rates for the specified species. - The background rates are set to a constant value of 0.01 counts/s for all bins. + Background rates and statistical uncertainties are read from the + ``imap_lo_l1b_bgrates`` dataset in ``sci_dependencies``. Each species + provides a 1-D array of shape ``(N_ESA_ENERGY_STEPS,)`` that is broadcast + uniformly across all spin-angle and off-angle bins. If the bgrates dataset + is absent, all arrays default to zero. Parameters ---------- - pointing_start_met : float - The start MET time of the pointing. - pointing_end_met : float - The end MET time of the pointing. species : FilterType The species to set the background rates for. Can be "h" or "o". - anc_dependencies : list - Ancillary files needed for L1C data product creation. + sci_dependencies : dict + Science dependency datasets. Expected to contain the key + ``"imap_lo_l1b_bgrates"`` with variables + ``"{species}_background_rates"`` and ``"{species}_background_variance"``. attr_mgr : ImapCdfAttributes Attribute manager used to get the L1C attributes. @@ -1094,43 +1076,30 @@ def set_background_rates( dtype=np.float16, ) - # read in the background rates from ancillary file - if species == FilterType.HYDROGEN: - background_df = lo_ancillary.read_ancillary_file( - next(str(s) for s in anc_dependencies if "hydrogen-background" in str(s)) - ) - else: - background_df = lo_ancillary.read_ancillary_file( - next(str(s) for s in anc_dependencies if "oxygen-background" in str(s)) - ) - - # find to the rows for the current pointing - # TODO: This assumes that the backgrounds will never change mid-pointing. - # Is that a safe assumption? - pointing_bg_df = background_df[ - (background_df["GoodTime_start"] < pointing_end_met) - & (background_df["GoodTime_end"] > pointing_start_met) - ] + bgrates_ds = sci_dependencies.get("imap_lo_l1b_bgrates") + if bgrates_ds is not None: + species_key = species.value + rate_field = f"{species_key}_background_rates" + variance_field = f"{species_key}_background_variance" + + if rate_field in bgrates_ds: + rates_per_esa = bgrates_ds[ + rate_field + ].values # shape: (N_ESA_ENERGY_STEPS,) + bg_rates = np.broadcast_to( + rates_per_esa[:, np.newaxis, np.newaxis], + (N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS), + ).astype(np.float16) + + if variance_field in bgrates_ds: + var_per_esa = bgrates_ds[ + variance_field + ].values # shape: (N_ESA_ENERGY_STEPS,) + bg_stat_uncert = np.broadcast_to( + var_per_esa[:, np.newaxis, np.newaxis], + (N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS), + ).astype(np.float16) - # convert the bin start and end resolution from 6 degrees to 0.1 degrees - pointing_bg_df["bin_start"] = pointing_bg_df["bin_start"] * 60 - # The bin_end index is inclusive, so add one and convert to 0.1 - # degree resolution - pointing_bg_df["bin_end"] = (pointing_bg_df["bin_end"] + 1) * 60 - # for each row in the bg ancillary file for this pointing - for _, row in pointing_bg_df.iterrows(): - bin_start = int(row["bin_start"]) - bin_end = int(row["bin_end"]) - # for each energy step, set the background rate and uncertainty - for esa_step in range(0, 7): - value = row[f"E-Step{esa_step + 1}"] - if row["rate/sigma"] == "rate": - bg_rates[esa_step, bin_start:bin_end, :] = value - elif row["rate/sigma"] == "sigma": - bg_sys_err[esa_step, bin_start:bin_end, :] = value - else: - raise ValueError("Unknown background type in ancillary file.") - # set the background rates, uncertainties, and systematic errors bg_rates_data = xr.DataArray( data=bg_rates[np.newaxis, :, :, :], dims=["epoch", "esa_energy_step", "spin_angle", "off_angle"], diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index dac8da5519..d487b3b6a7 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -34,8 +34,6 @@ lo_l1b, resweep_histogram_data, set_avg_spin_durations_per_event, - set_bad_or_goodtimes, - set_bad_times, set_coincidence_type, set_each_event_epoch, set_esa_mode, @@ -46,7 +44,6 @@ set_spin_cycle_from_spin_data, split_backgrounds_and_goodtimes_dataset, ) -from imap_processing.lo.lo_ancillary import read_ancillary_file from imap_processing.spice.spin import get_spin_data from imap_processing.spice.time import ( et_to_met, @@ -764,45 +761,6 @@ def test_identify_species(attr_mgr_l1b): np.testing.assert_array_equal(l1b_de["species"], expected_species) -def test_set_bad_times(anc_dependencies): - # Arrange - l1b_de = xr.Dataset( - { - "esa_step": ("epoch", [1, 1, 3, 1]), - "spin_bin": ("epoch", [1900, 2000, 3000, 2]), - }, - coords={ - "epoch": met_to_ttj2000ns([473385599, 473385600, 473385601, 473385602]), - }, - ) - - expected_bad_times = np.array([0, 1, 0, 0]) - - # Act - l1b_de = set_bad_times(l1b_de, anc_dependencies) - - # Assert - np.testing.assert_array_equal(l1b_de["badtimes"], expected_bad_times) - - -def test_set_bad_or_goodtimes(anc_dependencies): - # Arrange - # badtimes ancillary - df = read_ancillary_file(anc_dependencies[1]) - - epoch = met_to_ttj2000ns([473385599, 473385600, 473385601, 473385602]) - esa_step = np.array([1, 1, 3, 1]) - spin_bin = np.array([1900, 2000, 3000, 2]) - - expected_bad_times = np.array([0, 1, 0, 0]) - - # Act - badtimes = set_bad_or_goodtimes(df, epoch, esa_step, spin_bin) - - # Assert - np.testing.assert_array_equal(badtimes, expected_bad_times) - - @patch( "imap_processing.lo.l1b.lo_l1b.lo_instrument_pointing", return_value=np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]), diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 079dad205c..51d5be0e24 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -1,7 +1,6 @@ from unittest.mock import patch import numpy as np -import pandas as pd import pytest import xarray as xr @@ -179,41 +178,25 @@ def doubles_counts(counts): @pytest.fixture -def expected_bg(): - expected_rates = np.array( - [ - [ - np.full((3600, 40), 0.0098), - np.full((3600, 40), 0.0089), - np.full((3600, 40), 0.0118), - np.full((3600, 40), 0.0113), - np.full((3600, 40), 0.0056), - np.full((3600, 40), 0.0008), - np.full((3600, 40), 0.0), - ] - ], - dtype=np.float16, +def l1b_bgrates_ds(): + h_rates = np.array([0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07], dtype=np.float32) + h_var = np.array( + [0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007], dtype=np.float32 ) - - expected_err = np.array( - [ - [ - np.full((3600, 40), 0.0025), - np.full((3600, 40), 0.002), - np.full((3600, 40), 0.0015), - np.full((3600, 40), 0.0015), - np.full((3600, 40), 0.001), - np.full((3600, 40), 0.0008), - np.full((3600, 40), 0.0), - ] - ], - dtype=np.float16, + o_rates = np.array( + [0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007], dtype=np.float32 + ) + o_var = np.array( + [0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007], dtype=np.float32 + ) + return xr.Dataset( + { + "h_background_rates": ("esa_step", h_rates), + "h_background_variance": ("esa_step", h_var), + "o_background_rates": ("esa_step", o_rates), + "o_background_variance": ("esa_step", o_var), + } ) - - expected_uncert = np.zeros((1, 7, 3600, 40), dtype=np.float16) - - expected_bg = (expected_rates, expected_uncert, expected_err) - return expected_bg @patch("imap_processing.lo.l1c.lo_l1c.calculate_exposure_times") @@ -234,7 +217,16 @@ def test_lo_l1c( repoint_met, ): # Arrange - data = {"imap_lo_l1b_de": l1b_de_spin} + data = { + "imap_lo_l1b_de": l1b_de_spin, + "imap_lo_l1b_goodtimes": xr.Dataset( + { + "gt_start_met": ("epoch", [511000000.0]), + "gt_end_met": ("epoch", [511100000.0]), + }, + coords={"epoch": met_to_ttj2000ns([511000000.0])}, + ), + } use_fake_spin_data_for_time(511000000) use_fake_repoint_data_for_time(np.arange(511000000, 511000000 + 86400 * 5, 86400)) mock_set_background_rates.return_value = (None, None, None) @@ -271,35 +263,37 @@ def mock_add_sc_pos_vel(pset): assert "sc_velocity" in output_dataset -def test_filter_goodtimes(l1b_de, anc_dependencies): +def test_filter_goodtimes(): # Arrange + event_mets = [473389199, 473389200, 473389201, 473389202, 473389203, 473407619] l1b_de_all = xr.Dataset( { "esa_step": ("epoch", [1, 2, 1, 4, 5, 2]), "spin_bin": ("epoch", [1900, 2000, 3000, 3000, 3000, 3000]), }, - coords={ - "epoch": met_to_ttj2000ns( - [ - 473389199, - 473389200, - 473389201, - 473389202, - 473389203, - 473407619, - ] - ) + coords={"epoch": met_to_ttj2000ns(event_mets)}, + ) + + # Two goodtime windows: [473389201, 473389203] and [473407619, 473407620] + gt_starts = np.array([473389201.0, 473407619.0]) + goodtimes_ds = xr.Dataset( + { + "gt_start_met": ("epoch", gt_starts), + "gt_end_met": ("epoch", [473389203.0, 473407620.0]), }, + coords={"epoch": met_to_ttj2000ns(gt_starts)}, ) - expected_goodtimes_mask = [False, False, True, False, True, False] - l1b_goodtimes_onl_expected = l1b_de_all.isel(epoch=expected_goodtimes_mask) + # Events at MET 473389201-473389203 and 473407619 fall inside goodtime windows; + # 473389199 and 473389200 are before the first window. + expected_mask = [False, False, True, True, True, True] + expected = l1b_de_all.isel(epoch=expected_mask) # Act - l1b_goodtimes_only = filter_goodtimes(l1b_de_all, anc_dependencies) + result = filter_goodtimes(l1b_de_all, goodtimes_ds) # Assert - xr.testing.assert_equal(l1b_goodtimes_only, l1b_goodtimes_onl_expected) + xr.testing.assert_equal(result, expected) def test_lo_l1c_no_goodtimes( @@ -310,7 +304,17 @@ def test_lo_l1c_no_goodtimes( repoint_met, ): # Arrange - data = {"imap_lo_l1b_de": l1b_de_spin} + # Goodtime window [0, 1] does not cover any event + data = { + "imap_lo_l1b_de": l1b_de_spin, + "imap_lo_l1b_goodtimes": xr.Dataset( + { + "gt_start_met": ("epoch", [0.0]), + "gt_end_met": ("epoch", [1.0]), + }, + coords={"epoch": met_to_ttj2000ns([0.0])}, + ), + } use_fake_spin_data_for_time(511000000) use_fake_repoint_data_for_time(np.arange(511000000, 511000000 + 86400 * 5, 86400)) expected_logical_source = "imap_lo_l1c_pset" @@ -573,23 +577,24 @@ def test_calculate_bin_weights_distributed(): def test_create_goodtimes_fraction(): """Test good-times fractional coverage calculation from ancillary data.""" - # Arrange - create a simple goodtimes DataFrame + # Arrange - create a simple goodtimes Dataset # Good-times cover the full pointing duration for all spin bins # bin_start and bin_end are inclusive, 0-indexed (0-59 for 6-degree bins) - goodtimes_df = pd.DataFrame( + goodtimes_ds = xr.Dataset( { - "GoodTime_start": [500000000.0, 500000000.0], - "GoodTime_end": [500001000.0, 500001000.0], - "bin_start": [0, 30], # 6-degree bins: 0-29 and 30-59 - "bin_end": [29, 59], # inclusive: first half bins 0-29, second half 30-59 - "E-Step1": [1, 1], - "E-Step2": [1, 0], # ESA step 2 only good for first half - "E-Step3": [1, 1], - "E-Step4": [1, 1], - "E-Step5": [1, 1], - "E-Step6": [1, 1], - "E-Step7": [1, 1], - } + "gt_start_met": ("epoch", [500000000.0, 500000000.0]), + "gt_end_met": ("epoch", [500001000.0, 500001000.0]), + "bin_start": ("epoch", [0, 30]), # 6-degree bins: 0-29 and 30-59 + "bin_end": ("epoch", [29, 59]), # inclusive: first half 0-29, second 30-59 + "E-Step1": ("epoch", [1, 1]), + "E-Step2": ("epoch", [1, 0]), # ESA step 2 only good for first half + "E-Step3": ("epoch", [1, 1]), + "E-Step4": ("epoch", [1, 1]), + "E-Step5": ("epoch", [1, 1]), + "E-Step6": ("epoch", [1, 1]), + "E-Step7": ("epoch", [1, 1]), + }, + coords={"epoch": [0, 1]}, ) pointing_start_met = 500000000.0 @@ -597,7 +602,7 @@ def test_create_goodtimes_fraction(): # Act fraction = create_goodtimes_fraction( - goodtimes_df, pointing_start_met, pointing_end_met + goodtimes_ds, pointing_start_met, pointing_end_met ) # Assert @@ -618,20 +623,21 @@ def test_create_goodtimes_fraction(): def test_create_goodtimes_fraction_partial_coverage(): """Test good-times with partial time coverage of pointing period.""" # Arrange - good-times cover only half of the pointing duration - goodtimes_df = pd.DataFrame( + goodtimes_ds = xr.Dataset( { - "GoodTime_start": [500000000.0], - "GoodTime_end": [500000500.0], # Only first 500s of 1000s pointing - "bin_start": [0], - "bin_end": [59], # All spin bins (0-59 inclusive) - "E-Step1": [1], - "E-Step2": [1], - "E-Step3": [1], - "E-Step4": [1], - "E-Step5": [1], - "E-Step6": [1], - "E-Step7": [1], - } + "gt_start_met": ("epoch", [500000000.0]), + "gt_end_met": ("epoch", [500000500.0]), # Only first 500s of 1000s pointing + "bin_start": ("epoch", [0]), + "bin_end": ("epoch", [59]), # All spin bins (0-59 inclusive) + "E-Step1": ("epoch", [1]), + "E-Step2": ("epoch", [1]), + "E-Step3": ("epoch", [1]), + "E-Step4": ("epoch", [1]), + "E-Step5": ("epoch", [1]), + "E-Step6": ("epoch", [1]), + "E-Step7": ("epoch", [1]), + }, + coords={"epoch": [0]}, ) pointing_start_met = 500000000.0 @@ -639,7 +645,7 @@ def test_create_goodtimes_fraction_partial_coverage(): # Act fraction = create_goodtimes_fraction( - goodtimes_df, pointing_start_met, pointing_end_met + goodtimes_ds, pointing_start_met, pointing_end_met ) # Assert - all bins should have 50% coverage @@ -649,20 +655,21 @@ def test_create_goodtimes_fraction_partial_coverage(): def test_create_goodtimes_fraction_no_overlap(): """Test good-times fraction when no good-times overlap with pointing.""" # Arrange - goodtimes outside pointing period - goodtimes_df = pd.DataFrame( + goodtimes_ds = xr.Dataset( { - "GoodTime_start": [400000000.0], - "GoodTime_end": [400001000.0], - "bin_start": [0], - "bin_end": [59], # All spin bins (0-59 inclusive) - "E-Step1": [1], - "E-Step2": [1], - "E-Step3": [1], - "E-Step4": [1], - "E-Step5": [1], - "E-Step6": [1], - "E-Step7": [1], - } + "gt_start_met": ("epoch", [400000000.0]), + "gt_end_met": ("epoch", [400001000.0]), + "bin_start": ("epoch", [0]), + "bin_end": ("epoch", [59]), # All spin bins (0-59 inclusive) + "E-Step1": ("epoch", [1]), + "E-Step2": ("epoch", [1]), + "E-Step3": ("epoch", [1]), + "E-Step4": ("epoch", [1]), + "E-Step5": ("epoch", [1]), + "E-Step6": ("epoch", [1]), + "E-Step7": ("epoch", [1]), + }, + coords={"epoch": [0]}, ) pointing_start_met = 500000000.0 @@ -670,7 +677,7 @@ def test_create_goodtimes_fraction_no_overlap(): # Act fraction = create_goodtimes_fraction( - goodtimes_df, pointing_start_met, pointing_end_met + goodtimes_ds, pointing_start_met, pointing_end_met ) # Assert - all zeros since no overlap @@ -678,46 +685,57 @@ def test_create_goodtimes_fraction_no_overlap(): @pytest.mark.parametrize("species", [FilterType.HYDROGEN, FilterType.OXYGEN]) -def test_set_background_rates( - l1b_de_spin, anc_dependencies, attr_mgr, species, expected_bg -): +def test_set_background_rates(l1b_bgrates_ds, attr_mgr, species): # Arrange - pointing_start_met = 473389200.0 - pointing_end_met = 473472000.0 + sci_deps = {"imap_lo_l1b_bgrates": l1b_bgrates_ds} + species_key = species.value + expected_rates = l1b_bgrates_ds[f"{species_key}_background_rates"].values + expected_var = l1b_bgrates_ds[f"{species_key}_background_variance"].values # Act - rates, uncert, err = set_background_rates( - pointing_start_met, pointing_end_met, species, anc_dependencies, attr_mgr - ) + rates, uncert, err = set_background_rates(species, sci_deps, attr_mgr) + + # Assert shape + assert rates.shape == (1, N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS) + assert uncert.shape == (1, N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS) + assert err.shape == (1, N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS) + + # Rates and uncertainties must be uniform across spatial bins for each ESA step + for i in range(N_ESA_ENERGY_STEPS): + np.testing.assert_array_equal( + rates.values[0, i, :, :], + np.full( + (N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS), + expected_rates[i], + dtype=np.float16, + ), + ) + np.testing.assert_array_equal( + uncert.values[0, i, :, :], + np.full( + (N_SPIN_ANGLE_BINS, N_OFF_ANGLE_BINS), expected_var[i], dtype=np.float16 + ), + ) - # Assert - np.testing.assert_array_equal( - rates.values, - expected_bg[0], - ) - np.testing.assert_array_equal( - uncert.values, - expected_bg[1], - ) - np.testing.assert_array_equal( - err.values, - expected_bg[2], - ) + # Systematic error is always zero + np.testing.assert_array_equal(err.values, 0) -def test_set_background_rates_species_error(anc_dependencies, attr_mgr): - # Arrange - pointing_start_met = 473389100.0 - pointing_end_met = 473472100.0 - species = FilterType.DOUBLES +def test_set_background_rates_no_bgrates(attr_mgr): + """Returns zeros when imap_lo_l1b_bgrates is absent from sci_dependencies.""" + rates, uncert, err = set_background_rates(FilterType.HYDROGEN, {}, attr_mgr) + + np.testing.assert_array_equal(rates.values, 0) + np.testing.assert_array_equal(uncert.values, 0) + np.testing.assert_array_equal(err.values, 0) + +def test_set_background_rates_species_error(attr_mgr): # Act with pytest.raises( ValueError, match="Species must be 'h' or 'o', but got doubles." ): - rates, uncert, err = set_background_rates( - pointing_start_met, pointing_end_met, species, anc_dependencies, attr_mgr - ) + set_background_rates(FilterType.DOUBLES, {}, attr_mgr) def test_set_pointing_directions(attr_mgr): From 41a89882fe07347f264632867b50fb256bf1cdd9 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 13 May 2026 17:45:16 -0400 Subject: [PATCH 2/5] fixed failing test and simplified others --- imap_processing/tests/lo/test_lo_l1c.py | 40 ++----------------------- 1 file changed, 3 insertions(+), 37 deletions(-) diff --git a/imap_processing/tests/lo/test_lo_l1c.py b/imap_processing/tests/lo/test_lo_l1c.py index 51d5be0e24..dd59a4f188 100644 --- a/imap_processing/tests/lo/test_lo_l1c.py +++ b/imap_processing/tests/lo/test_lo_l1c.py @@ -584,15 +584,6 @@ def test_create_goodtimes_fraction(): { "gt_start_met": ("epoch", [500000000.0, 500000000.0]), "gt_end_met": ("epoch", [500001000.0, 500001000.0]), - "bin_start": ("epoch", [0, 30]), # 6-degree bins: 0-29 and 30-59 - "bin_end": ("epoch", [29, 59]), # inclusive: first half 0-29, second 30-59 - "E-Step1": ("epoch", [1, 1]), - "E-Step2": ("epoch", [1, 0]), # ESA step 2 only good for first half - "E-Step3": ("epoch", [1, 1]), - "E-Step4": ("epoch", [1, 1]), - "E-Step5": ("epoch", [1, 1]), - "E-Step6": ("epoch", [1, 1]), - "E-Step7": ("epoch", [1, 1]), }, coords={"epoch": [0, 1]}, ) @@ -608,16 +599,9 @@ def test_create_goodtimes_fraction(): # Assert assert fraction.shape == (N_ESA_ENERGY_STEPS, N_SPIN_ANGLE_BINS) - # ESA step 1 (index 0) should have 100% coverage (fraction = 1.0) - np.testing.assert_allclose(fraction[0, :], 1.0) - - # ESA step 2 (index 1) should only have 100% for first half, 0% for second half - np.testing.assert_allclose(fraction[1, :1800], 1.0) - np.testing.assert_allclose(fraction[1, 1800:], 0.0) - - # ESA steps 3-7 (indices 2-6) should have 100% coverage - for i in range(2, 7): - np.testing.assert_allclose(fraction[i, :], 1.0) + # The current implementation does not filter by bin range or E-Step flags, + # so coverage is 1.0 everywhere + np.testing.assert_allclose(fraction, 1.0) def test_create_goodtimes_fraction_partial_coverage(): @@ -627,15 +611,6 @@ def test_create_goodtimes_fraction_partial_coverage(): { "gt_start_met": ("epoch", [500000000.0]), "gt_end_met": ("epoch", [500000500.0]), # Only first 500s of 1000s pointing - "bin_start": ("epoch", [0]), - "bin_end": ("epoch", [59]), # All spin bins (0-59 inclusive) - "E-Step1": ("epoch", [1]), - "E-Step2": ("epoch", [1]), - "E-Step3": ("epoch", [1]), - "E-Step4": ("epoch", [1]), - "E-Step5": ("epoch", [1]), - "E-Step6": ("epoch", [1]), - "E-Step7": ("epoch", [1]), }, coords={"epoch": [0]}, ) @@ -659,15 +634,6 @@ def test_create_goodtimes_fraction_no_overlap(): { "gt_start_met": ("epoch", [400000000.0]), "gt_end_met": ("epoch", [400001000.0]), - "bin_start": ("epoch", [0]), - "bin_end": ("epoch", [59]), # All spin bins (0-59 inclusive) - "E-Step1": ("epoch", [1]), - "E-Step2": ("epoch", [1]), - "E-Step3": ("epoch", [1]), - "E-Step4": ("epoch", [1]), - "E-Step5": ("epoch", [1]), - "E-Step6": ("epoch", [1]), - "E-Step7": ("epoch", [1]), }, coords={"epoch": [0]}, ) From eac463173d8473b89f920e5b1c2cfcfe4429cda8 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 14 May 2026 10:03:25 -0400 Subject: [PATCH 3/5] changes suggested by CoPilot --- imap_processing/lo/l1c/lo_l1c.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 544b8b373a..50a6280a78 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -675,7 +675,7 @@ def create_goodtimes_fraction( ---------- goodtimes_ds : xarray.Dataset Dataset containing the good-times data with variables: - gt_start_met, gt_end_met, bin_start, bin_end, E-Step1 - E-Step7. + gt_start_met, gt_end_met. pointing_start_met : float The start MET time of the pointing. pointing_end_met : float @@ -730,15 +730,11 @@ def create_goodtimes_fraction( # Calculate fraction of pointing duration covered by this good-time time_fraction = overlap_duration / total_pointing_duration - # spin_bin_start and spin_bin_end cover the entire 360 degree range - # Convert to 0.1-degree resolution (multiply by 60) - spin_bin_start = 0 - spin_bin_end = 3600 - - # For each ESA step, accumulate the fractional coverage - for esa_idx in range(N_ESA_ENERGY_STEPS): - # Add this time fraction to the affected bins - goodtimes_fraction[esa_idx, spin_bin_start:spin_bin_end] += time_fraction + # For each ESA step, accumulate the fractional coverage. + # Add this time fraction to the affected bins. + # Note that all ESA Levels and all N_SPIN_ANGLE_BINS currently get the + # same increment, pending algorithmic changes in the future. + goodtimes_fraction += time_fraction # Clip to [0, 1] in case of overlapping good-time periods goodtimes_fraction = np.clip(goodtimes_fraction, 0.0, 1.0) From 5945bc04fdbab9b26654669580604c94b597d35c Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 15 May 2026 12:01:18 -0400 Subject: [PATCH 4/5] some change based on PR feedback --- imap_processing/lo/l1c/lo_l1c.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 50a6280a78..9cb98ad915 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -22,7 +22,6 @@ from imap_processing.spice.time import ( met_to_ttj2000ns, ttj2000ns_to_et, - ttj2000ns_to_met, ) N_ESA_ENERGY_STEPS = 7 @@ -103,9 +102,10 @@ def lo_l1c(sci_dependencies: dict, anc_dependencies: list) -> list[xr.Dataset]: pointing_start_met = 0.0 pointing_end_met = 0.0 else: - # Set the pointing start and end times based on the first epoch pointing_start_met, pointing_end_met = get_pointing_times( - float(ttj2000ns_to_met(l1b_goodtimes_only["epoch"][0].item()).item()) + float( + sci_dependencies["imap_lo_l1b_goodtimes"]["gt_start_met"].values[0] + ) ) pset = xr.Dataset( @@ -264,7 +264,7 @@ def filter_goodtimes(l1b_de: xr.Dataset, goodtimes_ds: xr.Dataset) -> xr.Dataset axis=1, ) - return l1b_de.sel(epoch=in_goodtime) + return l1b_de.isel(epoch=in_goodtime) def get_triple_coincidences(de: xr.Dataset) -> xr.Dataset: From 970bc0329c6e49ecb7c62da1c011b51a0dcc9765 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 15 May 2026 16:06:15 -0400 Subject: [PATCH 5/5] Update lo_l1c.py Co-authored-by: Tim Plummer --- imap_processing/lo/l1c/lo_l1c.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imap_processing/lo/l1c/lo_l1c.py b/imap_processing/lo/l1c/lo_l1c.py index 9cb98ad915..2eb43c99da 100644 --- a/imap_processing/lo/l1c/lo_l1c.py +++ b/imap_processing/lo/l1c/lo_l1c.py @@ -260,7 +260,7 @@ def filter_goodtimes(l1b_de: xr.Dataset, goodtimes_ds: xr.Dataset) -> xr.Dataset # Keep events that fall within any goodtime window in_goodtime = np.any( - (epochs[:, None] >= gt_starts) & (epochs[:, None] <= gt_ends), + (epochs[:, np.newaxis] >= gt_starts) & (epochs[:, np.newaxis] <= gt_ends), axis=1, )