sixsutils.py 5.19 KB
Newer Older
Philipp Oleynik's avatar
Philipp Oleynik committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#!/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']

Philipp Oleynik's avatar
bugfix    
Philipp Oleynik committed
21
22
23
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
Philipp Oleynik's avatar
Philipp Oleynik committed
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59


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
60
61
62
63
    if energy > 100:
        return energy
    else:
        return energy * (1 - scs.hyp2f1(1, 1 / 0.678, 1 + 1 / 0.678, -np.power(energy, 0.678) / (95 * kbirks)))
Philipp Oleynik's avatar
Philipp Oleynik committed
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126


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