diff --git a/imap_processing/tests/external_test_data_config.py b/imap_processing/tests/external_test_data_config.py index 8ab7c04af..d6295a213 100644 --- a/imap_processing/tests/external_test_data_config.py +++ b/imap_processing/tests/external_test_data_config.py @@ -153,6 +153,15 @@ ("validate_upstream_ion_1_culling_results_repoint00047_v1.csv", "ultra/data/l1/"), ("validate_spectral_culling_results_repoint00047_v1.csv", "ultra/data/l1/"), ("de_test_data_repoint00047.csv", "ultra/data/l1/"), + ("extendedspin_test_data_repoint00383.csv", "ultra/data/l1/"), + ("p1de_test_data_repoint00383.csv", "ultra/data/l1/"), + ("status_test_data_repoint00383.csv", "ultra/data/l1/"), + ("voltage_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), + ("validate_high_energy_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), + ("validate_stat_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), + ("validate_spectral_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), + ("validate_upstream_ion_1_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), + ("validate_upstream_ion_2_culling_results_repoint00383_v1.csv", "ultra/data/l1/"), ("FM45_UltraFM45Extra_TV_Tests_2024-01-22T0930_20240122T093008.CCSDS", "ultra/data/l0/"), ("ultra45_raw_sc_rawnrgevnt_19840122_00.csv", "ultra/data/l0/"), ("ultra45_raw_sc_enaphxtofhnrgimg_FM45_UltraFM45Extra_TV_Tests_2024-01-22T0930_20240122T093008.csv", diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py index 56a5bf456..bc6202b90 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1b_culling.py @@ -55,13 +55,21 @@ ) +# Get the ULTRA 90 extended spin config for this pointing +repoint = "repoint00383" +REPOINT_383_SPIN_CONFIG = ExtendedSpinConfig.from_csv( + TEST_PATH / "imap_ultra_l1b-90sensor-extendedspin-config_20251001_v001.csv", repoint +) + + @pytest.fixture -def setup_repoint_47_data(): - """Fixture to set up data for validation test using repoint 47.""" - de_df = pd.read_csv(TEST_PATH / "de_test_data_repoint00047.csv") +def setup_repoint_383_data(): + """Fixture to set up data for validation test using repoint 383 (ULTRA 90).""" + # The priority config for repoint 383 is p1, so only the priority 1 de is used. + de_df = pd.read_csv(TEST_PATH / "p1de_test_data_repoint00383.csv") de_ds = xr.Dataset( { - "de_event_met": ("epoch", de_df.event_times.values), + "de_event_met": ("epoch", de_df.de_event_met.values), "event_times": ("epoch", de_df.event_times.values), "energy_spacecraft": ("epoch", de_df.energy_spacecraft.values), "quality_outliers": ("epoch", de_df.quality_outliers.values), @@ -69,7 +77,7 @@ def setup_repoint_47_data(): "ebin": ("epoch", de_df.ebin.values), } ) - xspin = pd.read_csv(TEST_PATH / "extendedspin_test_data_repoint00047.csv") + xspin = pd.read_csv(TEST_PATH / "extendedspin_test_data_repoint00383.csv") spin_bin_size = UltraConstants.SPIN_BIN_SIZE spin_tbin_edges = get_binned_spins_edges( xspin.spin_number.values, @@ -87,7 +95,7 @@ def setup_repoint_47_data(): de_datasets, energy_ranges, spin_tbin_edges, - REPOINT_47_SPIN_CONFIG, + REPOINT_383_SPIN_CONFIG, 90, ) @@ -393,7 +401,8 @@ def test_flag_low_voltage(test_data): spin_tbin_edges = get_binned_spins_edges( spins, spin_period, spin_starttime, spin_bin_size ) - quality_flags = flag_low_voltage(spin_tbin_edges, mock_status_dataset) + v_threshold = REPOINT_383_SPIN_CONFIG.voltage_threshold + quality_flags = flag_low_voltage(spin_tbin_edges, mock_status_dataset, v_threshold) # There should be an extra bin edge for the last bin to indicate the end of the last # spin bin @@ -406,7 +415,7 @@ def test_flag_low_voltage(test_data): # Set only the first spin to be below threshold mock_status_dataset["rightdeflection_v"].data[1:] += 5000 mock_status_dataset["leftdeflection_v"].data[1:] += 5000 - quality_flags = flag_low_voltage(spin_tbin_edges, mock_status_dataset) + quality_flags = flag_low_voltage(spin_tbin_edges, mock_status_dataset, v_threshold) # Check that only the first spin is flagged for low voltage assert np.all(quality_flags[0]) # The rest should not be flagged @@ -432,7 +441,9 @@ def test_flag_low_voltage_incomplete_bins(test_data): spin_tbin_edges = get_binned_spins_edges( spins, spin_period, spin_starttime, spin_bin_size ) - quality_flags = flag_low_voltage(spin_tbin_edges, mock_status_dataset) + quality_flags = flag_low_voltage( + spin_tbin_edges, mock_status_dataset, REPOINT_383_SPIN_CONFIG.voltage_threshold + ) # check quality flag assert quality_flags.shape == (n_spins // spin_bin_size,) @@ -444,15 +455,16 @@ def test_expand_bin_flags_to_spins(caplog): spin_bin_size = 5 n_spins = 12 # Mock the shape of binned quality flags for 12 spins and a bin size of 5 - binned_qf = np.full((n_spins // spin_bin_size), 1) + binned_qf = np.array([1, 2]) quality_flags = expand_bin_flags_to_spins(n_spins, binned_qf, spin_bin_size) # Check the size assert quality_flags.shape == (n_spins,) - # The first 10 spins should be flagged since they fall into the first two bins - assert np.all(quality_flags[:10] == 1) - # The last 2 spins should not be flagged since they fall into the last incomplete - # bin - assert np.all(quality_flags[10:] == 0) + # The first 10 spins should have the flags of the two complete bins + assert np.all(quality_flags[:5] == 1) + assert np.all(quality_flags[5:10] == 2) + # The last 2 spins fall into the incomplete bin and should inherit the flag of + # the last complete bin + assert np.all(quality_flags[10:] == 2) binned_qf = np.full((n_spins // spin_bin_size) + 1, 1) # test that a warning is logged when there are incomplete bins found expand_bin_flags_to_spins(n_spins, binned_qf, spin_bin_size) @@ -511,13 +523,13 @@ def test_get_energy_and_spin_dependent_rejection_mask(): def test_validate_voltage_cull(): """Validate that low voltage spins are correctly flagged""" # read test data from csv files - xspin = pd.read_csv(TEST_PATH / "extendedspin_test_data_repoint00047.csv") + xspin = pd.read_csv(TEST_PATH / "extendedspin_test_data_repoint00383.csv") validation_low_voltage_qf = np.loadtxt( - TEST_PATH / "voltage_culling_results_repoint00047.csv", + TEST_PATH / "voltage_culling_results_repoint00383_v1.csv", delimiter=",", dtype=np.uint16, ) - status_df = pd.read_csv(TEST_PATH / "status_test_data_repoint00047.csv") + status_df = pd.read_csv(TEST_PATH / "status_test_data_repoint00383.csv") # build the status dataset including the variables needed for the low voltage flag status_ds = xr.Dataset( { @@ -528,14 +540,15 @@ def test_validate_voltage_cull(): ) # Use constants from the code to ensure consistency with the actual culling code spin_bin_size = UltraConstants.SPIN_BIN_SIZE - lv_threshold = 3000 spin_tbin_edges = get_binned_spins_edges( xspin.spin_number.values, xspin.spin_period.values, xspin.spin_start_time.values, spin_bin_size, ) - lv_flags = flag_low_voltage(spin_tbin_edges, status_ds, lv_threshold) + lv_flags = flag_low_voltage( + spin_tbin_edges, status_ds, REPOINT_47_SPIN_CONFIG.voltage_threshold + ) assert np.array_equal(lv_flags, validation_low_voltage_qf) @@ -569,7 +582,7 @@ def test_get_valid_earth_angle_events(mock_spkezr): vde = np.sqrt(np.sum(de_dps_velocity**2, 1)) uv = de_dps_velocity / vde[:, np.newaxis] local_uv = np.array([upos, yax, zax]) @ np.transpose(-uv) - expected_flags = local_uv[0, :] < np.cos(earth_angle_threshold) + expected_flags = np.abs(local_uv[0, :]) < np.cos(earth_angle_threshold) actual_flags = get_valid_earth_angle_events(de_dataset, earth_angle_threshold) np.testing.assert_array_equal(actual_flags, expected_flags) @@ -613,7 +626,7 @@ def test_get_valid_de_count_summary_valid_events(): {"p0": de_dataset, "p1": de_dataset}, energy_range_edges, spin_tbin_edges, - REPOINT_47_SPIN_CONFIG, + REPOINT_383_SPIN_CONFIG, 90, ) @@ -655,10 +668,9 @@ def test_get_valid_de_count_summary_ultra45(mock_spkezr): "ebin": ("epoch", np.full(len(energy), 10)), } ) - # ensure that all events fail the earth angle check by setting a very large # keepout angle - keepout_angle = np.radians(360) + keepout_angle = np.radians(90) counts = get_valid_de_count_summary( {"p0": de_dataset, "p1": de_dataset}, energy_range_edges, @@ -710,6 +722,7 @@ def test_flag_high_energy(): start=0, stop=len(energy) + 1, step=4 ) # create spin bins of 4 seconds de_datasets = {"p0": de_dataset, "p1": de_dataset} + de_counts_summary = get_valid_de_count_summary( de_datasets, energy_range_edges, @@ -721,8 +734,8 @@ def test_flag_high_energy(): de_counts_summary, spin_tbin_edges, energy_range_edges, - None, cull_thresholds, + None, combine_spin_bin_radius=0, ) @@ -744,34 +757,19 @@ def test_flag_high_energy(): assert not np.any(quality_flags[:, 2]) -@mock.patch( - "imap_processing.ultra.l1b.ultra_l1b_culling.UltraConstants.HIGH_ENERGY_CULL_CHANNEL", - 4, -) @pytest.mark.external_test_data -def test_validate_high_energy_cull(setup_repoint_47_data): +def test_validate_high_energy_cull(setup_repoint_383_data): """Validate that high energy spins are correctly flagged""" - # Mock thresholds to match the test data (I used fake ones to create more - # complexity) - mock_thresholds = np.array([0.05, 1.5, 0.6, 119.2, 0.2]) * 20 expected_qf = pd.read_csv( - TEST_PATH / "validate_high_energy_culling_results_repoint00047_v2.csv" + TEST_PATH / "validate_high_energy_culling_results_repoint00383_v1.csv" ).to_numpy() - de_datasets, _, spin_tbin_edges, _, _ = setup_repoint_47_data - - # Get the energy ranges - energy_ranges = np.array([4.2, 9.4425, 21.2116, 47.2388, 105.202, 316.335]) - de_counts_summary = get_valid_de_count_summary( - de_datasets, energy_ranges, spin_tbin_edges, REPOINT_47_SPIN_CONFIG - ) - high_energy_combined_spin_bin_radius = 3 + _, _, spin_tbin_edges, energy_ranges, de_counts_summary = setup_repoint_383_data e_flags = flag_high_energy( de_counts_summary, spin_tbin_edges, energy_ranges, + REPOINT_383_SPIN_CONFIG.energy_thresholds, None, - mock_thresholds, - combine_spin_bin_radius=high_energy_combined_spin_bin_radius, ) np.testing.assert_array_equal(e_flags, ~expected_qf.astype(bool)) @@ -805,6 +803,7 @@ def test_flag_statistical_outliers(): start=0, stop=len(energy) + 1, step=spin_step ) # create spin bins of 7 seconds de_datasets = {"p0": de_dataset, "p1": de_dataset} + de_counts_summary = get_valid_de_count_summary( de_datasets, energy_range_edges, @@ -907,23 +906,15 @@ def test_get_poisson_stats(): @pytest.mark.external_test_data -def test_validate_stat_cull(setup_repoint_47_data): +def test_validate_stat_cull(setup_repoint_383_data): """Validate that statistical-outlier quality flags match expected results.""" # read test data from csv files results_df = pd.read_csv( - TEST_PATH / "validate_stat_culling_results_repoint00047_v3.csv" - ) - de_datasets, _, spin_tbin_edges, energy_ranges, de_counts_summary = ( - setup_repoint_47_data + TEST_PATH / "validate_stat_culling_results_repoint00383_v1.csv" ) - - # Create a mask of flagged events to test that the stat cull algorithm - # properly ignores these. The test data was created using this exact mask as well. + _, _, spin_tbin_edges, energy_ranges, de_counts_summary = setup_repoint_383_data + # No spin bins flagged by previous steps mask = np.zeros((len(energy_ranges) - 1, len(spin_tbin_edges) - 1), dtype=bool) - mask[0:2, 0:2] = ( - True # This will mark the first 2 energy bins and first 2 spin bins as flagged - ) - # ignored in the statistics calculation and flagging. flags, con, it, std = flag_statistical_outliers( de_counts_summary, spin_tbin_edges, energy_ranges, mask ) @@ -960,38 +951,34 @@ def test_get_binned_energy_ranges(): @pytest.mark.external_test_data -def test_validate_upstream_ion_cull(setup_repoint_47_data): +@pytest.mark.parametrize( + "version, channels", + [ + (1, UltraConstants.UPSTREAM_ION_ENERGY_CHANNELS_1), + (2, UltraConstants.UPSTREAM_ION_ENERGY_CHANNELS_2), + ], +) +def test_validate_upstream_ion_cull(setup_repoint_383_data, version, channels): """Validate that upstream ion quality flags match expected results.""" # read test data from csv files expected_results = pd.read_csv( - TEST_PATH / "validate_upstream_ion_1_culling_results_repoint00047_v1.csv" + TEST_PATH + / f"validate_upstream_ion_{version}_culling_results_repoint00383_v1.csv" ).to_numpy() - de_datasets, _, spin_tbin_edges, energy_ranges, de_counts_summary = ( - setup_repoint_47_data - ) + _, _, spin_tbin_edges, energy_ranges, de_counts_summary = setup_repoint_383_data + # No spin bins flagged by previous steps mask = np.zeros((len(energy_ranges) - 1, len(spin_tbin_edges) - 1), dtype=bool) - mask[0:2, 0:2] = ( - True # This will mark the first 2 energy bins and first 2 spin bins as flagged - ) - - flags = flag_upstream_ion( - de_counts_summary, - energy_ranges, - mask, - UltraConstants.UPSTREAM_ION_ENERGY_CHANNELS_1, - ) - # Combine the flags with the mask to get the final expected results since the - # masked bins should be flagged as well. - results = flags | mask + flags = flag_upstream_ion(de_counts_summary, energy_ranges, mask, channels) + # The upstream ion flags are energy independent, so they apply to every energy + # range. + results = np.broadcast_to(flags, mask.shape) np.testing.assert_array_equal(results, ~expected_results.astype(bool)) @pytest.mark.external_test_data -def test_upstream_ion_cull_invalid_channels(setup_repoint_47_data): +def test_upstream_ion_cull_invalid_channels(setup_repoint_383_data): """Validate upstream ion error handling.""" - de_datasets, _, spin_tbin_edges, energy_ranges, de_counts_summary = ( - setup_repoint_47_data - ) + _, _, spin_tbin_edges, energy_ranges, de_counts_summary = setup_repoint_383_data mask = np.zeros((len(energy_ranges) - 1, len(spin_tbin_edges) - 1), dtype=bool) with pytest.raises( ValueError, @@ -1007,23 +994,18 @@ def test_upstream_ion_cull_invalid_channels(setup_repoint_47_data): @pytest.mark.external_test_data -def test_validate_spectral_cull(setup_repoint_47_data): +def test_validate_spectral_cull(setup_repoint_383_data): """Validate that spectral flags match expected results.""" # read test data from csv files expected_results = pd.read_csv( - TEST_PATH / "validate_spectral_culling_results_repoint00047_v1.csv" + TEST_PATH / "validate_spectral_culling_results_repoint00383_v1.csv" ).to_numpy() - de_datasets, xspin, spin_tbin_edges, energy_ranges, de_counts_summary = ( - setup_repoint_47_data - ) - mask = np.zeros((len(energy_ranges) - 1, len(spin_tbin_edges) - 1), dtype=bool) - mask[0:2, 0:2] = ( - True # This will mark the first 2 energy bins and first 2 spin bins as flagged - ) + _, _, _, energy_ranges, de_counts_summary = setup_repoint_383_data flags = flag_spectral_events( de_counts_summary, energy_ranges, UltraConstants.SPECTRAL_ENERGY_CHANNELS, ) - results = flags | mask + # The spectral flags are energy independent, so they apply to every energy range. + results = np.broadcast_to(flags, expected_results.shape) np.testing.assert_array_equal(results, ~expected_results.astype(bool)) diff --git a/imap_processing/ultra/constants.py b/imap_processing/ultra/constants.py index af6814e0d..925113349 100644 --- a/imap_processing/ultra/constants.py +++ b/imap_processing/ultra/constants.py @@ -207,13 +207,6 @@ class UltraConstants: # Angle threshold in radians for ULTRA 45 degree culling. # This is only needed for ULTRA 45 since Earth may be in the FOV. EARTH_ANGLE_45_THRESHOLD = np.radians(15) - # An array of energy thresholds to use for culling. Each one corresponds to - # the number of energy bins used. - # n_bins=len(PSET_ENERGY_BIN_EDGES)[BASE_CULL_EBIN:] // N_CULL_EBINS - # an error will be raised if this does not match n_bins - HIGH_ENERGY_CULL_THRESHOLDS = ( - np.array([4.0, 2.0, 1.20, 0.45, 0.1, 0.1]) * SPIN_BIN_SIZE - ) # Use the channel defined below to determine which spins are contaminated HIGH_ENERGY_CULL_CHANNEL = 5 # For the high energy cull, we want to combine spin bins because an SEP event is @@ -254,7 +247,3 @@ class UltraConstants: # When True, applies the scattering rejection mask based on the FWHM thresholds # to the L1C fine energy bin maps. APPLY_SCATTERING_REJECTION_L1C: bool = False - - # When true, only use raw direct events for culling in extendedspin.py. Otherwise, - # refer to the priority specified in the extendedspin-config ancillary file. - L1B_USE_RAW_DE_ONLY = False diff --git a/imap_processing/ultra/l1b/extendedspin.py b/imap_processing/ultra/l1b/extendedspin.py index b69348f42..b03b88e0c 100644 --- a/imap_processing/ultra/l1b/extendedspin.py +++ b/imap_processing/ultra/l1b/extendedspin.py @@ -64,29 +64,30 @@ def calculate_extendedspin( aux_dataset = dict_datasets[f"imap_ultra_l1a_{instrument_id}sensor-aux"] rates_dataset = dict_datasets[f"imap_ultra_l1a_{instrument_id}sensor-rates"] status_dataset = dict_datasets[f"imap_ultra_l1b_{instrument_id}sensor-status"] - # Use the raw de for the spin level quantities since it contains all events. + # Use the priority 1 de for the spin level quantities since it contains all events. # The energy dependent culling selects its de dataset per energy range. - de_dataset = de_datasets["p0"] + priority_1_de_dataset = de_datasets["p1"] extendedspin_dict = {} rates_qf, spin, energy_bin_geometric_mean, n_sigma_per_energy = flag_rates( - de_dataset["spin"].values, - de_dataset["energy"].values, + priority_1_de_dataset["spin"].values, + priority_1_de_dataset["energy"].values, ) count_rates, _, _counts, _ = get_energy_histogram( - de_dataset["spin"].values, de_dataset["energy"].values + priority_1_de_dataset["spin"].values, priority_1_de_dataset["energy"].values ) attitude_qf, spin_rates, spin_period, spin_starttime = flag_attitude( - de_dataset["spin"].values, aux_dataset + priority_1_de_dataset["spin"].values, aux_dataset ) # TODO: We will add to this later - hk_qf = flag_hk(de_dataset["spin"].values) - inst_qf = flag_imap_instruments(de_dataset["spin"].values) + hk_qf = flag_hk(priority_1_de_dataset["spin"].values) + inst_qf = flag_imap_instruments(priority_1_de_dataset["spin"].values) spin_bin_size = UltraConstants.SPIN_BIN_SIZE spin_tbin_edges = get_binned_spins_edges( spin, spin_period, spin_starttime, spin_bin_size ) + # Calculate goodtime quality flags. # The culling algorithms should be called in the following order # 1. Low voltage @@ -96,7 +97,9 @@ def calculate_extendedspin( # 5. Spectral cull # 6. Statistical outliers (energy dependent) - voltage_qf = flag_low_voltage(spin_tbin_edges, status_dataset) + voltage_qf = flag_low_voltage( + spin_tbin_edges, status_dataset, extendedspin_conf.voltage_threshold + ) # Get energy bins used at l1c intervals, _, _ = build_energy_bins() # Get the energy ranges @@ -112,14 +115,12 @@ def calculate_extendedspin( instrument_id, ) - # Calculate the high energy quality flags - energy_thresholds = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS high_energy_qf = flag_high_energy( de_counts_summary, spin_tbin_edges, energy_ranges, + extendedspin_conf.energy_thresholds, voltage_qf, - energy_thresholds, ) # For the following culls, mask the spins that have already been flagged for # low voltage @@ -155,9 +156,9 @@ def calculate_extendedspin( # Track rejected events in each spin based on # quality flags in de l1b data. rejected_counts = count_rejected_events_per_spin( - de_dataset["spin"].values, - de_dataset["quality_scattering"].values, - de_dataset["quality_outliers"].values, + priority_1_de_dataset["spin"].values, + priority_1_de_dataset["quality_scattering"].values, + priority_1_de_dataset["quality_outliers"].values, ) # These will be the coordinates. diff --git a/imap_processing/ultra/l1b/lookup_utils.py b/imap_processing/ultra/l1b/lookup_utils.py index 278e98408..a60eff5e1 100644 --- a/imap_processing/ultra/l1b/lookup_utils.py +++ b/imap_processing/ultra/l1b/lookup_utils.py @@ -692,6 +692,7 @@ class ExtendedSpinConfig: """Pointing dependent l1b culling configurations.""" energy_thresholds: np.ndarray # energy thresholds for culling + # (counts per 20-spin bin) voltage_threshold: float # voltage threshold for culling date: datetime.datetime # Date when configuration changed priority: str # Which de product to use priority 1-4 de or raw de. e.g. p0-p4 @@ -736,11 +737,14 @@ def from_csv(cls, config_file_path: str, repointing: str) -> "ExtendedSpinConfig f" per pointing." ) config = filtered_df.iloc[0] - # Get the column names that contain "cullThresh" + # Get the column names that contain "cullThresh" and sort them ascending. thresh_colnames = sorted( [col for col in df.columns if "cullThresh" in col], key=lambda threshold: int(threshold.split("_")[-1]), ) + # The thresholds are counts per spin bin of UltraConstants.SPIN_BIN_SIZE (20) + # spins, so they are used as is. If SPIN_BIN_SIZE changes, the thresholds in + # the config file need to be updated for the new bin size. e_thresholds = np.array([config[col] for col in thresh_colnames]) # The config provides thresholds for bins 0-4; bin 5 (> # UltraConstands.MAX_ENERGY_THRESHOLD keV) reuses the bin 4 threshold diff --git a/imap_processing/ultra/l1b/ultra_l1b_culling.py b/imap_processing/ultra/l1b/ultra_l1b_culling.py index e0b2d31d6..119e2055c 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_culling.py +++ b/imap_processing/ultra/l1b/ultra_l1b_culling.py @@ -631,7 +631,7 @@ def count_rejected_events_per_spin( def flag_low_voltage( spin_tbin_edges: NDArray, status_dataset: xr.Dataset, - voltage_threshold: float = UltraConstants.LOW_VOLTAGE_CULL_THRESHOLD, + voltage_threshold: float, ) -> NDArray: """ Flag low voltage events. @@ -666,8 +666,11 @@ def flag_low_voltage( low_voltage_times = status_dataset["shcoarse"].data[low_voltage_inds] # For each low voltage time, find the corresponding spin time + # Remove the last bin edge of spin_tbin_edges so any trailing spins + # (not assigned a bin) end up in the last bin. This ensures that + # a low voltage event during a trailing spin is flagged. lv_spin_inds = np.atleast_1d( - np.searchsorted(spin_tbin_edges, low_voltage_times, side="right") - 1 + np.searchsorted(spin_tbin_edges[:-1], low_voltage_times, side="right") - 1 ) # Ensure that the indices are within the valid range of spin groups valid_bin_inds = (lv_spin_inds >= 0) & (lv_spin_inds < spin_bin_size) @@ -688,8 +691,8 @@ def flag_high_energy( de_counts_summary: np.ndarray, spin_tbin_edges: NDArray, energy_ranges: NDArray, + energy_thresholds: NDArray, mask: NDArray = None, - energy_thresholds: np.ndarray = UltraConstants.HIGH_ENERGY_CULL_THRESHOLDS, combine_spin_bin_radius: int = UltraConstants.HIGH_ENERGY_COMBINED_SPIN_BIN_RADIUS, ) -> NDArray: """ @@ -704,13 +707,13 @@ def flag_high_energy( Edges of the spin time bins. energy_ranges : numpy.ndarray Array of energy range edges. + energy_thresholds : numpy.ndarray + Array of count thresholds for flagging high energy events corresponding to + each energy range. mask : numpy.ndarray, optional Mask indicating which events to consider for high energy flagging (e.g., after low voltage culling). True indicates the spin bins that should NOT be considered for high energy flagging. - energy_thresholds : numpy.ndarray - Array of count thresholds for flagging high energy events corresponding to - each energy range. combine_spin_bin_radius : int Number of spin bins to combine on either side of the current bin to get a smoother estimate of the counts per bin (see ``combine_de_counts_summary``). @@ -1102,10 +1105,9 @@ def get_valid_de_count_summary( # By default, use priority 1 de unless it's the last energy bin and the # config is p0 or if the config is set to use raw de only. de = de_datasets["p1"] - if UltraConstants.L1B_USE_RAW_DE_ONLY or ( - i == n_energy_ranges - 1 and priority_conf == "p0" - ): + if i == n_energy_ranges - 1 and priority_conf == "p0": de = de_datasets["p0"] + valid_outliers = de["quality_outliers"].values == 0 valid_scattering = de["quality_scattering"].values == 0 # TODO what about species non-proton? For those psets dont cull based on @@ -1172,16 +1174,18 @@ def combine_de_counts_summary( A 2D array of counts per energy range and spin bin for valid events, with counts combined across spin bins. """ - # Pad array along the spin bin axis to ensure sliding_window_view returns - # an array of the correct shape. + # Pad array along the spin bin axis with NaN so that sliding_window_view returns + # an array of the correct shape. To be consistent with the ULTRA IT + # implementation, the padded values are ignored in the mean, so edge bins are + # averaged over only the bins that exist. counts_padded = np.pad( - de_counts_summary, + de_counts_summary.astype(float), ((0, 0), (combine_spin_bin_radius, combine_spin_bin_radius)), - mode="edge", + constant_values=np.nan, ) window_size = combine_spin_bin_radius * 2 + 1 windows = sliding_window_view(counts_padded, window_shape=window_size, axis=1) - return np.mean(windows, axis=-1) + return np.nanmean(windows, axis=-1) def get_valid_earth_angle_events( @@ -1202,8 +1206,8 @@ def get_valid_earth_angle_events( Returns ------- valid_earth_angle_events : NDArray - A boolean array indicating which events have Earth angle greater than the - specified threshold. + A boolean array indicating which events have an angle greater than the + specified threshold from both the Earth and anti-Earth directions. """ velocity_dps_sc = de_dataset_subset["velocity_dps_sc"].values # Use the mean event time to compute the Earth unit vector since the spacecraft @@ -1230,13 +1234,9 @@ def get_valid_earth_angle_events( ) # shape (n_events, 3) # Get cos(theta) between each particle look direction and Earth direction cos_sep = np.dot(unit_look_dirs, earth_unit_vector) # shape (n_events,) - # Clip cos_sep to the valid range of [-1, 1] to avoid numerical issues with arccos - cos_sep = np.clip(cos_sep, -1.0, 1.0) - sep_angle = np.arccos(cos_sep) - # An event is valid if the separation angle between the particle look - # direction and Earth direction is greater than the Earth angle limit - # (i.e., the Earth is outside the field of view). - return sep_angle > earth_ang_45 + # An event is valid if the particle look direction is more than the Earth angle + # limit away from both the Earth direction and the anti-Earth direction. + return np.abs(cos_sep) < np.cos(earth_ang_45) def get_energy_range_flags(energy_ranges_edges: NDArray) -> NDArray: @@ -1404,5 +1404,8 @@ def expand_bin_flags_to_spins( ) repeated_flags = repeated_flags[:n_spins] quality_flags[: len(repeated_flags)] = repeated_flags + # Spins in an incomplete trailing bin inherit the flag of the last complete bin. + if len(repeated_flags) < n_spins: + quality_flags[len(repeated_flags) :] = binned_quality_flags[-1] return quality_flags