diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py b/imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py index 25f0b7f76..daccc1377 100644 --- a/imap_processing/tests/ultra/unit/test_ultra_l1b_extended.py +++ b/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 + ) diff --git a/imap_processing/ultra/l1b/ultra_l1b_extended.py b/imap_processing/ultra/l1b/ultra_l1b_extended.py index 9bc1dab4e..24861ca93 100644 --- a/imap_processing/ultra/l1b/ultra_l1b_extended.py +++ b/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) diff --git a/imap_processing/ultra/lookup_tables/FM45_Startup1_ULTRA_IMGPARAMS_20240719.csv b/imap_processing/ultra/lookup_tables/FM45_Startup1_ULTRA_IMGPARAMS_20240719.csv index ffae6b5c7..fc7a1ed63 100644 --- a/imap_processing/ultra/lookup_tables/FM45_Startup1_ULTRA_IMGPARAMS_20240719.csv +++ b/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,, \ No newline at end of file +,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,, \ No newline at end of file