Commit 6f58b2d4 authored by Philipp Oleynik's avatar Philipp Oleynik
Browse files

Workfromhome addition.

fold_spectrum is added.
parent 0b4444a7
......@@ -3,11 +3,12 @@ import numpy as np
import scipy.special as scs
import socket
import scipy.constants as const
import math
if __name__ == "__main__":
print("This is a module! Bye.")
else:
if socket.gethostname() == 'spicule':
if socket.gethostname() == 'spicule' or socket.gethostname() == 'philyra':
loglut = np.load('/home/phil/git/philipp-packages/loglut.npz')
lut = loglut['lut']
npzfile = np.load('/home/phil/git/philipp-packages/sixs_lut_n.npz')
......@@ -23,8 +24,8 @@ core_threshold = 0.520 # in MeV TODO: Another magic number. Chosen as an initia
radiation_area = 5 * 5 * 4 * const.pi # radius 5 cm sphere area in cm2
def get_basepath():
if socket.gethostname() == 'spicule':
def get_basepath( ):
if socket.gethostname() == 'spicule' or socket.gethostname() == 'philyra':
return "/home/phil/Work/Calibration/"
else:
return "/home/pholey/Geant4/"
......@@ -39,7 +40,7 @@ def saturate2048(x):
return 0 if x < 0 else 2047 if x > 2047 else x
def logbit11to8(*, adc11=0, clut=lut):
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.
......@@ -77,7 +78,7 @@ def core_energy_to_64(radiated_energy):
x = max(radiated_energy * core_cal_slope_high + core_cal_offset_high,
radiated_energy * core_cal_slope_low + core_cal_offset_low)
adc_11 = saturate2048(int(x)) - core_instr_offset
adc_8 = logbit11to8(adc11=int(adc_11))
adc_8 = logbit11to8(adc11 = int(adc_11))
mute = 0 if radiated_energy < core_threshold else 1 # A correction for the threshold
y = 1 if radiated_energy > core_threshold else 0 # A correction for the algorithm in the SIXS firmware
return mute * (int(adc_8) // 4) if adc_8 >= 4 else y
......@@ -95,15 +96,15 @@ def side_energy_to_64(energy, side): # MeV to 6 bit side
side_cal_intercept_high = np.array([151.416, 138.757, 153.346, 147.083, 154.024])
side_cal_slope_high = np.array([284.069, 286.196, 284.543, 286.994, 288.152])
side_instr_offset = np.array([168, 159, 174, 167, 175])
x = max(side_cal_intercept_high[side] + side_cal_slope_high[side] * energy,
side_cal_intercept_low[side] + side_cal_slope_low[side] * energy)
adc_11 = saturate2048(int(x)) - side_instr_offset[side]
adc_8 = logbit11to8(adc11=int(adc_11))
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):
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.
......@@ -115,18 +116,44 @@ def get_energy_grid(*, channels_per_decade=256, min_energy=0.01, max_energy=1.0E
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
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
:param spectrum: intensities defined at the midpoint of each energy bin
:param response: geometric factor curve defined at the midpoint of each energy bin
:return: countrate in the channel described by the response.
"""
if bin_width is None:
return math.nan
if spectrum is None or response is None:
return 0
if (len(spectrum) == len(response)) and (len(spectrum) == len(bin_width)):
spectrum = np.multiply(spectrum, bin_width)
return np.dot(spectrum, response)
else:
return 0
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)
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