Adds livetime correction to get data

stixpy/calibration/detector.py

@@ -22,7 +22,9 @@ def get_srm():
     -------
 
     """
-    drm_save = read_genx("/Users/shane/Projects/STIX/git/stix_drm_20220713.genx")
+    # drm_save = read_genx("/Users/shane/Projects/STIX/git/stix_drm_20220713.genx")
+    drm_save = np.load('/home/jmitchell/software/stixpy-dev/stixpy/config/data/detector/')
+
     drm = drm_save["SAVEGEN0"]["SMATRIX"] * u.count / u.keV / u.photon
     energies_in = drm_save["SAVEGEN0"]["EDGES_IN"] * u.keV
     energies_in_width = np.diff(energies_in)

stixpy/calibration/elut.py

@@ -0,0 +1,87 @@
+from pathlib import Path
+
+import numpy as np
+
+from stixpy.calibration.detector import get_sci_channels
+from stixpy.io.readers import read_elut, read_elut_index
+
+__all__ = ["get_elut", "get_elut_correction"]
+
+def get_elut(date):
+    r"""
+    Get the energy lookup table (ELUT) for the given date
+
+    Combines the ELUT with the science energy channels for the same date.
+
+    Parameters
+    ----------
+    date `astropy.time.Time`
+        Date to look up the ELUT.
+
+    Returns
+    -------
+
+    """
+    root = Path(__file__).parent.parent
+    elut_index_file = Path(root, *["config", "data", "elut", "elut_index.csv"])
+
+    elut_index = read_elut_index(elut_index_file)
+    elut_info = elut_index.at(date)
+    if len(elut_info) == 0:
+        raise ValueError(f"No ELUT for for date {date}")
+    elif len(elut_info) > 1:
+        raise ValueError(f"Multiple ELUTs for for date {date}")
+    start_date, end_date, elut_file = list(elut_info)[0]
+    sci_channels = get_sci_channels(date)
+    elut_table = read_elut(elut_file, sci_channels)
+
+
+    return elut_table
+
+def get_elut_correction(e_ind, pixel_data):
+    r"""
+    Get ELUT correction factors
+
+    Only correct the low and high energy edges as internal bins are contiguous.
+
+    Parameters
+    ----------
+    e_ind : `np.ndarray`
+        Energy channel indices
+    pixel_data : `~stixpy.products.sources.CompressedPixelData`
+        Pixel data
+
+    Returns
+    -------
+
+    """
+
+    energy_mask = pixel_data.energy_masks.energy_mask.astype(bool)
+    elut = get_elut(pixel_data.time_range.center)
+    ebin_edges_low = np.zeros((32, 12, 32), dtype=float)
+    ebin_edges_low[..., 1:] = elut.e_actual
+    ebin_edges_low = ebin_edges_low[..., energy_mask]
+    ebin_edges_high = np.zeros((32, 12, 32), dtype=float)
+    ebin_edges_high[..., 0:-1] = elut.e_actual
+    ebin_edges_high[..., -1] = np.nan
+    ebin_edges_high = ebin_edges_high[..., energy_mask]
+    ebin_widths = ebin_edges_high - ebin_edges_low
+    ebin_sci_edges_low = elut.e[..., 0:-1].value
+    ebin_sci_edges_low = ebin_sci_edges_low[..., energy_mask]
+    ebin_sci_edges_high = elut.e[..., 1:].value
+    ebin_sci_edges_high = ebin_sci_edges_high[..., energy_mask]
+    e_cor_low = (ebin_edges_high[..., e_ind[0]] - ebin_sci_edges_low[..., e_ind[0]]) / ebin_widths[..., e_ind[0]]
+    e_cor_high = (ebin_sci_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[-1]]) / ebin_widths[..., e_ind[-1]]
+
+    numbers = np.array([
+        0, 1, 2, 3, 4, 5, 6, 7, 13, 14, 15, 19, 20, 21, 22, 23, 24, 25,
+        26, 27, 28, 29, 30, 31
+    ])
+
+    bins =  ebin_sci_edges_high - ebin_sci_edges_low
+
+    bins_actual = ebin_widths[numbers, :8, :].mean(axis=1).mean(axis=0)
+
+    cor = bins / bins_actual
+
+    return e_cor_high, e_cor_low, cor

stixpy/calibration/energy.py

@@ -38,6 +38,10 @@ def get_elut(date):
     sci_channels = get_sci_channels(date)
     elut_table = read_elut(elut_file, sci_channels)
 
+
+    print('science_channels = ',np.shape(sci_channels))
+    print('elut_table_channels = ',np.shape(elut_table))
+
     return elut_table
 
 

stixpy/calibration/grid.py

@@ -7,13 +7,14 @@
 import astropy.units as u
 import numpy as np
 from astropy.table import Table
+import xraydb
 
 __all__ = ["get_grid_transmission", "_calculate_grid_transmission"]
 
 from stixpy.coordinates.frames import STIXImaging
 
 
-def get_grid_transmission(flare_location: STIXImaging):
+def get_grid_transmission(ph_energy, flare_location: STIXImaging):
     r"""
     Return the grid transmission for the 32 sub-collimators corrected for internal shadowing.
 
@@ -32,18 +33,141 @@ def get_grid_transmission(flare_location: STIXImaging):
     front = Table.read(grid_info / "grid_param_front.txt", format="ascii", names=column_names)
     rear = Table.read(grid_info / "grid_param_rear.txt", format="ascii", names=column_names)
 
-    transmission_front = _calculate_grid_transmission(front, flare_location)
-    transmission_rear = _calculate_grid_transmission(rear, flare_location)
-    total_transmission = transmission_front * transmission_rear
+    # ;; Orientation of the slits of the grid as seen from the detector side
+    grid_orient_front_all = front['o']
+    grid_orient_rear_all = rear['o']
+
+    pitch_front_all = front['p']
+    pitch_rear_all = rear['p']
+
+    thickness_front_all = front['thick']
+    thickness_rear_all = rear['thick']
+
+    sc = front['sc']
+
+
+    # fpath = loc_file( 'CFL_subcoll_transmission.txt', path = getenv('STX_GRID') )
+    column_names_cfl = ["subc_n", "subc_label", "intercept", "slope[1/deg]"]
+
+    subcol_transmission = Table.read(grid_info / "CFL_subcoll_transmission.txt", format="ascii", names=column_names_cfl)
+
+    subc_n_all = subcol_transmission['subc_n']
+    subc_label = subcol_transmission['subc_label']
+    intercept_all = subcol_transmission['intercept']
+    slope_all = subcol_transmission['slope[1/deg]']
+
+
+    muvals = xraydb.material_mu('W', ph_energy * 1e3, density=19.30, kind='total') / 10 # in units of mm^-1
+    L = 1 / muvals
+    print('L  = ', L)
+    # trans = np.exp(-0.4 / L)
+    subc_transm=L
+
+    det_indices_top24 =  np.array([0, 1, 2, 3, 4, 5, 6, 7, 13, 14, 15, 19, 
+                                     20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31])
+    # det_all = np.arange(0,32,1)
+
+    # det_indices_top24 =  [0, 1, 2, 3, 4, 5, 6, 7, 13, 14, 15, 19, 
+    #                                  20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]
+    idx_full = det_indices_top24
+    idx = [i for i,x in enumerate(sc-1) if x in  det_indices_top24]
+    print('idx = ', idx)
+    print('lenidx = ',len(idx))
+    # idx = det_all
+
+    # for i,idx in enumerate(idx_full):
+    print(grid_orient_front_all)
+
+    grid_orient_front = grid_orient_front_all[idx]
+    pitch_front = pitch_front_all[idx]
+    thickness_front = thickness_front_all[idx]
+
+    # ;;------ Rear grid
+    grid_orient_rear = grid_orient_rear_all[idx]
+    pitch_rear = pitch_rear_all[idx]
+    thickness_rear = thickness_rear_all[idx]
+
+    grid_orient_avg = (grid_orient_front + grid_orient_rear) / 2
+
+    flare_loc_deg = flare_location / 3600  #. ;; Convert coordinates to deg
+    theta = flare_loc_deg[0] * np.cos(np.deg2rad(grid_orient_avg)) + flare_loc_deg[1] * np.sin(np.deg2rad(grid_orient_avg)) 
+    print('THETA = ', theta)
+
+# ;;------ Subcollimator tranmsission at low energies
+# idx = np.where(subc_n_all eq (subc_n+1))
+
+    intercept = intercept_all[det_indices_top24]
+    slope = slope_all[det_indices_top24]
+    subc_transm_low_e = intercept + slope * theta
+
+    # ;;------ Transmission of front and rear grid
+    slit_to_pitch = np.sqrt(subc_transm_low_e)
+
+    slit_front = slit_to_pitch*pitch_front
+    slit_rear = slit_to_pitch*pitch_rear
+
+    transm_front = stx_grid_transmission(pitch_front, slit_front, thickness_front, L)
+    transm_rear = stx_grid_transmission(pitch_rear, slit_rear, thickness_rear, L)
+
+        # subc_transm.append(transm_front * transm_rear)
+
+    subc_transm = transm_front * transm_rear
+
+    # transmission_front = _calculate_grid_transmission(front, flare_location)
+    # transmission_rear = _calculate_grid_transmission(rear, flare_location)
+    # total_transmission = transmission_front * transmission_rear
 
     # The finest grids are made from multiple layers for the moment remove these and set 1
-    final_transmission = np.ones(32)
-    sc = front["sc"]
-    finest_scs = [11, 12, 13, 17, 18, 19]  # 1a, 2a, 1b, 2c, 1c, 2b
-    idx = np.argwhere(np.isin(sc, finest_scs, invert=True)).ravel()
-    final_transmission[sc[idx] - 1] = total_transmission[idx]
-
-    return final_transmission
+    # final_transmission = np.ones(32)
+    # sc = front["sc"]
+    # finest_scs = [11, 12, 13, 17, 18, 19]  # 1a, 2a, 1b, 2c, 1c, 2b
+    # idx = np.argwhere(np.isin(sc, finest_scs, invert=True)).ravel()
+    # final_transmission[sc[idx] - 1] = total_transmission[idx]
+
+    print('subc = ',subc_transm)
+    return subc_transm
+
+
+def stx_grid_transmission(pitch, slit, thickness, L):
+
+    ds = 5e-3
+    dh = 5e-2
+
+    # n_energies = np.shape(L)[0]
+    # n_subc = np.shape(pitch)
+
+    slit_rep = slit.reshape(1, len(slit))
+    pitch_rep = pitch.reshape(1, len(pitch))
+    H_rep = thickness.reshape(1, len(thickness))
+    L_rep = L.reshape(len(L), 1)
+
+    print('slit = ',slit[0])
+    print('pitch = ',pitch[0])
+    print('h_rep = ',thickness[0])
+
+    # slit_rep = np.tile(slit, (n_energies, 1))
+    # pitch_rep = np.tile(pitch, (n_energies, 1))
+    # H_rep = np.tile(thickness, (n_energies, 1))
+    # L_rep = np.tile(L, (n_subc, 1)).T
+
+    print(np.shape(slit_rep))
+    print(np.shape(pitch_rep))
+    print(np.shape(H_rep))
+    print(np.shape(L_rep))
+
+    # ;; Transmission for a wedge shape model for grid imperfections
+    g0 = slit_rep / pitch_rep + (pitch_rep - slit_rep) / pitch_rep * np.exp( - H_rep / L_rep )
+    ttt = L_rep / dh * ( 1. - np.exp(- dh / L_rep ) )
+    g1 = 2. * ds / pitch_rep * (ttt - np.exp( - H_rep / L_rep ))
+
+    print('g0 = ',g0[:,0])
+    print('g1 = ',g1[:,0])
+
+    g_transmission = g0 + g1
+
+    print('transmission = ',g_transmission)
+
+    return g_transmission
 
 
 def _calculate_grid_transmission(grid_params, flare_location):

stixpy/calibration/livetime.py

@@ -62,7 +62,7 @@
 import astropy.units as u
 import numpy as np
 
-__all__ = ["pileup_correction_factor", "get_livetime_fraction"]
+__all__ = ["pileup_correction_factor", "get_livetime_fraction","livetime_counts_corr"]
 
 from stixpy.io.readers import read_subc_params
 
@@ -93,7 +93,7 @@ def pileup_correction_factor():
     return prob_diff_pix
 
 
-def get_livetime_fraction(trigger_rate, *, eta=1.10 * u.us, tau=10.1 * u.us):
+def get_livetime_fraction(trigger_rate,time_del, eta=1.1e-6 *u.s, tau=10.1e-6 * u.s):
     """
     Return the live time fraction for the given trigger rate.
 
@@ -113,7 +113,16 @@ def get_livetime_fraction(trigger_rate, *, eta=1.10 * u.us, tau=10.1 * u.us):
     """
     beta = 0.94059104  # pileup_correction_factor()
 
+    tau = tau / time_del
+    eta = eta / time_del
+
     photons_in = trigger_rate / (1.0 - trigger_rate * (tau + eta))
     livetime_fraction = 1 / (1.0 + (tau + eta) * photons_in)
     two_photon = np.exp(-eta * beta * photons_in) * livetime_fraction
+
+
     return livetime_fraction, two_photon, photons_in
+
+
+
+

stixpy/calibration/visibility.py

@@ -177,6 +177,8 @@ def create_meta_pixels(
     t_ind = np.argwhere(t_mask).ravel()
     e_ind = np.argwhere(e_mask).ravel()
 
+    print('e_ind = ', e_ind)
+
     time_range = TimeRange(
         pixel_data.times[t_ind[0]] - pixel_data.duration[t_ind[0]] / 2,
         pixel_data.times[t_ind[-1]] + pixel_data.duration[t_ind[-1]] / 2,
@@ -201,21 +203,36 @@ def create_meta_pixels(
 
     pixel_data.data["livefrac"] = livefrac
 
+    print('e_ind = ',e_ind)
+
     e_cor_high, e_cor_low = get_elut_correction(e_ind, pixel_data)
 
+    print('e_cor_high_shape = ',np.shape(e_cor_high))
+    print('counts shape = ',np.shape(pixel_data.data["counts"]))
+
     # Get counts and other data.
     idx_pix = _PIXEL_SLICES.get(pixels.lower(), None)
+
+    print('pix = ',np.shape(e_cor_low[..., idx_pix]))
+
     if idx_pix is None:
         raise ValueError(f"Unrecognised input for 'pixels': {pixels}. Supported values: {list(_PIXEL_SLICES.keys())}")
     counts = pixel_data.data["counts"].astype(float)
     count_errors = np.sqrt(pixel_data.data["counts_comp_err"].astype(float).value ** 2 + counts.value) * u.ct
     ct = counts[t_ind][..., idx_pix, e_ind]
+    # print('ct = ',np.shape(ct))
+    ct_or = ct
     ct[..., 0] = ct[..., 0] * e_cor_low[..., idx_pix]
     ct[..., -1] = ct[..., -1] * e_cor_high[..., idx_pix]
     ct_error = count_errors[t_ind][..., idx_pix, e_ind]
     ct_error[..., 0] = ct_error[..., 0] * e_cor_low[..., idx_pix]
     ct_error[..., -1] = ct_error[..., -1] * e_cor_high[..., idx_pix]
 
+    # print(np.where( ((ct / ct_or) != 1) & (np.isnan(ct / ct_or) == False) )[0])
+
+    # indices = np.where(ct[...,0] != ct_or[...,0])
+    # print(indices)
+
     lt = (livefrac * pixel_data.data["timedel"].reshape(-1, 1).to("s"))[t_ind].sum(axis=0)
 
     ct_summed = ct.sum(axis=(0, 3))  # .astype(float)
@@ -290,6 +307,8 @@ def get_elut_correction(e_ind, pixel_data):
     ebin_sci_edges_high = ebin_sci_edges_high[..., energy_mask]
     e_cor_low = (ebin_edges_high[..., e_ind[0]] - ebin_sci_edges_low[..., e_ind[0]]) / ebin_widths[..., e_ind[0]]
     e_cor_high = (ebin_sci_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[-1]]) / ebin_widths[..., e_ind[-1]]
+
+
     return e_cor_high, e_cor_low