Ultra l1c pset sensitivity calculation and minor organization (IMAP-Science-Operations-Center#1604)

Author: laspsandoval

imap_processing/cdf/config/imap_ultra_l1b_variable_attrs.yaml

@@ -250,13 +250,43 @@ theta:
   LABLAXIS: theta
   UNITS: degrees
 
+phi_fwhm:
+  <<: *default_float32
+  CATDESC: FWHM of the phi distribution.
+  FIELDNAM: phi_fwhm
+  LABLAXIS: phi fwhm
+  UNITS: degrees
+
+theta_fwhm:
+  <<: *default_float32
+  CATDESC: FWHM of the theta distribution.
+  FIELDNAM: phi_fwhm
+  LABLAXIS: phi fwhm
+  UNITS: degrees
+
 event_efficiency:
   <<: *default_float32
   CATDESC: Estimated event efficiency for this path through the instrument.
   FIELDNAM: event_efficiency
   LABLAXIS: event efficiency
   UNITS: " "
 
+direct_event_unit_velocity:
+  <<: *default_float32
+  CATDESC: Unit vector of the direct event velocity in the instrument frame.
+  DEPEND_1: component
+  FIELDNAM: direct_event_unit_velocity
+  LABLAXIS: direct event unit velocity
+  UNITS: " "
+
+direct_event_unit_position:
+  <<: *default_float32
+  CATDESC: Unit vector of the direct event position in the instrument frame.
+  DEPEND_1: component
+  FIELDNAM: direct_event_unit_position
+  LABLAXIS: direct event unit position
+  UNITS: " "
+
 velocity_sc:
   <<: *default_float32
   CATDESC: x,y,z-components of the velocity vector of the ena in the spacecraft frame.

imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml

@@ -31,22 +31,6 @@ longitude:
   LABLAXIS: longitude
   UNITS: degrees
 
-spin_angle:
-  <<: *default
-  CATDESC: Angular bin defined by central location of the bin
-  FIELDNAM: spin_angle
-  LABLAXIS: spin angle
-  # TODO: come back to format
-  UNITS: degrees
-
-time_bin:
-  <<: *default
-  CATDESC: Time bin defined by central location of the bin.
-  FIELDNAM: time_bin
-  LABLAXIS: time bin
-  # TODO: come back to format
-  UNITS: seconds
-
 species_bin:
   <<: *default
   CATDESC: Species bin defined by central location.
@@ -95,30 +79,6 @@ sid:
   # TODO: come back to format
   UNITS: " "
 
-row:
-  <<: *default
-  CATDESC: Row number.
-  FIELDNAM: row
-  LABLAXIS: row
-  # TODO: come back to format
-  UNITS: " "
-
-column:
-  <<: *default
-  CATDESC: Column number.
-  FIELDNAM: column
-  LABLAXIS: column
-  # TODO: come back to format
-  UNITS: " "
-
-spin:
-  <<: *default
-  CATDESC: Spin number at integration start.
-  FIELDNAM: spin
-  LABLAXIS: spin
-  # TODO: come back to format
-  UNITS: " "
-
 shcoarse:
   <<: *default
   CATDESC: CCSDS packet timestamp.
@@ -127,14 +87,6 @@ shcoarse:
   # TODO: come back to format
   UNITS: seconds
 
-packetdata:
-  <<: *default
-  CATDESC: A 2D numpy array representing pixel values. Each pixel is stored as an unsigned 16-bit integer.
-  FIELDNAM: packetdata
-  LABLAXIS: packet data
-  # TODO: come back to format
-  UNITS: " "
-
 energy_bin_delta:
   <<: *default_float32
   CATDESC: Difference between the energy bin edges.

imap_processing/tests/ultra/unit/test_de.py

@@ -113,10 +113,3 @@ def test_calculate_de(l1b_de_dataset, df_filt):
         len(l1b_de_dataset["epoch"]),
         3,
     )
-
-    # Event efficiency
-    assert np.allclose(
-        l1b_de_dataset["event_efficiency"],
-        np.full(len(l1b_de_dataset["epoch"]), np.nan),
-        equal_nan=True,
-    )

imap_processing/tests/ultra/unit/test_lookup_utils.py

@@ -66,7 +66,7 @@ def test_get_egy_norm():
 def test_get_image_params():
     """Tests function get_image_params."""
 
-    image_params = get_image_params("XFTLTOFF")
+    image_params = get_image_params("XFTLTOFF", "ultra45")
 
     assert image_params == 49.3
 

imap_processing/tests/ultra/unit/test_spacecraft_pset.py

@@ -30,7 +30,7 @@ def test_pset():
     # This is just setting up the data so that it is in the format of l1b_de_dataset.
     test_path = TEST_PATH / "ultra-90_raw_event_data_shortened.csv"
     df = pd.read_csv(test_path)
-    instrument_velocity = get_de_velocity(
+    instrument_velocity, _, _ = get_de_velocity(
         (df["Xf"], df["Yf"]), (df["Xb"], df["Yb"]), df["d"], df["TOF"].values
     )
 

imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py

@@ -64,6 +64,7 @@ def test_get_front_x_position(
     xf = get_front_x_position(
         de_dataset["start_type"].data,
         de_dataset["start_pos_tdc"].data,
+        "ultra45",
     )
 
     assert xf == pytest.approx(df_filt["Xf"].astype("float"), 1e-5)
@@ -121,7 +122,9 @@ def test_get_ssd_back_position_and_tof_offset(
 ):
     """Tests get_ssd_back_position function."""
     _, _, _, de_dataset = test_fixture
-    yb, tof_offset, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset)
+    yb, tof_offset, ssd_number = get_ssd_back_position_and_tof_offset(
+        de_dataset, "ultra45"
+    )
 
     df = pd.read_csv(events_fsw_comparison_theta_0)
     df_filt = df[(df["StartType"] != -1) & (df["StopType"] >= 8)]
@@ -238,7 +241,7 @@ def test_get_de_velocity(test_fixture):
         for col in ["Xf", "Yf", "Xb", "Yb", "d", "TOF"]
     )
 
-    v = get_de_velocity(
+    v, vhat, r = get_de_velocity(
         (test_xf, test_yf),
         (test_xb, test_yb),
         test_d,
@@ -265,6 +268,42 @@ def test_get_de_velocity(test_fixture):
         atol=1e-01,
         rtol=0,
     )
+    np.testing.assert_allclose(
+        vhat[test_tof > 0][:, 0],
+        df_ph["vhatX"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
+    np.testing.assert_allclose(
+        vhat[test_tof > 0][:, 1],
+        df_ph["vhatY"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
+    np.testing.assert_allclose(
+        vhat[test_tof > 0][:, 2],
+        df_ph["vhatZ"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
+    np.testing.assert_allclose(
+        r[test_tof > 0][:, 0],
+        -df_ph["vhatX"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
+    np.testing.assert_allclose(
+        r[test_tof > 0][:, 1],
+        -df_ph["vhatY"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
+    np.testing.assert_allclose(
+        r[test_tof > 0][:, 2],
+        -df_ph["vhatZ"].astype("float").values[test_tof > 0],
+        atol=1e-01,
+        rtol=0,
+    )
 
 
 def test_get_ssd_tof(test_fixture):
@@ -273,7 +312,7 @@ def test_get_ssd_tof(test_fixture):
     df_ssd = df_filt[np.isin(df_filt["StopType"], [StopType.SSD.value])]
     test_xf = df_filt["Xf"].astype("float").values
 
-    ssd_tof = get_ssd_tof(de_dataset, test_xf)
+    ssd_tof = get_ssd_tof(de_dataset, test_xf, "ultra45")
 
     np.testing.assert_allclose(
         ssd_tof, df_ssd["TOF"].astype("float"), atol=1e-05, rtol=0
@@ -296,7 +335,7 @@ def test_get_de_energy_kev(test_fixture):
         for col in ["Xf", "Yf", "Xb", "Yb", "d", "TOF"]
     )
 
-    v = get_de_velocity(
+    v, v_hat, r_hat = get_de_velocity(
         (test_xf, test_yf),
         (test_xb, test_yb),
         test_d,
@@ -315,7 +354,7 @@ def test_get_energy_ssd(test_fixture):
     """Tests get_energy_ssd function."""
     df_filt, _, _, de_dataset = test_fixture
     df_ssd = df_filt[np.isin(df_filt["StopType"], [StopType.SSD.value])]
-    _, _, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset)
+    _, _, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset, "ultra45")
     energy = get_energy_ssd(de_dataset, ssd_number)
     test_energy = df_ssd["Energy"].astype("float")
 
@@ -334,7 +373,11 @@ def test_get_energy_pulse_height(test_fixture):
     test_yb = df_filt["Yb"].astype("float").values
 
     energy = get_energy_pulse_height(
-        de_dataset["stop_type"].data, de_dataset["energy_ph"].data, test_xb, test_yb
+        de_dataset["stop_type"].data,
+        de_dataset["energy_ph"].data,
+        test_xb,
+        test_yb,
+        "ultra45",
     )
     test_energy = df_ph["Energy"].astype("float")
 

imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py

@@ -16,6 +16,7 @@
     get_spacecraft_exposure_times,
     get_spacecraft_histogram,
     get_spacecraft_sensitivity,
+    grid_sensitivity,
 )
 
 BASE_PATH = imap_module_directory / "ultra" / "lookup_tables"
@@ -171,9 +172,12 @@ def test_get_spacecraft_sensitivity():
     df_efficiencies = pd.read_csv(efficiences)
     df_geometric_function = pd.read_csv(geometric_function)
 
-    sensitivity = get_spacecraft_sensitivity(df_efficiencies, df_geometric_function)
+    sensitivity, energy_vals, right_ascension, declination = get_spacecraft_sensitivity(
+        df_efficiencies, df_geometric_function
+    )
 
-    assert sensitivity.shape == df_efficiencies.shape
+    assert sensitivity.shape == (df_efficiencies.shape[0], df_efficiencies.shape[1] - 2)
+    assert np.array_equal(energy_vals, np.arange(3.0, 80.5, 0.5))
 
     df_efficiencies_test = pd.DataFrame(
         {"3.0keV": [1.0, 2.0], "3.5keV": [3.0, 4.0], "4.0keV": [5.0, 6.0]}
@@ -192,3 +196,26 @@ def test_get_spacecraft_sensitivity():
     assert np.allclose(
         df_sensitivity_test.to_numpy(), expected_sensitivity.to_numpy(), atol=1e-6
     )
+
+
+@pytest.mark.external_test_data
+def test_grid_sensitivity():
+    """Tests grid_sensitivity function."""
+    efficiencies_path = TEST_PATH / "Ultra_90_DPS_efficiencies_all.csv"
+    geometric_function_path = TEST_PATH / "ultra_90_dps_gf.csv"
+
+    df_efficiencies = pd.read_csv(efficiencies_path)
+    df_geometric_function = pd.read_csv(geometric_function_path)
+
+    sensitivity, energy_vals, ra, dec = get_spacecraft_sensitivity(
+        df_efficiencies, df_geometric_function
+    )
+
+    expected_result = sensitivity["3.0keV"].values
+    result = grid_sensitivity(df_efficiencies, df_geometric_function, 3.0)
+
+    assert np.allclose(result, expected_result, atol=1e-5)
+
+    # Check that out-of-bounds energy returns all NaNs
+    result = grid_sensitivity(df_efficiencies, df_geometric_function, 2.5)
+    assert np.isnan(result).all()

imap_processing/ultra/l1b/de.py

@@ -20,6 +20,7 @@
     get_eventtimes,
     get_front_x_position,
     get_front_y_position,
+    get_fwhm,
     get_path_length,
     get_ph_tof_and_back_positions,
     get_phi_theta,
@@ -104,11 +105,14 @@ def calculate_de(de_dataset: xr.Dataset, name: str) -> xr.Dataset:
     event_times = np.full(len(de_dataset["epoch"]), np.nan, dtype=np.float64)
     spin_starts = np.full(len(de_dataset["epoch"]), np.nan, dtype=np.float64)
     spin_period_sec = np.full(len(de_dataset["epoch"]), np.nan, dtype=np.float64)
+    start_type = np.full(len(de_dataset["epoch"]), np.nan, dtype=np.uint8)
 
     xf[valid_indices] = get_front_x_position(
         de_dataset["start_type"].data[valid_indices],
         de_dataset["start_pos_tdc"].data[valid_indices],
+        f"ultra{sensor}",
     )
+    start_type[valid_indices] = de_dataset["start_type"].data[valid_indices]
 
     (
         event_times[valid_indices],
@@ -131,6 +135,7 @@ def calculate_de(de_dataset: xr.Dataset, name: str) -> xr.Dataset:
         de_dataset["energy_ph"].data[ph_indices],
         xb[ph_indices],
         yb[ph_indices],
+        f"ultra{sensor}",
     )
     r[ph_indices] = get_path_length(
         (xf[ph_indices], yf[ph_indices]),
@@ -151,8 +156,10 @@ def calculate_de(de_dataset: xr.Dataset, name: str) -> xr.Dataset:
     )
 
     # SSD
-    tof[ssd_indices] = get_ssd_tof(de_dataset, xf)
-    yb[ssd_indices], _, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset)
+    tof[ssd_indices] = get_ssd_tof(de_dataset, xf, f"ultra{sensor}")
+    yb[ssd_indices], _, ssd_number = get_ssd_back_position_and_tof_offset(
+        de_dataset, f"ultra{sensor}"
+    )
     xc[ssd_indices] = np.zeros(len(ssd_indices))
     xb[ssd_indices] = np.zeros(len(ssd_indices))
     etof[ssd_indices] = np.zeros(len(ssd_indices))
@@ -195,13 +202,15 @@ def calculate_de(de_dataset: xr.Dataset, name: str) -> xr.Dataset:
     de_dict["phi"] = phi
     de_dict["theta"] = theta
 
-    v = get_de_velocity(
+    v, vhat, r = get_de_velocity(
         (de_dict["x_front"], de_dict["y_front"]),
         (de_dict["x_back"], de_dict["y_back"]),
         de_dict["front_back_distance"],
         de_dict["tof_start_stop"],
     )
     de_dict["direct_event_velocity"] = v.astype(np.float32)
+    de_dict["direct_event_unit_velocity"] = vhat.astype(np.float32)
+    de_dict["direct_event_unit_position"] = r.astype(np.float32)
 
     de_dict["tof_energy"] = get_de_energy_kev(v, species_bin)
     de_dict["energy"] = energy
@@ -224,9 +233,12 @@ def calculate_de(de_dataset: xr.Dataset, name: str) -> xr.Dataset:
     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)
 
-    # TODO: TBD.
-    de_dict["event_efficiency"] = np.full(
-        len(de_dataset["epoch"]), np.nan, dtype=np.float32
+    de_dict["phi_fwhm"], de_dict["theta_fwhm"] = get_fwhm(
+        start_type,
+        f"ultra{sensor}",
+        de_dict["tof_energy"],
+        de_dict["phi"],
+        de_dict["theta"],
     )
 
     dataset = create_dataset(de_dict, name, "l1b")