diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 902eeefad1..72fd96a683 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,7 +3,7 @@ version: 2 build: os: "ubuntu-22.04" tools: - python: "3.10" + python: "3.12" sphinx: configuration: docs/source/conf.py diff --git a/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml b/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml index 8d5a108237..ec6ef36a2b 100644 --- a/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_enamaps_l2-common_variable_attrs.yaml @@ -239,6 +239,30 @@ bg_rate_sys_err: VAR_NOTES: Systematic error in the background count rate from non-ENA (non-heliospheric) sources. VAR_TYPE: support_data +bg_rate_sys_err_minus: + <<: *default_float32 + DEPEND_0: epoch + DICT_KEY: SPASE>Particle>ParticleType:Atom,ParticleQuantity:NumberFlux,Qualifier:Uncertainty,CoordinateSystemName:HAE,CoordinateRepresentation:Spherical + DISPLAY_TYPE: map_image + FIELDNAM: Background Rate Sys Error Lower + FORMAT: F6.3 + LABLAXIS: Rate Error Lower + UNITS: count s-1 + VAR_NOTES: Lower bound of the systematic error in the background count rate from non-ENA (non-heliospheric) sources. + VAR_TYPE: support_data + +bg_rate_sys_err_plus: + <<: *default_float32 + DEPEND_0: epoch + DICT_KEY: SPASE>Particle>ParticleType:Atom,ParticleQuantity:NumberFlux,Qualifier:Uncertainty,CoordinateSystemName:HAE,CoordinateRepresentation:Spherical + DISPLAY_TYPE: map_image + FIELDNAM: Background Rate Sys Error Upper + FORMAT: F6.3 + LABLAXIS: Rate Error Upper + UNITS: count s-1 + VAR_NOTES: Upper bound of the systematic error in the background count rate from non-ENA (non-heliospheric) sources. + VAR_TYPE: support_data + bg_intensity: <<: *default_float32 DEPEND_0: epoch @@ -275,6 +299,30 @@ bg_intensity_sys_err: VAR_NOTES: Non-statistical error in the background intensity from non-ENA (non-heliospheric) sources. VAR_TYPE: support_data +bg_intensity_sys_err_minus: + <<: *default_float32 + DEPEND_0: epoch + DICT_KEY: SPASE>Particle>ParticleType:Atom,ParticleQuantity:NumberFlux,Qualifier:Uncertainty,CoordinateSystemName:HAE,CoordinateRepresentation:Spherical + DISPLAY_TYPE: map_image + FIELDNAM: Background Intensity Non-statistical Error Lower + FORMAT: F8.3 + LABLAXIS: Intensity Non-statistical Error Lower + UNITS: cm -2 s -1 sr -1 keV -1 + VAR_NOTES: Lower bound of the non-statistical error in the background intensity from non-ENA (non-heliospheric) sources. + VAR_TYPE: support_data + +bg_intensity_sys_err_plus: + <<: *default_float32 + DEPEND_0: epoch + DICT_KEY: SPASE>Particle>ParticleType:Atom,ParticleQuantity:NumberFlux,Qualifier:Uncertainty,CoordinateSystemName:HAE,CoordinateRepresentation:Spherical + DISPLAY_TYPE: map_image + FIELDNAM: Background Intensity Non-statistical Error Upper + FORMAT: F8.3 + LABLAXIS: Intensity Non-statistical Error Upper + UNITS: cm -2 s -1 sr -1 keV -1 + VAR_NOTES: Upper bound of the non-statistical error in the background intensity from non-ENA (non-heliospheric) sources. + VAR_TYPE: support_data + ena_count: <<: *default_float32 DEPEND_0: epoch diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py index 4eddbafda9..191239d2cf 100644 --- a/imap_processing/ena_maps/utils/naming.py +++ b/imap_processing/ena_maps/utils/naming.py @@ -203,9 +203,13 @@ def build_map_var_catdesc(self, support_var_name: str) -> str | None: "bg_intensity": "Background Inten", "bg_intensity_stat_uncert": "Background Inten Stat. Unc.", "bg_intensity_sys_err": "Background Inten Sys. Err.", + "bg_intensity_sys_err_minus": "Background Inten Sys. Err. Lower", + "bg_intensity_sys_err_plus": "Background Inten Sys. Err. Upper", "bg_rate": "Background Count Rate", "bg_rate_stat_uncert": "Background Count Rate Stat. Unc.", "bg_rate_sys_err": "Background Count Rate Sys. Err.", + "bg_rate_sys_err_minus": "Background Count Rate Sys. Err. Lower", + "bg_rate_sys_err_plus": "Background Count Rate Sys. Err. Upper", "ena_count_rate": "ENA Count Rate", "ena_count_rate_stat_uncert": "ENA Count Rate Stat. Unc.", "obs_date": "Mean Observation Date", diff --git a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv new file mode 100644 index 0000000000..bcfc4cf1d9 --- /dev/null +++ b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv @@ -0,0 +1,16 @@ +esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_H,GF_Trpl_H_unc_minus,GF_Trpl_H_unc_plus +# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV] +0,1,1,0.01633,0.28,4.45E-05,4.03E-05,4.20E-05 +0,2,2,0.03047,0.43,5.02E-05,4.53E-05,4.72E-05 +0,3,3,0.05576,0.89,6.16E-05,5.58E-05,5.82E-05 +0,4,4,0.10626,1.9,7.12E-05,4.51E-05,4.89E-05 +0,5,5,0.20004,3.4,8.89E-05,5.85E-05,6.31E-05 +0,6,6,0.40496,7.3,1.12E-04,6.58E-05,7.23E-05 +0,7,7,0.78729,22,1.43E-04,8.25E-05,9.09E-05 +1,1,1,0.01719,0.22,1.05E-04,8.82E-05,9.26E-05 +1,2,2,0.03236,0.36,1.22E-04,1.02E-04,1.07E-04 +1,3,3,0.05948,0.77,1.44E-04,1.21E-04,1.27E-04 +1,4,4,0.11441,1.3,1.73E-04,1.17E-04,1.26E-04 +1,5,5,0.2137,3,2.15E-04,1.37E-04,1.49E-04 +1,6,6,0.43736,5.3,2.82E-04,1.65E-04,1.81E-04 +1,7,7,0.83888,11.7,3.61E-04,2.08E-04,2.29E-04 diff --git a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv new file mode 100644 index 0000000000..f494346c19 --- /dev/null +++ b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv @@ -0,0 +1,16 @@ +esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_O,GF_Trpl_O_unc_minus,GF_Trpl_O_unc_plus +# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV] +0,1,1,0.01919,1.372344105,1.98E-05,1.74E-05,1.81E-05 +0,2,2,0.03675,2.537963082,3.18E-05,2.39E-05,2.53E-05 +0,3,3,0.07121,6.591339559,5.34E-05,4.05E-05,4.29E-05 +0,4,4,0.14141,17.43990907,8.64E-05,6.14E-05,6.56E-05 +0,5,5,0.274,35.83900763,1.32E-04,9.07E-05,9.72E-05 +0,6,6,0.58503,98.06681395,2.46E-04,1.53E-04,1.67E-04 +0,7,7,1.13506,238.1076888,3.81E-04,2.23E-04,2.45E-04 +1,1,1,0.02043303552,1.575924599,5.02E-05,4.41E-05,4.61E-05 +1,2,2,0.03972,2.953020557,8.13E-05,6.09E-05,6.47E-05 +1,3,3,0.07648,7.476494389,1.33E-04,1.01E-04,1.07E-04 +1,4,4,0.15353,14.87953071,2.38E-04,1.69E-04,1.80E-04 +1,5,5,0.29846,41.48950349,3.87E-04,2.67E-04,2.86E-04 +1,6,6,0.61524,83.57629594,6.61E-04,4.11E-04,4.47E-04 +1,7,7,1.24282,212.0806857,1.02E-03,6.10E-04,6.67E-04 diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 25cea7b722..aa0e6e4a68 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -58,9 +58,9 @@ class LoConstants: # The first matching open interval (low < pivot < high) is used; if none matches, # THRESHOLD_BG_RATE_RAM_DEFAULT / THRESHOLD_BG_RATE_ANTI_RAM_DEFAULT apply. PIVOT_ANGLE_THRESHOLDS: ClassVar[dict[tuple[float, float], tuple[float, float]]] = { - (88.0, 92.0): (0.014, 0.007), - (73.0, 77.0): (0.0175, 0.00875), - (103.0, 107.0): (0.0112, 0.0056), + (88.0, 92.0): (0.028, 0.014), + (73.0, 77.0): (0.035, 0.0175), + (103.0, 107.0): (0.0224, 0.0112), } # Default background-rate thresholds [counts/s] when no pivot range matches. @@ -75,3 +75,17 @@ class LoConstants: # Padding [s] added to begin/end of each goodtime interval to ensure complete # cycles are covered at interval edges. GOODTIME_PADDING: float = 2.0 + + # Star-sensor spin-angle binning offset (fractional bin-index shift used when + # computing sample centers), keyed by the IFB star-sync housekeeping state + # (ifb_ctrl_star_sync). Flight software 4.8 enabled star sync ("EN"), + # switching from binning to the bin center (+0.5) to the left edge (+0.0). + STAR_BIN_OFFSET_BY_SYNC: ClassVar[dict[str | None, float]] = { + "DS": 0.5, # star sync disabled (pre FSW 4.8) + "EN": 0.0, # star sync enabled (FSW 4.8+) + } + + # Number of ending bins to exclude from each star-sensor profile average. + STAR_END_BINS_TO_EXCLUDE: int = 2 + # Minimum COUNT value for a star-sensor record to be considered valid. + STAR_MIN_COUNT_THRESHOLD: int = 700 diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index a085a5b361..f375388ad7 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -1986,7 +1986,7 @@ def split_rate_dataset( def filter_valid_star_records( l1a_star: xr.Dataset, - min_count: int = 700, + min_count: int = c.STAR_MIN_COUNT_THRESHOLD, time_window_offset: float = 0.0, time_window_duration: float | None = None, ) -> np.ndarray: @@ -2014,7 +2014,7 @@ def filter_valid_star_records( valid_mask : np.ndarray Boolean array indicating valid records. """ - # Section 5: Acceptance Criteria - COUNT >= 700 + # Section 5: Acceptance Criteria - COUNT >= min_count count_mask = l1a_star["count"].values >= min_count # shcoarse is already in MET seconds @@ -2049,7 +2049,7 @@ def filter_valid_star_records( def calculate_star_sensor_profile_for_group( data: np.ndarray, counts: np.ndarray, - end_bins_to_exclude: int = 2, + end_bins_to_exclude: int = c.STAR_END_BINS_TO_EXCLUDE, ) -> tuple[np.ndarray, np.ndarray]: """ Calculate averaged star sensor amplitude profile for a group of records. @@ -2098,14 +2098,63 @@ def calculate_star_sensor_profile_for_group( return avg_amplitude, count_array +def get_star_bin_offset(l1b_nhk: xr.Dataset, reference_epoch: int) -> float: + """ + Determine the star-sensor binning offset from the IFB star-sync state. + + Reads ``ifb_ctrl_star_sync`` from the NHK housekeeping at the record nearest + on or before ``reference_epoch`` and maps it to a binning offset via + ``LoConstants.STAR_BIN_OFFSET_BY_SYNC``. + + Parameters + ---------- + l1b_nhk : xr.Dataset + L1B NHK dataset containing ``ifb_ctrl_star_sync`` and an ``epoch`` + coordinate (TT2000 nanoseconds since J2000). + reference_epoch : int + Epoch at which to evaluate the sync state, in TT2000 nanoseconds since + J2000. The NHK record in effect at or before this time is used. + + Returns + ------- + bin_offset : float + Fractional bin-index offset to use when computing sample spin-angle + centers. + + Raises + ------ + KeyError + If ``ifb_ctrl_star_sync`` is not present in ``l1b_nhk``. + """ + if "ifb_ctrl_star_sync" not in l1b_nhk: + raise KeyError( + "ifb_ctrl_star_sync field not found in L1B NHK dataset. " + "Cannot determine star-sensor binning offset." + ) + + nhk_epoch = l1b_nhk["epoch"].values + sync_state = l1b_nhk["ifb_ctrl_star_sync"].values + + # Use the housekeeping record in effect at the reference epoch (the last + # NHK sample at or before it), clamping to the first sample if the reference + # epoch falls before NHK coverage. + idx = max(int(np.searchsorted(nhk_epoch, reference_epoch, side="right")) - 1, 0) + state = str(sync_state[idx]) + + offset = c.STAR_BIN_OFFSET_BY_SYNC[state] + logger.info(f"Star sync state '{state}' -> bin offset {offset}") + return offset + + def calculate_star_sensor_profiles_by_group( l1a_star: xr.Dataset, sampling_cadence: float, spin_period: float, group_size: int = 64, start_angle_offset: float = 62.0, - end_bins_to_exclude: int = 2, - min_count_threshold: int = 700, + end_bins_to_exclude: int = c.STAR_END_BINS_TO_EXCLUDE, + min_count_threshold: int = c.STAR_MIN_COUNT_THRESHOLD, + bin_offset: float = 0.5, ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Calculate averaged star sensor amplitude profiles for groups of records. @@ -2129,6 +2178,10 @@ def calculate_star_sensor_profiles_by_group( Number of ending bins to exclude from each average (default: 2). min_count_threshold : int Minimum COUNT value for valid record (default: 700). + bin_offset : float + Fractional offset applied to bin indices when computing sample + spin-angle centers (default: 0.5). Use 0.5 to bin to the bin center + and 0.0 to bin to the left edge. Returns ------- @@ -2152,7 +2205,7 @@ def calculate_star_sensor_profiles_by_group( # Calculate spin angles (same for all groups) deg_per_bin = 360.0 * (sampling_cadence / 1000.0) / spin_period bin_indices = np.arange(720) - sample_centers = (bin_indices + 0.5) * deg_per_bin + sample_centers = (bin_indices + bin_offset) * deg_per_bin spin_angle = (start_angle_offset + sample_centers) % 360.0 if n_valid == 0: @@ -2282,13 +2335,20 @@ def l1b_star( logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") # TODO: Read from ancillary config file when available - lo_angle_offset = 2.0 - sc_to_inst_angle_offset = ( - 360 * get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) - + lo_angle_offset + sc_to_inst_angle_offset = 360 * get_spacecraft_to_instrument_spin_phase_offset( + SpiceFrame.IMAP_LO ) - end_bins_to_exclude = 2 - min_count_threshold = 700 + end_bins_to_exclude = c.STAR_END_BINS_TO_EXCLUDE + min_count_threshold = c.STAR_MIN_COUNT_THRESHOLD + + # Global epoch times from L1A data (used for start_doy/end_doy below). + global_start_epoch = l1a_star["epoch"].values[0] + global_end_epoch = l1a_star["epoch"].values[-1] + + # Select the star-sensor binning convention from the IFB star-sync state in + # housekeeping. Evaluate at the earliest star record's epoch so a pointing that + # spans the `EN` event uses the value corresponding to the state at its start. + bin_offset = get_star_bin_offset(l1b_nhk, int(global_start_epoch)) # Calculate profiles for each 64-spin group ( @@ -2304,12 +2364,9 @@ def l1b_star( start_angle_offset=sc_to_inst_angle_offset, end_bins_to_exclude=end_bins_to_exclude, min_count_threshold=min_count_threshold, + bin_offset=bin_offset, ) - # Get global epoch times from L1A data for start_doy and end_doy - global_start_epoch = l1a_star["epoch"].values[0] - global_end_epoch = l1a_star["epoch"].values[-1] - # Create dataset with spin_angle as coordinate and multiple epochs group_epochs = met_to_ttj2000ns(group_mets) l1b_star_ds = xr.Dataset( @@ -2375,7 +2432,6 @@ def l1b_star( l1b_star_ds.attrs["pointing_mid_met"] = pointing_mid_met l1b_star_ds.attrs["sampling_cadence_ms"] = sampling_cadence l1b_star_ds.attrs["spin_duration_sec"] = spin_duration - l1b_star_ds.attrs["lo_angle_offset_deg"] = lo_angle_offset l1b_star_ds.attrs["end_bins_excluded"] = end_bins_to_exclude l1b_star_ds.attrs["min_count_threshold"] = min_count_threshold l1b_star_ds.attrs["group_size"] = group_size diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index b6af437d98..2160a3eddd 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -710,9 +710,9 @@ def load_geometric_factor_data(species: str) -> pd.DataFrame: anc_path = Path(__file__).parent.parent / "ancillary_data" if species == "h": - gf_file = anc_path / "imap_lo_hydrogen-geometric-factor_v001.csv" + gf_file = sorted(anc_path.glob("*hydrogen-geometric-factor*"))[-1] else: # species == "o" - gf_file = anc_path / "imap_lo_oxygen-geometric-factor_v001.csv" + gf_file = sorted(anc_path.glob("*oxygen-geometric-factor*"))[-1] return lo_ancillary.read_ancillary_file(gf_file) @@ -776,7 +776,8 @@ def initialize_geometric_factor_variables( "energy_delta_minus", "energy_delta_plus", "geometric_factor", - "geometric_factor_stat_uncert", + "geometric_factor_stat_uncert_minus", + "geometric_factor_stat_uncert_plus", ] # Initialize variables with proper dimensions (energy only) @@ -817,7 +818,8 @@ def populate_geometric_factors( gf_coords = {"energy": "Cntr_E"} gf_vars = { "geometric_factor": f"GF_Trpl_{species.upper()}", - "geometric_factor_stat_uncert": f"GF_Trpl_{species.upper()}_unc", + "geometric_factor_stat_uncert_minus": f"GF_Trpl_{species.upper()}_unc_minus", + "geometric_factor_stat_uncert_plus": f"GF_Trpl_{species.upper()}_unc_plus", } if species == "h": # NOTE: From an e-mail from Nathan on 2025-09-11 (values converted to keV) @@ -1007,11 +1009,16 @@ def calculate_intensities(dataset: xr.Dataset) -> xr.Dataset: dataset["counts_over_eff_squared"] ) / (dataset["geometric_factor"] * dataset["energy"] * dataset["exposure_factor"]) - # Equation 5 from mapping document (systematic uncertainty) - dataset["ena_intensity_sys_err"] = ( - dataset["ena_intensity"] - * dataset["geometric_factor_stat_uncert"] - / dataset["geometric_factor"] + for suffix in ("minus", "plus"): + dataset[f"ena_intensity_sys_err_{suffix}"] = ( + dataset["ena_intensity"] + * dataset[f"geometric_factor_stat_uncert_{suffix}"] + / dataset["geometric_factor"] + ) + + # Symmetric systematic error (mean of the asymmetric minus/plus bounds) + dataset["ena_intensity_sys_err"] = 0.5 * ( + dataset["ena_intensity_sys_err_minus"] + dataset["ena_intensity_sys_err_plus"] ) return dataset @@ -1040,11 +1047,17 @@ def calculate_backgrounds(dataset: xr.Dataset) -> xr.Dataset: dataset["bg_rate_stat_uncert_exposure_factor2"] / dataset["exposure_factor"] ** 2 ) - # Equation 8 from mapping document (background systematic uncertainty) - dataset["bg_rate_sys_err"] = ( - dataset["bg_rate"] - * dataset["geometric_factor_stat_uncert"] - / dataset["geometric_factor"] + + for suffix in ("minus", "plus"): + dataset[f"bg_rate_sys_err_{suffix}"] = ( + dataset["bg_rate"] + * dataset[f"geometric_factor_stat_uncert_{suffix}"] + / dataset["geometric_factor"] + ) + + # Symmetric systematic error (mean of the asymmetric minus/plus bounds) + dataset["bg_rate_sys_err"] = 0.5 * ( + dataset["bg_rate_sys_err_minus"] + dataset["bg_rate_sys_err_plus"] ) # Background intensity @@ -1054,10 +1067,17 @@ def calculate_backgrounds(dataset: xr.Dataset) -> xr.Dataset: dataset["bg_intensity_stat_uncert"] = dataset["bg_rate_stat_uncert"] / ( dataset["geometric_factor"] * dataset["energy"] ) - dataset["bg_intensity_sys_err"] = ( - dataset["bg_intensity"] - * dataset["geometric_factor_stat_uncert"] - / dataset["geometric_factor"] + + for suffix in ("minus", "plus"): + dataset[f"bg_intensity_sys_err_{suffix}"] = ( + dataset["bg_intensity"] + * dataset[f"geometric_factor_stat_uncert_{suffix}"] + / dataset["geometric_factor"] + ) + + # Symmetric systematic error (mean of the asymmetric minus/plus bounds) + dataset["bg_intensity_sys_err"] = 0.5 * ( + dataset["bg_intensity_sys_err_minus"] + dataset["bg_intensity_sys_err_plus"] ) return dataset @@ -1401,6 +1421,8 @@ def cleanup_intermediate_variables(dataset: xr.Dataset) -> xr.Dataset: potential_vars = [ "geometric_factor", "geometric_factor_stat_uncert", + "geometric_factor_stat_uncert_minus", + "geometric_factor_stat_uncert_plus", "counts_over_eff", "counts_over_eff_squared", "bg_rate_exposure_factor", diff --git a/imap_processing/lo/lo_ancillary.py b/imap_processing/lo/lo_ancillary.py index 994b03f30c..0a90ddc174 100644 --- a/imap_processing/lo/lo_ancillary.py +++ b/imap_processing/lo/lo_ancillary.py @@ -1,6 +1,7 @@ """Ancillary file reading for IMAP-Lo processing.""" from pathlib import Path +from typing import Any import pandas as pd @@ -36,20 +37,27 @@ def read_ancillary_file(ancillary_file: str | Path) -> pd.DataFrame: pd.DataFrame DataFrame containing the ancillary data. """ - skiprows = None + legacy_format = False + read_csv_kwargs: dict[str, Any] = {} if "esa-mode-lut" in str(ancillary_file): # skip the first row which is a comment - skiprows = [0] + read_csv_kwargs["skiprows"] = [0] elif "geometric-factor" in str(ancillary_file): - # skip the rows with comment headers indicating Hi_Res and Hi_Thr - skiprows = [1, 38] - df = pd.read_csv(ancillary_file, converters=_CONVERTERS, skiprows=skiprows) + # legacy format - rows with comment headers indicating Hi_Res and Hi_Thr + legacy_format = "Hi_Thr,,," in Path(ancillary_file).read_text() + if legacy_format: + read_csv_kwargs["skiprows"] = [1, 38] + else: + read_csv_kwargs["comment"] = "#" + df = pd.read_csv(ancillary_file, converters=_CONVERTERS, **read_csv_kwargs) df = df.rename(columns=_RENAME_COLUMNS) if "geometric-factor" in str(ancillary_file): - # Add an ESA mode column based on the known structure of the file. - # The first 36 rows are ESA mode 0 (HiRes), the second 36 are ESA mode 1 (HiThr) - df["esa_mode"] = 0 - df.loc[36:, "esa_mode"] = 1 + if legacy_format and "esa_mode" not in df.columns: + # Add an ESA mode column based on the known structure of the file. + # The first 36 rows are ESA mode 0 (HiRes), the second 36 are ESA mode 1 + # (HiThr) + df["esa_mode"] = 0 + df.loc[36:, "esa_mode"] = 1 return df diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 17d0ae3db7..4291dbdd80 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -27,6 +27,7 @@ get_pivot_angle_from_nhk, get_sampling_cadence_from_nhk, get_spin_start_times, + get_star_bin_offset, identify_species, initialize_l1b_de, l1b_bgrates_and_goodtimes, @@ -1836,6 +1837,58 @@ def test_profiles_by_group_angle_wrapping(self, mock_repoint): assert np.any(spin_angle < 100) # Some angles wrapped to lower range +class TestStarBinOffset: + """Tests for the star-sensor binning offset derived from HK star-sync state.""" + + @staticmethod + def _nhk(states): + """Build a minimal NHK dataset with ifb_ctrl_star_sync at 1 Hz.""" + n = len(states) + return xr.Dataset( + {"ifb_ctrl_star_sync": ("epoch", list(states))}, + coords={"epoch": met_to_ttj2000ns(np.arange(n, dtype=np.float64))}, + ) + + @pytest.mark.parametrize("state, expected", [("DS", 0.5), ("EN", 0.0)]) + def test_offset_maps_sync_state(self, state, expected): + """DS (sync disabled) -> 0.5; EN (sync enabled) -> 0.0.""" + nhk = self._nhk([state] * 3) + star_end = int(nhk["epoch"].values[1]) + assert get_star_bin_offset(nhk, star_end) == expected + + def test_uses_state_at_or_before_reference(self): + """Offset reflects the NHK record in effect at the reference epoch.""" + # Sync flips DS->EN at index 2. + nhk = self._nhk(["DS", "DS", "EN", "EN"]) + before = int(nhk["epoch"].values[1]) + at = int(nhk["epoch"].values[2]) + assert get_star_bin_offset(nhk, before) == 0.5 + assert get_star_bin_offset(nhk, at) == 0.0 + + def test_straddling_pointing_uses_start_state(self): + """A pointing that begins before the enable event uses DS (bin center).""" + # Star data starts under DS but ends under EN: evaluating at the start + # (as l1b_star does) yields DS. + nhk = self._nhk(["DS", "DS", "EN", "EN"]) + star_start = int(nhk["epoch"].values[0]) + assert get_star_bin_offset(nhk, star_start) == 0.5 + + def test_missing_field_raises(self): + """A clear error is raised if the star-sync field is absent.""" + nhk = xr.Dataset( + {"ifb_data_interval": ("epoch", [21.0])}, + coords={"epoch": met_to_ttj2000ns([0.0])}, + ) + with pytest.raises(KeyError, match="ifb_ctrl_star_sync"): + get_star_bin_offset(nhk, int(nhk["epoch"].values[0])) + + def test_unexpected_state(self): + """An unrecognized star-sync state raises a KeyError.""" + nhk = self._nhk(["??"]) + with pytest.raises(KeyError): + assert get_star_bin_offset(nhk, int(nhk["epoch"].values[0])) == 0.5 + + class TestL1bStar: """Tests for l1b_star function.""" @@ -1873,6 +1926,7 @@ def test_initializes_with_spin_data( l1b_nhk = xr.Dataset( { "ifb_data_interval": ("epoch", [21.0] * n_records), + "ifb_ctrl_star_sync": ("epoch", ["DS"] * n_records), }, coords={"epoch": list(range(n_records))}, ) @@ -1953,6 +2007,7 @@ def test_dataset_structure_and_attributes( l1b_nhk = xr.Dataset( { "ifb_data_interval": ("epoch", [21.0]), + "ifb_ctrl_star_sync": ("epoch", ["DS"]), }, coords={"epoch": [0]}, ) @@ -1998,10 +2053,8 @@ def test_dataset_structure_and_attributes( assert l1b_star_ds["count_per_bin"].attrs["VALIDMAX"] == 100000 # Assert - Check processing parameter attributes - assert "lo_angle_offset_deg" in l1b_star_ds.attrs assert "end_bins_excluded" in l1b_star_ds.attrs assert "min_count_threshold" in l1b_star_ds.attrs - assert l1b_star_ds.attrs["lo_angle_offset_deg"] == 2.0 assert l1b_star_ds.attrs["end_bins_excluded"] == 2 assert l1b_star_ds.attrs["min_count_threshold"] == 700 @@ -2035,6 +2088,7 @@ def test_start_and_end_doy_variables( l1b_nhk = xr.Dataset( { "ifb_data_interval": ("epoch", [21.0, 21.0, 21.0]), + "ifb_ctrl_star_sync": ("epoch", ["DS", "DS", "DS"]), }, coords={"epoch": [0, 1, 2]}, ) @@ -2103,6 +2157,7 @@ def test_multiple_groups_created( l1b_nhk = xr.Dataset( { "ifb_data_interval": ("epoch", [21.0] * n_records), + "ifb_ctrl_star_sync": ("epoch", ["DS"] * n_records), }, coords={"epoch": list(range(n_records))}, ) diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 3422f20704..f4868791b8 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -278,7 +278,8 @@ def sample_geometric_factor_data(): "Cntr_E": 0.01 * (i + 1), # Simple energy values "Cntr_E_unc": 0.001 * (i + 1), "GF_Trpl_H": 1e-4 * (i + 1), - "GF_Trpl_H_unc": 1e-5 * (i + 1), + "GF_Trpl_H_unc_minus": 1e-5 * (i + 1), + "GF_Trpl_H_unc_plus": 2e-5 * (i + 1), "GF_Dbl_all": 2e-4 * (i + 1), "GF_Dbl_all_unc": 2e-5 * (i + 1), "GF_Trpl_all": 3e-4 * (i + 1), @@ -294,7 +295,8 @@ def sample_geometric_factor_data(): "Cntr_E": 0.015 * (i + 1), # Slightly different for oxygen "Cntr_E_unc": 0.0015 * (i + 1), "GF_Trpl_O": 1.5e-4 * (i + 1), - "GF_Trpl_O_unc": 1.5e-5 * (i + 1), + "GF_Trpl_O_unc_minus": 1.5e-5 * (i + 1), + "GF_Trpl_O_unc_plus": 3e-5 * (i + 1), } ) @@ -358,7 +360,8 @@ def sample_dataset_with_geometric_factors(): dataset["exposure_factor"] = (("epoch", "energy"), np.ones((1, 7)) * 1.0) dataset["geometric_factor"] = (("energy",), np.ones(7) * 1e-4) dataset["energy"] = (("energy",), np.ones(7) * 0.1) # Energy values - dataset["geometric_factor_stat_uncert"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_minus"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_plus"] = (("energy",), np.ones(7) * 3e-5) return dataset @@ -393,7 +396,14 @@ def sample_dataset_with_background_intermediates(): # Add geometric factors for systematic uncertainty calculation dataset["geometric_factor"] = (("energy",), np.ones(n_energy) * 1e-4) - dataset["geometric_factor_stat_uncert"] = (("energy",), np.ones(n_energy) * 1e-5) + dataset["geometric_factor_stat_uncert_minus"] = ( + ("energy",), + np.ones(n_energy) * 1e-5, + ) + dataset["geometric_factor_stat_uncert_plus"] = ( + ("energy",), + np.ones(n_energy) * 3e-5, + ) return dataset @@ -705,7 +715,8 @@ def test_reduce_geometric_factor_dataset_hydrogen( # Check that hydrogen-specific columns are present assert "Cntr_E" in result.data_vars assert "GF_Trpl_H" in result.data_vars - assert "GF_Trpl_H_unc" in result.data_vars + assert "GF_Trpl_H_unc_minus" in result.data_vars + assert "GF_Trpl_H_unc_plus" in result.data_vars # Verify energy values match expected for hydrogen expected_energies = [0.01 * (i + 1) for i in range(7)] @@ -724,7 +735,8 @@ def test_reduce_geometric_factor_dataset_oxygen( # Check that oxygen-specific columns are present assert "Cntr_E" in result.data_vars assert "GF_Trpl_O" in result.data_vars - assert "GF_Trpl_O_unc" in result.data_vars + assert "GF_Trpl_O_unc_minus" in result.data_vars + assert "GF_Trpl_O_unc_plus" in result.data_vars # Verify energy values match expected for oxygen expected_energies = [0.015 * (i + 1) for i in range(7)] @@ -854,9 +866,9 @@ def test_normalize_coordinates_basic(self, species): # Check that energy coordinate is present assert "energy" in result.coords expected_energies = ( - np.array([0.01633, 0.03047, 0.05576, 0.1063, 0.2, 0.405, 0.7873]) + np.array([0.01633, 0.03047, 0.05576, 0.10626, 0.20004, 0.40496, 0.78729]) if species == "h" - else np.array([0.016, 0.032, 0.065, 0.135, 0.279, 0.601, 1.206]) + else np.array([0.01919, 0.03675, 0.07121, 0.14141, 0.274, 0.58503, 1.13506]) ) np.testing.assert_array_equal(result.coords["energy"], expected_energies) @@ -1122,10 +1134,15 @@ def test_calculate_intensities_basic(self, sample_dataset_with_geometric_factors result["ena_intensity_stat_uncert"], expected_stat_uncert ) - # Check systematic uncertainty calculation + # Check systematic uncertainty calculation. The single `_sys_err` is the + # mean of the asymmetric minus/plus bounds. + mean_gf_stat_uncert = 0.5 * ( + sample_dataset_with_geometric_factors["geometric_factor_stat_uncert_minus"] + + sample_dataset_with_geometric_factors["geometric_factor_stat_uncert_plus"] + ) expected_sys_err = ( result["ena_intensity"] - * sample_dataset_with_geometric_factors["geometric_factor_stat_uncert"] + * mean_gf_stat_uncert / sample_dataset_with_geometric_factors["geometric_factor"] ) xr.testing.assert_allclose(result["ena_intensity_sys_err"], expected_sys_err) @@ -1180,11 +1197,13 @@ def test_calculate_backgrounds_basic( xr.testing.assert_allclose(result["bg_rate_stat_uncert"], expected_stat_uncert) # Check systematic uncertainty calculation - # (geometric_factor_stat_uncert / geometric_factor) * bg_rate + # (mean(geometric_factor_stat_uncert bounds) / geometric_factor) * bg_rate + mean_gf_stat_uncert = 0.5 * ( + dataset["geometric_factor_stat_uncert_minus"] + + dataset["geometric_factor_stat_uncert_plus"] + ) expected_sys_err = ( - result["bg_rate"] - * dataset["geometric_factor_stat_uncert"] - / dataset["geometric_factor"] + result["bg_rate"] * mean_gf_stat_uncert / dataset["geometric_factor"] ) xr.testing.assert_allclose(result["bg_rate_sys_err"], expected_sys_err) @@ -1202,7 +1221,8 @@ def test_calculate_backgrounds_zero_exposure(self): ), "exposure_factor": (("epoch", "energy"), np.zeros((1, 7))), "geometric_factor": (("energy",), np.ones(7) * 1e-4), - "geometric_factor_stat_uncert": (("energy",), np.ones(7) * 1e-5), + "geometric_factor_stat_uncert_minus": (("energy",), np.ones(7) * 1e-5), + "geometric_factor_stat_uncert_plus": (("energy",), np.ones(7) * 3e-5), }, coords={"epoch": [8.1794907049e17], "energy": list(range(7))}, ) @@ -1733,7 +1753,8 @@ def test_initialize_geometric_factor_variables(self): "energy_delta_minus", "energy_delta_plus", "geometric_factor", - "geometric_factor_stat_uncert", + "geometric_factor_stat_uncert_minus", + "geometric_factor_stat_uncert_plus", ] for var in expected_vars: @@ -1773,15 +1794,21 @@ def test_populate_geometric_factors( # Check hydrogen values assert result["energy"].values[i] == 0.01 * (i + 1) assert result["geometric_factor"].values[i] == 1e-4 * (i + 1) - assert result["geometric_factor_stat_uncert"].values[i] == ( + assert result["geometric_factor_stat_uncert_minus"].values[i] == ( 1e-5 * (i + 1) ) + assert result["geometric_factor_stat_uncert_plus"].values[i] == ( + 2e-5 * (i + 1) + ) else: # oxygen assert result["energy"].values[i] == 0.015 * (i + 1) assert result["geometric_factor"].values[i] == 1.5e-4 * (i + 1) - assert result["geometric_factor_stat_uncert"].values[i] == ( + assert result["geometric_factor_stat_uncert_minus"].values[i] == ( 1.5e-5 * (i + 1) ) + assert result["geometric_factor_stat_uncert_plus"].values[i] == ( + 3e-5 * (i + 1) + ) # Ensure that energy_deltas are in units of keV assert np.all(result["energy_delta_plus"].values < 1) assert np.all(result["energy_delta_minus"].values < 1) @@ -2292,7 +2319,8 @@ def test_calculate_all_rates_and_intensities_complete(self): "exposure_factor": (("energy",), np.ones(7) * 1.0), "geometric_factor": (("energy",), np.ones(7) * 1e-4), "energy": (("energy",), np.ones(7) * 0.1), - "geometric_factor_stat_uncert": (("energy",), np.ones(7) * 1e-5), + "geometric_factor_stat_uncert_minus": (("energy",), np.ones(7) * 1e-5), + "geometric_factor_stat_uncert_plus": (("energy",), np.ones(7) * 3e-5), # Background intermediate data "bg_rate_exposure_factor": (("energy",), np.ones(7) * 0.3), "bg_rate_stat_uncert_exposure_factor2": ( @@ -2335,7 +2363,8 @@ def test_calculate_all_rates_with_cg_correction( dataset["counts_over_eff_squared"] = (("epoch", "energy"), np.ones((1, 7)) * 12) dataset["exposure_factor"] = (("epoch", "energy"), np.ones((1, 7)) * 1.0) dataset["geometric_factor"] = (("energy",), np.ones(7) * 1e-4) - dataset["geometric_factor_stat_uncert"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_minus"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_plus"] = (("energy",), np.ones(7) * 3e-5) dataset["bg_rate_exposure_factor"] = ( ("epoch", "energy"), np.ones((1, 7)) * 0.3, @@ -2386,7 +2415,8 @@ def test_calculate_all_rates_cg_with_other_corrections( dataset["counts_over_eff_squared"] = (("epoch", "energy"), np.ones((1, 7)) * 12) dataset["exposure_factor"] = (("epoch", "energy"), np.ones((1, 7)) * 1.0) dataset["geometric_factor"] = (("energy",), np.ones(7) * 1e-4) - dataset["geometric_factor_stat_uncert"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_minus"] = (("energy",), np.ones(7) * 1e-5) + dataset["geometric_factor_stat_uncert_plus"] = (("energy",), np.ones(7) * 3e-5) dataset["bg_rate_exposure_factor"] = ( ("epoch", "energy"), np.ones((1, 7)) * 0.3,