diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index 8c457edd37..3ff06f7d3d 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -171,3 +171,135 @@ count_per_bin: LABLAXIS: Sample Count DEPEND_0: epoch DEPEND_1: spin_angle + +pivot: + CATDESC: Pivot angle derived from housekeeping data + FIELDNAM: Pivot Angle + FORMAT: F8.3 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 180.0 + VAR_TYPE: support_data + UNITS: deg + LABLAXIS: Pivot Angle + +pivot_de: + CATDESC: Pivot angle derived from direct event data + FIELDNAM: Pivot Angle (DE) + FORMAT: F8.3 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 180.0 + VAR_TYPE: support_data + UNITS: deg + LABLAXIS: Pivot Angle DE + +gt_start_met: + <<: *default + CATDESC: Good time interval start in Mission Elapsed Time + FIELDNAM: Good Time Start MET + FORMAT: F19.3 + FILLVAL: -1.0000000E+31 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: s + LABLAXIS: GT Start + +gt_end_met: + <<: *default + CATDESC: Good time interval end in Mission Elapsed Time + FIELDNAM: Good Time End MET + FORMAT: F19.3 + FILLVAL: -1.0000000E+31 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: s + LABLAXIS: GT End + +h_background_rates: + CATDESC: Hydrogen background rates per ESA level + FIELDNAM: H Background Rates + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: data + UNITS: counts/s + LABLAXIS: H BG Rate + +h_background_variance: + CATDESC: Hydrogen background rate variance per ESA level + FIELDNAM: H Background Variance + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: data + UNITS: (counts/s)^2 + LABLAXIS: H BG Variance + +h_synthetic_floor: + CATDESC: Hydrogen synthetic background floor + FIELDNAM: H Synthetic Floor + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: counts/s + LABLAXIS: H Synthetic Floor + +h_proxy_floor: + CATDESC: Hydrogen proxy background floor + FIELDNAM: H Proxy Floor + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: counts/s + LABLAXIS: H Proxy Floor + +o_background_rates: + CATDESC: Oxygen background rates per ESA level + FIELDNAM: O Background Rates + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: data + UNITS: counts/s + LABLAXIS: O BG Rate + +o_background_variance: + CATDESC: Oxygen background rate variance per ESA level + FIELDNAM: O Background Variance + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: data + UNITS: (counts/s)^2 + LABLAXIS: O BG Variance + +o_synthetic_floor: + CATDESC: Oxygen synthetic background floor + FIELDNAM: O Synthetic Floor + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: counts/s + LABLAXIS: O Synthetic Floor + +o_proxy_floor: + CATDESC: Oxygen proxy background floor + FIELDNAM: O Proxy Floor + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0000000E+31 + VAR_TYPE: support_data + UNITS: counts/s + LABLAXIS: O Proxy Floor diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py new file mode 100644 index 0000000000..5745799ad9 --- /dev/null +++ b/imap_processing/lo/constants.py @@ -0,0 +1,59 @@ +"""Contains constants variables to support IMAP-Lo processing.""" + +# ------- Constants for L1B processing ------- + +# Ion species tracked. "H" is mandatory; any others for which we have histrates +# may be added here. +ELEMS = ("H", "O") + +# Hours into the day (UTC) for HK data to calculate median for pivot angle estimation. +PIVOT_HK_HOUR_RANGE: tuple[int, int] = (3, 15) + +# Per-day overrides for the anti-RAM background rate. Keyed by (year, day-of-year) +BG_RATE_ANTI_RAM_OVERRIDES: dict[tuple[int, int], float] = { + (2026, 62): 0.0014, + (2026, 64): 0.0, + (2026, 65): 0.0, + (2026, 91): 0.03, +} + +# One histogram accumulation cycle duration [s] +HISTOGRAM_CYCLE_EPOCHS: int = 420 + +N_CYCLE_SUM: int = 1 # Granularity of goodtime boundaries +N_CYCLE_AVE: int = 7 # Cycles to average over when estimating background rates +N_ESA_LEVELS: int = 7 # Total number of ESA levels +RAM_ESA_LEVELS: tuple[int, ...] = ( + 6, + 7, +) # ESA levels for RAM estimation (1-indexed) + +# Histogram angular bins (0-indexed) corresponding to the RAM and anti-RAM look +# directions +RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(0, 20), slice(50, 60)) +ANTI_RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(20, 50),) + +# Nominal background rates [counts/s] for each species +BG_RATES = {"H": 0.0014925, "O": 0.000136635} +# When no exposure is available, scale the nominal rate down as a conservative estimate. +BG_RATE_FALLBACK_SCALE: dict[str, float] = {"H": 1.0, "O": 0.3} +# Minimum non-zero background rate floor = nominal / divisor +BG_RATE_FLOOR_DIVISOR: dict[str, float] = {"H": 50.0, "O": 150.0} + +# Maximum acceptable background count rates [counts/s]. There are separate thresholds +# for RAM vs. anti-RAM, and for pivot near 90 deg vs. others. +THRESHOLD_BG_RATE_RAM_90: float = 0.014 +THRESHOLD_BG_RATE_ANTI_RAM_90: float = 0.007 +THRESHOLD_BG_RATE_RAM_NON_90: float = 0.0175 +THRESHOLD_BG_RATE_ANTI_RAM_NON_90: float = 0.00875 + +# Maximum time gap [s] between consecutive histogram epochs before treating them as +# separate intervals. +DELAY_MAX: int = 100 +# Pivot angles within this range [degrees] are treated as "near 90". +PIVOT_90_RANGE: tuple[float, float] = 88.0, 92.0 +# Fraction of each cycle duration that contributes actual exposure. +EXPOSURE_FACTOR: float = 0.5 +# Padding [s] added to begin/end of each goodtime interval to ensure complete +# cycles are covered at interval edges. +GOODTIME_PADDING: float = 2.0 diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index b6ce01d517..9b834f5c90 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -2,6 +2,7 @@ import logging from dataclasses import Field +from datetime import datetime, timedelta from pathlib import Path import numpy as np @@ -9,6 +10,7 @@ import xarray as xr from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes +from imap_processing.lo import constants as c from imap_processing.lo import lo_ancillary from imap_processing.lo.l1b.tof_conversions import ( TOF0_CONV, @@ -34,6 +36,7 @@ interpolate_spin_data, ) from imap_processing.spice.time import ( + epoch_to_doy, epoch_to_fractional_doy, et_to_utc, met_to_ttj2000ns, @@ -172,22 +175,17 @@ # Fields to include in the split background rates/goodtimes datasets BACKGROUND_RATE_FIELDS = [ - "start_met", - "end_met", - "bin_start", - "bin_end", "h_background_rates", "h_background_variance", "o_background_rates", "o_background_variance", + "h_synthetic_floor", + "h_proxy_floor", + "o_synthetic_floor", + "o_proxy_floor", ] -GOODTIMES_FIELDS = [ - "gt_start_met", - "gt_end_met", - "bin_start", - "bin_end", - "esa_goodtime_flags", -] + +GOODTIMES_FIELDS = ["gt_start_met", "gt_end_met", "pivot", "pivot_de"] # ------------------------------------------------------------------- DE_CLOCK_TICK_S = 4.096e-3 # seconds per DE clock tick @@ -2501,8 +2499,7 @@ def l1b_star( def l1b_bgrates_and_goodtimes( # noqa: PLR0912 sci_dependencies: dict, attr_mgr_l1b: ImapCdfAttributes, - cycle_count: int = 10, - delay_max: int = 840, + delay_max: int | None = None, ) -> xr.Dataset: """ Create the IMAP-Lo L1B Background dataset. @@ -2515,10 +2512,9 @@ def l1b_bgrates_and_goodtimes( # noqa: PLR0912 Dictionary of datasets needed for L1B data product creation in xarray Datasets. attr_mgr_l1b : ImapCdfAttributes Attribute manager for L1B dataset metadata. - cycle_count : int - Maximum number of ASCs to group together (default: 10). - delay_max : int - Maximum allowed delay between entries in seconds (default: 840). + delay_max : int | None + Maximum time gap [s] between consecutive histogram epochs before treating them + as separate intervals. If None, the default value of 100 is used. Returns ------- @@ -2526,308 +2522,300 @@ def l1b_bgrates_and_goodtimes( # noqa: PLR0912 L1B bgrates dataset with ESA flags per epoch and bin. Each dataset also includes a background rate. """ - l1b_histrates = sci_dependencies["imap_lo_l1b_histrates"] - # l1b_nhk = sci_dependencies["imap_lo_l1b_nhk"] + if delay_max is None: + delay_max = c.DELAY_MAX - # Initialize the dataset - l1b_backgrounds_and_goodtimes_ds = xr.Dataset() - datasets_to_return = [] + elems = c.ELEMS # shortcut - # Set the expected background rate based on the pivot angle - # This assumes a static pivot_angle for the entire pointing - # pivot_angle = _get_nearest_pivot_angle(l1b_histrates["epoch"].values[0], l1b_nhk) - # if (pivot_angle < 95.0) & (pivot_angle > 85.0): - # h_bg_rate_nom = 0.0028 - # else: - # h_bg_rate_nom = 0.0033 - h_bg_rate_nom = 0.0028 - o_bg_rate_nom = h_bg_rate_nom / 100 - - interval_nom = 420 * cycle_count # seconds - exposure = interval_nom * 0.5 # 50% duty cycle - - h_intensity = np.sum( - l1b_histrates["h_counts"][:, 0:NUM_ESA_STEPS, 20:50], axis=(1, 2) - ) - o_intensity = np.sum( - l1b_histrates["o_counts"][:, 0:NUM_ESA_STEPS, 20:50], axis=(1, 2) - ) - - # Use proper SPICE-based time conversion with current kernels - # Note: The reference script adds +9 seconds because they use an - # "older time kernel (pre 2012)" - # We use current SPICE kernels, so we should NOT add that offset - met = ttj2000ns_to_met(l1b_histrates["epoch"].values) - - max_row_count = np.shape(h_intensity)[0] - bg_start_met = xr.DataArray([0.0]) - bg_end_met = xr.DataArray([0.0]) - epochs = l1b_histrates["epoch"].values.copy() - epochs = xr.DataArray(epochs, dims=["epoch"]) - goodtimes = xr.DataArray(np.zeros((max_row_count, 2), dtype=np.int64)) - h_background_rate = xr.DataArray(np.zeros((1, NUM_ESA_STEPS), dtype=np.float32)) - h_background_rate_variance = xr.DataArray( - np.zeros((1, NUM_ESA_STEPS), dtype=np.float32) - ) - o_background_rate = xr.DataArray(np.zeros((1, NUM_ESA_STEPS), dtype=np.float32)) - o_background_rate_variance = xr.DataArray( - np.zeros((1, NUM_ESA_STEPS), dtype=np.float32) - ) - - # Walk through the histrate data in chunks of cycle_count (10) - # and identify goodtime intervals and calculate background rates - row_count = 0 - sum_h_bg_counts = 0.0 - sum_h_bg_exposure = 0.0 - sum_o_bg_counts = 0.0 - begin = 0.0 - end = 0.0 - logger.debug( - f"Starting goodtimes calculation with {max_row_count} epochs, " - f"cycle_count={cycle_count}, delay_max={delay_max}" - ) - logger.debug(f"h_bg_rate_nom={h_bg_rate_nom}, exposure={exposure}") - for index in range(0, max_row_count, cycle_count): - # Calculate the interval for this chunk - if (index + cycle_count - 1) < max_row_count: - interval = met[index + cycle_count - 1] - met[index] - else: - interval = interval_nom * max_row_count + pivot_de: float = 0.0 + cdf_de = sci_dependencies.get("imap_lo_l1b_de") + if cdf_de is not None: + pivot_de = cdf_de["pivot_angle"].item() if "pivot_angle" in cdf_de else 0.0 - logger.debug( - f"\n Index {index}: met[{index}]=" - f"{met[index] if index < max_row_count else 'N/A'}, " - f"interval={interval}, begin={begin}" + pivot: float = 90.0 + cdf_hk = sci_dependencies.get("imap_lo_l1b_nhk") + if cdf_hk is not None and "pcc_coarse_pot_pri" in cdf_hk: + hk_epoch_ets = ttj2000ns_to_et(cdf_hk["epoch"]) + start_et_hk = ( + hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[0]).total_seconds() + ) + end_et_hk = ( + hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[1]).total_seconds() ) - # Skip this chunk if the interval is too large (indicates a gap) - if interval > (interval_nom + delay_max): - logger.debug( - f" Skipping chunk due to large interval ({interval} > " - f"{interval_nom + delay_max})" - ) - # If we were tracking a goodtime interval, close it before the gap - if begin > 0.0: - end = met[index - 1] - logger.debug(f" Closing interval before gap: {begin} -> {end}") - - epochs[row_count] = l1b_histrates["epoch"][index - 1].values.item() - goodtimes[row_count, :] = [int(begin - 620), int(end + 320)] - logger.debug( - f" STORED interval {row_count} (large interval): " - f"{int(begin - 620)} -> {int(end + 320)} (raw: {begin} -> {end})" - ) + coarse_pot_pri = cdf_hk["pcc_coarse_pot_pri"].values + pivot = np.nanmedian( # type: ignore + coarse_pot_pri[(hk_epoch_ets >= start_et_hk) & (hk_epoch_ets <= end_et_hk)] + ) + if np.isnan(pivot): + pivot = 90.0 + + cdf_hist = sci_dependencies["imap_lo_l1b_histrates"] + epoch_ttj2000 = cdf_hist["epoch"].values + n_epochs = epoch_ttj2000.shape[0] + met = ttj2000ns_to_met(epoch_ttj2000) + + # Get year and day-of-year for the anti-RAM threshold override lookup + epoch_utc_str = et_to_utc(ttj2000ns_to_et(epoch_ttj2000[0])) + epoch_start_dt = datetime.strptime(epoch_utc_str.split("T")[0], "%Y-%m-%d") + epoch_year = epoch_start_dt.year + epoch_doy = epoch_to_doy(epoch_ttj2000[:1])[0] + + # Choose background rate thresholds based on pivot orientation. + if c.PIVOT_90_RANGE[0] < pivot < c.PIVOT_90_RANGE[1]: + bg_rate_ram_nominal = c.THRESHOLD_BG_RATE_RAM_90 + bg_rate_anti_ram_nominal = c.THRESHOLD_BG_RATE_ANTI_RAM_90 + else: + bg_rate_ram_nominal = c.THRESHOLD_BG_RATE_RAM_NON_90 + bg_rate_anti_ram_nominal = c.THRESHOLD_BG_RATE_ANTI_RAM_NON_90 - row_count += 1 - begin = 0.0 - end = 0.0 + # Manual overrides of the anti-RAM threshold for anomalous days. + bg_rate_anti_ram_nominal = c.BG_RATE_ANTI_RAM_OVERRIDES.get( + (epoch_year, epoch_doy), bg_rate_anti_ram_nominal + ) - # Skip this chunk after closing interval - continue + ram_esa_indices = [i - 1 for i in c.RAM_ESA_LEVELS] # Convert to 0-indexed - # Check for time gap from previous chunk - delta_time = 0.0 - if index > 0: - delta_time = met[index] - (met[index - 1] + 420) - logger.debug( - f" Delta time from previous: {delta_time} (max: {delay_max})" - ) + # Sum histogram counts over the relevant angular bins for each species and + # direction. RAM counts use only certain ESA steps; anti-RAM counts use all. + elem_ram_counts = {} + elem_anti_ram_counts = {} + for elem in elems: + elem_counts = cdf_hist[f"{elem.lower()}_counts"].values + elem_ram_counts[elem] = sum( + np.sum(elem_counts[:, ram_esa_indices, b], axis=(1, 2)) + for b in c.RAM_HISTOGRAM_BINS + ) + elem_anti_ram_counts[elem] = sum( + np.sum(elem_counts[:, :, b], axis=(1, 2)) for b in c.ANTI_RAM_HISTOGRAM_BINS + ) - # If there's a gap and we have an active interval, close it - if (delta_time > delay_max) & (begin > 0.0): - end = met[index - 1] - logger.debug(f" Closing interval due to time gap: {begin} -> {end}") + # Pre-compute expected exposure times [s] for the averaging and summing windows. + exposure = c.HISTOGRAM_CYCLE_EPOCHS * c.N_CYCLE_AVE * c.EXPOSURE_FACTOR + exposure_ram = exposure * len(c.RAM_ESA_LEVELS) / c.N_ESA_LEVELS + exposure_sum = c.HISTOGRAM_CYCLE_EPOCHS * c.N_CYCLE_SUM * c.EXPOSURE_FACTOR + + # Walk through histogram epochs one N_CYCLE_SUM block at a time. + begin = end = 0.0 + interval = c.HISTOGRAM_CYCLE_EPOCHS * c.N_CYCLE_SUM + synthetic_floors = {e: 0.0 for e in elems} # Accumulated model-predicted BG counts + proxy_floors = { + e: 0.0 for e in elems + } # Accumulated measured anti-RAM counts (BG proxy) + goodtime_exposure_avg = goodtime_exposure_sum = 0.0 + goodtime_rows = [] + + for i in range(0, n_epochs, c.N_CYCLE_SUM): + measured_interval = interval + if i + c.N_CYCLE_SUM < n_epochs: + measured_interval = met[i + c.N_CYCLE_SUM] - met[i] + + if measured_interval > (interval + delay_max): + if begin > 0.0: + end = met[i - 1] + goodtime_rows.append( + ( + begin, + end, + bg_rate_anti_ram_nominal, + goodtime_exposure_avg, + goodtime_exposure_sum, + ) + ) + begin = end = 0.0 + continue - epochs[row_count] = l1b_histrates["epoch"][index - 1].values.item() - goodtimes[row_count, :] = [int(begin - 620), int(end + 320)] - logger.debug( - f" STORED interval {row_count} (time gap): " - f"{int(begin - 620)} -> {int(end + 320)} (raw: {begin} -> {end})" + # A large gap (missing data) forces the current good-time interval to close. + delta_time = 0.0 + if i > 0: + delta_time = met[i] - (met[i - 1] + c.HISTOGRAM_CYCLE_EPOCHS) + + if (delta_time > c.DELAY_MAX) and (begin > 0.0): + end = met[i - 1] + goodtime_rows.append( + ( + begin, + end, + bg_rate_anti_ram_nominal, + goodtime_exposure_avg, + goodtime_exposure_sum, + ) ) - - row_count += 1 - begin = 0.0 - end = 0.0 - - # Calculate counts and rate for this chunk - antiram_h_counts = float(np.sum(h_intensity[index : index + cycle_count])) - antiram_o_counts = float(np.sum(o_intensity[index : index + cycle_count])) - antiram_h_rate = antiram_h_counts / exposure - - logger.debug( - f" Rate: {antiram_h_rate:.6f}, threshold: {h_bg_rate_nom:.6f}, " - f"counts: {antiram_h_counts}" + begin = end = 0.0 + + # Sliding window centered on epoch i for rate averaging + window_avg_start = max(int(i - c.N_CYCLE_AVE // 2), 0) + window_avg_end = min(n_epochs, window_avg_start + c.N_CYCLE_AVE) + if (window_avg_end - window_avg_start) < c.N_CYCLE_AVE: + window_avg_start = max(window_avg_end - c.N_CYCLE_AVE, 0) + + # Sliding window centered on epoch i for accumulating counts + window_sum_start = max(int(i - c.N_CYCLE_SUM // 2), 0) + window_sum_end = min(n_epochs, window_sum_start + c.N_CYCLE_SUM) + if (window_sum_end - window_sum_start) < c.N_CYCLE_SUM: + window_sum_start = max(window_avg_end - c.N_CYCLE_SUM, 0) + + # Estimate background rates from the averaged H counts + ram_rate = ( + np.sum(elem_ram_counts["H"][window_avg_start:window_avg_end]) / exposure_ram + ) + anti_ram_rate = ( + np.sum(elem_anti_ram_counts["H"][window_avg_start:window_avg_end]) + / exposure ) - # If rate is below threshold, accumulate for background - if antiram_h_rate < h_bg_rate_nom: + # good-time = intervals where background rates are below threshold + if (ram_rate < bg_rate_ram_nominal) and ( + anti_ram_rate < bg_rate_anti_ram_nominal + ): if begin == 0.0: - begin = met[index] - logger.debug(f" Starting new interval at {begin}") + begin = met[i] # Start a new good-time interval - sum_h_bg_counts = sum_h_bg_counts + antiram_h_counts - sum_o_bg_counts = sum_o_bg_counts + antiram_o_counts - sum_h_bg_exposure = sum_h_bg_exposure + exposure - - # If rate exceeds threshold, close the interval if one is active - if antiram_h_rate >= h_bg_rate_nom: - if begin > 0.0: - end = met[index - 1] - logger.debug( - f" Closing interval due to rate threshold: {begin} -> {end}" - ) - print(" antiram_h_rate: ", antiram_h_rate, " at index ", index) - print("l1b_histrates epoch: ", l1b_histrates["epoch"][index - 1].values) - epochs[row_count] = l1b_histrates["epoch"][index - 1].values.item() - goodtimes[row_count, :] = [int(begin - 620), int(end + 320)] - logger.debug( - f" STORED interval {row_count} (rate threshold): " - f"{int(begin - 620)} -> {int(end + 320)} (raw: {begin} -> {end})" + for elem in elems: + synthetic_floors[elem] += c.BG_RATES[elem] * exposure + proxy_floors[elem] += np.sum( + elem_anti_ram_counts["H"][window_sum_start:window_sum_end] ) - row_count += 1 - begin = 0.0 - end = 0.0 + goodtime_exposure_avg += exposure + goodtime_exposure_sum += exposure_sum + + elif begin > 0.0: + # Background exceeded threshold; close the current good-time interval. + end = met[i - 1] + goodtime_rows.append( + ( + begin, + end, + bg_rate_anti_ram_nominal, + goodtime_exposure_avg, + goodtime_exposure_sum, + ) + ) + begin = end = 0.0 - # Handle the final interval if one is still open - if (end == 0.0) & (begin > 0.0): - end = met[max_row_count - 1] + if (end == 0.0) and (begin > 0.0): + end = met[n_epochs - 1] if end > begin: - epochs[row_count] = l1b_histrates["epoch"][max_row_count - 1] - goodtimes[row_count, :] = [int(begin - 620), int(end + 320)] - logger.debug( - f" STORED interval {row_count} (final): " - f"{int(begin - 620)} -> {int(end + 320)} (raw: {begin} -> {end})" + goodtime_rows.append( + ( + begin, + end, + bg_rate_anti_ram_nominal, + goodtime_exposure_avg, + goodtime_exposure_sum, + ) ) - row_count += 1 - begin = 0.0 - end = 0.0 - - # Record the background rates for the entire pointing - if sum_h_bg_exposure > 0.0: - h_bg_rate = sum_h_bg_counts / sum_h_bg_exposure - h_bg_rate_variance = np.sqrt(sum_h_bg_counts) / sum_h_bg_exposure - o_bg_rate = sum_o_bg_counts / sum_h_bg_exposure - o_bg_rate_variance = np.sqrt(sum_o_bg_counts) / sum_h_bg_exposure - - if h_bg_rate_variance <= 0.0: - h_bg_rate_variance = h_bg_rate - - if o_bg_rate_variance <= 0.0: - o_bg_rate_variance = o_bg_rate - - if h_bg_rate <= 0.0: - h_bg_rate = h_bg_rate_nom / 50.0 - h_bg_rate_variance = h_bg_rate - - if o_bg_rate <= 0.0: - o_bg_rate = o_bg_rate_nom * 0.3 - o_bg_rate_variance = o_bg_rate - - h_background_rate[0, :] = np.full(NUM_ESA_STEPS, h_bg_rate) - h_background_rate_variance[0, :] = np.full(NUM_ESA_STEPS, h_bg_rate_variance) - o_background_rate[0, :] = np.full(NUM_ESA_STEPS, o_bg_rate) - o_background_rate_variance[0, :] = np.full(NUM_ESA_STEPS, o_bg_rate_variance) - bg_start_met[0] = met[0] - bg_end_met[0] = met[max_row_count - 1] - - # Handle case where no goodtimes were found -- produce a - # single record with invalid times (the defaults above) - if row_count == 0: - row_count = 1 - - # Trim arrays to actual size - epoch = epochs.isel(epoch=slice(0, row_count)) - goodtimes = goodtimes.isel(dim_0=slice(0, row_count)) - - l1b_backgrounds_and_goodtimes_ds["epoch"] = xr.DataArray( - data=epoch, + # Compute background rates per species + bg_rates_out = {} + sigma_bg_rates_out = {} + for elem in elems: + if goodtime_exposure_avg == 0: + bg_rate = bg_rate_anti_ram_nominal * c.BG_RATE_FALLBACK_SCALE[elem] + sigma_bg_rate = bg_rate + else: + bg_rate = synthetic_floors[elem] / goodtime_exposure_avg + sigma_bg_rate = np.sqrt(synthetic_floors[elem]) / goodtime_exposure_avg + + if bg_rate == 0.0: + bg_rate = bg_rate_anti_ram_nominal / c.BG_RATE_FLOOR_DIVISOR[elem] + sigma_bg_rate = bg_rate + if sigma_bg_rate == 0.0: + sigma_bg_rate = bg_rate + + bg_rates_out[elem] = bg_rate + sigma_bg_rates_out[elem] = sigma_bg_rate + + # Final adjustment - add padding to each goodtime interval + for i, (begin, end, *other) in enumerate(goodtime_rows): + goodtime_rows[i] = ( + begin - c.GOODTIME_PADDING, + end + c.GOODTIME_PADDING, + *other, + ) + + if len(goodtime_rows) == 0: + goodtime_rows = [(0, 0, 0, 0, 0)] + + # Initialize the dataset + datasets_to_return = [] + l1b_combined_ds = xr.Dataset() + + epoch_values = met_to_ttj2000ns(np.array([r[0] for r in goodtime_rows])) + + l1b_combined_ds["epoch"] = xr.DataArray( + data=epoch_values, name="epoch", dims=["epoch"], attrs=attr_mgr_l1b.get_variable_attributes("epoch"), ) - l1b_backgrounds_and_goodtimes_ds["epoch"].attrs["DEPEND_0"] = "epoch" - l1b_backgrounds_and_goodtimes_ds["start_met"] = xr.DataArray( - data=bg_start_met, - name="start_met", - dims=["met"], - attrs=attr_mgr_l1b.get_variable_attributes("met"), + l1b_combined_ds["epoch"].attrs["DEPEND_0"] = "epoch" + + l1b_combined_ds["pivot"] = xr.DataArray( + data=np.float32(pivot), + name="pivot", + attrs=attr_mgr_l1b.get_variable_attributes("pivot"), ) - l1b_backgrounds_and_goodtimes_ds["end_met"] = xr.DataArray( - data=bg_end_met, - name="end_met", - dims=["met"], - attrs=attr_mgr_l1b.get_variable_attributes("met"), + l1b_combined_ds["pivot_de"] = xr.DataArray( + data=np.float32(pivot_de), + name="pivot_de", + attrs=attr_mgr_l1b.get_variable_attributes("pivot_de"), ) - l1b_backgrounds_and_goodtimes_ds["gt_start_met"] = xr.DataArray( - data=goodtimes[:, 0], + + l1b_combined_ds["gt_start_met"] = xr.DataArray( + data=np.array([r[0] for r in goodtime_rows], dtype=np.float32), name="Goodtime_start", dims=["epoch"], - # attrs=attr_mgr_l1b.get_variable_attributes("epoch"), + attrs=attr_mgr_l1b.get_variable_attributes("gt_start_met"), ) - l1b_backgrounds_and_goodtimes_ds["gt_end_met"] = xr.DataArray( - data=goodtimes[:, 1], + l1b_combined_ds["gt_end_met"] = xr.DataArray( + data=np.array([r[1] for r in goodtime_rows], dtype=np.float32), name="Goodtime_end", dims=["epoch"], - # attrs=attr_mgr_l1b.get_variable_attributes("epoch"), - ) - l1b_backgrounds_and_goodtimes_ds["h_background_rates"] = xr.DataArray( - data=h_background_rate, - name="h_bg_rate", - dims=["met", "esa_step"], - # attrs=attr_mgr_l1b.get_variable_attributes("esa_background_rates"), - ) - l1b_backgrounds_and_goodtimes_ds["h_background_variance"] = xr.DataArray( - data=h_background_rate_variance, - name="h_bg_rate_variance", - dims=["met", "esa_step"], - ) - l1b_backgrounds_and_goodtimes_ds["o_background_rates"] = xr.DataArray( - data=o_background_rate, - name="o_bg_rate", - dims=["met", "esa_step"], - # attrs=attr_mgr_l1b.get_variable_attributes("esa_background_rates"), - ) - l1b_backgrounds_and_goodtimes_ds["o_background_variance"] = xr.DataArray( - data=o_background_rate_variance, - name="o_bg_rate_variance", - dims=["met", "esa_step"], - ) - - # We're only creating one record for all bins for now - # Note that this is true for both GoodTimes and background rates, - # so we cheat here by just using one record. - l1b_backgrounds_and_goodtimes_ds["bin_start"] = xr.DataArray( - data=np.zeros(row_count, dtype=int), - name="bin_start", - dims=["epoch"], - # attrs=attr_mgr_l1b.get_variable_attributes("bin_start"), - ) - l1b_backgrounds_and_goodtimes_ds["bin_end"] = xr.DataArray( - data=np.zeros(row_count, dtype=int) + 59, - name="bin_end", - dims=["epoch"], - # attrs=attr_mgr_l1b.get_variable_attributes("bin_end"), - ) - - # For now, set all ESA flags to 1 (good) since we don't have - # an algorithm for this yet - l1b_backgrounds_and_goodtimes_ds["esa_goodtime_flags"] = xr.DataArray( - data=np.zeros((row_count, NUM_ESA_STEPS), dtype=int) + 1, - name="E-step", - dims=["epoch", "esa_step"], - # attrs=attr_mgr_l1b.get_variable_attributes("esa_goodtime_flags"), - ) + attrs=attr_mgr_l1b.get_variable_attributes("gt_end_met"), + ) + + # Per-species scalar variables + for elem in elems: + elem_lower = elem.lower() + # For *_background_rates, and *_background_variance for each species, + # we return a (N_ESA_LEVELS) array of identical values to be backward + # compatible with an old implementation of the algorithm. + l1b_combined_ds[f"{elem_lower}_background_rates"] = xr.DataArray( + data=np.full(c.N_ESA_LEVELS, bg_rates_out[elem]), + name=f"{elem_lower}_background_rates", + attrs=attr_mgr_l1b.get_variable_attributes( + f"{elem_lower}_background_rates" + ), + dims=["esa_step"], + ) + l1b_combined_ds[f"{elem_lower}_background_variance"] = xr.DataArray( + data=np.full(c.N_ESA_LEVELS, sigma_bg_rates_out[elem]), + name=f"{elem_lower}_background_variance", + attrs=attr_mgr_l1b.get_variable_attributes( + f"{elem_lower}_background_variance" + ), + dims=["esa_step"], + ) + l1b_combined_ds[f"{elem_lower}_synthetic_floor"] = xr.DataArray( + data=np.float32(synthetic_floors[elem]), + name=f"{elem_lower}_synthetic_floor", + attrs=attr_mgr_l1b.get_variable_attributes(f"{elem_lower}_synthetic_floor"), + ) + l1b_combined_ds[f"{elem_lower}_proxy_floor"] = xr.DataArray( + data=np.float32(proxy_floors[elem]), + name=f"{elem_lower}_proxy_floor", + attrs=attr_mgr_l1b.get_variable_attributes(f"{elem_lower}_proxy_floor"), + ) - logger.info("L1B Background Rates and Goodtimes created successfully") + logger.info("L1B Background Rates and Bettertimes created successfully") l1b_bgrates_ds, l1b_goodtimes_ds = split_backgrounds_and_goodtimes_dataset( - l1b_backgrounds_and_goodtimes_ds, attr_mgr_l1b + l1b_combined_ds, attr_mgr_l1b ) datasets_to_return.extend([l1b_bgrates_ds, l1b_goodtimes_ds]) - print("epoch bgrates meta", l1b_bgrates_ds["epoch"].attrs) - print("epoch goodtimes meta", l1b_goodtimes_ds["epoch"].attrs) + return datasets_to_return @@ -2853,7 +2841,6 @@ def split_backgrounds_and_goodtimes_dataset( l1b_goodtimes_rates : xr.Dataset The L1B goodtimes rates dataset. """ - # Use centralized lists for fields to include in split datasets l1b_goodtimes_ds = l1b_backgrounds_and_goodtimes_ds[GOODTIMES_FIELDS] l1b_goodtimes_ds.attrs = attr_mgr_l1b.get_global_attributes("imap_lo_l1b_goodtimes") lib_bgrates_ds = l1b_backgrounds_and_goodtimes_ds[BACKGROUND_RATE_FIELDS] diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index a6fc24b845..3e269202fd 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -2195,25 +2195,16 @@ def test_l1b_bgrates_and_goodtimes_basic(attr_mgr_l1b): """Test basic functionality of l1b_bgrates_and_goodtimes.""" # Arrange - Create a simple L1B histogram rates dataset # with enough data points to create goodtime intervals - num_epochs = 100 # 10 cycles of 10 epochs each - met_start = 473389200 # Start MET time - met_spacing = 42 # seconds between epochs + num_epochs = 100 + met_start = 473389200 + met_spacing = 42 - # Create evenly spaced MET times met_times = np.arange(met_start, met_start + num_epochs * met_spacing, met_spacing) epoch_times = met_to_ttj2000ns(met_times) - # Create counts data with low background rates (below threshold) - # h_bg_rate_nom = 0.0028, exposure = 420*10*0.5 = 2100 seconds - # To be below threshold: rate = counts / exposure < 0.0028 - # Summed over 7 ESA steps * 30 azimuth bins * 10 epochs = 2100 values - # Max total counts per chunk: 0.0028 * 2100 = 5.88 counts - # Use 10% of max for safety: 5.88 / 2100 / 10 = 0.00028 per element - h_counts_per_epoch = 0.00028 # Low counts to ensure below threshold - o_counts_per_epoch = 0.000028 # 10x smaller for oxygen - - h_counts = np.ones((num_epochs, 7, 60)) * h_counts_per_epoch - o_counts = np.ones((num_epochs, 7, 60)) * o_counts_per_epoch + # Low counts to keep background rates below threshold + h_counts = np.ones((num_epochs, 7, 60)) * 0.00028 + o_counts = np.ones((num_epochs, 7, 60)) * 0.000028 l1b_histrates = xr.Dataset( { @@ -2227,12 +2218,14 @@ def test_l1b_bgrates_and_goodtimes_basic(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert - Should return a list with two datasets assert isinstance(result, list) @@ -2240,23 +2233,21 @@ def test_l1b_bgrates_and_goodtimes_basic(attr_mgr_l1b): l1b_bgrates_ds, l1b_goodtimes_ds = result - # Check bgrates dataset structure + # Check bgrates dataset structure (BACKGROUND_RATE_FIELDS) assert "h_background_rates" in l1b_bgrates_ds.data_vars assert "h_background_variance" in l1b_bgrates_ds.data_vars + assert "h_synthetic_floor" in l1b_bgrates_ds.data_vars + assert "h_proxy_floor" in l1b_bgrates_ds.data_vars assert "o_background_rates" in l1b_bgrates_ds.data_vars assert "o_background_variance" in l1b_bgrates_ds.data_vars - # Note: bgrates uses 'met' dimension, goodtimes has epoch in data vars + assert "o_synthetic_floor" in l1b_bgrates_ds.data_vars + assert "o_proxy_floor" in l1b_bgrates_ds.data_vars - # Check goodtimes dataset structure + # Check goodtimes dataset structure (GOODTIMES_FIELDS) assert "gt_start_met" in l1b_goodtimes_ds.data_vars assert "gt_end_met" in l1b_goodtimes_ds.data_vars - assert "bin_start" in l1b_goodtimes_ds.data_vars - assert "bin_end" in l1b_goodtimes_ds.data_vars - assert "esa_goodtime_flags" in l1b_goodtimes_ds.data_vars - - # Check dimensions - assert l1b_bgrates_ds["h_background_rates"].dims == ("met", "esa_step") - assert l1b_bgrates_ds["h_background_rates"].shape[1] == 7 # 7 ESA steps + assert "pivot" in l1b_goodtimes_ds.data_vars + assert "pivot_de" in l1b_goodtimes_ds.data_vars # Check that goodtime intervals were created assert len(l1b_goodtimes_ds["gt_start_met"]) > 0 @@ -2267,13 +2258,6 @@ def test_l1b_bgrates_and_goodtimes_basic(attr_mgr_l1b): l1b_goodtimes_ds["gt_start_met"].values <= l1b_goodtimes_ds["gt_end_met"].values ) - # Check bin_start and bin_end values - assert np.all(l1b_goodtimes_ds["bin_start"].values == 0) - assert np.all(l1b_goodtimes_ds["bin_end"].values == 59) - - # Check ESA goodtime flags are all 1 (good) - assert np.all(l1b_goodtimes_ds["esa_goodtime_flags"].values == 1) - def test_l1b_bgrates_and_goodtimes_with_gap(attr_mgr_l1b): """Test l1b_bgrates_and_goodtimes handles data gaps correctly.""" @@ -2317,12 +2301,14 @@ def test_l1b_bgrates_and_goodtimes_with_gap(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2351,24 +2337,19 @@ def test_l1b_bgrates_and_goodtimes_high_rate(attr_mgr_l1b): epoch_times = met_to_ttj2000ns(met_times) # Create high counts (above threshold) - # h_bg_rate_nom = 0.0028, exposure = 420*10*0.5 = 2100 seconds - # To be above threshold: rate > 0.0028 - # Use 10x threshold for high rate periods: 0.028 counts/sec - # That's 0.028 * 2100 / 2100_values = 0.028 per element - h_counts = np.ones((num_epochs, 7, 60)) * 0.028 # High rate (10x threshold) - o_counts = np.ones((num_epochs, 7, 60)) * 0.0028 + # h_bg_rate_nom = 0.0014925, exposure = 420*7*0.5 = 1470 seconds + # To be above threshold: rate > 0.0014925 + # Use 10x threshold for high rate periods: 0.014925 counts/sec + h_counts = np.ones((num_epochs, 7, 60)) * 0.014925 # High rate (10x threshold) + o_counts = np.ones((num_epochs, 7, 60)) * 0.0014925 # Make first 20 epochs low (below threshold) - h_counts[:20, :, :] = 0.00028 - o_counts[:20, :, :] = 0.000028 - - # Make middle 60 epochs high (above threshold) - h_counts[20:80, :, :] = 0.028 - o_counts[20:80, :, :] = 0.0028 + h_counts[:20, :, :] = 0.00014925 + o_counts[:20, :, :] = 0.000014925 - # Make last 20 epochs low again - h_counts[80:, :, :] = 0.00028 - o_counts[80:, :, :] = 0.000028 + # Make last 20 epochs low + h_counts[80:, :, :] = 0.00014925 + o_counts[80:, :, :] = 0.000014925 l1b_histrates = xr.Dataset( { @@ -2382,12 +2363,14 @@ def test_l1b_bgrates_and_goodtimes_high_rate(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2425,66 +2408,20 @@ def test_l1b_bgrates_and_goodtimes_no_goodtimes(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } - bgrates_ds, goodtimes_ds = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 + _, goodtimes_ds = l1b_bgrates_and_goodtimes( + sci_dependencies, attr_mgr_l1b, delay_max=840 ) - # Function should return two datasets - assert "h_background_rates" in bgrates_ds.data_vars - # Goodtimes dataset should exist and contain the gt_* fields - # (defaults when none found) - assert "gt_start_met" in goodtimes_ds.data_vars - assert "gt_end_met" in goodtimes_ds.data_vars - # When no goodtimes were detected the default invalid times are used (zeros) + # When no goodtimes are detected a single fallback row (0, 0) is used. + # The padding loop runs before the fallback is inserted, so the zeros are unchanged. assert int(goodtimes_ds["gt_start_met"].values[0]) == 0 assert int(goodtimes_ds["gt_end_met"].values[0]) == 0 - assert int(bgrates_ds["start_met"].values[0]) == 0 - assert int(bgrates_ds["end_met"].values[0]) == 0 - - -def test_l1b_bgrates_and_goodtimes_custom_cycle_count(attr_mgr_l1b): - """Test l1b_bgrates_and_goodtimes with custom cycle_count parameter.""" - # Arrange - num_epochs = 50 - met_start = 473389200 - met_spacing = 42 - - met_times = np.arange(met_start, met_start + num_epochs * met_spacing, met_spacing) - epoch_times = met_to_ttj2000ns(met_times) - - # Low counts (below threshold) - h_counts = np.ones((num_epochs, 7, 60)) * 0.00028 - o_counts = np.ones((num_epochs, 7, 60)) * 0.000028 - - l1b_histrates = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), h_counts), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), o_counts), - }, - coords={ - "epoch": epoch_times, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - }, - ) - - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} - - # Act - Use different cycle_count - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=5, delay_max=420 - ) - - # Assert - l1b_bgrates_ds, l1b_goodtimes_ds = result - - # Should successfully create datasets with custom parameters - assert len(l1b_goodtimes_ds["gt_start_met"]) > 0 - # Background rates should be calculated from the low-count period - assert np.all(l1b_bgrates_ds["h_background_rates"].values > 0) - assert np.all(l1b_bgrates_ds["o_background_rates"].values > 0) def test_l1b_bgrates_and_goodtimes_empty_dataset(attr_mgr_l1b): @@ -2513,12 +2450,14 @@ def test_l1b_bgrates_and_goodtimes_empty_dataset(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert - Should still create valid datasets even with minimal data l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2529,102 +2468,71 @@ def test_l1b_bgrates_and_goodtimes_empty_dataset(attr_mgr_l1b): def test_split_backgrounds_and_goodtimes_dataset(attr_mgr_l1b): """Test split_backgrounds_and_goodtimes_dataset separates fields correctly.""" - # Arrange - Create a combined dataset with both background and goodtime fields - num_records = 5 - epoch_times = met_to_ttj2000ns( - np.arange(473389200, 473389200 + num_records * 420, 420) - ) + # Arrange - Create a combined dataset matching the structure produced by + # l1b_bgrates_and_goodtimes: scalar background rate fields and epoch-indexed + # goodtime interval fields. + num_records = 3 + met_starts = np.arange(473389200, 473389200 + num_records * 420, 420) + epoch_times = met_to_ttj2000ns(met_starts) combined_ds = xr.Dataset( - { - # Background rate fields - "epoch": ("epoch", epoch_times), - "h_background_rates": (("met", "esa_step"), np.random.rand(num_records, 7)), - "h_background_variance": ( - ("met", "esa_step"), - np.random.rand(num_records, 7), - ), - "o_background_rates": (("met", "esa_step"), np.random.rand(num_records, 7)), - "o_background_variance": ( - ("met", "esa_step"), - np.random.rand(num_records, 7), - ), - # Goodtime fields - "gt_start_met": ( - "met", - np.arange(473389200, 473389200 + num_records * 420, 420), - ), - "gt_end_met": ( - "met", - np.arange(473389200 + 400, 473389200 + num_records * 420 + 400, 420), - ), - # Also include non-prefixed background start/end fields so - # split_backgrounds_and_goodtimes_dataset can select - "start_met": ( - "met", - np.arange(473389200, 473389200 + num_records * 420, 420), - ), - "end_met": ( - "met", - np.arange(473389200 + 400, 473389200 + num_records * 420 + 400, 420), - ), - "bin_start": ("met", np.zeros(num_records, dtype=int)), - "bin_end": ("met", np.zeros(num_records, dtype=int) + 59), - "esa_goodtime_flags": ( - ("met", "esa_step"), - np.ones((num_records, 7), dtype=int), - ), - }, - coords={ - "met": np.arange(num_records), - "esa_step": np.arange(1, 8), - }, + coords={"epoch": epoch_times}, ) + combined_ds["gt_start_met"] = xr.DataArray( + met_starts.astype(np.int64), dims=["epoch"] + ) + combined_ds["gt_end_met"] = xr.DataArray( + (met_starts + 400).astype(np.int64), dims=["epoch"] + ) + combined_ds["pivot"] = xr.DataArray(np.float32(90.0)) + combined_ds["pivot_de"] = xr.DataArray(np.float32(89.5)) + combined_ds["h_background_rates"] = xr.DataArray(np.float32(0.01)) + combined_ds["h_background_variance"] = xr.DataArray(np.float32(0.001)) + combined_ds["o_background_rates"] = xr.DataArray(np.float32(0.002)) + combined_ds["o_background_variance"] = xr.DataArray(np.float32(0.0002)) + combined_ds["h_synthetic_floor"] = xr.DataArray(np.float32(5.0)) + combined_ds["h_proxy_floor"] = xr.DataArray(np.float32(4.0)) + combined_ds["o_synthetic_floor"] = xr.DataArray(np.float32(0.5)) + combined_ds["o_proxy_floor"] = xr.DataArray(np.float32(0.4)) # Act bgrates_ds, goodtimes_ds = split_backgrounds_and_goodtimes_dataset( combined_ds, attr_mgr_l1b ) - # Assert - Check bgrates dataset has background fields - # Note: bgrates includes 'start_met', 'end_met', 'bin_start', 'bin_end' per - # BACKGROUND_RATE_FIELDS - assert "h_background_rates" in bgrates_ds.data_vars - assert "h_background_variance" in bgrates_ds.data_vars - assert "o_background_rates" in bgrates_ds.data_vars - assert "o_background_variance" in bgrates_ds.data_vars - # Note: bgrates uses 'met' dimension, goodtimes has epoch in data vars - - # Check goodtimes dataset structure + # Assert - bgrates dataset contains all background rate fields (scalar) + for field in [ + "h_background_rates", + "h_background_variance", + "o_background_rates", + "o_background_variance", + "h_synthetic_floor", + "h_proxy_floor", + "o_synthetic_floor", + "o_proxy_floor", + ]: + assert field in bgrates_ds.data_vars + assert bgrates_ds[field].dims == () # scalar + + # Assert - goodtimes dataset contains the expected fields assert "gt_start_met" in goodtimes_ds.data_vars assert "gt_end_met" in goodtimes_ds.data_vars - assert "bin_start" in goodtimes_ds.data_vars - assert "bin_end" in goodtimes_ds.data_vars - assert "esa_goodtime_flags" in goodtimes_ds.data_vars + assert "pivot" in goodtimes_ds.data_vars + assert "pivot_de" in goodtimes_ds.data_vars - # Check dimensions - assert bgrates_ds["h_background_rates"].dims == ("met", "esa_step") - assert bgrates_ds["h_background_rates"].shape[1] == 7 # 7 ESA steps - - # Check that goodtime intervals were created + # Assert - goodtime intervals were created and are valid assert len(goodtimes_ds["gt_start_met"]) > 0 - assert len(goodtimes_ds["gt_end_met"]) > 0 - - # Check that start times are before end times assert np.all( goodtimes_ds["gt_start_met"].values <= goodtimes_ds["gt_end_met"].values ) - # Check bin_start and bin_end values - assert np.all(goodtimes_ds["bin_start"].values == 0) - assert np.all(goodtimes_ds["bin_end"].values == 59) - - # Check ESA goodtime flags are all 1 (good) - assert np.all(goodtimes_ds["esa_goodtime_flags"].values == 1) + # Assert - pivot fields are scalar + assert goodtimes_ds["pivot"].dims == () + assert goodtimes_ds["pivot_de"].dims == () -def test_l1b_bgrates_and_goodtimes_azimuth_bins(attr_mgr_l1b): - """Test that the function correctly uses azimuth bins 20-50 for calculations.""" +def test_l1b_bgrates_and_goodtimes_ram_and_anti_ram_bins(attr_mgr_l1b): + """Test that the function correctly uses bins anti-RAM 20-50 and RAM 0-20/50-60.""" # Arrange - Create dataset with specific counts in different azimuth bins num_epochs = 30 met_start = 473389200 @@ -2633,16 +2541,21 @@ def test_l1b_bgrates_and_goodtimes_azimuth_bins(attr_mgr_l1b): met_times = np.arange(met_start, met_start + num_epochs * met_spacing, met_spacing) epoch_times = met_to_ttj2000ns(met_times) - # Set high counts outside bins 20-50, low counts inside bins 20-50 - h_counts = ( - np.ones((num_epochs, 7, 60)) * 0.028 - ) # High counts (10x threshold) everywhere + # High counts everywhere by default + h_counts = np.ones((num_epochs, 7, 60)) * 0.028 o_counts = np.ones((num_epochs, 7, 60)) * 0.0028 - # Set low counts in the bins that are actually used (20-50) - h_counts[:, :, 20:50] = 0.00028 # Low counts in used bins (below threshold) + # Anti-RAM bins (20-50): set low across all ESA steps (below anti-RAM threshold) + h_counts[:, :, 20:50] = 0.00028 o_counts[:, :, 20:50] = 0.000028 + # RAM bins (0-20, 50-60) for RAM ESA steps (0-indexed 5, 6): + # set low (below RAM threshold) + h_counts[:, 5:7, 0:20] = 0.00028 + h_counts[:, 5:7, 50:60] = 0.00028 + o_counts[:, 5:7, 0:20] = 0.000028 + o_counts[:, 5:7, 50:60] = 0.000028 + l1b_histrates = xr.Dataset( { "h_counts": (("epoch", "esa_step", "spin_bin_6"), h_counts), @@ -2655,19 +2568,23 @@ def test_l1b_bgrates_and_goodtimes_azimuth_bins(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + # Required dependencies added in the updated function signature + cdf_de = xr.Dataset({"pivot_angle": xr.DataArray(90.0)}) + cdf_hk = xr.Dataset() # No pcc_coarse_pot_pri; pivot defaults to 90.0 - # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": cdf_de, + "imap_lo_l1b_nhk": cdf_hk, + } - # Assert - Should create goodtime intervals because bins 20-50 have low counts + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) l1b_bgrates_ds, l1b_goodtimes_ds = result + # Should create goodtime intervals because RAM and anti-RAM bins have low counts assert len(l1b_goodtimes_ds["gt_start_met"]) > 0 - # Background rates should be calculated from the low-count bins - assert np.all(l1b_bgrates_ds["h_background_rates"].values < 1.0) + # h_synthetic_floor accumulates the modeled background during good times + assert l1b_bgrates_ds["h_synthetic_floor"].values > 0 def test_l1b_bgrates_and_goodtimes_variance_calculation(attr_mgr_l1b): @@ -2700,12 +2617,14 @@ def test_l1b_bgrates_and_goodtimes_variance_calculation(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2720,7 +2639,7 @@ def test_l1b_bgrates_and_goodtimes_variance_calculation(attr_mgr_l1b): def test_l1b_bgrates_and_goodtimes_offset_application(attr_mgr_l1b): - """Test that goodtime start/end offsets (-620, +320) are applied correctly.""" + """Test that padding is applied to goodtime intervals.""" # Arrange num_epochs = 30 met_start = 473389200 @@ -2745,26 +2664,27 @@ def test_l1b_bgrates_and_goodtimes_offset_application(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result - # Check that gt_start_met is earlier than gt_end_met (accounting for offsets) - for i in range(len(l1b_goodtimes_ds["gt_start_met"])): - start = l1b_goodtimes_ds["gt_start_met"].values[i] - end = l1b_goodtimes_ds["gt_end_met"].values[i] + # All epochs are below threshold, so one goodtime interval spanning the full + # dataset is expected. GOODTIME_PADDING is subtracted from begin and added to end. + assert len(l1b_goodtimes_ds["gt_start_met"]) == 1 - # Start should be before end - assert start < end + raw_begin = met_times[0] + raw_end = met_times[-1] - # The difference should be reasonable (not negative due to offset) - assert (end - start) > 0 + assert l1b_goodtimes_ds["gt_start_met"].values[0] <= raw_begin + assert l1b_goodtimes_ds["gt_end_met"].values[0] >= raw_end def test_l1b_bgrates_and_goodtimes_rate_transition_low_to_high(attr_mgr_l1b): @@ -2800,12 +2720,14 @@ def test_l1b_bgrates_and_goodtimes_rate_transition_low_to_high(attr_mgr_l1b): }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2866,12 +2788,14 @@ def test_l1b_bgrates_and_goodtimes_rate_transition_high_to_low_to_high(attr_mgr_ }, ) - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} + sci_dependencies = { + "imap_lo_l1b_histrates": l1b_histrates, + "imap_lo_l1b_de": xr.Dataset(), + "imap_lo_l1b_nhk": xr.Dataset(), + } # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=10, delay_max=840 - ) + result = l1b_bgrates_and_goodtimes(sci_dependencies, attr_mgr_l1b, delay_max=840) # Assert l1b_bgrates_ds, l1b_goodtimes_ds = result @@ -2889,101 +2813,3 @@ def test_l1b_bgrates_and_goodtimes_rate_transition_high_to_low_to_high(attr_mgr_ # Background rates should be positive for all intervals assert np.all(l1b_bgrates_ds["h_background_rates"].values > 0) assert np.all(l1b_bgrates_ds["o_background_rates"].values > 0) - - -def test_l1b_bgrates_and_goodtimes_large_interval_with_active_tracking(attr_mgr_l1b): - """ - Test that an active goodtime interval is properly closed when a large interval - gap is encountered. - - This test ensures the code path where: - 1. We're actively tracking an interval (begin > 0.0) - 2. A chunk with interval > (interval_nom + delay_max) is encountered - 3. The active interval is closed before skipping the gap - """ - # Arrange - Create dataset where we start tracking, then hit a large interval - cycle_count = 10 - delay_max = 840 - met_spacing = 42 - - # First: Create enough low-rate chunks to start tracking (begin > 0.0) - num_chunks_before_gap = 2 # 2 chunks of 10 epochs each = 20 epochs - epochs_per_chunk = 10 - num_epochs_first = num_chunks_before_gap * epochs_per_chunk - - met_start = 473389200 - met_times_first = np.arange( - met_start, met_start + num_epochs_first * met_spacing, met_spacing - ) - - # Second: Create a chunk where the interval is too large - # The interval is measured from the first epoch of the chunk to the last - # We need interval > (interval_nom + delay_max) = 4200 + 840 = 5040 - large_gap = 6000 # Larger than threshold - met_times_gap_chunk_start = met_times_first[-1] + met_spacing - - # Create the problematic chunk (10 more epochs) - met_times_gap_chunk = np.arange( - met_times_gap_chunk_start, - met_times_gap_chunk_start + epochs_per_chunk * met_spacing, - met_spacing, - ) - - met_times_gap_chunk_adjusted = met_times_gap_chunk.copy() - met_times_gap_chunk_adjusted[-1] = met_times_gap_chunk[0] + large_gap - - # Third: Add more normal data after the gap - met_times_after = np.arange( - met_times_gap_chunk_adjusted[-1] + met_spacing, - met_times_gap_chunk_adjusted[-1] + met_spacing + 200 * met_spacing, - met_spacing, - ) - - met_times = np.concatenate( - [met_times_first, met_times_gap_chunk_adjusted, met_times_after] - ) - epoch_times = met_to_ttj2000ns(met_times) - - # All counts are low (below h_bg_rate_nom = 0.0028) to ensure we start tracking - h_counts = np.ones((len(met_times), 7, 60)) * 0.00025 - o_counts = np.ones((len(met_times), 7, 60)) * 0.000025 - - l1b_histrates = xr.Dataset( - { - "h_counts": (("epoch", "esa_step", "spin_bin_6"), h_counts), - "o_counts": (("epoch", "esa_step", "spin_bin_6"), o_counts), - }, - coords={ - "epoch": epoch_times, - "esa_step": np.arange(1, 8), - "spin_bin_6": np.arange(60), - }, - ) - - sci_dependencies = {"imap_lo_l1b_histrates": l1b_histrates} - - # Act - result = l1b_bgrates_and_goodtimes( - sci_dependencies, attr_mgr_l1b, cycle_count=cycle_count, delay_max=delay_max - ) - - # Assert - l1b_bgrates_ds, l1b_goodtimes_ds = result - - # Should have created at least 2 intervals: - # 1. The interval that was closed before the gap - # 2. The interval after the gap - assert len(l1b_goodtimes_ds["gt_start_met"]) >= 2 - - # The first interval should end before the gap chunk - # (it should be closed when we detect the large interval) - first_interval_end = l1b_goodtimes_ds["gt_end_met"].values[0] - gap_chunk_start = met_times_gap_chunk_adjusted[0] - - # The first interval should end before the gap chunk starts - # (with the +320 offset applied in the code) - assert first_interval_end < gap_chunk_start + 320 - - # Verify background rates are valid - assert np.all(l1b_bgrates_ds["h_background_rates"].values > 0) - assert np.all(l1b_bgrates_ds["o_background_rates"].values > 0)