ULTRA l1b use energy bin's geometric mean for scattering culling (#2373)

* use the energy bins mean instead of the event energy for calculating scattering flags

Author: lacoak21

imap_processing/quality_flags.py

@@ -43,7 +43,8 @@ class ImapDEOutliersUltraFlags(FlagNameMixin):
     NONE = CommonFlags.NONE
     FOV = 2**0  # bit 0
     PHCORR = 2**1  # bit 1
-    COINPH = 2**2  # bit 4 # Event validity
+    COINPH = 2**2  # bit 2 # Event validity
+    INVALID_ENERGY = 2**3  # bit 3
 
 
 class ImapHkUltraFlags(FlagNameMixin):

imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py

@@ -354,12 +354,13 @@ def test_get_de_energy_kev(test_fixture):
         test_d,
         test_tof,
     )
-
-    energy = get_de_energy_kev(v, species_bin_ph)
+    quality_flags = np.zeros_like(species_bin_ph)
+    energy = get_de_energy_kev(v, species_bin_ph, quality_flags)
     index_hydrogen = np.where(species_bin_ph == 1)
     actual_energy = energy[index_hydrogen[0]]
     expected_energy = df_ph["energy_revised"].astype("float")
-
+    # All flags should be zero
+    assert np.sum(quality_flags) == 0
     np.testing.assert_allclose(actual_energy, expected_energy, atol=1e-01, rtol=0)
 
 

imap_processing/ultra/l1b/de.py

@@ -304,7 +304,9 @@ def calculate_de(
     de_dict["direct_event_unit_position"] = r_hat.astype(np.float32)
 
     tof_energy[valid_indices] = get_de_energy_kev(
-        velocities[valid_indices], species_bin[valid_indices]
+        velocities[valid_indices],
+        species_bin[valid_indices],
+        quality_flags[valid_indices],
     )
     de_dict["tof_energy"] = tof_energy
     de_dict["energy"] = energy
@@ -348,8 +350,12 @@ def calculate_de(
     de_dict["velocity_dps_sc"] = sc_dps_velocity
     de_dict["velocity_dps_helio"] = helio_velocity
 
-    de_dict["energy_spacecraft"] = get_de_energy_kev(sc_dps_velocity, species_bin)
-    de_dict["energy_heliosphere"] = get_de_energy_kev(helio_velocity, species_bin)
+    de_dict["energy_spacecraft"] = get_de_energy_kev(
+        sc_dps_velocity, species_bin, quality_flags
+    )
+    de_dict["energy_heliosphere"] = get_de_energy_kev(
+        helio_velocity, species_bin, quality_flags
+    )
 
     de_dict["phi_fwhm"], de_dict["theta_fwhm"] = get_fwhm(
         start_type,
@@ -376,6 +382,7 @@ def calculate_de(
         ancillary_files,
         "l1b-sensor-gf-noblades",
     )
+
     de_dict["quality_outliers"] = quality_flags
     flag_scattering(
         de_dict["tof_energy"],

imap_processing/ultra/l1b/ultra_l1b_culling.py

@@ -22,6 +22,7 @@
     get_scattering_thresholds,
 )
 from imap_processing.ultra.l1b.quality_flag_filters import DE_QUALITY_FLAG_FILTERS
+from imap_processing.ultra.l1c.l1c_lookup_utils import build_energy_bins
 
 logging.basicConfig(level=logging.INFO)
 logger = logging.getLogger(__name__)
@@ -514,7 +515,13 @@ def flag_scattering(
         Quality flags.
     """
     scattering_thresholds = get_scattering_thresholds(ancillary_files)
-
+    _, _, energy_bin_geometric_means = build_energy_bins()
+    energy_bin_inds = np.digitize(tof_energy, UltraConstants.PSET_ENERGY_BIN_EDGES)
+    # Clip indices to valid range (events outside the energy bins get assigned
+    # to the nearest bin. These events have already been flagged and
+    # will be ignored in l1c)
+    energy_bin_inds = np.clip(energy_bin_inds, 1, len(energy_bin_geometric_means))
+    energy_geom_means = energy_bin_geometric_means[energy_bin_inds - 1]
     for (e_min, e_max), threshold in scattering_thresholds.items():
         event_mask = (tof_energy >= e_min) & (tof_energy < e_max)
         # Input the theta and phi values for the current energy range.
@@ -528,8 +535,11 @@ def flag_scattering(
         )
         # FWHM_PHI = A_PHI * E^G_PHI
         # FWHM_THETA = A_THETA * E^G_THETA
-        fwhm_theta = theta_coeffs[:, 0] * tof_energy[event_mask] ** theta_coeffs[:, 1]
-        fwhm_phi = phi_coeffs[:, 0] * tof_energy[event_mask] ** phi_coeffs[:, 1]
+        # Use the geometric mean of the energy bin for the scattering check
+        fwhm_theta = (
+            theta_coeffs[:, 0] * energy_geom_means[event_mask] ** theta_coeffs[:, 1]
+        )
+        fwhm_phi = phi_coeffs[:, 0] * energy_geom_means[event_mask] ** phi_coeffs[:, 1]
         is_nan = np.isnan(fwhm_theta) | np.isnan(fwhm_phi)
         quality_flags[np.where(event_mask)[0][is_nan]] |= (
             ImapDEScatteringUltraFlags.NAN_PHI_OR_THETA.value

imap_processing/ultra/l1b/ultra_l1b_extended.py

@@ -616,7 +616,9 @@ def get_ssd_tof(
     return np.asarray(tof, dtype=np.float64)
 
 
-def get_de_energy_kev(v: np.ndarray, species: np.ndarray) -> NDArray:
+def get_de_energy_kev(
+    v: np.ndarray, species: np.ndarray, quality_flags: np.ndarray | None = None
+) -> NDArray:
     """
     Calculate the direct event energy.
 
@@ -626,6 +628,8 @@ def get_de_energy_kev(v: np.ndarray, species: np.ndarray) -> NDArray:
         N x 3 array of velocity components (vx, vy, vz) in km/s.
     species : np.ndarray
         Species of the particle.
+    quality_flags : np.ndarray, optional
+        Quality flags to set when there is an outlier.
 
     Returns
     -------
@@ -648,6 +652,14 @@ def get_de_energy_kev(v: np.ndarray, species: np.ndarray) -> NDArray:
     energy[valid_mask] = (
         0.5 * UltraConstants.MASS_H * v2[valid_mask] * UltraConstants.J_KEV
     )
+    # Flag out of range energies
+    if quality_flags is not None:
+        energy_out_of_range = (energy < UltraConstants.PSET_ENERGY_BIN_EDGES[0]) | (
+            energy > UltraConstants.PSET_ENERGY_BIN_EDGES[-1]
+        )
+        quality_flags[energy_out_of_range] |= (
+            ImapDEOutliersUltraFlags.INVALID_ENERGY.value
+        )
 
     return energy