diff --git a/imap_processing/tests/ultra/unit/test_l1c_lookup_utils.py b/imap_processing/tests/ultra/unit/test_l1c_lookup_utils.py index ae1e1bbe4..d2149d4ea 100644 --- a/imap_processing/tests/ultra/unit/test_l1c_lookup_utils.py +++ b/imap_processing/tests/ultra/unit/test_l1c_lookup_utils.py @@ -4,10 +4,11 @@ import xarray as xr from imap_processing.ultra.l1c.l1c_lookup_utils import ( - calculate_fwhm_spun_scattering, + calculate_accepted_pixels, get_scattering_thresholds_for_energy, get_spacecraft_pointing_lookup_tables, get_static_deadtime_ratios, + in_restricted_fov, mask_below_fwhm_scattering_threshold, ) @@ -31,6 +32,7 @@ def test_get_spacecraft_pointing_lookup_tables(ancillary_files): assert theta_vals.shape == (cols, npix) assert phi_vals.shape == (cols, npix) assert ra_and_dec.shape == (2, npix) + assert boundary_scale_factors.shape == (cols, npix) # Value tests assert for_indices_by_spin_phase.dtype == bool @@ -115,21 +117,31 @@ def test_get_static_deadtime_ratios(ancillary_files): assert np.all((dt_ratio >= 0.0) & (dt_ratio <= 1.0)) -def test_calculate_fwhm_spun_scattering(ancillary_files): - """Test calculate_fwhm_spun_scattering function.""" +def test_in_restricted_fov(): + """Test in_restricted_fov function.""" + # Create mock theta and phi values + # First two values are outside the restricted FOV, the rest are inside + theta_vals = np.array([[-50, 49, 10, 40, 35]]) + # The last value is outside the restricted FOV + phi_vals = np.array([[20, 30, 40, 50, 70]]) # shape (1, 5) + accepted_pixels = in_restricted_fov(theta_vals, phi_vals, 45) + expected_accepted_pixels = np.array([[False, False, True, True, False]]) + np.testing.assert_array_equal(accepted_pixels, expected_accepted_pixels) + + +def test_calculate_accepted_pixels(ancillary_files): + """Test calculate_accepted_pixels function.""" # Make array with ones (we are only testing the shape here) for_pixels = np.ones((50, 10)) theta_vals = np.ones((50, 10)) * 20 # All theta values are 20 phi_vals = np.ones((50, 5)) * 15 # All phi with pytest.raises(ValueError, match="Shape mismatch"): - calculate_fwhm_spun_scattering( - for_pixels, theta_vals, phi_vals, ancillary_files, 45 - ) + calculate_accepted_pixels(for_pixels, theta_vals, phi_vals, ancillary_files, 45) @pytest.mark.external_test_data -def test_calculate_fwhm_spun_scattering_reject(ancillary_files): - """Test calculate_fwhm_spun_scattering function.""" +def test_calculate_accepted_pixels_reject(ancillary_files): + """Test calculate_accepted_pixels function.""" nside = 8 pix = hp.nside2npix(nside) steps = 5 # Reduced for testing @@ -144,16 +156,49 @@ def test_calculate_fwhm_spun_scattering_reject(ancillary_files): # Simulate first 100 pixels are in the FOR for all spin phases inside_inds = 100 for_pixels[:, :, :inside_inds] = True - valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = ( - calculate_fwhm_spun_scattering( - for_pixels, - mock_theta, - mock_phi, - ancillary_files, - 45, - reject_scattering=True, - ) + valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = calculate_accepted_pixels( + for_pixels, + mock_theta, + mock_phi, + ancillary_files, + 45, + reject_scattering=True, ) assert valid_spun_pixels.shape == (steps, energy_dim, pix) # Check that some pixels are rejected assert not np.array_equal(valid_spun_pixels, for_pixels) + + +@pytest.mark.external_test_data +def test_calculate_accepted_pixels_restrict_fov(ancillary_files): + """Test calculate_accepted_pixels function for FOV restrictions.""" + nside = 8 + pix = hp.nside2npix(nside) + steps = 5 # Reduced for testing + energy_dim = 46 + np.random.seed(42) + mock_theta = np.random.uniform(-60, 60, (steps, energy_dim, pix)) + mock_phi = np.random.uniform(-60, 60, (steps, energy_dim, pix)) + # Create for_pixels with all True values to isolate the effect of FOV restrictions + for_pixels = xr.DataArray( + np.ones((steps, energy_dim, pix)).astype(bool), + dims=("spin_phase_step", "energy", "pixel"), + ) + # Set first 30 pixels to be outside of the FOR for all spin phases and energies + outside_inds = 30 + for_pixels[:, :, :outside_inds] = False + valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = calculate_accepted_pixels( + for_pixels, + mock_theta, + mock_phi, + ancillary_files, + 45, + apply_fov_restriction=True, + ) + assert valid_spun_pixels.shape == (steps, energy_dim, pix) + # In this case, the valid spun pixels should be the same as the pixels that are + # in the restricted FOV except the first 30 pixels which are outside the FOR and + # should be False for all spin phases and energies. + expected_accepted_px = in_restricted_fov(mock_theta, mock_phi, 45) + expected_accepted_px[:, :, :outside_inds] = False + assert np.array_equal(expected_accepted_px, valid_spun_pixels) diff --git a/imap_processing/tests/ultra/unit/test_spacecraft_pset.py b/imap_processing/tests/ultra/unit/test_spacecraft_pset.py index 5c57b5458..2a0be5d50 100644 --- a/imap_processing/tests/ultra/unit/test_spacecraft_pset.py +++ b/imap_processing/tests/ultra/unit/test_spacecraft_pset.py @@ -79,6 +79,8 @@ def test_calculate_spacecraft_pset( ["epoch", "component"], particle_velocity_dps_spacecraft, ), + "theta": (["epoch"], np.zeros(len(species), dtype=np.float32)), + "phi": (["epoch"], np.zeros(len(species), dtype=np.float32)), "energy_spacecraft": (["epoch"], energy_dps_spacecraft), "spin": (["epoch"], df["Spin"].values), "quality_scattering": ( @@ -180,6 +182,8 @@ def test_calculate_spacecraft_pset_with_cdf( de_dict["event_times"] = 817561854.185627 + ( df_subset["tdb"].values - df_subset["tdb"].values[0] ) + de_dict["theta"] = np.zeros(len(df_subset), dtype=np.float32) + de_dict["phi"] = np.zeros(len(df_subset), dtype=np.float32) name = "imap_ultra_l1b_45sensor-de" dataset = create_dataset(de_dict, name, "l1b") diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1c.py b/imap_processing/tests/ultra/unit/test_ultra_l1c.py index fbe6903e2..bb8b8d1d4 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1c.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1c.py @@ -179,6 +179,8 @@ def test_calculate_spacecraft_pset_with_cdf( de_dict["spin"] = np.full(len(sc_dps_velocity), 0) de_dict["species"] = np.ones(len(sc_dps_velocity), dtype=np.uint8) de_dict["ebin"] = np.ones(len(sc_dps_velocity), dtype=np.uint8) + de_dict["theta"] = np.zeros(len(df_subset), dtype=np.float32) + de_dict["phi"] = np.zeros(len(df_subset), dtype=np.float32) de_dict["event_times"] = df_subset["tdb"].values name = "imap_ultra_l1b_45sensor-de" @@ -245,6 +247,8 @@ def test_calculate_helio_pset_with_cdf( # Fake SCLK in seconds that matches SPICE. de_dict["event_times"] = np.full(len(df_subset), 2.41187e13) de_dict["ebin"] = np.ones(len(df_subset), dtype=np.uint8) + de_dict["theta"] = np.zeros(len(df_subset), dtype=np.float32) + de_dict["phi"] = np.zeros(len(df_subset), dtype=np.float32) species_bin = np.full(len(df_subset), 1, dtype=np.uint8) # PosYSlit is True for left (start_type = 1) diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py b/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py index 375b0e237..1751dd6af 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py +++ b/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py @@ -14,7 +14,7 @@ from imap_processing.ultra.constants import UltraConstants from imap_processing.ultra.l1c import ultra_l1c_pset_bins from imap_processing.ultra.l1c.spacecraft_pset import ( - calculate_fwhm_spun_scattering, + calculate_accepted_pixels, ) from imap_processing.ultra.l1c.ultra_l1c_pset_bins import ( build_energy_bins, @@ -269,15 +269,13 @@ def test_apply_deadtime_correction( """Tests apply_deadtime_correction function.""" mock_theta, mock_phi, spin_phase_steps, inside_inds, pix, steps = spun_index_data deadtime_ratios = xr.DataArray(np.ones(steps), dims="spin_phase_step") - valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = ( - calculate_fwhm_spun_scattering( - spin_phase_steps, - mock_theta, - mock_phi, - ancillary_files, - 45, - reject_scattering=False, - ) + valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = calculate_accepted_pixels( + spin_phase_steps, + mock_theta, + mock_phi, + ancillary_files, + 45, + reject_scattering=False, ) boundary_sf = xr.DataArray(np.ones((pix, steps)), dims=("pixel", "spin_phase_step")) exposure_pointing_adjusted = calculate_exposure_time( @@ -303,15 +301,13 @@ def test_apply_deadtime_correction_energy_dep( deadtime_ratios = xr.DataArray(np.ones(steps), dims="spin_phase_step") boundary_sf = xr.DataArray(np.ones((steps, pix)), dims=("spin_phase_step", "pixel")) - valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = ( - calculate_fwhm_spun_scattering( - spin_phase_steps, - mock_theta, - mock_phi, - ancillary_files, - 45, - reject_scattering=True, - ) + valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = calculate_accepted_pixels( + spin_phase_steps, + mock_theta, + mock_phi, + ancillary_files, + 45, + reject_scattering=True, ) exposure_pointing_adjusted = calculate_exposure_time( @@ -349,15 +345,13 @@ def test_get_eff_and_gf(imap_ena_sim_metakernel, ancillary_files, spun_index_dat # Simulate first 100 pixels are in the FOR for all spin phases inside_inds = 100 spin_phase_steps[:, :, :inside_inds] = True - valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = ( - calculate_fwhm_spun_scattering( - spin_phase_steps, - mock_theta, - mock_phi, - ancillary_files, - 45, - reject_scattering=False, - ) + valid_spun_pixels, fwhm_theta, fwhm_phi, thresholds = calculate_accepted_pixels( + spin_phase_steps, + mock_theta, + mock_phi, + ancillary_files, + 45, + reject_scattering=False, ) boundary_sf = xr.DataArray( np.ones((steps, energy_dim, pix)), dims=("spin_phase_step", "energy", "pixel") @@ -409,7 +403,7 @@ def test_get_spacecraft_exposure_times( ) # Spin phase steps, random 0 or 1 pixels_below_threshold, fwhm_theta, fwhm_phi, thresholds = ( - calculate_fwhm_spun_scattering( + calculate_accepted_pixels( spin_phase_steps, mock_theta, mock_phi, ancillary_files, 45 ) ) diff --git a/imap_processing/ultra/constants.py b/imap_processing/ultra/constants.py index f3a7e76f1..ac315a138 100644 --- a/imap_processing/ultra/constants.py +++ b/imap_processing/ultra/constants.py @@ -181,6 +181,13 @@ class UltraConstants: FOV_THETA_OFFSET_DEG = 0.0 FOV_PHI_LIMIT_DEG = 60.0 + # Restricted FOV theta/phi acceptance limits (degrees). + # Samples outside these bounds are excluded from GF, efficiency, exposure, + # and counts maps at L1C (fine energy bin maps only). + RESTRICTED_FOV_THETA_LOW_DEG_45: float = -43.0 + RESTRICTED_FOV_THETA_HIGH_DEG_45: float = 43.0 + RESTRICTED_FOV_THETA_LOW_DEG_90: float = -43.0 + RESTRICTED_FOV_THETA_HIGH_DEG_90: float = 43.0 # For spatiotemporal culling EARTH_RADIUS_KM: float = 6378.1 @@ -224,3 +231,19 @@ class UltraConstants: # and pad with fill values if we use fewer bins. MAX_ENERGY_RANGES = 16 MAX_ENERGY_RANGE_EDGES = MAX_ENERGY_RANGES + 1 + + # L1C PSET constants + + # When True, applies the FOV restrictions defined above to the L1C fine energy bin + # maps (GF, efficiency, exposure, counts). This culls regions of the instrument + # field of view with poor efficiency calibration from inclusion into the map making + # process. + APPLY_FOV_RESTRICTIONS_L1C: bool = True + + # When True, applies the boundary scale factors from the ancillary file to exposure + # time, efficiency, and geometric factor maps. + APPLY_BOUNDARY_SCALE_FACTORS_L1C: bool = False + + # 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 diff --git a/imap_processing/ultra/l1c/helio_pset.py b/imap_processing/ultra/l1c/helio_pset.py index cf98368e8..9d562aa5f 100644 --- a/imap_processing/ultra/l1c/helio_pset.py +++ b/imap_processing/ultra/l1c/helio_pset.py @@ -21,7 +21,8 @@ ) from imap_processing.ultra.l1c.l1c_lookup_utils import ( build_energy_bins, - calculate_fwhm_spun_scattering, + calculate_accepted_pixels, + in_restricted_fov, ) from imap_processing.ultra.l1c.make_helio_index_maps import ( make_helio_index_maps_with_nominal_kernels, @@ -76,10 +77,6 @@ def calculate_helio_pset( dataset : xarray.Dataset Dataset containing the data. """ - # Do not cull events based on scattering thresholds - reject_scattering = False - # Do not apply boundary scale factor corrections - apply_bsf = False nside = 32 num_spin_steps = 720 sensor_id = int(parse_filename_like(name)["sensor"][0:2]) @@ -95,7 +92,7 @@ def calculate_helio_pset( rejected = get_de_rejection_mask( species_dataset["quality_scattering"].values, species_dataset["quality_outliers"].values, - reject_scattering, + UltraConstants.APPLY_SCATTERING_REJECTION_L1C, ) species_dataset = species_dataset.isel(epoch=~rejected) # Check if spin_number is in the goodtimes dataset, if not then we can @@ -114,13 +111,6 @@ def calculate_helio_pset( species_dataset["spin"].values, ) species_dataset = species_dataset.isel(epoch=~energy_dependent_rejected) - v_mag_helio_spacecraft = np.linalg.norm( - species_dataset["velocity_dps_helio"].values, axis=1 - ) - vhat_dps_helio = ( - species_dataset["velocity_dps_helio"].values - / v_mag_helio_spacecraft[:, np.newaxis] - ) # Get the start and stop times of the pointing period repoint_id = species_dataset.attrs.get("Repointing", None) if repoint_id is None: @@ -138,7 +128,7 @@ def calculate_helio_pset( spin_duration=15.0, num_steps=num_spin_steps, instrument_frame=instrument_frame, - compute_bsf=apply_bsf, + compute_bsf=UltraConstants.APPLY_BOUNDARY_SCALE_FACTORS_L1C, ) boundary_scale_factors = helio_pointing_ds.bsf theta_vals = helio_pointing_ds.theta @@ -147,16 +137,37 @@ def calculate_helio_pset( logger.info("calculating spun FWHM scattering values.") pixels_below_scattering, scattering_theta, scattering_phi, scattering_thresholds = ( - calculate_fwhm_spun_scattering( + calculate_accepted_pixels( fov_index, theta_vals, phi_vals, ancillary_files, instrument_id, - reject_scattering, + reject_scattering=UltraConstants.APPLY_SCATTERING_REJECTION_L1C, + apply_fov_restriction=UltraConstants.APPLY_FOV_RESTRICTIONS_L1C, ) ) - + # Apply restricted FOV filter to counts. + # If a valid event's instrument frame theta/phi falls outside the restricted + # acceptance window it is excluded from the L1C fine energy-bin counts map. + if UltraConstants.APPLY_FOV_RESTRICTIONS_L1C: + fov_accepted = in_restricted_fov( + species_dataset["theta"].values, + species_dataset["phi"].values, + instrument_id, + ) + logger.info( + f"Restricted FOV counts filter: keeping {fov_accepted.sum()} / " + f"{len(fov_accepted)} events." + ) + species_dataset = species_dataset.isel(epoch=fov_accepted) + v_mag_helio_spacecraft = np.linalg.norm( + species_dataset["velocity_dps_helio"].values, axis=1 + ) + vhat_dps_helio = ( + species_dataset["velocity_dps_helio"].values + / v_mag_helio_spacecraft[:, np.newaxis] + ) counts, counts_n_pix = get_spacecraft_histogram( vhat_dps_helio, species_dataset["energy_heliosphere"].values, @@ -184,7 +195,7 @@ def calculate_helio_pset( energy_bins=energy_bin_geometric_means, sensor_id=sensor_id, ancillary_files=ancillary_files, - apply_bsf=apply_bsf, + apply_bsf=UltraConstants.APPLY_BOUNDARY_SCALE_FACTORS_L1C, goodtimes_dataset=goodtimes_dataset, ) logger.info("Calculating spun efficiencies and geometric function.") @@ -197,7 +208,7 @@ def calculate_helio_pset( n_pix, sensor_id, ancillary_files, - apply_bsf, + apply_bsf=UltraConstants.APPLY_BOUNDARY_SCALE_FACTORS_L1C, ) logger.info("Calculating background rates.") diff --git a/imap_processing/ultra/l1c/l1c_lookup_utils.py b/imap_processing/ultra/l1c/l1c_lookup_utils.py index 1d9d87fa5..15d5f541a 100644 --- a/imap_processing/ultra/l1c/l1c_lookup_utils.py +++ b/imap_processing/ultra/l1c/l1c_lookup_utils.py @@ -17,6 +17,46 @@ logger = logging.getLogger(__name__) +def in_restricted_fov( + theta: np.ndarray, + phi: np.ndarray, + instrument_id: int, +) -> np.ndarray: + """ + Determine whether the theta/phi pairs are inside the restricted FOV. + + Parameters + ---------- + theta : np.ndarray + Array of theta values in degrees. + phi : np.ndarray + Array of phi values in degrees. + instrument_id : int + Instrument ID, either 45 or 90. + + Returns + ------- + np.ndarray + Boolean array indicating whether each theta/phi pair is within the restricted + FOV. + """ + # Theta limits are dependent on the instrument id + if instrument_id == 90: + low_theta_lim = UltraConstants.RESTRICTED_FOV_THETA_LOW_DEG_90 + high_theta_lim = UltraConstants.RESTRICTED_FOV_THETA_HIGH_DEG_90 + elif instrument_id == 45: + low_theta_lim = UltraConstants.RESTRICTED_FOV_THETA_LOW_DEG_45 + high_theta_lim = UltraConstants.RESTRICTED_FOV_THETA_HIGH_DEG_45 + else: + raise ValueError(f"Invalid instrument ID: {instrument_id}. Must be 45 or 90.") + + return ( + (theta >= low_theta_lim) + & (theta <= high_theta_lim) + & (np.abs(phi) < UltraConstants.FOV_PHI_LIMIT_DEG) + ) + + def mask_below_fwhm_scattering_threshold( theta_coeffs: np.ndarray, phi_coeffs: np.ndarray, @@ -67,18 +107,20 @@ def mask_below_fwhm_scattering_threshold( return scattering_mask, fwhm_theta, fwhm_phi -def calculate_fwhm_spun_scattering( +def calculate_accepted_pixels( # noqa: PLR0912 for_indices_by_spin_phase: xr.DataArray, theta_vals: np.ndarray, phi_vals: np.ndarray, ancillary_files: dict, instrument_id: int, reject_scattering: bool = False, + apply_fov_restriction: bool = False, ) -> tuple[xr.DataArray, NDArray, NDArray, NDArray]: """ - Calculate FWHM scattering values for each pixel, energy bin, and spin phase step. + Calculate the accepted pixels based scattering and FOV restrictions. - This function also calculates a mask for pixels that are below the FWHM threshold. + This function also calculates a mask for pixels that are below the FWHM threshold + and FWHM scattering values for each pixel, energy bin, and spin phase step. Parameters ---------- @@ -99,15 +141,23 @@ def calculate_fwhm_spun_scattering( Instrument ID, either 45 or 90. reject_scattering : bool Whether to reject pixels based on scattering thresholds. + apply_fov_restriction : bool + Whether to apply the restricted FOV theta/phi acceptance test. When + ``True``, any spin-step/pixel combination whose instrument-frame theta + and phi do not satisfy :func:`in_restricted_fov` is skipped entirely: + it is not added to the averaging numerator *or* denominator, and it + does not contribute exposure time. This gates GF, efficiency, and + exposure for the fine L1C energy-bin maps. Returns ------- valid_spun_pixels : xarray.DataArray - Boolean array indicating, for each spin phase step, energy_bin, pixel, - the pixel is inside the Field Of Regard (FOR) and whether the pixel is inside the - FOR at that spin phase and its computed FWHM at that energy is below the - scattering threshold. If reject_scattering is False, this will just reflect - the FOR mask (for_indices_by_spin_phase). + Boolean array of shape ``(spin_phase_step, energy_bin, pixel)`` indicating + which pixel samples are accepted for L1C accumulation at each spin step. + Acceptance always requires the sample to be inside the Field of Regard (FOR), + can optionally require passing the restricted theta/phi FOV criteria when + ``apply_fov_restriction=True``, and can optionally require passing the FWHM + scattering threshold when ``reject_scattering=True``. scattering_fwhm_theta : NDArray Calculated FWHM scatting values for theta at each energy bin and averaged over spin phase. @@ -143,9 +193,10 @@ def calculate_fwhm_spun_scattering( steps = for_indices_by_spin_phase.sizes["spin_phase_step"] energies = energy_bin_geometric_means[np.newaxis, :] # Initialize DataArray to hold boolean of valid pixels at each spin phase step - # If reject_scattering if false, this will just be the FOR mask. + # If reject_scattering and apply_fov_restriction are both False, this will just + # be the FOR mask. spun_dims = ("spin_phase_step", "energy", "pixel") - if reject_scattering: + if reject_scattering or apply_fov_restriction: valid_pixels = xr.DataArray( np.zeros((steps, len(energy_bin_geometric_means), n_pix), dtype=bool), dims=spun_dims, @@ -156,6 +207,8 @@ def calculate_fwhm_spun_scattering( ).transpose(*spun_dims) else: valid_pixels = for_indices_by_spin_phase + # TODO refactor loop below and combine the energy dependent and independent cases + # if possible to avoid code duplication. # The "for_indices_by_spin_phase" lookup table contains the boolean values of each # pixel at each spin phase step, indicating whether the pixel is inside the FOR. # It starts at Spin-phase = 0, and increments in fine steps (1 ms), spinning the @@ -167,14 +220,27 @@ def calculate_fwhm_spun_scattering( if for_inds.ndim > 1: # Energy dependent FOR indices for e_ind in range(len(energy_bin_geometric_means)): - for_inds_energy = for_inds[e_ind, :] - # Skip if no pixels in FOR - if not np.any(for_inds_energy): + if not np.any(for_inds[e_ind, :]): + continue + accepted_pix = np.flatnonzero(for_inds[e_ind, :]) + theta = theta_vals[i, e_ind, accepted_pix] + phi = phi_vals[i, e_ind, accepted_pix] + + # Check if we need to restrict the fov further using theta/phi + # acceptance limits + if apply_fov_restriction: + fov_mask = in_restricted_fov(theta, phi, instrument_id) + # update accepted pixel indices to reflect the FOV restriction + accepted_pix = accepted_pix[fov_mask] + theta = theta[fov_mask] + phi = phi[fov_mask] + # update valid pixels to reflect the FOV restriction + valid_pixels[i, e_ind, accepted_pix] = True + + if len(accepted_pix) == 0: continue - theta = theta_vals[i, e_ind, for_inds_energy] - phi = phi_vals[i, e_ind, for_inds_energy] theta_coeffs, phi_coeffs = get_scattering_coefficients( theta.data, phi.data, lookup_tables=scattering_luts ) @@ -192,19 +258,37 @@ def calculate_fwhm_spun_scattering( ) # If rejecting scattering, store the mask if reject_scattering: - valid_pixels[i, e_ind, for_inds_energy] = scattering_mask.flatten() + valid_pixels[i, e_ind, accepted_pix] = scattering_mask.flatten() # Accumulate FWHM values - fwhm_theta_sum[e_ind, for_inds_energy] += fwhm_theta.flatten() - fwhm_phi_sum[e_ind, for_inds_energy] += fwhm_phi.flatten() - sample_count[e_ind, for_inds_energy] += 1 + fwhm_theta_sum[e_ind, accepted_pix] += fwhm_theta.flatten() + fwhm_phi_sum[e_ind, accepted_pix] += fwhm_phi.flatten() + sample_count[e_ind, accepted_pix] += 1 else: # Energy independent FOR indices if not np.any(for_inds): continue - theta = theta_vals[i, for_inds] - phi = phi_vals[i, for_inds] + accepted_pix = np.flatnonzero(for_inds) # Get indices of pixels in FOR for + # this spin phase step + + theta = theta_vals[i, accepted_pix] + phi = phi_vals[i, accepted_pix] + + # Check if we need to apply the restricted FOV mask before calculating + # scattering coefficients + if apply_fov_restriction: + fov_mask = in_restricted_fov(theta, phi, instrument_id) + # update accepted pixel indices to reflect the FOV restriction + accepted_pix = accepted_pix[fov_mask] + theta = theta[fov_mask] + phi = phi[fov_mask] + # update valid pixels to reflect the FOV restriction + valid_pixels[i, :, accepted_pix] = True + + if len(accepted_pix) == 0: + continue + theta_coeffs, phi_coeffs = get_scattering_coefficients( theta, phi, lookup_tables=scattering_luts ) @@ -216,14 +300,13 @@ def calculate_fwhm_spun_scattering( scattering_thresholds=scattering_thresholds_for_energy_mean, ) ) - if reject_scattering: - valid_pixels[i, :, for_inds] = scattering_mask.T + valid_pixels[i, :, accepted_pix] = scattering_mask.T # Accumulate FWHM values - fwhm_theta_sum[:, for_inds] += fwhm_theta.T - fwhm_phi_sum[:, for_inds] += fwhm_phi.T - sample_count[:, for_inds] += 1 + fwhm_theta_sum[:, accepted_pix] += fwhm_theta.T + fwhm_phi_sum[:, accepted_pix] += fwhm_phi.T + sample_count[:, accepted_pix] += 1 fwhm_phi_avg = np.zeros_like(fwhm_phi_sum) fwhm_theta_avg = np.zeros_like(fwhm_theta_sum) diff --git a/imap_processing/ultra/l1c/spacecraft_pset.py b/imap_processing/ultra/l1c/spacecraft_pset.py index 80324a6c8..313b39df1 100644 --- a/imap_processing/ultra/l1c/spacecraft_pset.py +++ b/imap_processing/ultra/l1c/spacecraft_pset.py @@ -1,4 +1,4 @@ -"""Calculate Pointing Set Grids.""" +"""Calculate Spacecraft Pointing Set Grids.""" import logging @@ -20,8 +20,9 @@ ) from imap_processing.ultra.l1c.l1c_lookup_utils import ( build_energy_bins, - calculate_fwhm_spun_scattering, + calculate_accepted_pixels, get_spacecraft_pointing_lookup_tables, + in_restricted_fov, ) from imap_processing.ultra.l1c.ultra_l1c_culling import compute_culling_mask from imap_processing.ultra.l1c.ultra_l1c_pset_bins import ( @@ -73,10 +74,6 @@ def calculate_spacecraft_pset( dataset : xarray.Dataset Dataset containing the data. """ - # Do not cull events based on scattering thresholds - reject_scattering = False - # Do not apply boundary scale factor corrections - apply_bsf = False pset_dict: dict[str, np.ndarray] = {} sensor_id = int(parse_filename_like(name)["sensor"][0:2]) @@ -93,7 +90,7 @@ def calculate_spacecraft_pset( de_rejected = get_de_rejection_mask( species_dataset["quality_scattering"].values, species_dataset["quality_outliers"].values, - reject_scattering, + UltraConstants.APPLY_SCATTERING_REJECTION_L1C, ) species_dataset = species_dataset.isel(epoch=~de_rejected) # Check if spin_number is in the goodtimes dataset, if not then we can @@ -112,12 +109,6 @@ def calculate_spacecraft_pset( species_dataset["spin"].values, ) species_dataset = species_dataset.isel(epoch=~energy_dependent_rejected) - v_mag_dps_spacecraft = np.linalg.norm( - species_dataset["velocity_dps_sc"].values, axis=1 - ) - vhat_dps_spacecraft = ( - species_dataset["velocity_dps_sc"].values / v_mag_dps_spacecraft[:, np.newaxis] - ) # Get lookup table for FOR indices by spin phase step ( @@ -130,15 +121,38 @@ def calculate_spacecraft_pset( logger.info("calculating spun FWHM scattering values.") valid_spun_pixels, scattering_theta, scattering_phi, scattering_thresholds = ( - calculate_fwhm_spun_scattering( + calculate_accepted_pixels( for_indices_by_spin_phase, theta_vals, phi_vals, ancillary_files, instrument_id, - reject_scattering, + reject_scattering=UltraConstants.APPLY_SCATTERING_REJECTION_L1C, + apply_fov_restriction=UltraConstants.APPLY_FOV_RESTRICTIONS_L1C, ) ) + + # Apply restricted FOV filter to counts. + # If a valid event's instrument frame theta/phi falls outside the restricted + # acceptance window it is excluded from the L1C fine energy-bin counts map. + if UltraConstants.APPLY_FOV_RESTRICTIONS_L1C: + fov_accepted = in_restricted_fov( + species_dataset["theta"].values, + species_dataset["phi"].values, + instrument_id, + ) + logger.info( + f"Restricted FOV counts filter: keeping {fov_accepted.sum()} / " + f"{len(fov_accepted)} events." + ) + species_dataset = species_dataset.isel(epoch=fov_accepted) + + v_mag_dps_spacecraft = np.linalg.norm( + species_dataset["velocity_dps_sc"].values, axis=1 + ) + vhat_dps_spacecraft = ( + species_dataset["velocity_dps_sc"].values / v_mag_dps_spacecraft[:, np.newaxis] + ) counts, counts_n_pix = get_spacecraft_histogram( vhat_dps_spacecraft, species_dataset["energy_spacecraft"].values, @@ -172,7 +186,7 @@ def calculate_spacecraft_pset( energy_bins=energy_bin_geometric_means, sensor_id=sensor_id, ancillary_files=ancillary_files, - apply_bsf=apply_bsf, + apply_bsf=UltraConstants.APPLY_BOUNDARY_SCALE_FACTORS_L1C, goodtimes_dataset=goodtimes_dataset, ) logger.info("Calculating spun efficiencies and geometric function.") @@ -185,7 +199,7 @@ def calculate_spacecraft_pset( n_pix, sensor_id, ancillary_files, - apply_bsf, + apply_bsf=UltraConstants.APPLY_BOUNDARY_SCALE_FACTORS_L1C, ) sensitivity = efficiencies * geometric_function diff --git a/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py b/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py index 9097d1926..29ab305c9 100644 --- a/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py +++ b/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py @@ -661,7 +661,7 @@ def get_efficiencies_and_geometric_function( # Accumulate and sum eff and gf values bsfs = ( - boundary_scale_factors[pixel_inds, i] + boundary_scale_factors[i, pixel_inds] if apply_bsf else np.ones(len(pixel_inds)) )