Skip to content
Snippets Groups Projects
Commit a8312e7b authored by Philipp Schmidt's avatar Philipp Schmidt
Browse files

Make REMI plots more robust against empty signal

parent 3077e2ad
Branches fix/REMI-robuster-plots
No related tags found
1 merge request!1005[REMI] Make plots robuster against missing signal
This commit is part of merge request !1005. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
notebooks/REMI/REMI_Digitize_and_Transform.ipynb
+ 40
− 25
Edit
View file @ a8312e7b
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Data selection parameters. # Data selection parameters.
run = 104 # Run ID. run = 104 # Run ID.
in_folder = '/gpfs/exfel/exp/SQS/202101/p002535/raw' # Partial input path appended with run ID. in_folder = '/gpfs/exfel/exp/SQS/202101/p002535/raw' # Partial input path appended with run ID.
out_folder = '/gpfs/exfel/exp/SQS/202101/p002535/scratch/cal_test' # Full path to output folder. out_folder = '/gpfs/exfel/exp/SQS/202101/p002535/scratch/cal_test' # Full path to output folder.
calib_config_path = '/gpfs/exfel/exp/SQS/202101/p002535/usr/config_board2+4.yaml' # Path to correction and transform configuration calib_config_path = '/gpfs/exfel/exp/SQS/202101/p002535/usr/config_board2+4.yaml' # Path to correction and transform configuration
# Output parameters. # Output parameters.
karabo_id = 'SQS_REMI_DLD6' # Karabo device ID root for virtual output device. karabo_id = 'SQS_REMI_DLD6' # Karabo device ID root for virtual output device.
proposal = '' # Proposal, leave empty for auto detection based on in_folder proposal = '' # Proposal, leave empty for auto detection based on in_folder
out_aggregator = 'REMI01' # Aggregator name for output files. out_aggregator = 'REMI01' # Aggregator name for output files.
out_seq_len = 5000 # Number of trains per sequence file in output. out_seq_len = 5000 # Number of trains per sequence file in output.
det_device_id = '{karabo_id}/DET/{det_name}' # Karabo device ID for virtual output device. det_device_id = '{karabo_id}/DET/{det_name}' # Karabo device ID for virtual output device.
det_output_key = 'output' # Pipeline name for fast data output. det_output_key = 'output' # Pipeline name for fast data output.
chunks_triggers = [500] # HDF chunk size for triggers. chunks_triggers = [500] # HDF chunk size for triggers.
chunks_edges = [500, 7, 50] # HDF chunk size for edges. chunks_edges = [500, 7, 50] # HDF chunk size for edges.
chunks_amplitudes = [500, 7, 50] # HDF chunk size for amplitudes. chunks_amplitudes = [500, 7, 50] # HDF chunk size for amplitudes.
chunks_hits = [50, 50] # HDF chunk size for hits. chunks_hits = [50, 50] # HDF chunk size for hits.
chunks_signals = [50, 50] # HDF chunk size for signals. chunks_signals = [50, 50] # HDF chunk size for signals.
dataset_compression = 'gzip' # HDF compression method. dataset_compression = 'gzip' # HDF compression method.
dataset_compression_opts = 3 # HDF GZIP compression level. dataset_compression_opts = 3 # HDF GZIP compression level.
# Detector parameters. # Detector parameters.
quad_anode = False # Reconstruction assumes a hex anode by default, change for quad anodes. quad_anode = False # Reconstruction assumes a hex anode by default, change for quad anodes.
# Trigger parameters. # Trigger parameters.
ppt_source = 'SQS_RR_UTC/TSYS/TIMESERVER:outputBunchPattern' ppt_source = 'SQS_RR_UTC/TSYS/TIMESERVER:outputBunchPattern'
ignore_fel = False # Ignore any FEL entries in the PPT. ignore_fel = False # Ignore any FEL entries in the PPT.
ignore_ppl = False # Ignore any PPL entries in the PPT. ignore_ppl = False # Ignore any PPL entries in the PPT.
trailing_trigger = False # Add a trigger after all regular pulses with the remaining trace. trailing_trigger = False # Add a trigger after all regular pulses with the remaining trace.
ppl_offset = 0 # In units of the PPT. ppl_offset = 0 # In units of the PPT.
laser_ppt_mask = -1 # Bit mask for used laser, negative to auto-detect from instrument. laser_ppt_mask = -1 # Bit mask for used laser, negative to auto-detect from instrument.
instrument_sase = 3 # Which SASE we're running at for PPT decoding. instrument_sase = 3 # Which SASE we're running at for PPT decoding.
first_pulse_offset = 10000 # Sample position where the first pulse begins, ignored when PPT is reconstructed. first_pulse_offset = 10000 # Sample position where the first pulse begins, ignored when PPT is reconstructed.
single_pulse_length = 25000 # How many samples if there's only one pulse. single_pulse_length = 25000 # How many samples if there's only one pulse.
pulse_start_offset = 0 # Signal offset at the start of each pulse. pulse_start_offset = 0 # Signal offset at the start of each pulse.
pulse_end_offset = 0 # Signal offset at the end of each pulse. pulse_end_offset = 0 # Signal offset at the end of each pulse.
# PPT reconstruction parameters. # PPT reconstruction parameters.
reconstruct_ppt = False # Reconstruct PPT from some trigger edges. reconstruct_ppt = False # Reconstruct PPT from some trigger edges.
trigger_edge_channel = '4_D' # Channel to use for triggering. trigger_edge_channel = '4_D' # Channel to use for triggering.
trigger_edge_offset = 0 # Offset to apply to the first trigger edge position to compute first pulse offset. trigger_edge_offset = 0 # Offset to apply to the first trigger edge position to compute first pulse offset.
fake_ppt_offset = 0 # Offset in reconstructed PPT for pulses. fake_ppt_offset = 0 # Offset in reconstructed PPT for pulses.
# Parallelization parameters. # Parallelization parameters.
mp_find_triggers = 0.5 # Parallelization for finding triggers. mp_find_triggers = 0.5 # Parallelization for finding triggers.
mp_find_edges = 0.5 # Parallelization for digitizing analog signal. mp_find_edges = 0.5 # Parallelization for digitizing analog signal.
mt_avg_trace = 2 # Parallelization for trace averaging. mt_avg_trace = 2 # Parallelization for trace averaging.
mp_rec_hits = 1.0 # Parallelization for hit reconstruction. mp_rec_hits = 1.0 # Parallelization for hit reconstruction.
# DEPRECATED AND IGNORED # DEPRECATED AND IGNORED
# Left for compatibility with webservice or legacy configuration. # Left for compatibility with webservice or legacy configuration.
cycle = '' # Proposal cycle, passed by webservice but not used. cycle = '' # Proposal cycle, passed by webservice but not used.
karabo_da = 'bar' # Karabo data aggregator name, passed by webservice but not used karabo_da = 'bar' # Karabo data aggregator name, passed by webservice but not used
cal_db_timeout = 0 # Calibration DB timeout, passed by webservice but not used. cal_db_timeout = 0 # Calibration DB timeout, passed by webservice but not used.
cal_db_interface = 'foo' # Calibration DB interface, passed by webservice but not used. cal_db_interface = 'foo' # Calibration DB interface, passed by webservice but not used.
save_raw_triggers = True # Whether to save trigger position in files, ignored and always enabled. save_raw_triggers = True # Whether to save trigger position in files, ignored and always enabled.
save_raw_edges = True # Whether to save digitized edge positions in files, ignored and always enabled. save_raw_edges = True # Whether to save digitized edge positions in files, ignored and always enabled.
save_raw_amplitudes = True # Whether to save analog pulse amplitudes in files, ignored and always enabled. save_raw_amplitudes = True # Whether to save analog pulse amplitudes in files, ignored and always enabled.
save_rec_signals = True # Whether to save reconstructed signals (u1-w2, mcp) in files, ignored and always enabled. save_rec_signals = True # Whether to save reconstructed signals (u1-w2, mcp) in files, ignored and always enabled.
save_rec_hits = True # Whether to save reoncstructed hits (x,y,t,m) in files, ignored and always enabled. save_rec_hits = True # Whether to save reoncstructed hits (x,y,t,m) in files, ignored and always enabled.
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from datetime import datetime from datetime import datetime
from logging import warning from logging import warning
from pathlib import Path from pathlib import Path
import re import re
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm from matplotlib.colors import LogNorm
from matplotlib.patches import Circle from matplotlib.patches import Circle
from threadpoolctl import threadpool_limits from threadpoolctl import threadpool_limits
import tabulate import tabulate
from IPython.display import Latex, Markdown, display from IPython.display import Latex, Markdown, display
import h5py import h5py
import pasha as psh import pasha as psh
from euxfel_bunch_pattern import indices_at_sase, indices_at_laser from euxfel_bunch_pattern import indices_at_sase, indices_at_laser
from extra_data import RunDirectory, by_id from extra_data import RunDirectory, by_id
from extra_remi import Analysis, trigger_dt from extra_remi import Analysis, trigger_dt
from extra_remi.util import timing from extra_remi.util import timing
from extra_remi.rd_resort import signal_dt, hit_dt from extra_remi.rd_resort import signal_dt, hit_dt
from extra_remi.files import DataFile, sequence_pulses from extra_remi.files import DataFile, sequence_pulses
if quad_anode: if quad_anode:
from extra_remi.plots import plot_detector_diagnostics_quad as plot_detector_diagnostics from extra_remi.plots import plot_detector_diagnostics_quad as plot_detector_diagnostics
else: else:
from extra_remi.plots import plot_detector_diagnostics_hex as plot_detector_diagnostics from extra_remi.plots import plot_detector_diagnostics_hex as plot_detector_diagnostics
%matplotlib inline %matplotlib inline
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def finite_flattened_slice(array, slice_=np.s_[:]): def finite_flattened_slice(array, slice_=np.s_[:]):
"""Return flattened and finite values for a given slice.""" """Return flattened and finite values for a given slice."""
sliced_array = array[slice_] sliced_array = array[slice_]
return sliced_array[np.isfinite(sliced_array)] return sliced_array[np.isfinite(sliced_array)]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
calib_config_path = Path(calib_config_path) calib_config_path = Path(calib_config_path)
if not calib_config_path.is_file(): if not calib_config_path.is_file():
# If the path cannot be resolved right now, try the same path relative to in_folder. # If the path cannot be resolved right now, try the same path relative to in_folder.
calib_config_path = Path(in_folder) / calib_config_path calib_config_path = Path(in_folder) / calib_config_path
if not calib_config_path.is_file(): if not calib_config_path.is_file():
# Disallow implicit config file creation. # Disallow implicit config file creation.
raise ValueError('calib_config_path not found - neither absolute nor relative to in_folder') raise ValueError('calib_config_path not found - neither absolute nor relative to in_folder')
remi = Analysis(calib_config_path, use_hex=not quad_anode) remi = Analysis(calib_config_path, use_hex=not quad_anode)
# Collect required sources and keys required. # Collect required sources and keys required.
sourcekeys = set() sourcekeys = set()
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
sourcekeys |= remi.get_detector_sourcekeys(det_name) sourcekeys |= remi.get_detector_sourcekeys(det_name)
if not reconstruct_ppt: if not reconstruct_ppt:
sourcekeys.add((ppt_source, 'data.bunchPatternTable')) sourcekeys.add((ppt_source, 'data.bunchPatternTable'))
with timing('open_run'): with timing('open_run'):
# Initial opening of input data. # Initial opening of input data.
base_dc = RunDirectory(Path(in_folder) / f'r{run:04d}', inc_suspect_trains=True) base_dc = RunDirectory(Path(in_folder) / f'r{run:04d}', inc_suspect_trains=True)
with timing('select_data'): with timing('select_data'):
# Filter down to those trains with data for all required sources. # Filter down to those trains with data for all required sources.
filter_run = base_dc.select(sourcekeys, require_all=True) filter_run = base_dc.select(sourcekeys, require_all=True)
# Re-select entire data collection to the trains with data. # Re-select entire data collection to the trains with data.
dc = base_dc.select_trains(by_id[filter_run.train_ids]) dc = base_dc.select_trains(by_id[filter_run.train_ids])
base_dc = None base_dc = None
filter_run = None filter_run = None
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Transformation parameters # Transformation parameters
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Additional parameters through the user-side configuration file for analog channels and detector settings. Parameters that are deprecated and ignored, but present in the file, are excluded. Additional parameters through the user-side configuration file for analog channels and detector settings. Parameters that are deprecated and ignored, but present in the file, are excluded.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def print_leaf(leaf, indent=0, ignored_keys={}): def print_leaf(leaf, indent=0, ignored_keys={}):
for key, value in leaf.items(): for key, value in leaf.items():
if key in ignored_keys: if key in ignored_keys:
continue continue
if isinstance(value, dict): if isinstance(value, dict):
print(indent * 4 * ' ' + key) print(indent * 4 * ' ' + key)
print_leaf(value, indent=indent+1, ignored_keys=ignored_keys) print_leaf(value, indent=indent+1, ignored_keys=ignored_keys)
else: else:
print(indent * 4 * ' ' + f'{key}: {value}') print(indent * 4 * ' ' + f'{key}: {value}')
print(calib_config_path.resolve()) print(calib_config_path.resolve())
print_leaf(remi.tree, ignored_keys={'instrument', 'trigger'}) print_leaf(remi.tree, ignored_keys={'instrument', 'trigger'})
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
discr_table = [] discr_table = []
if quad_anode: if quad_anode:
signals = ['u1', 'u2', 'v1', 'v2', 'mcp'] signals = ['u1', 'u2', 'v1', 'v2', 'mcp']
wire_angles = [np.pi*(3/4), np.pi*(1/4)] wire_angles = [np.pi*(3/4), np.pi*(1/4)]
else: else:
signals = ['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp'] signals = ['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']
wire_angles = [np.pi*(3/4), np.pi*(3/4+1/3), np.pi*(3/4+2/3)] wire_angles = [np.pi*(3/4), np.pi*(3/4+1/3), np.pi*(3/4+2/3)]
N = 15 N = 15
shifts = np.linspace(-0.4, 0.4, N) shifts = np.linspace(-0.4, 0.4, N)
for det_name, cur_det in remi['detector'].items(): for det_name, cur_det in remi['detector'].items():
fig = plt.figure(num=f'wiring_{det_name}', figsize=(9, 5)) fig = plt.figure(num=f'wiring_{det_name}', figsize=(9, 5))
fig.text(0.5, 1.0, det_name, ha='center', va='top', size='xx-large') fig.text(0.5, 1.0, det_name, ha='center', va='top', size='xx-large')
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0]) ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
ax.set_axis_off() ax.set_axis_off()
ax.add_patch(Circle((0,0), 1, ec='black', fc='none', lw=2)) ax.add_patch(Circle((0,0), 1, ec='black', fc='none', lw=2))
ax.set_xlim(-1.5*(9/5), 1.5*(9/5)) ax.set_xlim(-1.5*(9/5), 1.5*(9/5))
ax.set_ylim(-1.5, 1.5) ax.set_ylim(-1.5, 1.5)
_, params = remi.get_discriminator(cur_det['channels']) _, params = remi.get_discriminator(cur_det['channels'])
discr_header = params[0].keys() discr_header = params[0].keys()
for channel_idx in range(len(signals)): for channel_idx in range(len(signals)):
index = cur_det['indices'].index(channel_idx) index = cur_det['indices'].index(channel_idx)
discr_table.append((det_name, signals[channel_idx], discr_table.append((det_name, signals[channel_idx],
cur_det['channels'][index], cur_det['channels'][index],
remi['digitizer']['discriminator'], remi['digitizer']['discriminator'],
*params[index].values())) *params[index].values()))
for j, start_angle in enumerate(wire_angles): for j, start_angle in enumerate(wire_angles):
x1 = np.cos(start_angle+np.pi/4) x1 = np.cos(start_angle+np.pi/4)
x2 = np.cos(start_angle+5*np.pi/4) x2 = np.cos(start_angle+5*np.pi/4)
y1 = np.sin(start_angle+np.pi/4) y1 = np.sin(start_angle+np.pi/4)
y2 = np.sin(start_angle+5*np.pi/4) y2 = np.sin(start_angle+5*np.pi/4)
channel = cur_det['channels'][cur_det['indices'].index(2*j)] channel = cur_det['channels'][cur_det['indices'].index(2*j)]
ax.text(x1*1.2, y1*1.2, f'{signals[2*j]}\n{channel}', ax.text(x1*1.2, y1*1.2, f'{signals[2*j]}\n{channel}',
c=f'C{j}', fontsize='xx-large', va='center', ha='center') c=f'C{j}', fontsize='xx-large', va='center', ha='center')
channel = cur_det['channels'][cur_det['indices'].index(2*j+1)] channel = cur_det['channels'][cur_det['indices'].index(2*j+1)]
ax.text(x2*1.2, y2*1.2, f'{signals[2*j+1]}\n{channel}', ax.text(x2*1.2, y2*1.2, f'{signals[2*j+1]}\n{channel}',
c=f'C{j}', fontsize='xx-large', va='center', ha='center') c=f'C{j}', fontsize='xx-large', va='center', ha='center')
for k, shift in enumerate(shifts): for k, shift in enumerate(shifts):
x1 = np.cos(start_angle+np.pi/4+shifts[k]) x1 = np.cos(start_angle+np.pi/4+shifts[k])
x2 = np.cos(start_angle+5*np.pi/4+shifts[N-k-1]) x2 = np.cos(start_angle+5*np.pi/4+shifts[N-k-1])
y1 = np.sin(start_angle+np.pi/4+shifts[k]) y1 = np.sin(start_angle+np.pi/4+shifts[k])
y2 = np.sin(start_angle+5*np.pi/4+shifts[N-k-1]) y2 = np.sin(start_angle+5*np.pi/4+shifts[N-k-1])
ax.plot([x1, x2], [y1, y2], c=f'C{j}') ax.plot([x1, x2], [y1, y2], c=f'C{j}')
mcp_angle = np.pi/6 mcp_angle = np.pi/6
channel = cur_det['channels'][cur_det['indices'].index(6)] channel = cur_det['channels'][cur_det['indices'].index(6)]
ax.text(1.4*np.cos(mcp_angle), 1.2*np.sin(mcp_angle), f'mcp\n{channel}', ax.text(1.4*np.cos(mcp_angle), 1.2*np.sin(mcp_angle), f'mcp\n{channel}',
c='k', fontsize='xx-large', va='center', ha='center') c='k', fontsize='xx-large', va='center', ha='center')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
display(Latex(tabulate.tabulate( display(Latex(tabulate.tabulate(
discr_table, tablefmt='latex', headers=['', '', '', 'discriminator', *discr_header]))) discr_table, tablefmt='latex', headers=['', '', '', 'discriminator', *discr_header])))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Pulse and trigger information # Pulse and trigger information
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Read PPT from file or reconstruct PPT for older data ### Read PPT from file or reconstruct PPT for older data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
if reconstruct_ppt: if reconstruct_ppt:
# Take up to the first hundred trains for now. # Take up to the first hundred trains for now.
# Could be done for each train individually, but likely not necessary for now. # Could be done for each train individually, but likely not necessary for now.
trigger_trace = dc[remi['digitizer']['source'], remi['digitizer']['key_pattern'].format(trigger_edge_channel)] \ trigger_trace = dc[remi['digitizer']['source'], remi['digitizer']['key_pattern'].format(trigger_edge_channel)] \
[:100].ndarray().mean(axis=0).astype(np.float64) [:100].ndarray().mean(axis=0).astype(np.float64)
trigger_trace -= trigger_trace[0] # Use simple offset correction. trigger_trace -= trigger_trace[0] # Use simple offset correction.
fake_ppt = np.zeros(2700, dtype=np.uint32) fake_ppt = np.zeros(2700, dtype=np.uint32)
discr_func, discr_params = remi.get_discriminator([trigger_edge_channel]) discr_func, discr_params = remi.get_discriminator([trigger_edge_channel])
edges = np.zeros(1000, dtype=np.float64) edges = np.zeros(1000, dtype=np.float64)
num_pulses = discr_func(trigger_trace, edges=edges, **discr_params[0]) num_pulses = discr_func(trigger_trace, edges=edges, **discr_params[0])
edges = edges[:num_pulses] edges = edges[:num_pulses]
first_edge = edges[0] first_edge = edges[0]
rel_edges = np.round(edges - first_edge) rel_edges = np.round(edges - first_edge)
edge_diff = rel_edges[1] - rel_edges[0] edge_diff = rel_edges[1] - rel_edges[0]
if not np.allclose(rel_edges[1:] - rel_edges[:-1], edge_diff): if not np.allclose(rel_edges[1:] - rel_edges[:-1], edge_diff):
raise ValueError('PPT reconstruction for unstable edge intervals not supported') raise ValueError('PPT reconstruction for unstable edge intervals not supported')
pulse_spacing = edge_diff / (2 * remi['digitizer']['clock_factor']) # In units of PPT pulse_spacing = edge_diff / (2 * remi['digitizer']['clock_factor']) # In units of PPT
if not float.is_integer(pulse_spacing): if not float.is_integer(pulse_spacing):
raise ValueError('PPT reconstruction encountered non-integer pulse spacing') raise ValueError('PPT reconstruction encountered non-integer pulse spacing')
pulse_spacing = int(pulse_spacing) pulse_spacing = int(pulse_spacing)
# Taken from euxfel_bunch_pattern/__init__.py # Taken from euxfel_bunch_pattern/__init__.py
from euxfel_bunch_pattern import DESTINATION_T4D, DESTINATION_T5D, PHOTON_LINE_DEFLECTION from euxfel_bunch_pattern import DESTINATION_T4D, DESTINATION_T5D, PHOTON_LINE_DEFLECTION
if instrument_sase == 1: if instrument_sase == 1:
flag = DESTINATION_T4D flag = DESTINATION_T4D
elif instrument_sase == 2: elif instrument_sase == 2:
flag = DESTINATION_T5D flag = DESTINATION_T5D
elif instrument_sase == 3: elif instrument_sase == 3:
flag = DESTINATION_T4D | PHOTON_LINE_DEFLECTION flag = DESTINATION_T4D | PHOTON_LINE_DEFLECTION
first_pulse_offset = int(first_edge + trigger_edge_offset) # Overwrite notebook argument. first_pulse_offset = int(first_edge + trigger_edge_offset) # Overwrite notebook argument.
fake_ppt[fake_ppt_offset:fake_ppt_offset + (pulse_spacing * num_pulses):pulse_spacing] = flag fake_ppt[fake_ppt_offset:fake_ppt_offset + (pulse_spacing * num_pulses):pulse_spacing] = flag
from pasha.functor import Functor, gen_split_slices from pasha.functor import Functor, gen_split_slices
class FakeKeyDataFunctor(Functor): class FakeKeyDataFunctor(Functor):
"""Functor appearing KeyData-like with constant data. """Functor appearing KeyData-like with constant data.
This functor serves a constant data row for a given number This functor serves a constant data row for a given number
of train IDs the same way a KeyData object would. of train IDs the same way a KeyData object would.
""" """
def __init__(self, row, train_ids): def __init__(self, row, train_ids):
self.row = row self.row = row
self.train_ids = train_ids self.train_ids = train_ids
def split(self, num_workers): def split(self, num_workers):
return gen_split_slices(len(self.train_ids), n_parts=num_workers) return gen_split_slices(len(self.train_ids), n_parts=num_workers)
def iterate(self, share): def iterate(self, share):
it = zip(range(*share.indices(len(self.train_ids))), self.train_ids) it = zip(range(*share.indices(len(self.train_ids))), self.train_ids)
for index, train_id in it: for index, train_id in it:
yield index, train_id, self.row yield index, train_id, self.row
ppt_data = FakeKeyDataFunctor(fake_ppt, dc.train_ids) ppt_data = FakeKeyDataFunctor(fake_ppt, dc.train_ids)
fig, ax = plt.subplots(num='reconstructed_ppt_triggers', figsize=(9, 6), clear=True, ncols=1, nrows=1) fig, ax = plt.subplots(num='reconstructed_ppt_triggers', figsize=(9, 6), clear=True, ncols=1, nrows=1)
ax.set_title('Edge trigger signal') ax.set_title('Edge trigger signal')
ax.plot(trigger_trace, lw=1, label=f'Mean {trigger_edge_channel} trace') ax.plot(trigger_trace, lw=1, label=f'Mean {trigger_edge_channel} trace')
ax.vlines(edges, trigger_trace.min()*1.1, trigger_trace.max()*1.1, ax.vlines(edges, trigger_trace.min()*1.1, trigger_trace.max()*1.1,
color='red', linewidth=3, alpha=0.3, label='Edge positions') color='red', linewidth=3, alpha=0.3, label='Edge positions')
ax.set_xlabel('Samples') ax.set_xlabel('Samples')
ax.set_ylabel('Intensity / ADU') ax.set_ylabel('Intensity / ADU')
ax.legend() ax.legend()
else: else:
ppt_data = dc[ppt_source, 'data.bunchPatternTable'] ppt_data = dc[ppt_source, 'data.bunchPatternTable']
print(f'Pulse pattern entries for {(ppt_data.data_counts() > 0).sum()} trains') print(f'Pulse pattern entries for {(ppt_data.data_counts() > 0).sum()} trains')
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Count pulses per train and compute offsets ### Count pulses per train and compute offsets
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Based on the pulse pattern tables, three global variables are obtained: # Based on the pulse pattern tables, three global variables are obtained:
# #
# * `pulse_counts [int32: len(dc.train_ids)]` containing the number of pulses per train. # * `pulse_counts [int32: len(dc.train_ids)]` containing the number of pulses per train.
# * `pulse_offsets [int32: len(dc.train_ids)]` containing the global offset for the first pulse of each train. # * `pulse_offsets [int32: len(dc.train_ids)]` containing the global offset for the first pulse of each train.
# * `num_pulses = pulse_counts.sum(axis=0)` # * `num_pulses = pulse_counts.sum(axis=0)`
def get_pulse_positions(ppt, sase, laser, ppl_offset): def get_pulse_positions(ppt, sase, laser, ppl_offset):
# Combine FEL and PPL positions. # Combine FEL and PPL positions.
fel_pos = indices_at_sase(ppt, sase) if not ignore_fel else np.array([]) fel_pos = indices_at_sase(ppt, sase) if not ignore_fel else np.array([])
ppl_pos = indices_at_laser(ppt, laser) if not ignore_ppl else np.array([]) ppl_pos = indices_at_laser(ppt, laser) if not ignore_ppl else np.array([])
if len(fel_pos) > 0: if len(fel_pos) > 0:
# Move PPL up to the FEL position. # Move PPL up to the FEL position.
ppl_pos += fel_pos[0] + ppl_offset ppl_pos += fel_pos[0] + ppl_offset
return np.union1d(fel_pos, ppl_pos), fel_pos, ppl_pos return np.union1d(fel_pos, ppl_pos), fel_pos, ppl_pos
if laser_ppt_mask < 0: if laser_ppt_mask < 0:
# If laser PPT mask is not specified, try to figure it out from device IDs. # If laser PPT mask is not specified, try to figure it out from device IDs.
from euxfel_bunch_pattern import PPL_BITS from euxfel_bunch_pattern import PPL_BITS
instrument = karabo_id[:karabo_id.index('_')] instrument = karabo_id[:karabo_id.index('_')]
try: try:
laser_ppt_mask = PPL_BITS[f'LP_{instrument}'] laser_ppt_mask = PPL_BITS[f'LP_{instrument}']
except KeyError: except KeyError:
raise ValueError(f'Laser PPT mask unknown for instrument `{instrument}`') raise ValueError(f'Laser PPT mask unknown for instrument `{instrument}`')
with timing('pulse_info'): with timing('pulse_info'):
psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_triggers)) psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_triggers))
# Build the pulse index # Build the pulse index
pulse_counts = psh.alloc(shape=len(dc.train_ids), dtype=np.uint64) pulse_counts = psh.alloc(shape=len(dc.train_ids), dtype=np.uint64)
has_ppt = psh.alloc(shape=len(dc.train_ids), dtype=bool, fill=False) has_ppt = psh.alloc(shape=len(dc.train_ids), dtype=bool, fill=False)
def count_pulses(wid, index, tid, ppt): def count_pulses(wid, index, tid, ppt):
pulse_counts[index] = len(get_pulse_positions(ppt, instrument_sase, laser_ppt_mask, ppl_offset)[0]) pulse_counts[index] = len(get_pulse_positions(ppt, instrument_sase, laser_ppt_mask, ppl_offset)[0])
has_ppt[index] = True has_ppt[index] = True
psh.map(count_pulses, ppt_data) psh.map(count_pulses, ppt_data)
# Fill any missing values with the highest. # Fill any missing values with the highest.
pulse_counts[has_ppt == False] = pulse_counts.max() pulse_counts[has_ppt == False] = pulse_counts.max()
if trailing_trigger: if trailing_trigger:
# Add a single count to every train for trailing trigger. # Add a single count to every train for trailing trigger.
warning('Trailing trigger active, all pulse counts are one higher than expected') warning('Trailing trigger active, all pulse counts are one higher than expected')
pulse_counts += 1 pulse_counts += 1
# Compute offsets based on pulse counts. # Compute offsets based on pulse counts.
pulse_offsets = np.zeros_like(pulse_counts) pulse_offsets = np.zeros_like(pulse_counts)
np.cumsum(pulse_counts[:-1], out=pulse_offsets[1:]) np.cumsum(pulse_counts[:-1], out=pulse_offsets[1:])
# Total number of pulses. # Total number of pulses.
num_pulses = int(pulse_counts.sum()) num_pulses = int(pulse_counts.sum())
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, ax = plt.subplots(num='pulse_counts', ncols=1, nrows=1, figsize=(9, 4), clear=True) fig, ax = plt.subplots(num='pulse_counts', ncols=1, nrows=1, figsize=(9, 4), clear=True)
ax.set_title('Pulse count') ax.set_title('Pulse count')
ax.plot(dc.train_ids, pulse_counts, lw=1) ax.plot(dc.train_ids, pulse_counts, lw=1)
ax.set_xlabel('Train ID') ax.set_xlabel('Train ID')
ax.set_ylabel('Number of pulses') ax.set_ylabel('Number of pulses')
ax.set_ylim(0, max(300, pulse_counts.max() + 10)) ax.set_ylim(0, max(300, pulse_counts.max() + 10))
ax.ticklabel_format(style='plain') ax.ticklabel_format(style='plain')
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Find triggers ### Find triggers
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# A trigger defines the boundary of a pulse on the digitizer trace stored per train. This cell creates a # A trigger defines the boundary of a pulse on the digitizer trace stored per train. This cell creates a
# global variable: # global variable:
# * `triggers [(start: int32, stop: int32, offset: float64, fel: bool, ppl: bool): num_pulses]` # * `triggers [(start: int32, stop: int32, offset: float64, fel: bool, ppl: bool): num_pulses]`
# containing the triggers for each pulse. # containing the triggers for each pulse.
# #
# This uses the pulse puttern table to locate the pulse positions on the trace. Only number of pulses and # This uses the pulse puttern table to locate the pulse positions on the trace. Only number of pulses and
# their distance can be drawn this way, leaving the absolute offset for the very first pulse to be # their distance can be drawn this way, leaving the absolute offset for the very first pulse to be
# configured via `trigger/ppt/first_pulse_offset`. If a PPL is used, it will be included in the trigger # configured via `trigger/ppt/first_pulse_offset`. If a PPL is used, it will be included in the trigger
# pattern. The ppt_offset parameter allows taking into account an offset betwen PPL and FEL. # pattern. The ppt_offset parameter allows taking into account an offset betwen PPL and FEL.
psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_triggers)) psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_triggers))
triggers = psh.alloc(shape=(num_pulses,), dtype=trigger_dt, fill=(-1, -1, np.nan, -1, 0, 0)) triggers = psh.alloc(shape=(num_pulses,), dtype=trigger_dt, fill=(-1, -1, np.nan, -1, 0, 0))
clock_factor = remi['digitizer']['clock_factor'] clock_factor = remi['digitizer']['clock_factor']
min_trace_len = min([ min_trace_len = min([
dc[src, key].entry_shape[0] for det_name in remi['detector'].keys() dc[src, key].entry_shape[0] for det_name in remi['detector'].keys()
for src, key in remi.get_detector_sourcekeys(det_name) for src, key in remi.get_detector_sourcekeys(det_name)
]) ])
def trigger_by_ppt(worker_id, index, train_id, ppt): def trigger_by_ppt(worker_id, index, train_id, ppt):
all_pos, fel_pos, ppl_pos = get_pulse_positions(ppt, instrument_sase, laser_ppt_mask, ppl_offset) all_pos, fel_pos, ppl_pos = get_pulse_positions(ppt, instrument_sase, laser_ppt_mask, ppl_offset)
num_pulses = len(all_pos) num_pulses = len(all_pos)
if num_pulses > 0: if num_pulses > 0:
if len(ppl_pos) == 0 and ppl_offset < 0: if len(ppl_pos) == 0 and ppl_offset < 0:
# No PPL pulses, but a negative offset is configured. This will cause # No PPL pulses, but a negative offset is configured. This will cause
# first_pulse_offset to start early and most likely miss pulses at the # first_pulse_offset to start early and most likely miss pulses at the
# end, so we correct by adding the ppl_offset to relative positions # end, so we correct by adding the ppl_offset to relative positions
# when computing trace positions. # when computing trace positions.
pos_corr = abs(ppl_offset) pos_corr = abs(ppl_offset)
else: else:
pos_corr = 0 pos_corr = 0
rel_pos = all_pos - all_pos[0] rel_pos = all_pos - all_pos[0]
if num_pulses > 1: if num_pulses > 1:
pulse_len = np.unique(rel_pos[1:] - rel_pos[:-1]).min() pulse_len = np.unique(rel_pos[1:] - rel_pos[:-1]).min()
elif num_pulses == 1: elif num_pulses == 1:
pulse_len = single_pulse_length pulse_len = single_pulse_length
start_frac = first_pulse_offset + (rel_pos + pos_corr) * 2 * clock_factor start_frac = first_pulse_offset + (rel_pos + pos_corr) * 2 * clock_factor
start_int = start_frac.astype(int) start_int = start_frac.astype(int)
train_triggers = triggers[pulse_offsets[index]:int(pulse_offsets[index]+num_pulses)] train_triggers = triggers[pulse_offsets[index]:int(pulse_offsets[index]+num_pulses)]
train_triggers['start'] = start_int + pulse_start_offset train_triggers['start'] = start_int + pulse_start_offset
train_triggers['stop'] = start_int + int(pulse_len * 2 * clock_factor) - 1 + pulse_end_offset train_triggers['stop'] = start_int + int(pulse_len * 2 * clock_factor) - 1 + pulse_end_offset
train_triggers['offset'] = start_frac - start_int train_triggers['offset'] = start_frac - start_int
train_triggers['pulse'] = all_pos.astype(np.int16) train_triggers['pulse'] = all_pos.astype(np.int16)
train_triggers['fel'] = [pos in fel_pos for pos in all_pos] train_triggers['fel'] = [pos in fel_pos for pos in all_pos]
train_triggers['ppl'] = [pos in ppl_pos for pos in all_pos] train_triggers['ppl'] = [pos in ppl_pos for pos in all_pos]
last_sample = train_triggers['stop'].max() last_sample = train_triggers['stop'].max()
else: else:
last_sample = first_pulse_offset last_sample = first_pulse_offset
if trailing_trigger: if trailing_trigger:
# Add trailing trigger if required. # Add trailing trigger if required.
trigger = triggers[int(pulse_offsets[index]+pulse_counts[index]-1)] trigger = triggers[int(pulse_offsets[index]+pulse_counts[index]-1)]
trigger['start'] = last_sample trigger['start'] = last_sample
trigger['stop'] = min_trace_len trigger['stop'] = min_trace_len
trigger['offset'] = 0.0 trigger['offset'] = 0.0
trigger['pulse'] = -1 trigger['pulse'] = -1
trigger['fel'] = False trigger['fel'] = False
trigger['ppl'] = False trigger['ppl'] = False
with timing('find_triggers'): with timing('find_triggers'):
psh.map(trigger_by_ppt, ppt_data) psh.map(trigger_by_ppt, ppt_data)
if (np.unique(triggers['pulse'][1:] - triggers['pulse'][:-1]) > 0).sum() > 1: if (np.unique(triggers['pulse'][1:] - triggers['pulse'][:-1]) > 0).sum() > 1:
# There is more than one delta between pulse entries across all pulses. This is not # There is more than one delta between pulse entries across all pulses. This is not
# necessarily a problem, as the pattern could simply have changed in between trains # necessarily a problem, as the pattern could simply have changed in between trains
# with each train being split properly. # with each train being split properly.
# If there's more than one delta in a single train, this likely points to a mismatch # If there's more than one delta in a single train, this likely points to a mismatch
# of FEL and PPL repetition rate. This is most likely not intended. # of FEL and PPL repetition rate. This is most likely not intended.
one = np.uint64(1) # Because np.uint64 + int = np.float64 one = np.uint64(1) # Because np.uint64 + int = np.float64
pulse_deltas = set() pulse_deltas = set()
for pulse_id, (offset, count) in enumerate(zip( for pulse_id, (offset, count) in enumerate(zip(
pulse_offsets, pulse_counts - one if trailing_trigger else pulse_counts pulse_offsets, pulse_counts - one if trailing_trigger else pulse_counts
)): )):
deltas = triggers['pulse'][offset+one:offset+count] - triggers['pulse'][offset:offset+count-one] deltas = triggers['pulse'][offset+one:offset+count] - triggers['pulse'][offset:offset+count-one]
if len(np.unique(deltas)) > 1: if len(np.unique(deltas)) > 1:
for delta in deltas: for delta in deltas:
pulse_deltas.add(delta) pulse_deltas.add(delta)
if len(pulse_deltas) > 1: if len(pulse_deltas) > 1:
delta_str = ', '.join([str(x) for x in sorted(pulse_deltas)]) delta_str = ', '.join([str(x) for x in sorted(pulse_deltas)])
warning(f'Different pulse lengths (PPT: {delta_str}) encountered within single trains, ' warning(f'Different pulse lengths (PPT: {delta_str}) encountered within single trains, '
f'separated pulse spectra may split up signals!') f'separated pulse spectra may split up signals!')
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, (lx, rx) = plt.subplots(num='trigger_positions', ncols=2, nrows=1, figsize=(9, 4), clear=True, fig, (lx, rx) = plt.subplots(num='trigger_positions', ncols=2, nrows=1, figsize=(9, 4), clear=True,
gridspec_kw=dict(top=0.75)) gridspec_kw=dict(top=0.75))
# Display ~400 pulses or 10 trains, whatever is lower # Display ~400 pulses or 10 trains, whatever is lower
n_trains = max(abs(pulse_offsets - 200).argmin(), 5) n_trains = max(abs(pulse_offsets - 200).argmin(), 5)
visible_triggers = triggers[:pulse_offsets[n_trains]] visible_triggers = triggers[:pulse_offsets[n_trains]]
pulse_index = np.arange(len(visible_triggers)) pulse_index = np.arange(len(visible_triggers))
pumped = visible_triggers['fel'] & visible_triggers['ppl'] pumped = visible_triggers['fel'] & visible_triggers['ppl']
fel_only = visible_triggers['fel'] & ~pumped fel_only = visible_triggers['fel'] & ~pumped
ppl_only = visible_triggers['ppl'] & ~pumped ppl_only = visible_triggers['ppl'] & ~pumped
lx.plot(pulse_index[pumped], visible_triggers[pumped]['start'], ' .', ms=3, c='C0', label='FEL+PPL') lx.plot(pulse_index[pumped], visible_triggers[pumped]['start'], ' .', ms=3, c='C0', label='FEL+PPL')
lx.plot(pulse_index[fel_only], visible_triggers[fel_only]['start'], '.', ms=3, c='C1', label='FEL-only') lx.plot(pulse_index[fel_only], visible_triggers[fel_only]['start'], '.', ms=3, c='C1', label='FEL-only')
lx.plot(pulse_index[ppl_only], visible_triggers[ppl_only]['start'], '.', ms=2, c='C2', label='PPL-only') lx.plot(pulse_index[ppl_only], visible_triggers[ppl_only]['start'], '.', ms=2, c='C2', label='PPL-only')
max_start = visible_triggers['start'].max() max_start = visible_triggers['start'].max()
lx.vlines(pulse_offsets[:n_trains], 0, max_start, color='grey', linewidth=1, alpha=0.2) lx.vlines(pulse_offsets[:n_trains], 0, max_start, color='grey', linewidth=1, alpha=0.2)
lx.tick_params(right=True) lx.tick_params(right=True)
lx.set_xlabel('Pulse index') lx.set_xlabel('Pulse index')
lx.set_xlim(-15, pulse_offsets[n_trains]+15) lx.set_xlim(-15, pulse_offsets[n_trains]+15)
lx.set_ylabel('Trigger position') lx.set_ylabel('Trigger position')
lx.set_ylim(-max_start // 20, max_start + max_start // 20) lx.set_ylim(-max_start // 20, max_start + max_start // 20)
lx.legend(fontsize='small', loc='lower right') lx.legend(fontsize='small', loc='lower right')
train_lx = lx.twiny() train_lx = lx.twiny()
train_lx.set_xlabel('Train ID', labelpad=8) train_lx.set_xlabel('Train ID', labelpad=8)
train_lx.set_xlim(lx.get_xlim()) train_lx.set_xlim(lx.get_xlim())
train_lx.set_xticks(pulse_offsets[:n_trains]) train_lx.set_xticks(pulse_offsets[:n_trains])
train_lx.set_xticklabels([str(int(x)) for x in dc.train_ids[:n_trains]], train_lx.set_xticklabels([str(int(x)) for x in dc.train_ids[:n_trains]],
rotation=-45, fontsize='x-small') rotation=-45, fontsize='x-small')
rx.plot(triggers['start'], lw=0.2) rx.plot(triggers['start'], lw=0.2)
rx.set_xlabel('Pulse index') rx.set_xlabel('Pulse index')
rx.tick_params(left=False, labelleft=False, right=True, labelright=True) rx.tick_params(left=False, labelleft=False, right=True, labelright=True)
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The trigger defines the boundary of each pulse on the digitizer trace acquired by train. The starting position in samples of each found trigger is shown for the first few trains in detail on the left and all trains on the right. The trigger defines the boundary of each pulse on the digitizer trace acquired by train. The starting position in samples of each found trigger is shown for the first few trains in detail on the left and all trains on the right.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Analog signal to digital edges # Analog signal to digital edges
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Find edges in analog signal ### Find edges in analog signal
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_edges)) psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_edges))
threadpool_limits(limits=remi.get_num_workers(mt_avg_trace)) threadpool_limits(limits=remi.get_num_workers(mt_avg_trace))
det_data = {} det_data = {}
for det_name, det in remi['detector'].items(): for det_name, det in remi['detector'].items():
det_sourcekeys = remi.get_detector_sourcekeys(det_name) det_sourcekeys = remi.get_detector_sourcekeys(det_name)
det_get_traces = remi.get_traces_getter(det_name) det_get_traces = remi.get_traces_getter(det_name)
trace_len = dc[next(iter(det_sourcekeys))].entry_shape[0] trace_len = dc[next(iter(det_sourcekeys))].entry_shape[0]
edges = psh.alloc(shape=(num_pulses, 7, det['max_hits']), edges = psh.alloc(shape=(num_pulses, 7, det['max_hits']),
dtype=np.float64, fill=np.nan) dtype=np.float64, fill=np.nan)
amplitudes = psh.alloc(shape=(num_pulses, 7, det['max_hits']), amplitudes = psh.alloc(shape=(num_pulses, 7, det['max_hits']),
dtype=np.float64, fill=np.nan) dtype=np.float64, fill=np.nan)
avg_traces = psh.alloc_per_worker(shape=(7, trace_len), dtype=np.float64) avg_traces = psh.alloc_per_worker(shape=(7, trace_len), dtype=np.float64)
def prepare_edge_worker(worker_id): def prepare_edge_worker(worker_id):
correct_func = remi.get_baseline_corrector() correct_func = remi.get_baseline_corrector()
discr_func, discr_params = remi.get_discriminator(det['channels']) discr_func, discr_params = remi.get_discriminator(det['channels'])
source_name = remi['digitizer']['source'] source_name = remi['digitizer']['source']
bl_start, bl_stop, _ = remi.get_baseline_limits(trace_len) bl_start, bl_stop, _ = remi.get_baseline_limits(trace_len)
bl_sym = remi['digitizer']['baseline_symmetry'] bl_sym = remi['digitizer']['baseline_symmetry']
time_cal = 1e9 / (2 * remi['digitizer']['clock_factor'] * (1.3e9 / 288)) time_cal = 1e9 / (2 * remi['digitizer']['clock_factor'] * (1.3e9 / 288))
traces_corr = np.empty((7, trace_len), dtype=np.float64) traces_corr = np.empty((7, trace_len), dtype=np.float64)
baselines = np.empty(bl_sym, dtype=np.float64) baselines = np.empty(bl_sym, dtype=np.float64)
yield yield
@psh.with_init(prepare_edge_worker) @psh.with_init(prepare_edge_worker)
def find_edges(worker_id, index, train_id, data): def find_edges(worker_id, index, train_id, data):
try: try:
data = det_get_traces(data[source_name]) data = det_get_traces(data[source_name])
except KeyError: except KeyError:
return return
for channel_idx in range(7): for channel_idx in range(7):
correct_func(data[channel_idx], traces_corr[channel_idx], correct_func(data[channel_idx], traces_corr[channel_idx],
baselines, bl_start, bl_stop) baselines, bl_start, bl_stop)
avg_traces[worker_id] += traces_corr avg_traces[worker_id] += traces_corr
pulses_slice = np.s_[pulse_offsets[index]:pulse_offsets[index]+pulse_counts[index]] pulses_slice = np.s_[pulse_offsets[index]:pulse_offsets[index]+pulse_counts[index]]
for trigger, pulse_edges, pulse_amplitudes in zip( for trigger, pulse_edges, pulse_amplitudes in zip(
triggers[pulses_slice], edges[pulses_slice], amplitudes[pulses_slice] triggers[pulses_slice], edges[pulses_slice], amplitudes[pulses_slice]
): ):
trigger_slice = np.s_[trigger['start']:trigger['stop']] trigger_slice = np.s_[trigger['start']:trigger['stop']]
for trace, channel_params, channel_edges, channel_amplitudes in zip( for trace, channel_params, channel_edges, channel_amplitudes in zip(
traces_corr, discr_params, pulse_edges, pulse_amplitudes traces_corr, discr_params, pulse_edges, pulse_amplitudes
): ):
discr_func(trace[trigger_slice], edges=channel_edges, discr_func(trace[trigger_slice], edges=channel_edges,
amplitudes=channel_amplitudes, **channel_params) amplitudes=channel_amplitudes, **channel_params)
if np.isfinite(pulse_edges).sum(axis=1).max() == det['max_hits']: if np.isfinite(pulse_edges).sum(axis=1).max() == det['max_hits']:
warning(f'Maximum number of edges reached in train {train_id}, pulse: {trigger["pulse"]}') warning(f'Maximum number of edges reached in train {train_id}, pulse: {trigger["pulse"]}')
with timing(f'find_edges, {det_name}'): with timing(f'find_edges, {det_name}'):
psh.map(find_edges, dc.select(det_sourcekeys)) psh.map(find_edges, dc.select(det_sourcekeys))
# Properly offset edges to their trigger offset and convert to time.
# This is not done earlier to preserve the information for plotting.
edges += triggers['offset'][:, None, None]
edges *= remi.get_time_calibration()
det_data[det_name] = {
'edges': edges,
'amplitudes': amplitudes,
'avg_trace': avg_traces.sum(axis=0) / len(dc.train_ids)
}
if not np.isfinite(edges).any(): if not np.isfinite(edges).any():
warning(f'No edges found for {det_name}') warning(f'No edges found for {det_name}')
continue
fig, (ux, bx) = plt.subplots(num=f'digitize_result_{det_name}', ncols=1, nrows=2, figsize=(9.5, 8), clear=True, fig, (ux, bx) = plt.subplots(num=f'digitize_result_{det_name}', ncols=1, nrows=2, figsize=(9.5, 8), clear=True,
gridspec_kw=dict(left=0.1, right=0.98, top=0.98, bottom=0.1, hspace=0.25)) gridspec_kw=dict(left=0.1, right=0.98, top=0.98, bottom=0.1, hspace=0.25))
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
max_num = 0 max_num = 0
for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']): for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
n, _, _ = ux.hist(finite_flattened_slice(amplitudes, np.s_[:, edge_idx, :]), n, _, _ = ux.hist(finite_flattened_slice(amplitudes, np.s_[:, edge_idx, :]),
bins=1000, range=(0, 4096), histtype='step', lw=1, bins=1000, range=(0, 4096), histtype='step', lw=1,
color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name) color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)
max_num = max(max_num, n.max()) max_num = max(max_num, n.max())
cur_edges = finite_flattened_slice(edges, np.s_[:, edge_idx, :]) cur_edges = finite_flattened_slice(edges, np.s_[:, edge_idx, :])
bx.hist(cur_edges - np.floor(cur_edges), bins=500, range=(0, 1), histtype='step', bx.hist(cur_edges - np.floor(cur_edges), bins=500, range=(0, 1), histtype='step',
lw=1, color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name) lw=1, color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)
ux.legend() ux.legend()
ux.set_title('Pulse height distributions') ux.set_title('Pulse height distributions')
ux.set_xlabel('Pulse height') ux.set_xlabel('Pulse height')
ux.set_yscale('log') ux.set_yscale('log')
ux.set_xlim(0, 4096) ux.set_xlim(0, 4096)
ux.set_ylim(10, 1.5*max(max_num, 10)) ux.set_ylim(10, 1.5*max(max_num, 10))
if remi['digitizer']['discriminator'] == 'cfd': if remi['digitizer']['discriminator'] == 'cfd':
ux.text(1024, 12.5, 'No pulse height feedback for constant fraction discrimination', ux.text(1024, 12.5, 'No pulse height feedback for constant fraction discrimination',
ha='center', va='center') ha='center', va='center')
bx.set_title('Fractional edge distributions') bx.set_title('Fractional edge distributions')
bx.set_xlabel('Edge positions - ⌊edge positions⌋') bx.set_xlabel('Edge positions - ⌊edge positions⌋')
bx.set_yscale('log') bx.set_yscale('log')
bx.set_xlim(-0.05, 1.2) bx.set_xlim(-0.05, 1.2)
bx.legend() bx.legend()
# Properly offset edges to their trigger offset and convert to time.
# This is not done earlier to preserve the information for plotting.
edges += triggers['offset'][:, None, None]
edges *= remi.get_time_calibration()
det_data[det_name] = {
'edges': edges,
'amplitudes': amplitudes,
'avg_trace': avg_traces.sum(axis=0) / len(dc.train_ids)
}
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
The analog signal is digitized into discrete edges using a fast timing discriminator. The result of this operation is available in files in the `raw.triggers` dataset. The analog signal is digitized into discrete edges using a fast timing discriminator. The result of this operation is available in files in the `raw.triggers` dataset.
The pulse height distribution is an integral view about the chosen digitization thresholds. For more detail, please refer to the spectral pulse height distributions further below. The pulse height distribution is an integral view about the chosen digitization thresholds. For more detail, please refer to the spectral pulse height distributions further below.
The fractional edge distribution visualizes the interpolated component of edge positions, i.e. between discrete digitizer samples. This should in general be flat, in particular a convex shape indicates poor interpolation due to too fast rise times. The fractional edge distribution visualizes the interpolated component of edge positions, i.e. between discrete digitizer samples. This should in general be flat, in particular a convex shape indicates poor interpolation due to too fast rise times.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Global average of analog signals ### Global average of analog signals
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
fig, axs = plt.subplots(num=f'global_average_{det_name}', nrows=7, figsize=(9.5, 8), clear=True, fig, axs = plt.subplots(num=f'global_average_{det_name}', nrows=7, figsize=(9.5, 8), clear=True,
gridspec_kw=dict(left=0.1, right=0.98, top=0.98, bottom=0.1)) gridspec_kw=dict(left=0.1, right=0.98, top=0.98, bottom=0.1))
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']): for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
axs[edge_idx].plot(det_data[det_name]['avg_trace'][edge_idx], lw=1) axs[edge_idx].plot(det_data[det_name]['avg_trace'][edge_idx], lw=1)
axs[edge_idx].tick_params(labelbottom=False) axs[edge_idx].tick_params(labelbottom=False)
axs[edge_idx].set_ylabel(edge_name) axs[edge_idx].set_ylabel(edge_name)
axs[-1].tick_params(labelbottom=True) axs[-1].tick_params(labelbottom=True)
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Sample for found edges ### Sample for found edges
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
edges = det_data[det_name]['edges'] edges = det_data[det_name]['edges']
fig = plt.figure(num=f'edge_samples_{det_name}', figsize=(9.5, 8)) fig = plt.figure(num=f'edge_samples_{det_name}', figsize=(9.5, 8))
grid = fig.add_gridspec(ncols=2, nrows=4, left=0.1, right=0.98, top=0.98, bottom=0.1) grid = fig.add_gridspec(ncols=2, nrows=4, left=0.1, right=0.98, top=0.98, bottom=0.1)
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
for signal_idx, signal_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']): for signal_idx, signal_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
row = (1 + signal_idx // 2) if signal_idx < 6 else 0 row = (1 + signal_idx // 2) if signal_idx < 6 else 0
col = (signal_idx % 2) if signal_idx < 6 else np.s_[:] col = (signal_idx % 2) if signal_idx < 6 else np.s_[:]
ax = fig.add_subplot(grid[row, col]) ax = fig.add_subplot(grid[row, col])
finite_edges = np.isfinite(edges[:, signal_idx, 0]) finite_edges = np.isfinite(edges[:, signal_idx, 0])
if not finite_edges.any(): if not finite_edges.any():
warning(f'No edges found for {det_name}/{signal_name}') warning(f'No edges found for {det_name}/{signal_name}')
continue continue
pulse_idx = np.uint64(finite_edges.nonzero()[0][0]) # Is combined with other uint64 values below. pulse_idx = np.uint64(finite_edges.nonzero()[0][0]) # Is combined with other uint64 values below.
train_idx = (pulse_idx >= pulse_offsets).nonzero()[0][-1] train_idx = (pulse_idx >= pulse_offsets).nonzero()[0][-1]
trigger = triggers[pulse_idx] trigger = triggers[pulse_idx]
sourcekey = remi.get_channel_sourcekey( sourcekey = remi.get_channel_sourcekey(
remi['detector'][det_name]['channels'][signal_idx]) remi['detector'][det_name]['channels'][signal_idx])
train_trace = dc[sourcekey].select_trains(np.s_[train_idx:train_idx+1]).ndarray()[0] train_trace = dc[sourcekey].select_trains(np.s_[train_idx:train_idx+1]).ndarray()[0]
corr_trace = np.zeros_like(train_trace, dtype=np.float64) corr_trace = np.zeros_like(train_trace, dtype=np.float64)
remi.get_baseline_corrector()( remi.get_baseline_corrector()(
train_trace, corr_trace, train_trace, corr_trace,
np.empty(remi['digitizer']['baseline_symmetry'], dtype=np.float64), np.empty(remi['digitizer']['baseline_symmetry'], dtype=np.float64),
*remi.get_baseline_limits(len(train_trace))[:2]) *remi.get_baseline_limits(len(train_trace))[:2])
pulse_trace = corr_trace[np.s_[trigger['start']:trigger['stop']]] pulse_trace = corr_trace[np.s_[trigger['start']:trigger['stop']]]
x_time = remi.get_time_calibration() * (np.arange(len(pulse_trace)) + trigger['offset']) x_time = remi.get_time_calibration() * (np.arange(len(pulse_trace)) + trigger['offset'])
ax.plot(x_time, pulse_trace, lw=1) ax.plot(x_time, pulse_trace, lw=1)
ax.set_xlim(x_time[0], x_time[-1]) ax.set_xlim(x_time[0], x_time[-1])
ax.set_ylim(-200, pulse_trace.max()*1.1) ax.set_ylim(-200, pulse_trace.max()*1.1)
ax.text(x_time[-1], pulse_trace.max(), ax.text(x_time[-1], pulse_trace.max(),
f'T{train_idx} P{pulse_idx - pulse_offsets[train_idx]} ', f'T{train_idx} P{pulse_idx - pulse_offsets[train_idx]} ',
va='top', ha='right') va='top', ha='right')
ax.tick_params(labelbottom=False) ax.tick_params(labelbottom=False)
ax.set_ylabel(signal_name) ax.set_ylabel(signal_name)
ax.vlines(edges[pulse_idx, signal_idx, :], *ax.get_ylim(), color='red', linewidth=1) ax.vlines(edges[pulse_idx, signal_idx, :], *ax.get_ylim(), color='red', linewidth=1)
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Digitized channel spectra ### Digitized channel spectra
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
fig = plt.figure(num=f'digitized_spectra_{det_name}', figsize=(9.5, 6)) fig = plt.figure(num=f'digitized_spectra_{det_name}', figsize=(9.5, 6))
edges = det_data[det_name]['edges'] edges = det_data[det_name]['edges']
amplitudes = det_data[det_name]['amplitudes']
min_edge = np.nanmin(edges) min_edge = np.nanmin(edges)
max_edge = np.nanmax(edges) max_edge = np.nanmax(edges)
if np.isnan(min_edge):
warning(f'No edges found for {det_name}')
continue
grid = fig.add_gridspec(ncols=3, nrows=3, left=0.08, right=0.98, top=0.95, hspace=0.4) grid = fig.add_gridspec(ncols=3, nrows=3, left=0.08, right=0.98, top=0.95, hspace=0.4)
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
numx = fig.add_subplot(grid[0, 0]) numx = fig.add_subplot(grid[0, 0])
numx.set_title('Edges per pulse') numx.set_title('Edges per pulse')
agg_window = num_pulses // 60 agg_window = num_pulses // 60
max_num_edges = 0.0 max_num_edges = 0.0
max_spectral_intensity = 0 max_spectral_intensity = 0
hist_axs = [] hist_axs = []
for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']): for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
if edge_idx < 6: if edge_idx < 6:
row = 1 + edge_idx % 2 row = 1 + edge_idx % 2
col = edge_idx // 2 col = edge_idx // 2
else: else:
row = 0 row = 0
col = np.s_[1:3] col = np.s_[1:3]
ax = fig.add_subplot(grid[row, col]) ax = fig.add_subplot(grid[row, col])
ax.set_title(f'TOF spectrum: {edge_name}') ax.set_title(f'TOF spectrum: {edge_name}')
num_edges = np.isfinite(edges[:, edge_idx, :]).sum(axis=1) num_edges = np.isfinite(edges[:, edge_idx, :]).sum(axis=1)
num_edges = num_edges[:((len(num_edges) // agg_window) * agg_window)] num_edges = num_edges[:((len(num_edges) // agg_window) * agg_window)]
num_edges = num_edges.reshape(-1, agg_window).mean(axis=1) num_edges = num_edges.reshape(-1, agg_window).mean(axis=1)
if (num_edges == 0).all(): if (num_edges == 0).all():
warning(f'No edges found for {det_name}/{edge_name}') warning(f'No edges found for {det_name}/{edge_name}')
continue continue
if edge_idx < 6: if edge_idx < 6:
plot_kwargs = dict(c=f'C{edge_idx}', ls='solid', lw=1.0) plot_kwargs = dict(c=f'C{edge_idx}', ls='solid', lw=1.0)
else: else:
plot_kwargs = dict(c='k', ls='dashed', lw=1.0) plot_kwargs = dict(c='k', ls='dashed', lw=1.0)
numx.plot(np.arange(len(num_edges)) * agg_window, num_edges, label=edge_name, **plot_kwargs) numx.plot(np.arange(len(num_edges)) * agg_window, num_edges, label=edge_name, **plot_kwargs)
max_num_edges = max(max_num_edges, num_edges.max()) max_num_edges = max(max_num_edges, num_edges.max())
y, _, _ = ax.hist(finite_flattened_slice(edges, np.s_[:, edge_idx, :]), y, _, _ = ax.hist(finite_flattened_slice(edges, np.s_[:, edge_idx, :]),
bins=int((max_edge - min_edge) // 5), range=(min_edge, max_edge), bins=int((max_edge - min_edge) // 5), range=(min_edge, max_edge),
color=plot_kwargs['c'], histtype='step', linewidth=1) color=plot_kwargs['c'], histtype='step', linewidth=1)
hist_axs.append(ax) hist_axs.append(ax)
max_spectral_intensity = max(max_spectral_intensity, y.max()) max_spectral_intensity = max(max_spectral_intensity, y.max())
numx.tick_params(labelbottom=False) numx.tick_params(labelbottom=False)
numx.set_ylim(0, 1.2*max_num_edges) numx.set_ylim(0, 1.2*max_num_edges)
for ax in hist_axs: for ax in hist_axs:
ax.set_ylim(0, max_spectral_intensity*1.1) ax.set_ylim(0, max_spectral_intensity*1.1)
ax.ticklabel_format(axis='y', style='sci', scilimits=(0, 3)) ax.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Spectral pulse height distributions ### Spectral pulse height distributions
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
fig = plt.figure(num=f'spectral_pulse_heights_{det_name}', figsize=(9.5, 12.0)) fig = plt.figure(num=f'spectral_pulse_heights_{det_name}', figsize=(9.5, 12.0))
grid = fig.add_gridspec(ncols=2, nrows=4, left=0.08, right=0.98, top=0.95, hspace=0.3) grid = fig.add_gridspec(ncols=2, nrows=4, left=0.08, right=0.98, top=0.95, hspace=0.3)
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
edges = det_data[det_name]['edges'] edges = det_data[det_name]['edges']
amplitudes = det_data[det_name]['amplitudes'] amplitudes = det_data[det_name]['amplitudes']
min_edge = np.nanmin(edges) min_edge = np.nanmin(edges)
max_edge = np.nanmax(edges) max_edge = np.nanmax(edges)
if np.isnan(min_edge):
warning(f'No edges found for {det_name}')
continue
max_amplitude = np.nanmax(amplitudes) max_amplitude = np.nanmax(amplitudes)
for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']): for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
if edge_idx < 6: if edge_idx < 6:
row = 1 + edge_idx // 2 row = 1 + edge_idx // 2
col = edge_idx % 2 col = edge_idx % 2
tof_bins = int((max_edge - min_edge) // 20) tof_bins = int((max_edge - min_edge) // 20)
else: else:
row = 0 row = 0
col = np.s_[:] col = np.s_[:]
tof_bins = int((max_edge - min_edge) // 10) tof_bins = int((max_edge - min_edge) // 10)
ax = fig.add_subplot(grid[row, col]) ax = fig.add_subplot(grid[row, col])
ax.set_title(f'Spectral pulse amplitudes: {edge_name}') ax.set_title(f'Spectral pulse amplitudes: {edge_name}')
flat_edges = finite_flattened_slice(edges, np.s_[:, edge_idx, :]) flat_edges = finite_flattened_slice(edges, np.s_[:, edge_idx, :])
flat_amplitudes = finite_flattened_slice(amplitudes, np.s_[:, edge_idx, :]) flat_amplitudes = finite_flattened_slice(amplitudes, np.s_[:, edge_idx, :])
ax.hist2d(flat_edges, flat_amplitudes,
bins=[tof_bins, 512], norm=LogNorm(), if flat_edges.size > 1:
range=[[min_edge, max_edge], [0, max_amplitude]]) ax.hist2d(flat_edges, flat_amplitudes,
bins=[tof_bins, 512], norm=LogNorm(),
range=[[min_edge, max_edge], [0, max_amplitude]])
if edge_idx == 6: if edge_idx == 6:
ax.set_ylabel('Pulse height') ax.set_ylabel('Pulse height')
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
A more detailed view into the distribution of pulse heights as a function of TOF, e.g. to indicate whether the spectrometer transmission may depend on the kinetic energy and/or (in the case of ions) mass. A more detailed view into the distribution of pulse heights as a function of TOF, e.g. to indicate whether the spectrometer transmission may depend on the kinetic energy and/or (in the case of ions) mass.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Detector diagnostics # Detector diagnostics
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
edges = det_data[det_name]['edges'] edges = det_data[det_name]['edges']
sort = remi.get_dld_sorter(det_name) sort = remi.get_dld_sorter(det_name)
sum_shifts = sort.sum_shifts if sort.sum_shifts != (0.0, 0.0, 0.0) else None sum_shifts = sort.sum_shifts if sort.sum_shifts != (0.0, 0.0, 0.0) else None
is_valid = remi.get_presort_mask(edges, edge_idx=0, w=not quad_anode, is_valid = remi.get_presort_mask(edges, edge_idx=0, w=not quad_anode,
sum_limit=max(sort.uncorrected_time_sum_half_widths), sum_limit=max(sort.uncorrected_time_sum_half_widths),
sum_shifts=sum_shifts) sum_shifts=sum_shifts)
if not is_valid.any(): if not is_valid.any():
warning(f'No valid preliminary edge combinations found for {det_name}') warning(f'No preliminary edge combinations found for {det_name}')
continue
signals, sums = remi.get_signals_and_sums(edges, indices=sort.channel_indices, sum_shifts=sum_shifts, signals, sums = remi.get_signals_and_sums(edges, indices=sort.channel_indices, sum_shifts=sum_shifts,
mask=is_valid) mask=is_valid)
fig = plot_detector_diagnostics(signals=signals, sums=sums, fig = plot_detector_diagnostics(signals=signals, sums=sums,
fig_num=f'diagnostics_{det_name}', im_scale=1.5, fig_num=f'diagnostics_{det_name}', im_scale=1.5,
sum_range=max(sort.uncorrected_time_sum_half_widths), sum_range=max(sort.uncorrected_time_sum_half_widths),
sorter=sort) sorter=sort)
fig.text(0.02, 0.98, det_name.upper() + ' before corrections', rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper() + ' before corrections', rotation=90, ha='left', va='top', size='x-large')
if remi['detector'][det_name]['use_sum_correction'] or remi['detector'][det_name]['use_pos_correction']: if remi['detector'][det_name]['use_sum_correction'] or remi['detector'][det_name]['use_pos_correction']:
n_masked = is_valid.sum() n_masked = is_valid.sum()
signals = np.full((n_masked, 3), np.nan, dtype=np.float64) signals = np.full((n_masked, 3), np.nan, dtype=np.float64)
sums = np.full((n_masked, 3), np.nan, dtype=np.float64) sums = np.full((n_masked, 3), np.nan, dtype=np.float64)
sort.correct(edges[is_valid], signals, sums) sort.correct(edges[is_valid], signals, sums)
fig = plot_detector_diagnostics(signals=signals, sums=sums, fig = plot_detector_diagnostics(signals=signals, sums=sums,
fig_num=f'corr_diagnostics_{det_name}', im_scale=1.5, fig_num=f'corr_diagnostics_{det_name}', im_scale=1.5,
sum_range=max(sort.uncorrected_time_sum_half_widths), sum_range=max(sort.uncorrected_time_sum_half_widths),
sorter=sort) sorter=sort)
fig.text(0.02, 0.98, det_name.upper() + ' after corrections', rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper() + ' after corrections', rotation=90, ha='left', va='top', size='x-large')
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Overview of initial detector signal correlations before actual hit reconstruction takes place. Only the firsts edge on each channel occuring for each trigger is included, if their times are compatible with a rough time sum window. Overview of initial detector signal correlations before actual hit reconstruction takes place. Only the firsts edge on each channel occuring for each trigger is included, if their times are compatible with a rough time sum window.
* The top row contains the spectrum of time differences on each wire in temporal coordinates on the left and spatial coordinates on the right (according to configured scale factors). * The top row contains the spectrum of time differences on each wire in temporal coordinates on the left and spatial coordinates on the right (according to configured scale factors).
* The middle row depicts time sums, first integrated and then as a function of time difference. The time sum should generally be somewhat constant, a spectrum-like appearance indicates wire ends have been swapped entirely. * The middle row depicts time sums, first integrated and then as a function of time difference. The time sum should generally be somewhat constant, a spectrum-like appearance indicates wire ends have been swapped entirely.
* [HEX-only] The bottom row shows the detector image for each combination of wires based on this limited dataset. There should be no deformations or rotations in any of the wire pairs, else likely channels are misassigned. * [HEX-only] The bottom row shows the detector image for each combination of wires based on this limited dataset. There should be no deformations or rotations in any of the wire pairs, else likely channels are misassigned.
The plot occurs twice if signal-level corrections for time sum or position are enabled. The plot occurs twice if signal-level corrections for time sum or position are enabled.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Hit reconstruction # Hit reconstruction
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_rec_hits)) psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_rec_hits))
for det_name, det in remi['detector'].items(): for det_name, det in remi['detector'].items():
edges = det_data[det_name]['edges'] edges = det_data[det_name]['edges']
signals = psh.alloc(shape=(num_pulses, 50), dtype=signal_dt, fill=np.nan) signals = psh.alloc(shape=(num_pulses, 50), dtype=signal_dt, fill=np.nan)
hits = psh.alloc(shape=(num_pulses, 50), dtype=hit_dt, fill=(np.nan, np.nan, np.nan, -1)) hits = psh.alloc(shape=(num_pulses, 50), dtype=hit_dt, fill=(np.nan, np.nan, np.nan, -1))
hit_counts = psh.alloc(shape=len(dc.train_ids), dtype=np.uint32) hit_counts = psh.alloc(shape=len(dc.train_ids), dtype=np.uint32)
def prepare_hit_worker(worker_id): def prepare_hit_worker(worker_id):
sort = remi.get_dld_sorter(det_name) sort = remi.get_dld_sorter(det_name)
yield yield
@psh.with_init(prepare_hit_worker) @psh.with_init(prepare_hit_worker)
def reconstruct_hits(worker_id, index, train_id): def reconstruct_hits(worker_id, index, train_id):
hit_counts[index] += sort.run_on_train( hit_counts[index] += sort.run_on_train(
edges, signals, hits, pulse_offsets[index], pulse_offsets[index] + pulse_counts[index]) edges, signals, hits, pulse_offsets[index], pulse_offsets[index] + pulse_counts[index])
with timing(f'rec_hits, {det_name}'): with timing(f'rec_hits, {det_name}'):
psh.map(reconstruct_hits, dc.train_ids) psh.map(reconstruct_hits, dc.train_ids)
det_data[det_name].update(signals=signals, hits=hits, hit_counts=hit_counts) det_data[det_name].update(signals=signals, hits=hits, hit_counts=hit_counts)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig, ax = plt.subplots(num='hit_count_per_trigger', figsize=(9.5, 4), ncols=1, clear=True, fig, ax = plt.subplots(num='hit_count_per_trigger', figsize=(9.5, 4), ncols=1, clear=True,
gridspec_kw=dict(top=0.92, right=0.98, left=0.05, bottom=0.12)) gridspec_kw=dict(top=0.92, right=0.98, left=0.05, bottom=0.12))
max_num_hits = 0.0 max_num_hits = 0.0
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
agg_window = num_pulses // min(1000, num_pulses) agg_window = num_pulses // min(1000, num_pulses)
num_hits = np.isfinite(det_data[det_name]['hits']['x']).sum(axis=1) num_hits = np.isfinite(det_data[det_name]['hits']['x']).sum(axis=1)
num_hits = num_hits[:(len(num_hits) // agg_window) * agg_window] num_hits = num_hits[:(len(num_hits) // agg_window) * agg_window]
num_hits = num_hits.reshape(-1, agg_window).mean(axis=1) num_hits = num_hits.reshape(-1, agg_window).mean(axis=1)
max_num_hits = max(max_num_hits, num_hits.max()) max_num_hits = max(max_num_hits, num_hits.max())
if max_num_hits == 0:
warning(f'No hits found for {det_name}')
continue
ax.plot(np.arange(0, (num_pulses // agg_window) * agg_window, agg_window), num_hits, ax.plot(np.arange(0, (num_pulses // agg_window) * agg_window, agg_window), num_hits,
lw=1, label=det_name.upper()) lw=1, label=det_name.upper())
ax.set_title('Hits per pulse') ax.set_title('Hits per pulse')
ax.set_xlabel('Pulse index') ax.set_xlabel('Pulse index')
ax.set_ylim(0, max_num_hits*1.1) ax.set_ylim(0, max_num_hits*1.1)
ax.legend() ax.legend()
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Reconstruction methods ### Reconstruction methods
Each hit may be reconstructed by one of 19 different methods. These differ by the number of real signals across the channels, which could be combined to form the hit. Each of these methods is designed by a number between `0` and `19` (with empty hits using `-1`), which can be found in the `m` key of a hit, e.g.: Each hit may be reconstructed by one of 19 different methods. These differ by the number of real signals across the channels, which could be combined to form the hit. Each of these methods is designed by a number between `0` and `19` (with empty hits using `-1`), which can be found in the `m` key of a hit, e.g.:
* `0`: All six anode signals and the corresponding MCP signal were found. * `0`: All six anode signals and the corresponding MCP signal were found.
* `4`: One signal on layer `u` is missing, all other signals for this event were found. * `4`: One signal on layer `u` is missing, all other signals for this event were found.
* `18`: Only one anode signal on each layer was found and the MCP signal is missing. There is no way to check whether this combination of signals is actually valid based on the detector data alone. * `18`: Only one anode signal on each layer was found and the MCP signal is missing. There is no way to check whether this combination of signals is actually valid based on the detector data alone.
| Method | `u+v+w +mcp` | | Method | `u+v+w +mcp` |
| - | - | | - | - |
| 0 | `2+2+2 +1` | | 0 | `2+2+2 +1` |
| 1 | `0+2+2 +1` | | 1 | `0+2+2 +1` |
| 2 | `2+0+2 +1` | | 2 | `2+0+2 +1` |
| 3 | `2+2+0 +1` | | 3 | `2+2+0 +1` |
| 4 | `1+2+2 +1` (2 permutations) | | 4 | `1+2+2 +1` (2 permutations) |
| 5 | `2+1+2 +1` (2 permutations) | | 5 | `2+1+2 +1` (2 permutations) |
| 6 | `2+2+1 +1` (2 permutations) | | 6 | `2+2+1 +1` (2 permutations) |
| 7 | `2+2+2 +0` | | 7 | `2+2+2 +0` |
| 8 | `0+2+2 +0` | | 8 | `0+2+2 +0` |
| 9 | `2+0+2 +0` | | 9 | `2+0+2 +0` |
| 10 | `2+2+0 +0` | | 10 | `2+2+0 +0` |
| 11 | `1+2+2 +0` (2 permutations) | | 11 | `1+2+2 +0` (2 permutations) |
| 12 | `2+1+2 +0` (2 permutations) | | 12 | `2+1+2 +0` (2 permutations) |
| 13 | `2+2+1 +0` (2 permutations) | | 13 | `2+2+1 +0` (2 permutations) |
| 14 | `2+1+1 +1` `1+2+1 +1` `1+1+2 +1` (12 permutations) | | 14 | `2+1+1 +1` `1+2+1 +1` `1+1+2 +1` (12 permutations) |
| 15 | `2+1+0 +1` `2+0+1 +1` `1+2+0 +1` `1+0+2 +1` `0+2+1 +1` `0+1+2 +1` (12 permutations) | | 15 | `2+1+0 +1` `2+0+1 +1` `1+2+0 +1` `1+0+2 +1` `0+2+1 +1` `0+1+2 +1` (12 permutations) |
| 16 | `1+1+1 +1` (8 permutations) | | 16 | `1+1+1 +1` (8 permutations) |
| 17 | `2+1+1 +0` `1+2+1 +0` `1+1+2 +0` (12 permutations) | | 17 | `2+1+1 +0` `1+2+1 +0` `1+1+2 +0` (12 permutations) |
| 18 | `1+1+1 +0` (8 permutations) | | 18 | `1+1+1 +0` (8 permutations) |
| 19 | `2+1+0 +0` `2+0+1 +0` `1+2+0 +0` `1+0+2 +0` `0+1+2 +0` `0+2+1 +0` (12 permutations) | | 19 | `2+1+0 +0` `2+0+1 +0` `1+2+0 +0` `1+0+2 +0` `0+1+2 +0` `0+2+1 +0` (12 permutations) |
* For hits reconstructed with method `> 10`, extra attention should be given to ensure they add meaningful signal. * For hits reconstructed with method `> 10`, extra attention should be given to ensure they add meaningful signal.
* Any method `> 14` has to considered risky, because neither a time sum nor the position can be checked. If the scale factors and/or `w` shift are not correct, then the number of events reconstructed with the risky methods will increase. They will most likely be *ghost hits*, which do not correspond to actual impacts on the detector. * Any method `> 14` has to considered risky, because neither a time sum nor the position can be checked. If the scale factors and/or `w` shift are not correct, then the number of events reconstructed with the risky methods will increase. They will most likely be *ghost hits*, which do not correspond to actual impacts on the detector.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
hits = det_data[det_name]['hits'] hits = det_data[det_name]['hits']
if not (hits['m'] >= 0).any():
warning(f'No hits found for {det_name}')
continue
fig, ax = plt.subplots(num=f'reconstruction_methods_{det_name}', figsize=(9.5, 5), ncols=1, clear=True, fig, ax = plt.subplots(num=f'reconstruction_methods_{det_name}', figsize=(9.5, 5), ncols=1, clear=True,
gridspec_kw=dict(left=0.08, right=0.91, top=0.8)) gridspec_kw=dict(left=0.08, right=0.91, top=0.8))
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
if not (hits['m'] >= 0).any():
warning(f'No hits found for {det_name}')
continue
method_bins = np.bincount(hits['m'][hits['m'] >= 0], minlength=20) method_bins = np.bincount(hits['m'][hits['m'] >= 0], minlength=20)
ax.bar(np.arange(20), method_bins, width=0.5) ax.bar(np.arange(20), method_bins, width=0.5)
ax.set_xlabel('Reconstruction method') ax.set_xlabel('Reconstruction method')
ax.set_xlim(-0.5, 19.5) ax.set_xlim(-0.5, 19.5)
ax.set_xticks(np.arange(20)) ax.set_xticks(np.arange(20))
ax.set_ylabel('Number of hits') ax.set_ylabel('Number of hits')
ax.set_ylim(0, method_bins.max()*1.05) ax.set_ylim(0, method_bins.max()*1.05)
ylims = ax.get_ylim() ylims = ax.get_ylim()
ax.tick_params(which='both', right=True, labelright=True) ax.tick_params(which='both', right=True, labelright=True)
num_risky = method_bins[15:].sum() num_risky = method_bins[15:].sum()
num_total = method_bins.sum() num_total = method_bins.sum()
ax.text(14.2, method_bins.max(), f'{(100*(num_total-num_risky)/num_total):.2g}%', ax.text(14.2, method_bins.max(), f'{(100*(num_total-num_risky)/num_total):.2g}%',
va='top', ha='right', color='black') va='top', ha='right', color='black')
ax.text(14.8, method_bins.max(), f'{(100*num_risky/num_total):.2g}%', ax.text(14.8, method_bins.max(), f'{(100*num_risky/num_total):.2g}%',
va='top', ha='left', color='red') va='top', ha='left', color='red')
ax.fill([14.5, 19.5, 19.5, 14.5], [ylims[0], ylims[0], ylims[1], ylims[1]], c='r', alpha=0.2) ax.fill([14.5, 19.5, 19.5, 14.5], [ylims[0], ylims[0], ylims[1], ylims[1]], c='r', alpha=0.2)
labelx = ax.twiny() labelx = ax.twiny()
labelx.set_xlim(*ax.get_xlim()) labelx.set_xlim(*ax.get_xlim())
labelx.set_xticks(ax.get_xticks()) labelx.set_xticks(ax.get_xticks())
labelx.set_xticklabels([ labelx.set_xticklabels([
'2+2+2 +1', '2+2+2 +1',
'0+2+2 +1', '2+0+2 +1', '2+2+0 +1', '0+2+2 +1', '2+0+2 +1', '2+2+0 +1',
'1+2+2 +1', '2+1+2 +1', '2+2+1 +1', '1+2+2 +1', '2+1+2 +1', '2+2+1 +1',
'2+2+2 +0', '2+2+2 +0',
'0+2+2 +0', '2+0+2 +0', '2+2+0 +0', '1+2+2 +0', '2+1+2 +0', '2+2+1 +0', '0+2+2 +0', '2+0+2 +0', '2+2+0 +0', '1+2+2 +0', '2+1+2 +0', '2+2+1 +0',
'2+1+1 +1', '2+1+1 +1',
'2+1+0 +1', '2+1+0 +1',
'1+1+1 +1', '1+1+1 +1',
'2+1+1 +0', '2+1+1 +0',
'1+1+1 +0', '1+1+1 +0',
'2+1+0 +0', '2+1+0 +0',
], rotation=90) ], rotation=90)
min_rel_tick = np.ceil((ax.get_ylim()[0] / num_total) / 0.1) * 0.1 min_rel_tick = np.ceil((ax.get_ylim()[0] / num_total) / 0.1) * 0.1
max_rel_tick = np.floor((method_bins.max() / num_total) / 0.1) * 0.1 max_rel_tick = np.floor((method_bins.max() / num_total) / 0.1) * 0.1
rely = ax.twinx() rely = ax.twinx()
rely.set_ylim(*ax.get_ylim()) rely.set_ylim(*ax.get_ylim())
rely.set_yticks(np.arange(0.0, max_rel_tick+0.01, 0.1)*num_total) rely.set_yticks(np.arange(0.0, max_rel_tick+0.01, 0.1)*num_total)
rely.set_yticks(np.arange(0.0, ylims[1]/num_total, 0.02)*num_total, minor=True) rely.set_yticks(np.arange(0.0, ylims[1]/num_total, 0.02)*num_total, minor=True)
rely.set_yticklabels([f'{(y/num_total)*100:.0f}%' for y in rely.get_yticks()]) rely.set_yticklabels([f'{(y/num_total)*100:.0f}%' for y in rely.get_yticks()])
rely.set_ylabel('Percentage of total hits') rely.set_ylabel('Percentage of total hits')
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Detector image and fishes ### Detector image and fishes
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
for det_name in remi['detector'].keys(): for det_name in remi['detector'].keys():
flat_hits = det_data[det_name]['hits'].reshape(-1) flat_hits = det_data[det_name]['hits'].reshape(-1)
flat_hits = flat_hits[np.isfinite(flat_hits[:]['x'])] flat_hits = flat_hits[np.isfinite(flat_hits[:]['x'])]
flat_hits = flat_hits[flat_hits['m'] <= 10] flat_hits = flat_hits[flat_hits['m'] <= 10]
if flat_hits.size == 0:
warning(f'No hits found for {det_name}')
continue
fig = plt.figure(num=f'detector_results_{det_name}', figsize=(9, 10.5)) fig = plt.figure(num=f'detector_results_{det_name}', figsize=(9, 10.5))
fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large') fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')
fig.text(0.02, 0.02, det_name.upper(), rotation=90, ha='left', va='bottom', size='x-large') fig.text(0.02, 0.02, det_name.upper(), rotation=90, ha='left', va='bottom', size='x-large')
imp = fig.add_axes([0.1 + 0.25/2, 0.56, 0.5, 0.45]) imp = fig.add_axes([0.1 + 0.25/2, 0.56, 0.5, 0.45])
txp = fig.add_axes([0.1, 0.27, 0.85, 0.23]) txp = fig.add_axes([0.1, 0.27, 0.85, 0.23])
typ = fig.add_axes([0.1, 0.02, 0.85, 0.23]) typ = fig.add_axes([0.1, 0.02, 0.85, 0.23])
if flat_hits.size == 0:
warning(f'No hits found for {det_name}')
continue
mcp_radius = remi['detector'][det_name]['mcp_radius'] mcp_radius = remi['detector'][det_name]['mcp_radius']
im_radius = mcp_radius * 1.1 im_radius = mcp_radius * 1.1
imp.hist2d(flat_hits['x'], flat_hits['y'], bins=(256, 256), imp.hist2d(flat_hits['x'], flat_hits['y'], bins=(256, 256),
range=[[-im_radius, im_radius], [-im_radius, im_radius]], norm=LogNorm()) range=[[-im_radius, im_radius], [-im_radius, im_radius]], norm=LogNorm())
imp.add_patch(Circle( imp.add_patch(Circle(
(0, 0), mcp_radius, (0, 0), mcp_radius,
linestyle='dashed', edgecolor='red', facecolor='none', linewidth=1)) linestyle='dashed', edgecolor='red', facecolor='none', linewidth=1))
imp.xaxis.set_label_position('top') imp.xaxis.set_label_position('top')
imp.set_xlabel('X / mm') imp.set_xlabel('X / mm')
imp.set_ylabel('Y / mm') imp.set_ylabel('Y / mm')
imp.tick_params(right=True, labelright=True, top=True, labeltop=True) imp.tick_params(right=True, labelright=True, top=True, labeltop=True)
imp.grid() imp.grid()
text_pos = 1.05*mcp_radius*np.sin(np.pi/4) text_pos = 1.05*mcp_radius*np.sin(np.pi/4)
imp.text(text_pos, text_pos, 'MCP', c='red', ha='left', va='bottom') imp.text(text_pos, text_pos, 'MCP', c='red', ha='left', va='bottom')
min_tof = flat_hits['t'].min() min_tof = flat_hits['t'].min()
max_tof = flat_hits['t'].max() max_tof = flat_hits['t'].max()
num_tof_bins = min(int((max_tof - min_tof) // 10), 500) num_tof_bins = min(int((max_tof - min_tof) // 10), 500)
if num_tof_bins == 0: if num_tof_bins == 0:
warning(f'All TOFs limited to single bin for {det_name}') warning(f'All TOFs limited to single bin for {det_name}')
continue continue
for ax, dim_label in zip([txp, typ], ['x', 'y']): for ax, dim_label in zip([txp, typ], ['x', 'y']):
ax.hist2d(flat_hits['t'], flat_hits[dim_label], bins=(num_tof_bins, 256), ax.hist2d(flat_hits['t'], flat_hits[dim_label], bins=(num_tof_bins, 256),
range=[[min_tof, max_tof], [-im_radius, im_radius]], norm=LogNorm()) range=[[min_tof, max_tof], [-im_radius, im_radius]], norm=LogNorm())
ax.set_ylabel(f'{dim_label.upper()} / mm') ax.set_ylabel(f'{dim_label.upper()} / mm')
typ.set_xlabel('Time-of-flight / ns') typ.set_xlabel('Time-of-flight / ns')
txp.tick_params(bottom=True, labelbottom=False, top=True, labeltop=True, right=True, labelright=True) txp.tick_params(bottom=True, labelbottom=False, top=True, labeltop=True, right=True, labelright=True)
typ.tick_params(right=True, labelright=True, top=True) typ.tick_params(right=True, labelright=True, top=True)
pass pass
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Transformed data files # Transformed data files
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Try to figure out proposal number from in_folder to work with older files. # Try to figure out proposal number from in_folder to work with older files.
m = re.match(r'p(\d{6})', Path(in_folder).parts[-2]) m = re.match(r'p(\d{6})', Path(in_folder).parts[-2])
if not proposal and m is not None: if not proposal and m is not None:
proposal = int(m[1]) proposal = int(m[1])
seq_len = out_seq_len if out_seq_len > 0 else len(dc.files[0].train_ids) seq_len = out_seq_len if out_seq_len > 0 else len(dc.files[0].train_ids)
dataset_kwargs = {k[8:]: v for k, v in locals().items() if k.startswith('dataset_compression')} dataset_kwargs = {k[8:]: v for k, v in locals().items() if k.startswith('dataset_compression')}
control_sources = [det_device_id.format(karabo_id=karabo_id, det_name=det_name.upper()) control_sources = [det_device_id.format(karabo_id=karabo_id, det_name=det_name.upper())
for det_name in remi['detector']] for det_name in remi['detector']]
instrument_channels = [ instrument_channels = [
f'{device_id}:{det_output_key}/{channel}' f'{device_id}:{det_output_key}/{channel}'
for device_id in control_sources for device_id in control_sources
for channel in ['raw', 'rec'] for channel in ['raw', 'rec']
] ]
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
Path(out_folder).mkdir(parents=True, exist_ok=True) Path(out_folder).mkdir(parents=True, exist_ok=True)
print('Writing sequence files', flush=True, end='') print('Writing sequence files', flush=True, end='')
t_write = timing('write_files') t_write = timing('write_files')
t_write.__enter__() t_write.__enter__()
for seq_id, train_mask, pulse_mask in sequence_pulses(dc.train_ids, pulse_counts, pulse_offsets, seq_len): for seq_id, train_mask, pulse_mask in sequence_pulses(dc.train_ids, pulse_counts, pulse_offsets, seq_len):
seq_train_ids = dc.train_ids[train_mask] seq_train_ids = dc.train_ids[train_mask]
with DataFile.from_details(out_folder, out_aggregator, run, seq_id) as outp: with DataFile.from_details(out_folder, out_aggregator, run, seq_id) as outp:
outp.create_metadata(like=dc, proposal=proposal, run=run, sequence=seq_id, outp.create_metadata(like=dc, proposal=proposal, run=run, sequence=seq_id,
control_sources=control_sources, instrument_channels=instrument_channels) control_sources=control_sources, instrument_channels=instrument_channels)
outp.create_index( outp.create_index(
seq_train_ids, seq_train_ids,
timestamps=dc.select_trains(by_id[seq_train_ids]).train_timestamps().astype(np.uint64) timestamps=dc.select_trains(by_id[seq_train_ids]).train_timestamps().astype(np.uint64)
) )
for det_name in remi['detector']: for det_name in remi['detector']:
cur_device_id = det_device_id.format(karabo_id=karabo_id, det_name=det_name.upper()) cur_device_id = det_device_id.format(karabo_id=karabo_id, det_name=det_name.upper())
cur_max_hits = remi['detector'][det_name]['max_hits'] cur_max_hits = remi['detector'][det_name]['max_hits']
cur_control_data = outp.create_control_source(cur_device_id) cur_control_data = outp.create_control_source(cur_device_id)
# Manually manipulate the file here, still creates the index properly. # Manually manipulate the file here, still creates the index properly.
remi.attach_detector_config(det_name, cur_control_data.get_run_group()) remi.attach_detector_config(det_name, cur_control_data.get_run_group())
cur_control_data.create_index(len(seq_train_ids)) cur_control_data.create_index(len(seq_train_ids))
cur_fast_data = outp.create_instrument_source(f'{cur_device_id}:{det_output_key}') cur_fast_data = outp.create_instrument_source(f'{cur_device_id}:{det_output_key}')
cur_data = det_data[det_name] cur_data = det_data[det_name]
cur_fast_data.create_key('raw.triggers', triggers[pulse_mask], cur_fast_data.create_key('raw.triggers', triggers[pulse_mask],
maxshape=(None,) + triggers.shape[1:], maxshape=(None,) + triggers.shape[1:],
chunks=tuple(chunks_triggers), **dataset_kwargs) chunks=tuple(chunks_triggers), **dataset_kwargs)
cur_fast_data.create_key('raw.edges', cur_data['edges'][pulse_mask], cur_fast_data.create_key('raw.edges', cur_data['edges'][pulse_mask],
maxshape=(None,) + cur_data['edges'].shape[1:], maxshape=(None,) + cur_data['edges'].shape[1:],
chunks=tuple(chunks_edges if chunks_edges[-1] <= cur_max_hits chunks=tuple(chunks_edges if chunks_edges[-1] <= cur_max_hits
else chunks_edges[:-1] + [cur_max_hits]), else chunks_edges[:-1] + [cur_max_hits]),
**dataset_kwargs) **dataset_kwargs)
cur_fast_data.create_key('raw.amplitudes', cur_data['amplitudes'][pulse_mask], cur_fast_data.create_key('raw.amplitudes', cur_data['amplitudes'][pulse_mask],
maxshape=(None,) + cur_data['amplitudes'].shape[1:], maxshape=(None,) + cur_data['amplitudes'].shape[1:],
chunks=tuple(chunks_amplitudes if chunks_amplitudes[-1] <= cur_max_hits chunks=tuple(chunks_amplitudes if chunks_amplitudes[-1] <= cur_max_hits
else chunks_amplitudes[:-1] + [cur_max_hits]), else chunks_amplitudes[:-1] + [cur_max_hits]),
**dataset_kwargs) **dataset_kwargs)
cur_fast_data.create_key('rec.signals', cur_data['signals'][pulse_mask], cur_fast_data.create_key('rec.signals', cur_data['signals'][pulse_mask],
maxshape=(None,) + cur_data['signals'].shape[1:], maxshape=(None,) + cur_data['signals'].shape[1:],
chunks=tuple(chunks_signals if chunks_signals[-1] <= cur_max_hits chunks=tuple(chunks_signals if chunks_signals[-1] <= cur_max_hits
else chunks_signals[:-1] + [cur_max_hits]), else chunks_signals[:-1] + [cur_max_hits]),
**dataset_kwargs) **dataset_kwargs)
cur_fast_data.create_key('rec.hits', cur_data['hits'][pulse_mask], cur_fast_data.create_key('rec.hits', cur_data['hits'][pulse_mask],
maxshape=(None,) + hits.shape[1:], maxshape=(None,) + hits.shape[1:],
chunks=tuple(chunks_hits if chunks_hits[-1] <= cur_max_hits chunks=tuple(chunks_hits if chunks_hits[-1] <= cur_max_hits
else chunks_hits[:-1] + [cur_max_hits]), else chunks_hits[:-1] + [cur_max_hits]),
**dataset_kwargs) **dataset_kwargs)
cur_fast_data.create_index(raw=pulse_counts[train_mask], rec=pulse_counts[train_mask]) cur_fast_data.create_index(raw=pulse_counts[train_mask], rec=pulse_counts[train_mask])
print('.', flush=True, end='') print('.', flush=True, end='')
print('') print('')
t_write.__exit__() t_write.__exit__()
``` ```
......
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

Legal Notice | Contacts