LUT added

sixsutils.py

+#!/usr/bin/python3
+import numpy as np
+import scipy.special as scs
+import socket
+import scipy.constants as const
+
+if __name__ == "__main__":
+    print("This is a module! Bye.")
+else:
+    if socket.gethostname() == 'spicule':
+        loglut = np.load('/home/phil/git/philipp-packages/loglut.npz')
+        lut = loglut['lut']
+        npzfile = np.load('/home/phil/git/philipp-packages/sixs_lut.npz')
+        sixs_lut = npzfile['arr_0']
+    else:
+        loglut = np.load('/home/pholey/git/philipp-packages/loglut.npz')
+        lut = loglut['lut']
+        npzfile = np.load('/home/pholey/git/philipp-packages/sixs_lut.npz')
+        sixs_lut = npzfile['arr_0']
+
+    side_threshold = 0.035  # TODO: It is a magic number.
+    core_threshold = 0.50  # in MeV TODO: Another magic number. Chosen as an initial value.
+    radiation_area = 5 * 5 * 4 * const.pi  # radius 5 cm sphere area in cm2
+
+
+def get_basepath():
+    if socket.gethostname() == 'spicule':
+        return "/home/phil/Work/Calibration/"
+    else:
+        return "/home/pholey/Geant4/"
+
+
+def saturate2048(x):
+    """
+    Saturates value to 11-bit and clips negative values.
+    :type x: int, float
+    :rtype: int, float
+    """
+    return 0 if x < 0 else 2047 if x > 2047 else x
+
+
+def logbit11to8(*, adc11=0, clut=lut):
+    """
+    Transforms a 11-bit ADC value ot a 8-bit number using a logarithmic LUT
+    :param adc11: float, int. Floats are floored.
+    :param clut: a 2048 element array containing the LUT
+    :return: an 8-bit integer positive value
+    """
+    adc11 = int(saturate2048(adc11))
+    return clut[adc11]
+
+
+def birks_correction(energy):  # Birks coefficient and parameter estimates taken from the calibration of RADMON
+    """
+    Converts deposited energy in CsI to radiated photon energy using Birks correction by Oleynik et al.,ASR, 2019
+    :param energy: Input energy in MeV
+    :return: Corrected energy in Mev
+    """
+    kbirks = 0.68E-3
+    return energy * (1 - scs.hyp2f1(1, 1 / 0.678, 1 + 1 / 0.678, -np.power(energy, 0.678) / (95 * kbirks)))
+
+
+def core_energy_to_64(energy):
+    """
+    Converts energy deposited in the core detector to a 6-bit value by the algorithms of SIXS.
+    :param energy: Input energy in MeV, float.
+    :return: An output channel, 0..63 integer value.
+    """
+    # core_cal_intercept = 178.8581  # old ones, should not be used
+    # core_cal_slope = 15.5249
+    # core_cal_power = 2.5937
+    core_instr_offset = 166
+    core_cal_slope = 47.121
+    core_cal_offset = 157  # TODO: It is a magic number. Fitted only, not verified.
+    radiated_energy = birks_correction(energy)
+    adc_11 = saturate2048(int(radiated_energy * core_cal_slope + core_cal_offset)) - core_instr_offset
+    adc_8 = logbit11to8(adc11=int(adc_11))
+    x = 1 if radiated_energy > core_threshold else 0  # A correction for the algorithm in the SIXS firmware
+    return int(round(adc_8)) // 4 if adc_8 >= 4 else x
+
+
+def side_energy_to_64(energy, side):  # MeV to 6 bit side
+    """
+    Converts energy deposited in a side detector to a 6-bit value by the algorithms of SIXS.
+    :type side: Side number, integer 0..4. Side 0 is the top one.
+    :param energy: Input energy in MeV, float.
+    :return: An output channel, 0..63 integer value.
+    """
+    side_cal_intercept = np.array([167.513, 158.190, 172.329, 165.407, 173.593])
+    side_cal_slope = np.array([243.464, 250.992, 253.015, 252.194, 253.492])
+    side_instr_offset = np.array([168, 159, 174, 167, 175])
+    adc_11 = saturate2048(int(side_cal_intercept[side] + side_cal_slope[side] * energy)) - side_instr_offset[side]
+    adc_8 = logbit11to8(adc11=int(adc_11))
+    return int(round(adc_8)) // 4 if adc_8 >= 0 else 0
+
+
+def get_energy_grid(*, channels_per_decade=256, min_energy=0.01, max_energy=1.0E5):
+    """
+    Calculates a standard logarithmic energy grid.
+    :param channels_per_decade: Number of channels between 1 and 10, excluding exact value of 10.
+    :param min_energy: Starting energy. Can be anything, not necessarily powers of 10
+    :param max_energy: Upper limit of energy. Can be anything, not necessarily powers of 10
+    :return: An integer and three float arrays: number_of_energy_steps, energy_midpoint, energy_cut, energy_bin_width
+    """
+    emin_start = (np.floor(np.log10(min_energy) * channels_per_decade) / channels_per_decade)
+    emax_stop = (np.floor(np.log10(max_energy) * channels_per_decade) / channels_per_decade)
+    number_of_energy_steps = int((emax_stop - emin_start) * channels_per_decade + 1)
+    log_step = 1.0 / channels_per_decade
+
+    energy_midpoint = np.zeros(shape=(number_of_energy_steps,), dtype=float)
+    energy_cut = np.zeros(shape=(number_of_energy_steps,), dtype=float)
+    energy_bin_width = np.zeros(shape=(number_of_energy_steps,), dtype=float)
+
+    for i in range(0, number_of_energy_steps, 1):
+        midpoint = np.power(10, emin_start) * np.power(10, log_step * (i + 0.5))
+        energy_bin_low = np.power(10, emin_start) * np.power(10, log_step * i)
+        energy_bin_high = np.power(10, emin_start) * np.power(10, log_step * (i + 1))
+
+        energy_cut[i] = energy_bin_high
+        energy_midpoint[i] = midpoint
+        energy_bin_width[i] = energy_bin_high - energy_bin_low
+
+    return number_of_energy_steps, energy_midpoint, energy_cut, energy_bin_width