this notebook loads data acquired with a segmented silicon detector
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp
import copy
import re
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
shot_no = 37144
Download logfile and time of plasma start and plasma end
log_URL = f"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/StripDetector/PH32_{shot_no}.log"
plasma_URL = "http://golem.fjfi.cvut.cz/shots/{}/Diagnostics/BasicDiagnostics/Results/{}"
ds = np.DataSource(None)
logfile = ds.open(log_URL.format(shot_no))
#logfile = ds.open(f"PH32_{shot_no}.log")
try:
t_plasma_start = ds.open(plasma_URL.format(shot_no, "t_plasma_start"))
t_plasma_start = np.loadtxt(t_plasma_start)
t_plasma_end = ds.open(plasma_URL.format(shot_no, "t_plasma_end"))
t_plasma_end = np.loadtxt(t_plasma_end)
except IOSError:
pass
else:
t_plasma_start = -1.0
t_plasma_end = -1.0
Communication between rPi (or any other computer) and PH32 readout chip is achieved with SURE (Simple USB Readout Equipment).
SURE uses NMEA 0183 protocol for communication. Each message for configuring PH32 is acknowledged. These 'echos' are saved in the log file.
For example $ACKNL,SEATM -- acquisition time set to 100us*1b
is response to $SEATM,100
command which sets acquisition time.
$STARB,[wait-for-trigger],[delay]
sets detector for measurement and is fololowed by 200 frames of measured data eg.:
$FRAME,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0*5d
From log file we extract
acquisition mode
frames
frames = []
# defaults
acc_time = 0.15 #ms
delay = 0.10 #ms
start_delay = 0.50 #ms
acc_mode = ['hit count'] * 32
# read logfile
for line in logfile:
regex_seatm = re.match("\$ACKNL,SEATM -- acquisition time set to (\d+)us", line)
if regex_seatm:
acc_time = regex_seatm.group(1)
acc_time = int(acc_time)/1000.0 ## ms
regex_delay = re.match("\$ACKNL,DELAY -- delay set to (\d+)us", line)
if regex_delay:
delay = regex_delay.group(1)
delay = int(delay)/1000.0 ## ms
regex_selo1= re.match(
"\$ACKNL,SELO1 -- local configuration for channel (\d+) was set to value (\d+)", line)
if regex_selo1:
channel = int(regex_selo1.group(1))
if ( int(regex_selo1.group(2)) & 0x40):
acc_mode[channel] = "ToT"
elif ( int(regex_selo1.group(2)) & 0x20):
acc_mode[channel] = "first hit"
elif ( int(regex_selo1.group(2)) & 0x60):
acc_mode[channel] = "ToA"
else:
acc_mode[channel] = "hit count"
regex_startb = re.match("\$STARB,(\d+),(\d+)", line)
if regex_startb:
start_delay = regex_startb.group(2)
start_delay = int(start_delay)/1000.0 ## ms
regex_frame = re.match("\$FRAME,((?:\d+,)+)(\d+)\*[0-9a-f]{2}", line)
if regex_frame:
line = np.fromstring(''.join(regex_frame.groups()), sep=',')
### saturation or error in measeurement or chanel is disabled
line[line == 65535] = -1
frames.append(line)
frames = np.array(frames)
plt.pcolormesh(frames.T);
plt.axis('off')
plt.savefig('rawdata.png')
The order of channels on the PH32 chip does not correspond to the order of strips on the sensor; we need to rearrange them.
Strip n. 0 is on channel 15, strip n. 1 on channel 16, ...
channelOrder = [15,16,14,17,13,18,12,19,11,20,10,21,9,22,8,23,24,7,25,6,26,5,27,4,28,3,29,2,30,1,31,0]
permutation_matrix = np.zeros((len(channelOrder), len(channelOrder)))
for idx, i in enumerate(channelOrder):
permutation_matrix[i, idx] = 1
frames = np.dot(frames, permutation_matrix)
## rearange acc_mode
new_acc_mode = [''] *32
for new, old in enumerate(channelOrder):
new_acc_mode[new] = acc_mode[old]
acc_mode = new_acc_mode
For energy measeuremets convert number of clock cycles to deposited energy
energy_calibration = 1.1504 #keV/clock
mask_energy = [ch == "ToT" or ch == "first hit" for ch in acc_mode]
frames[:,mask_energy] *= energy_calibration
About half of the frames are usually empty because the detector records after plasma end, so we'll plot only the interesting part.
def all_eq(iterable, val):
return all([i == val for i in iterable])
time_btw_frames = acc_time + delay
y_tics = np.arange(0,32,1)
t_tics = np.arange(0, frames.shape[0]*time_btw_frames, time_btw_frames)
t_tics += start_delay
cmap = copy.copy(mp.cm.get_cmap("viridis"))
#bad measeurement in grey
cmap.set_under("grey")
fig,ax = plt.subplots()
im = ax.pcolormesh(t_tics, y_tics, frames.T, vmin = 0, vmax = np.amax(frames[5:80]),
cmap=cmap, shading = 'auto')
ax.set_ylabel("strip n.")
ax.set_xlabel("time [ms]")
if t_plasma_start != -1.0 and t_plasma_start != -1.0:
ax.set_xlim(start_delay, t_plasma_end + 2)
ax.axvline(t_plasma_start, color = "white", dashes=[1, 2], label='plasma start')
ax.axvline(t_plasma_end, color = "white", dashes=[1, 2], label='plasma end')
else:
ax.set_xlim(start_delay, 30.0)
cbar = fig.colorbar(im, ax=ax);
if all_eq(acc_mode, "hit count"):
cbar.set_label(acc_mode)
elif all_eq(acc_mode, "ToT"):
cbar.set_label("Deposited energy [keV/{}ms]".format(acc_time))
elif all_eq(acc_mode,"first hit"):
cbar.set_label("First hit energy [kev]")
else:
cbar.set_label("mix")
fig.savefig('icon-fig.png')
mask_hit = [ch == 'hit count' for ch in acc_mode]
mask_ToT = [ch == 'ToT' for ch in acc_mode]
mask_fToT = [ch == 'first hit' for ch in acc_mode]
hit_curves = [hv.Curve((t_tics,ch), 't [ms]', 'hit count') for ch in frames[:,mask_hit].T]
tot_curves = [hv.Curve((t_tics,ch), 't [ms]', 'E dep} [keV]') for ch in frames[:,mask_ToT].T]
first_hit_curves = [hv.Curve((t_tics,ch), 't [ms]', '$E first-dep [keV]') for ch in frames[:,mask_fToT].T]
hit_overlay = hv.Overlay(hit_curves) if hit_curves else hv.Empty()
tot_overlay = hv.Overlay(tot_curves) if tot_curves else hv.Empty()
first_hit_overlay = hv.Overlay(first_hit_curves) if first_hit_curves else hv.Empty()
layout = hv.Layout(hit_overlay + tot_overlay + first_hit_overlay).cols(1)
x_max = t_plasma_end + 2 if t_plasma_end != -1.0 else 30.0
layout.opts(opts.Overlay(height=200, width=900, xlim=(start_delay, x_max)))