vectorized e to 2048

plotutil.py

@@ -54,6 +54,38 @@ def set_log_axes_simple(axes):
     axes.grid(True, which = 'both', alpha = 0.3, zorder = 0)
 
 
+def set_log_axes_2048(axes):
+    axes.tick_params(direction = 'in', which = 'both', zorder = 4)
+    axes.set_xscale("log", nonposx = 'clip', basex = 2.0)
+    axes.set_yscale("log", nonposy = 'clip', basey = 2.0)
+    axes.grid(True, which = 'both', alpha = 0.3, zorder = 0)
+    axes.set_xticks([1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048], minor = False)
+    axes.set_yticks([1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048], minor = False)
+    axes.set_xticklabels([r'1', r'2', r'4', r'8', r'16', r'32', r'64', r'128', r'256', r'512', r'1024', r'2048'], minor = False)
+    axes.set_yticklabels([r'1', r'2', r'4', r'8', r'16', r'32', r'64', r'128', r'256', r'512', r'1024', r'2048'], minor = False)
+    axes.set_ylim(1, 2048)
+    axes.set_xlim(1, 2048)
+
+
+def set_log_axes_bin16(axes):
+    axes.tick_params(direction = 'in', which = 'both', zorder = 4)
+    axes.set_xscale("log", nonposx = 'clip', basex = 2.0)
+    axes.set_yscale("log", nonposy = 'clip', basey = 2.0)
+    axes.grid(True, which = 'both', alpha = 0.3, zorder = 0)
+    axes.set_xticks([1, 2, 4, 8, 16, 32, 64, 128, 256,
+                     512, 1024, 2048, 4096, 8192, 16384, 32768, 65536], minor = False)
+    axes.set_yticks([1, 2, 4, 8, 16, 32, 64, 128, 256,
+                     512, 1024, 2048, 4096, 8192, 16384, 32768, 65536], minor = False)
+    axes.set_xticklabels([r'1', r'2', r'4', r'8', r'16', r'32', r'64', r'128', r'256',
+                          r'512', r'1024', r'2048', r'4096', r'8192', r'16384', r'32768', r'65536'],
+                         minor = False)
+    axes.set_yticklabels([r'1', r'2', r'4', r'8', r'16', r'32', r'64', r'128', r'256',
+                          r'512', r'1024', r'2048', r'4096', r'8192', r'16384', r'32768', r'65536'],
+                         minor = False)
+    axes.set_ylim(1, 2048)
+    axes.set_xlim(1, 2048)
+
+
 def set_log_axes(axes):
     set_log_axes_noaspect(axes)
     axes.set_aspect(0.56, adjustable = 'box')

sixsutils.py

@@ -4,6 +4,8 @@ import scipy.special as scs
 import socket
 import scipy.constants as const
 import math
+import fitsio
+from fitsio import FITS, FITSHDR
 
 if __name__ == "__main__":
     print("This is a module! Bye.")
@@ -31,13 +33,22 @@ def get_basepath( ):
         return "/home/pholey/Geant4/"
 
 
+def saturate2048_(x):
+    """
+    Saturates value to 11-bit and clips negative values. Old implementation.
+    :type x: int, float
+    :rtype: int, float
+    """
+    return 0 if x < 0 else 2047 if x > 2047 else x
+
+
 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
+    return np.clip(x, 0, 2047)
 
 
 def logbit11to8(*, adc11 = 0, clut = lut):
@@ -84,6 +95,24 @@ def core_energy_to_64(radiated_energy):
     return mute * (int(adc_8) // 4) if adc_8 >= 4 else y
 
 
+def core_energy_to_2048(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.
+    """
+    core_instr_offset = 166
+    core_cal_slope_high = 47.363
+    core_cal_offset_high = 155
+    core_cal_slope_low = 33.1541
+    core_cal_offset_low = 169.3  # TODO: It is a magic number. Fitted only, not verified.
+    h = radiated_energy * core_cal_slope_high + core_cal_offset_high
+    l = radiated_energy * core_cal_slope_low + core_cal_offset_low
+    x = np.floor(np.fmax(h, l)).astype(int)
+    adc_11 = np.clip(x, 0, 2047)  # - core_instr_offset
+    return adc_11
+
+
 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.