Ultra particle velocity and SSD TOF (#823)

* added particle velocity calcs and ssd tof

imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py

@@ -12,9 +12,11 @@
     get_coincidence_positions,
     get_front_x_position,
     get_front_y_position,
+    get_particle_velocity,
     get_path_length,
     get_ph_tof_and_back_positions,
     get_ssd_back_position_and_tof_offset,
+    get_ssd_tof,
 )
 
 
@@ -204,3 +206,54 @@ def test_calculate_etof_xc(de_dataset, yf_fixture):
     np.testing.assert_allclose(
         etof_bottom, df_bottom["eTOF"].astype("float").values, atol=1e-06, rtol=0
     )
+
+
+def test_get_particle_velocity(de_dataset, yf_fixture):
+    """Tests get_particle_velocity function."""
+    df_filt, _, _ = yf_fixture
+
+    ph_indices = np.nonzero(
+        np.isin(de_dataset["STOP_TYPE"], [StopType.Top.value, StopType.Bottom.value])
+    )[0]
+
+    ph_rows = df_filt.iloc[ph_indices]
+    test_xf = ph_rows["Xf"].astype("float").values
+    test_yf = ph_rows["Yf"].astype("float").values
+    test_xb = ph_rows["Xb"].astype("float").values
+    test_yb = ph_rows["Yb"].astype("float").values
+    test_d = ph_rows["d"].astype("float").values
+    test_tof = ph_rows["TOF"].astype("float").values
+
+    vhat_x, vhat_y, vhat_z = get_particle_velocity(
+        (test_xf, test_yf),
+        (test_xb, test_yb),
+        test_d,
+        test_tof,
+    )
+    # FSW test data should be negative and not have an analysis
+    # for negative tof values.
+    assert vhat_x[test_tof > 0] == pytest.approx(
+        -df_filt["vhatX"].iloc[ph_indices].astype("float").values[test_tof > 0],
+        rel=1e-2,
+    )
+    assert vhat_y[test_tof > 0] == pytest.approx(
+        -df_filt["vhatY"].iloc[ph_indices].astype("float").values[test_tof > 0],
+        rel=1e-2,
+    )
+    assert vhat_z[test_tof > 0] == pytest.approx(
+        -df_filt["vhatZ"].iloc[ph_indices].astype("float").values[test_tof > 0],
+        rel=1e-2,
+    )
+
+
+def test_get_ssd_tof(de_dataset, yf_fixture):
+    """Tests get_ssd_tof function."""
+    df_filt, _, _ = yf_fixture
+    df_ssd = df_filt[df_filt["StopType"].isin(StopType.SSD.value)]
+    test_xf = df_filt["Xf"].astype("float").values
+
+    ssd_tof = get_ssd_tof(de_dataset, test_xf)
+
+    np.testing.assert_allclose(
+        ssd_tof, df_ssd["TOF"].astype("float"), atol=1e-05, rtol=0
+    )

imap_processing/ultra/l1b/ultra_l1b_extended.py

@@ -1,11 +1,13 @@
 """Calculates Extended Raw Events for ULTRA L1b."""
 
+import logging
 from enum import Enum
 from typing import ClassVar
 
 import numpy as np
 import xarray
 from numpy import ndarray
+from numpy.typing import NDArray
 
 from imap_processing.ultra.l1b.lookup_utils import (
     get_back_position,
@@ -14,6 +16,8 @@
     get_y_adjust,
 )
 
+logger = logging.getLogger(__name__)
+
 # Constants in IMAP-Ultra Flight Software Specification document.
 D_SLIT_FOIL = 3.39  # shortest distance from slit to foil (mm)
 SLIT_Z = 44.89  # position of slit on Z axis (mm)
@@ -22,6 +26,7 @@
 N_ELEMENTS = 256  # number of elements in lookup table
 TRIG_CONSTANT = 81.92  # trigonometric constant (mm)
 # TODO: make lookup tables into config files.
+# TODO: put logic from Ultra FSW in here.
 
 
 class StartType(Enum):
@@ -222,9 +227,10 @@ def get_ph_tof_and_back_positions(
     t2[stop_type_top] = get_image_params("TOFSC") * t1[
         stop_type_top
     ] + get_image_params("TOFTPOFF")
-    tof[stop_type_top] = t2[stop_type_top] + xf_ph[stop_type_top] * get_image_params(
-        "XFTTOF"
-    )
+    # Variable xf_ph divided by 10 to convert to mm.
+    tof[stop_type_top] = t2[stop_type_top] + xf_ph[
+        stop_type_top
+    ] / 10 * get_image_params("XFTTOF")
 
     stop_type_bottom = de_filtered["STOP_TYPE"].data == StopType.Bottom.value
     xb[stop_type_bottom] = get_back_position(
@@ -239,9 +245,10 @@ def get_ph_tof_and_back_positions(
         stop_type_bottom
     ] + get_image_params("TOFBTOFF")  # 10*ns
 
+    # Variable xf_ph divided by 10 to convert to mm.
     tof[stop_type_bottom] = t2[stop_type_bottom] + xf_ph[
         stop_type_bottom
-    ] * get_image_params("XFTTOF")
+    ] / 10 * get_image_params("XFTTOF")
 
     return tof, t2, xb, yb
 
@@ -426,3 +433,111 @@ def get_coincidence_positions(
 
     # Convert to hundredths of a millimeter by multiplying times 100
     return etof, xc_array * 100
+
+
+def get_particle_velocity(
+    front_position: tuple[float, float],
+    back_position: tuple[float, float],
+    d: np.ndarray,
+    tof: np.ndarray,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Determine the particle velocity.
+
+    The equation is: velocity = ((xf - xb), (yf - yb), d).
+
+    Further description is available on pages 39 of
+    IMAP-Ultra Flight Software Specification document
+    (7523-9009_Rev_-.pdf).
+
+    Parameters
+    ----------
+    front_position : tuple
+        Front position (xf,yf) (hundredths of a millimeter).
+    back_position : tuple
+        Back position (xb,yb) (hundredths of a millimeter).
+    d : np.array
+        Distance from slit to foil (hundredths of a millimeter).
+    tof : np.array
+        Time of flight (tenths of a nanosecond).
+
+    Returns
+    -------
+    vhat_x : np.array
+        Normalized component of the velocity vector in x direction.
+    vhat_y : np.array
+        Normalized component of the velocity vector in y direction.
+    vhat_z : np.array
+        Normalized component of the velocity vector in z direction.
+    """
+    if tof[tof < 0].any():
+        logger.info("Negative tof values found.")
+
+    delta_x = front_position[0] - back_position[0]
+    delta_y = front_position[1] - back_position[1]
+
+    v_x = delta_x / tof
+    v_y = delta_y / tof
+    v_z = d / tof
+
+    # Magnitude of the velocity vector
+    magnitude_v = np.sqrt(v_x**2 + v_y**2 + v_z**2)
+
+    vhat_x = -v_x / magnitude_v
+    vhat_y = -v_y / magnitude_v
+    vhat_z = -v_z / magnitude_v
+
+    vhat_x[tof < 0] = np.iinfo(np.int64).min  # used as fillvals
+    vhat_y[tof < 0] = np.iinfo(np.int64).min
+    vhat_z[tof < 0] = np.iinfo(np.int64).min
+
+    return vhat_x, vhat_y, vhat_z
+
+
+def get_ssd_tof(de_dataset: xarray.Dataset, xf: np.ndarray) -> NDArray[np.float64]:
+    """
+    Calculate back xb, yb position for the SSDs.
+
+    An incoming particle could miss the stop anodes and instead
+    hit one of the SSDs between the anodes. Which SSD is hit
+    gives a coarse measurement of the y back position;
+    the x back position will be fixed.
+
+    Before hitting the SSD, particles pass through the stop foil;
+    dislodged electrons are accelerated back towards the coincidence anode.
+    The Coincidence Discrete provides a measure of the TOF.
+    A scale factor and offsets, and a multiplier convert xf to a tof offset.
+
+    Further description is available on pages 36 of
+    IMAP-Ultra Flight Software Specification document
+    (7523-9009_Rev_-.pdf).
+
+    Parameters
+    ----------
+    de_dataset : xarray.Dataset
+        Data in xarray format.
+    xf : np.array
+        Front x position (hundredths of a millimeter).
+
+    Returns
+    -------
+    tof : np.ndarray
+        Time of flight (tenths of a nanosecond).
+    """
+    _, tof_offset, ssd_number = get_ssd_back_position_and_tof_offset(de_dataset)
+    indices = np.nonzero(np.isin(de_dataset["STOP_TYPE"], [StopType.SSD.value]))[0]
+
+    de_discrete = de_dataset.isel(epoch=indices)["COIN_DISCRETE_TDC"]
+
+    time = get_image_params("TOFSSDSC") * de_discrete.values + tof_offset
+
+    # The scale factor and offsets, and a multiplier to convert xf to a tof offset.
+    # Convert xf to mm by dividing by 100.
+    tof = (
+        time
+        + get_image_params("TOFSSDTOTOFF")
+        + xf[indices] / 100 * get_image_params("XFTTOF")
+    ) * 10
+
+    # Convert TOF to tenths of a nanosecond.
+    return np.asarray(tof, dtype=np.float64)

imap_processing/ultra/lookup_tables/FM45_Startup1_ULTRA_IMGPARAMS_20240719.csv

@@ -1,2 +1,2 @@
 SHCOARSE,XFTSC,XFTLTOFF,XFTRTOFF,TOFSC,TOFTPOFF,TOFBTOFF,XFTTOF,XCOINTPSC,XCOINTPOFF,XCOINBTSC,XCOINBTOFF,ETOFSC,ETOFTPOFF,ETOFBTOFF,TOFDIFFTPMIN,TOFDIFFTPMAX,TOFDIFFBTMIN,TOFDIFFBTMAX,ETOFMIN,ETOFMAX,ETOFSLOPE1,ETOFOFF1,ETOFSLOPE2,ETOFOFF2,SPTPPHOFF,SPBTPHOFF,YBKSSD0,YBKSSD1,YBKSSD2,YBKSSD3,YBKSSD4,YBKSSD5,YBKSSD6,YBKSSD7,TOFSSDSC,TOFSSDLTOFF0,TOFSSDLTOFF1,TOFSSDLTOFF2,TOFSSDLTOFF3,TOFSSDLTOFF4,TOFSSDLTOFF5,TOFSSDLTOFF6,TOFSSDLTOFF7,TOFSSDRTOFF0,TOFSSDRTOFF1,TOFSSDRTOFF2,TOFSSDRTOFF3,TOFSSDRTOFF4,TOFSSDRTOFF5,TOFSSDRTOFF6,TOFSSDRTOFF7,TOFSSDTOTOFF,PATHSTEEPTHRESH,PATHMEDIUMTHRESH
-,0.172998047,49.3,48.25,0.5,-528,-525,0.001831055,0.067929688,41.75,0.067929688,-40.75,0.1,-44.5,-44.5,,,,,,,,,,,,,29.3,37.3,7.1,15.1,-15.1,-7.1,-37.3,-29.3,0.196484375,-6,-7.3,-3.8,-4.2,-3.8,-3.7,-6.3,-5,-5,-6.3,-3.7,-3.8,-4,-4.2,-7.3,-6,5.9,,
+,0.172998047,49.3,48.25,0.5,-528,-525,0.01831055,0.067929688,41.75,0.067929688,-40.75,0.1,-44.5,-44.5,,,,,,,,,,,,,29.3,37.3,7.1,15.1,-15.1,-7.1,-37.3,-29.3,0.196484375,-6,-7.3,-3.8,-4.2,-3.8,-3.7,-6.3,-5,-5,-6.3,-3.7,-3.8,-4,-4.2,-7.3,-6,5.9,,