ULTRA l1c theta limit check

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)

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")

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)

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
         )
     )

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

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.")