Binning prototype

ultra_user/energy_check/de_extended_calcs.py

@@ -0,0 +1,127 @@
+import numpy as np
+from numpy.typing import NDArray
+import logging
+import xarray
+import imap_processing.ultra.l1b.ultra_l1b_extended as l1b_ext
+import imap_processing.ultra.constants as constants
+from numpy import ndarray
+from numpy.typing import NDArray
+
+logger = logging.getLogger(__name__)
+
+
+# pull 1a stuff
+def get_xf(l1a_dset: xarray.Dataset) -> NDArray:
+    return l1b_ext.get_front_x_position(l1a_dset["START_TYPE"].data, l1a_dset["START_POS_TDC"].data)
+
+
+def get_1bdict(l1a_dset_uf: xarray.Dataset, sensor="ultra45") -> dict:
+    result = dict()
+    indices = np.nonzero(
+        np.isin(l1a_dset_uf["STOP_TYPE"], [l1b_ext.StopType.Top.value,
+                                           l1b_ext.StopType.Bottom.value])
+    )[0]
+    l1a_dset = l1a_dset_uf.isel(epoch=indices)
+
+    xf = get_xf(l1a_dset)
+    result["xf"] = xf
+    tof, t2, xb, yb = l1b_ext.get_ph_tof_and_back_positions(l1a_dset, xf, sensor)
+    result["tof"] = tof
+    result["t2"] = t2
+    result["xb"] = xb
+    result["yb"] = yb
+    d, yf = l1b_ext.get_front_y_position(l1a_dset["START_TYPE"].data, yb)
+    result["d"] = d
+    result["yf"] = yf
+    r = l1b_ext.get_path_length((xf, yf), (xb, yb), d)
+    result["r"] = r
+    ctof = l1b_ext.get_ctof(tof, r,"PH")
+    result["ctof"] = ctof
+    dmin = ndarray(ctof.size)
+    dmin[:] = constants.UltraConstants.DMIN_PH_CTOF
+    species = species_from_ctof(ctof)
+    result["species"] = species
+    v = de_velocity((xf, yf), (xb, yb), d, tof)
+    result["v"] = np.asarray(v)
+    v_mag = velocity_magnitude(ctof, dmin)
+    result["v_mag"] = v_mag
+    energy = de_energy_kev(v, species)
+    result["energy"] = energy
+    ih = np.where(species == "H")
+    result["ih"] = ih
+    ssd_indices = np.where(l1a_dset["STOP_TYPE"].data >= 8)[0]
+    ph_indices = np.where(l1a_dset["STOP_TYPE"].data < 8)[0]
+    result["issd"] = ssd_indices
+    result["ph"] = ph_indices
+    theta,phi = event_az_el(v)
+    result["event_theta"]=theta
+    result["event_phi"]=phi
+
+    return result
+
+
+# new stuff
+def de_velocity(
+        front_position: tuple[NDArray, NDArray],
+        back_position: tuple[NDArray, NDArray],
+        d: np.ndarray,
+        tof: np.ndarray,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    if tof[tof < 0].any():
+        logger.info("Negative tof values found.")
+
+    # distances in .1 mm
+    delta_x = .1 * ((front_position[0] - back_position[0]))
+    delta_y = .1 * ((front_position[1] - back_position[1]))
+    delta_z = .1 * d
+
+    v_x = mmperns_kmpers * delta_x / tof
+    v_y = mmperns_kmpers * delta_y / tof
+    v_z = mmperns_kmpers * delta_z / tof
+
+    return v_x, v_y, v_z
+
+
+def velocity_magnitude(ctof: ndarray, dmin: ndarray) -> np.ndarray:
+    cctof = ctof
+    cctof[ctof < 1] = 1
+    # ctof ihas units of 0.1ns
+    return 10 * mmperns_kmpers * dmin / cctof
+
+
+mH_kg = 1.6735575e-27
+Joule_to_keV = 6.242e+15
+mmperns_kmpers = 1000
+
+
+def de_energy_kev(v: tuple[NDArray, NDArray, NDArray],
+                  species: NDArray) -> NDArray:
+    vv = np.asarray(v) * 1.e3  # convert km/s to m/s
+    v2 = np.sum((vv * vv), 0)
+
+    iH = np.where(species == "H")
+    result = np.full_like(v2, np.nan)
+
+    result[iH] = 0.5 * mH_kg * v2[iH] * Joule_to_keV
+    return result
+
+
+def species_from_ctof(ctof: ndarray) -> ndarray:
+    result = np.ndarray(ctof.size, dtype=object)
+    result[:] = "unknown"
+    result[np.where(np.logical_and(ctof > 50, ctof < 200))] = "H"
+    return result
+
+
+def event_az_el(v: tuple[NDArray, NDArray, NDArray]) -> tuple[NDArray, NDArray]:
+    vv = np.asarray(v)
+    vmag = np.sqrt(np.sum((vv * vv), 0))
+
+    ux = vv[0, :] / vmag
+    uy = vv[1, :] / vmag
+    uz = vv[2, :] / vmag
+
+    theta = np.arccos(uz / vmag)
+    phi = np.arctan2(uy, ux)
+
+    return theta, phi

ultra_user/energy_check/echeck_fragments.py

@@ -0,0 +1,96 @@
+import sys
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import importlib
+
+import imap_processing.ultra.constants as constants
+import ultra_user.energy_check.event_dataset as ed
+import ultra_user.pipeline.test_data as td
+import imap_processing.ultra.l1b.ultra_l1b_extended as l1b_ext
+import ultra_user.energy_check.de_extended_calcs as de_calcs
+plt.interactive(True)
+
+importlib.reload(ed)
+importlib.reload(td)
+importlib.reload(de_calcs)
+
+d = ed.EventDataset()
+dtest = td.de_dataset()
+
+l1b=de_calcs.get_1bdict(dtest)
+v = np.asarray(l1b["v"])
+v_mag = l1b["v_mag"]
+
+v_mag0=np.sqrt(np.sum((v*v),0))
+
+ehist,echan = np.histogram(l1b["energy"],bins=20,range = (20,60))
+vhist,vchan = np.histogram(v_mag,bins=20,range = (100,5000))
+
+plt.plot(vchan[1:],vhist)
+plt.xlabel("Velocity (km/s)")
+plt.title("Histogram of test dataset")
+plt.show()
+
+plt.plot(echan[1:],ehist)
+plt.xlabel("Energy (keV)")
+plt.title("Histogram of test dataset")
+plt.show()
+
+
+d.ctof_ph_plot()
+
+v0 = d.get_v_ctof()
+v1 = d.get_v_r()
+
+
+################# older fragments
+ctof = d.get_cTOF()
+ph = d.get_ph_or_e()
+bin = d.get_cTOF_bin()
+f = pd.read_csv(testcsv)
+df_filt = df[df["StartType"] != -1]
+
+cTOF = df_filt["cTOF"].astype("float").values
+bin = df_filt["ComputedBin"].astype("float").values
+e_xx = np.array(df_filt["Energy"].astype("float"))
+e_yy = np.array(df_filt["EnergyOrPH"].astype("float"))
+r = np.array(df_filt["r"].astype("float"))
+dmin = constants.UltraConstants.DMIN
+
+v = dmin/cTOF *1.e7  #(unit conversion mm/tenths of a ns)
+amu = 1.66054e-27
+
+Ej = .05*amu*v*v
+J_per_Kev = 1.60218e-16
+Ekev = Ej/J_per_Kev
+hbins = 10**(np.arange(20)/20)
+Ehist,edges = np.histogram(Ekev,bins=hbins)
+binhist,binedge = np.histogram(bin,bins=np.arange(64))
+
+bina_index = np.where(np.logical_and( bin >5 ,bin < 10))
+binb_index = np.where(np.logical_and( bin >21 ,bin < 27))
+binc_index = np.where(np.logical_and( bin >40 ,bin < 50))
+
+Ehista,edges = np.histogram(Ekev[bina_index],bins=hbins)
+Ehistb,edges = np.histogram(Ekev[binb_index],bins=hbins)
+Ehistc,edges = np.histogram(Ekev[binc_index],bins=hbins)
+
+
+plt.plot(np.sqrt(hbins[0:-1]*hbins[1:]),Ehista)
+plt.plot(np.sqrt(hbins[0:-1]*hbins[1:]),Ehistb)
+plt.plot(np.sqrt(hbins[0:-1]*hbins[1:]),Ehistc)
+plt.xscale('log')
+plt.yscale('log')
+plt.xlabel("Energy")
+plt.ylabel("Number of events")
+plt.show()
+
+plt.plot(Ekev[bina_index],e_xx[bina_index],"o")
+plt.xscale('log')
+plt.yscale('log')
+plt.xlabel("Ekev")
+plt.ylabel("E_ph")
+plt.show()

ultra_user/energy_check/event_dataset.py

@@ -0,0 +1,62 @@
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import imap_processing.ultra.constants as constants
+
+
+class EventDataset:
+    def __init__(self, clean=True,
+                 csv="imap_processing/tests/ultra/test_data/l0/"+
+                     "ultra45_raw_sc_ultrarawimg_withFSWcalcs_FM45_40P_Phi28p5_BeamCal_LinearScan_phi2850_theta-000_20240207T102740.csv"):
+        df_raw = pd.read_csv(csv)
+        if clean:
+            self.df = df_raw[df_raw["StartType"] != -1]
+        else:
+            self.df = df_raw
+
+
+    def get_cTOF(self):
+        return np.array(self.df["cTOF"].astype("float").values)
+
+    def get_ph_or_e(self):
+        return np.array(self.df["EnergyOrPH"].astype("float"))
+
+    def get_cTOF_bin(self):
+        return np.array(self.df["ComputedBin"].astype("float"))
+
+    def get_v_ctof(self):
+        cTOF = self.get_cTOF()
+        issd = self.ssd_indices()
+        dmin_mcp = constants.UltraConstants.DMIN # will change to DMIN_PH_CTOF
+        dmin_ssd = constants.UltraConstants.DMIN_SSD_CTOF
+
+        v = dmin_mcp / cTOF
+        v[issd] = dmin_ssd/cTOF[issd]
+        return v * 1.e7  # (unit conversion mm/tenths of a ns)
+
+    def get_v_r(self):
+        r = np.array(self.df["r"].astype("float"))
+        tof = np.array(self.df["TOF"].astype("float"))
+        return (r / tof) * 1.e5
+
+    def get_stop_type(self):
+        return np.array(self.df["StopType"].astype("int"))
+
+    def ssd_indices(self):
+        stop_type = self.get_stop_type()
+        return np.where(stop_type >2)
+    def mcp_indices(self):
+        stop_type = self.get_stop_type()
+        return np.where(stop_type <=2)
+
+    def ctof_ph_plot(self,binlim=50, xr=[0, 300], cmap='viridis', vmin=None, vmax=None):
+        ctof = self.get_cTOF()
+        ph = self.get_ph_or_e()
+        bin = self.get_cTOF_bin()
+        ibin = np.where(bin < binlim)
+        plt.scatter(ctof[ibin], ph[ibin], c=bin[ibin], cmap=cmap, vmin=vmin, vmax=vmax)
+        plt.colorbar()
+        plt.xlim(xr)
+        plt.xlabel("cTOF")
+        plt.ylabel("Pulse Height")
+        plt.show()

ultra_user/energy_check/extended_1b_picklescript.py

@@ -0,0 +1,25 @@
+import ultra_user.energy_check.event_dataset as ed
+import ultra_user.pipeline.test_data as td
+import ultra_user.energy_check.de_extended_calcs as de_calcs
+import pickle
+
+d_csv = ed.EventDataset()
+l1a = td.de_dataset()
+
+l1b=de_calcs.get_1bdict(l1a)
+#pkldir = '/Users/demajr1/tmp/'
+pkldir = 'data/pkls/'
+
+
+pkfile = open(pkldir+'dataset_from_csv.pkl',"wb")
+pickle.dump(d_csv,pkfile)
+pkfile.close()
+
+pkfile = open(pkldir+'dataset_l1a.pkl',"wb")
+pickle.dump(l1a,pkfile)
+pkfile.close()
+
+pkfile = open(pkldir+'dataset_l1b.pkl',"wb")
+pickle.dump(l1b,pkfile)
+pkfile.close()
+

ultra_user/mapbin/binfragments.py

@@ -0,0 +1,57 @@
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+import ultra_user.mapbin.lonLatscheme as lls
+import ultra_user.mapbin.hpscheme as hps
+import ultra_user.mapbin.mapBinner as mapBinner
+import importlib
+import healpy as hp
+import pickle
+
+#matplotlib.use('TkAgg')
+#%matplotlib inline
+plt.interactive(True)
+importlib.reload(hps)
+importlib.reload(lls)
+
+binner0 = lls.LonLatScheme()
+
+b00 = binner0.bin_number(0,0)
+
+p00 = binner0.center_position(b00)
+c00 = binner0.boundary_position(b00)
+
+binner1 = hps.HpScheme()
+b11 = binner1.bin_number(0,0)
+p00 = binner1.center_position(b11)
+
+
+# get L1B data in
+pkl = open('data/pkls/dataset_l1b.pkl','rb')
+l1b = pickle.load(pkl)
+pkl.close()
+
+v1b =l1b['v']
+energy = l1b['energy']
+vmag = l1b['v_mag']
+
+
+az = np.degrees(np.atan2(v1b[1,:],v1b[0,:]))
+el = np.degrees(np.arccos(v1b[2,:]/vmag))
+lat = 90-el
+
+mp0 = mapBinner.MapBinner(binner0)
+mp1 = mapBinner.MapBinner(binner1)
+for ic,oneaz in enumerate(az):
+    mp0.add_value(1,oneaz,lat[ic])
+    mp1.add_value(1,oneaz,lat[ic])
+
+llmap = binner0.asGrid(mp0.bins)
+
+plt.imshow(np.transpose(llmap))
+plt.colorbar()
+plt.show()
+
+hp.cartview(mp1.bins)
+hp.graticule()
+plt.show()