Skip to content
GitLab
Commit e3568bfc authored by Philipp Oleynik's avatar Philipp Oleynik
Browse files

Argh.

parent 50a84c74
Hide whitespace changes
Inline Side-by-side
sixsutils.py
View file @ e3568bfc
......@@ -8,7 +8,7 @@ from sixs_main_lut import sixs_main_lut_matrix
from sixs_log_lut import sixs_log_lut_array
def generate_powerlaw(*, energy_grid = None, power_index = -2, sp_norm = 1.0):
def generate_powerlaw(*, energy_grid=None, power_index=-2, sp_norm=1.0):
"""
:param energy_grid:
......@@ -19,11 +19,11 @@ def generate_powerlaw(*, energy_grid = None, power_index = -2, sp_norm = 1.0):
spectrum = []
for energy in energy_grid:
spectrum.append(sp_norm * np.power(energy, power_index))
return np.array(spectrum, dtype = float)
return np.array(spectrum, dtype=float)
def generate_integral_powerlaw(*, energy_grid = None, grid_bin_width = None,
power_index = -2, sp_norm = 1.0, fast = False):
def generate_integral_powerlaw(*, energy_grid=None, grid_bin_width=None,
power_index=-2, sp_norm=1.0, fast=False):
"""
:param energy_grid:
......@@ -48,11 +48,11 @@ def generate_integral_powerlaw(*, energy_grid = None, grid_bin_width = None,
for energy_j in range(energy_i, len(energy_grid)):
integral_sp += spectrum[energy_j] * grid_bin_width[energy_j]
spectrum[energy_i] = integral_sp
return np.array(spectrum, dtype = float)
return np.array(spectrum, dtype=float)
def fold_spectrum(*, bin_width = None, spectrum = None, response = None):
def fold_spectrum(*, bin_width=None, spectrum=None, response=None):
"""
Folds incident spectrum with an instrument response. Very simple function to start with. Int( spectrum * response * dE)
:param bin_width: energy grid intervals defined at midpoints of each energy bin
......@@ -71,7 +71,7 @@ def fold_spectrum(*, bin_width = None, spectrum = None, response = None):
return 0
def get_energy_grid(*, channels_per_decade = 256, min_energy = 0.01, max_energy = 1.0E5):
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.
......@@ -84,20 +84,20 @@ def get_energy_grid(*, channels_per_decade = 256, min_energy = 0.01, max_energy
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)
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
......@@ -111,34 +111,34 @@ def saturate2048(x):
class Sixs:
def __init__(self, *, old_lut = False, old_lut_file_name = ""):
def __init__(self, *, old_lut=False, old_lut_file_name=""):
self.lut = sixs_main_lut_matrix
self.loglut = sixs_log_lut_array
self.core_instr_offset = 166
self.core_cal_slope_high = 47.363
self.core_cal_offset_high = 155
self.core_cal_slope_low = 33.1541
self.core_cal_offset_low = 169.3 # TODO: It is a magic number. Fitted only, not verified.
self.side_threshold = 0.035 # TODO: It is a magic number.
self.core_threshold = 0.520 # in MeV TODO: Another magic number. Chosen as an initial value.
self.radiation_area = 5 * 5 * 4 * const.pi # radius 5 cm sphere area in cm2
self.kbirks = 0.68E-3
self.side_cal_intercept_low = np.array([167.513, 158.190, 172.329, 165.407, 173.593])
self.side_cal_slope_low = np.array([243.464, 250.992, 253.015, 252.194, 253.492])
self.side_cal_intercept_high = np.array([151.416, 138.757, 153.346, 147.083, 154.024])
self.side_cal_slope_high = np.array([284.069, 286.196, 284.543, 286.994, 288.152])
self.side_instr_offset = np.array([168, 159, 174, 167, 175])
if old_lut:
npzfile = np.load(old_lut_file_name)
self.lut = np.array(npzfile['arr_0'])
else:
pass
# WARNING! These are heavily dependant on the user name and the machine! Change it to something sane.
if socket.gethostname() == 'spicule' or socket.gethostname() == 'philyra':
self.basepath = "/home/phil/Work/Calibration/"
......@@ -146,14 +146,14 @@ class Sixs:
else:
self.basepath = "/home/pholey/Geant4/"
self.datapath = "/home/pholey/SIXS_Data/bc_mpo_sixs/data_calibrated/"
def get_basepath(self):
return self.basepath
def get_datapath(self):
return self.datapath
def logbit11to8(self, *, adc11 = 0):
def logbit11to8(self, *, adc11=0):
"""
Transforms a 11-bit ADC value ot a 8-bit number using a logarithmic LUT
:param adc11: float, int, ndarray. Floats are floored.
......@@ -164,7 +164,7 @@ class Sixs:
else:
adc11 = (list(map(int, np.clip(adc11, 0, 2047))))
return self.loglut[adc11]
def birks_correction(self, 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
......@@ -179,7 +179,7 @@ class Sixs:
else:
return np.array([(x * (1 - scs.hyp2f1(1, 1 / 0.678, 1 + 1 / 0.678, -np.power(x, 0.678) / (95 * self.kbirks)))
if x < 100 else x) for x in energy])
def core_energy_to_64(self, radiated_energy):
"""
Converts energy deposited in the core detector to a 6-bit value by the algorithms of SIXS.
......@@ -190,10 +190,10 @@ class Sixs:
radiated_energy * self.core_cal_slope_low + self.core_cal_offset_low)
adc_11 = np.clip(np.floor(x), 0, 2047)
adc_11 -= self.core_instr_offset
adc_6 = self.logbit11to8(adc11 = adc_11) // 4
adc_6 = self.logbit11to8(adc11=adc_11) // 4
return np.maximum(adc_6, 0 if radiated_energy < self.core_threshold else 1)
def core_energy_to_64_multi(self, radiated_energy):
"""
Converts energy deposited in the core detector to a 6-bit value by the algorithms of SIXS.
......@@ -203,11 +203,11 @@ class Sixs:
x = np.maximum(radiated_energy * self.core_cal_slope_high + self.core_cal_offset_high,
radiated_energy * self.core_cal_slope_low + self.core_cal_offset_low)
adc_11 = saturate2048(int(x)) - self.core_instr_offset
adc_6 = self.logbit11to8(adc11 = adc_11) // 4
adc_6 = self.logbit11to8(adc11=adc_11) // 4
core_floor = list(map(int, radiated_energy < self.core_threshold))
return np.maximum(adc_6, core_floor)
def core_energy_to_64_uni(self, radiated_energy):
"""
Converts energy deposited in the core detector to a 6-bit value by the algorithms of SIXS.
......@@ -218,21 +218,21 @@ class Sixs:
radiated_energy * self.core_cal_slope_low + self.core_cal_offset_low)
adc_11 = np.clip(np.floor(x), 0, 2047)
adc_11 -= self.core_instr_offset
adc_6 = self.logbit11to8(adc11 = adc_11) // 4
adc_6 = self.logbit11to8(adc11=adc_11) // 4
if np.isscalar(radiated_energy):
core_floor = int(radiated_energy < self.core_threshold)
else:
core_floor = list(map(int, radiated_energy < self.core_threshold))
core_floor = list(map(int, radiated_energy > self.core_threshold))
return np.maximum(adc_6, core_floor)
def core_energy_to_2048(self, radiated_energy):
"""
Experimental. Converts energy deposited in the core detector to a 11-bit value by the algorithms of SIXS.
:param radiated_energy: Input energy in MeV, float.
:return: An output channel, 0..2047 integer value.
"""
h = radiated_energy * self.core_cal_slope_high + self.core_cal_offset_high
ll = radiated_energy * self.core_cal_slope_low + self.core_cal_offset_low
ll *= 0.0
......@@ -240,7 +240,7 @@ class Sixs:
adc_11 = np.clip(x, 0, 2047) # - core_instr_offset
adc_11[adc_11 < self.core_instr_offset] = 0
return adc_11
def side_energy_to_64(self, energy, side): # MeV to 6 bit side
"""
Converts energy deposited in a side detector to a 6-bit value by the algorithms of SIXS.
......@@ -249,13 +249,13 @@ class Sixs:
:param energy: Input energy in MeV, float.
:return: An output channel, 0..63 integer value.
"""
x = max(self.side_cal_intercept_high[side] + self.side_cal_slope_high[side] * energy,
self.side_cal_intercept_low[side] + self.side_cal_slope_low[side] * energy)
adc_11 = saturate2048(int(x)) - self.side_instr_offset[side]
adc_8 = self.logbit11to8(adc11 = int(adc_11))
adc_8 = self.logbit11to8(adc11=int(adc_11))
return int(round(adc_8)) // 4 if adc_8 >= 0 else 0
def side_energy_to_64_uni(self, energy, side): # MeV to 6 bit side
"""
Converts energy deposited in side detectors to 6-bit values by the algorithms of SIXS.
......@@ -267,16 +267,16 @@ class Sixs:
x = max(self.side_cal_intercept_high[side] + self.side_cal_slope_high[side] * energy,
self.side_cal_intercept_low[side] + self.side_cal_slope_low[side] * energy)
adc_11 = saturate2048(int(x)) - self.side_instr_offset[side]
adc_8 = self.logbit11to8(adc11 = int(adc_11))
adc_8 = self.logbit11to8(adc11=int(adc_11))
return int(round(adc_8)) // 4 if adc_8 >= 0 else 0
else:
x = np.maximum(self.side_cal_intercept_high[side] + np.multiply(self.side_cal_slope_high[side], energy),
self.side_cal_intercept_low[side] + np.multiply(self.side_cal_slope_low[side], energy))
adc_11 = np.clip(np.floor(x), 0, 2047) - self.side_instr_offset[side]
adc_8 = self.logbit11to8(adc11 = adc_11)
adc_8 = self.logbit11to8(adc11=adc_11)
return adc_8 // 4
def classifier_64(self, side6bit, core6bit, sidenum, thlevel = 1):
def classifier_64(self, side6bit, core6bit, sidenum=0, thlevel=1):
"""
:param side6bit:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment