Improve automatic determination of TOF bins in DLD report plots

6e26a244 · Philipp Schmidt · 1bb17a9c · 6e26a244
Commit 6e26a244 authored 1 year ago by Philipp Schmidt
--- a/notebooks/REMI/REMI_Digitize_and_Transform.ipynb
+++ b/notebooks/REMI/REMI_Digitize_and_Transform.ipynb
@@ -1074,7 +1074,7 @@
    "    min_tof = flat_hits['t'].min()\n",
    "    max_tof = flat_hits['t'].max()\n",
    "    \n",
-    "    num_tof_bins = int((max_tof - min_tof) // 5)\n",
+    "    num_tof_bins = min(int((max_tof - min_tof) // 10), 500)\n",
    "    \n",
    "    if num_tof_bins == 0:\n",
    "        warning(f'All TOFs limited to single bin for {det_name}')\n",

 %% Cell type:code id: tags:

 ``` python
 # Data selection parameters.
 run = 104  # 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.

 calib_config_path = '/gpfs/exfel/exp/SQS/202101/p002535/usr/config_board2+4.yaml'  # Path to correction and transform configuration

 # These parameters are required by xfel-calibrate but ignored in this notebook.
 cycle = ''  # Proposal cycle, currently not used.
 cal_db_timeout = 0  # Calibration DB timeout, currently not used.
 cal_db_interface = 'foo'  # Calibration DB interface, currently not used.
 karabo_da = 'bar'  # Karabo data aggregator name, currently not used

 # Output parameters.
 karabo_id = 'SQS_REMI_DLD6'  # Karabo device ID root for virtual output device.
 proposal = ''  # Proposal, leave empty for auto detection based on in_folder
 out_aggregator = 'REMI01'  # Aggregator name for output files.
 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_output_key = 'output'  # Pipeline name for fast data output.
 save_raw_triggers = True  # Whether to save trigger position in files.
 save_raw_edges = True  # Whether to save digitized edge positions in files.
 save_raw_amplitudes = True  # Whether to save analog pulse amplitudes in files.
 save_rec_signals = True  # Whether to save reconstructed signals (u1-w2, mcp) in files.
 save_rec_hits = True  # Whether to save reoncstructed hits (x,y,t,m) in files.
 chunks_triggers = [500]  # HDF chunk size for triggers.
 chunks_edges = [500, 7, 50]  # HDF chunk size for edges.
 chunks_amplitudes = [500, 7, 50]  # HDF chunk size for amplitudes.
 chunks_hits = [50, 50]  # HDF chunk size for hits.
 chunks_signals = [50, 50]  # HDF chunk size for signals.
 dataset_compression = 'gzip'  # HDF compression method.
 dataset_compression_opts = 3  # HDF GZIP compression level.

 # Detector parameters.
 quad_anode = False  # Reconstruction assumes a hex anode by default, change for quad anodes.

 # Trigger parameters.
 ppt_source = 'SQS_RR_UTC/TSYS/TIMESERVER:outputBunchPattern'
 ignore_fel = False  # Ignore any FEL entries in the PPT.
 ignore_ppl = False  # Ignore any PPL entries in the PPT.
 trailing_trigger = True  # Add a trigger after all regular pulses with the remaining trace.
 ppl_offset = 0  # In units of the PPT.
 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.
 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.
 pulse_start_offset = 0  # Signal offset at the start of each pulse.
 pulse_end_offset = 0  # Signal offset at the end of each pulse.

 # PPT reconstruction parameters.
 reconstruct_ppt = False  # Reconstruct PPT from some trigger edges.
 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.
 fake_ppt_offset = 0  # Offset in reconstructed PPT for pulses.

 # Parallelization parameters.
 mp_find_triggers = 0.5  # Parallelization for finding triggers.
 mp_find_edges = 0.5  # Parallelization for digitizing analog signal.
 mt_avg_trace = 2  # Parallelization for trace averaging.
 mp_rec_hits = 1.0  # Parallelization for hit reconstruction.
 ```

 %% Cell type:code id: tags:

 ``` python
 from datetime import datetime
 from logging import warning
 from pathlib import Path
 import re

 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.colors import LogNorm
 from threadpoolctl import threadpool_limits

 import h5py

 import pasha as psh
 from euxfel_bunch_pattern import indices_at_sase, indices_at_laser
 from extra_data import RunDirectory, by_id
 from extra_remi import Analysis, trigger_dt
 from extra_remi.util import timing
 from extra_remi.rd_resort import signal_dt, hit_dt
 from extra_remi.files import DataFile, sequence_pulses

 if quad_anode:
    from extra_remi.plots import plot_detector_diagnostics_quad as plot_detector_diagnostics
 else:
    from extra_remi.plots import plot_detector_diagnostics_hex as plot_detector_diagnostics

 %matplotlib inline
 ```

 %% Cell type:code id: tags:

 ``` python
 def finite_flattened_slice(array, slice_=np.s_[:]):
    """Return flattened and finite values for a given slice."""
    sliced_array = array[slice_]
    return sliced_array[np.isfinite(sliced_array)]
 ```

 %% Cell type:code id: tags:

 ``` python
 calib_config_path = Path(calib_config_path)

 if not calib_config_path.is_file():
    # 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

    if not calib_config_path.is_file():
        # Disallow implicit config file creation.
        raise ValueError('calib_config_path not found - neither absolute nor relative to in_folder')

 remi = Analysis(calib_config_path, use_hex=not quad_anode)

 with timing('open_run'):
    dc = remi.prepare_dc(RunDirectory(Path(in_folder) / f'r{run:04d}', inc_suspect_trains=True),
                         require_ppt=not reconstruct_ppt)
 ```

 %% Cell type:markdown id: tags:

 # Transformation parameters

 %% Cell type:code id: tags:

 ``` python
 def print_leaf(leaf, indent=0):
    for key, value in leaf.items():
        if isinstance(value, dict):
            print(indent * 4 * ' ' + key)
            print_leaf(value, indent=indent+1)
        else:
            print(indent * 4 * ' ' + f'{key}: {value}')

 print_leaf(remi.tree)
 ```

 %% Cell type:markdown id: tags:

 # Pulse and trigger information

 %% Cell type:markdown id: tags:

 ### Read PPT from file or reconstruct PPT for older data

 %% Cell type:code id: tags:

 ``` python
 if reconstruct_ppt:
    # Take up to the first hundred trains 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)] \
        [:100].ndarray().mean(axis=0).astype(np.float64)
    trigger_trace -= trigger_trace[0]  # Use simple offset correction.

    fake_ppt = np.zeros(2700, dtype=np.uint32)

    discr_func, discr_params = remi.get_discriminator([trigger_edge_channel])

    edges = np.zeros(1000, dtype=np.float64)
    num_pulses = discr_func(trigger_trace, edges=edges, **discr_params[0])
    edges = edges[:num_pulses]

    first_edge = edges[0]
    rel_edges = np.round(edges - first_edge)
    edge_diff = rel_edges[1] - rel_edges[0]

    if not np.allclose(rel_edges[1:] - rel_edges[:-1], edge_diff):
        raise ValueError('PPT reconstruction for unstable edge intervals not supported')

    pulse_spacing = edge_diff / (2 * remi['digitizer']['clock_factor'])  # In units of PPT

    if not float.is_integer(pulse_spacing):
        raise ValueError('PPT reconstruction encountered non-integer pulse spacing')

    pulse_spacing = int(pulse_spacing)

    # Taken from euxfel_bunch_pattern/__init__.py
    from euxfel_bunch_pattern import DESTINATION_T4D, DESTINATION_T5D, PHOTON_LINE_DEFLECTION
    if instrument_sase == 1:
        flag = DESTINATION_T4D
    elif instrument_sase == 2:
        flag = DESTINATION_T5D
    elif instrument_sase == 3:
        flag = DESTINATION_T4D | PHOTON_LINE_DEFLECTION

    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

    from pasha.functor import Functor, gen_split_slices
    class FakeKeyDataFunctor(Functor):
        """Functor appearing KeyData-like with constant data.

        This functor serves a constant data row for a given number
        of train IDs the same way a KeyData object would.
        """

        def __init__(self, row, train_ids):
            self.row = row
            self.train_ids = train_ids

        def split(self, num_workers):
            return gen_split_slices(len(self.train_ids), n_parts=num_workers)

        def iterate(self, share):
            it = zip(range(*share.indices(len(self.train_ids))), self.train_ids)

            for index, train_id in it:
                yield index, train_id, self.row

    ppt_data = FakeKeyDataFunctor(fake_ppt, dc.train_ids)

    fig, ax = plt.subplots(num=99, figsize=(9, 6), clear=True, ncols=1, nrows=1)

    ax.set_title('Edge trigger signal')
    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,
              color='red', linewidth=3, alpha=0.3, label='Edge positions')

    ax.set_xlabel('Samples')
    ax.set_ylabel('Intensity / ADU')
    ax.legend()

 else:
    ppt_data = dc[ppt_source, 'data.bunchPatternTable']
 ```

 %% Cell type:markdown id: tags:

 ### Count pulses per train and compute offsets

 %% Cell type:code id: tags:

 ``` python
 # 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_offsets [int32: len(dc.train_ids)]` containing the global offset for the first pulse of each train.
 # * `num_pulses = pulse_counts.sum(axis=0)`

 def get_pulse_positions(ppt, sase, laser, ppl_offset):
    # Combine FEL and PPL positions.

    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([])

    if len(fel_pos) > 0:
        # Move PPL up to the FEL position.
        ppl_pos += fel_pos[0] + ppl_offset

    return np.union1d(fel_pos, ppl_pos), fel_pos, ppl_pos

 if laser_ppt_mask < 0:
    # If laser PPT mask is not specified, try to figure it out from device IDs.
    from euxfel_bunch_pattern import PPL_BITS

    instrument = karabo_id[:karabo_id.index('_')]

    try:
        laser_ppt_mask = PPL_BITS[f'LP_{instrument}']
    except KeyError:
        raise ValueError(f'Laser PPT mask unknown for instrument `{instrument}`')

 with timing('pulse_info'):
    psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_triggers))

    # Build the pulse index
    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)

    def count_pulses(wid, index, tid, ppt):
        pulse_counts[index] = len(get_pulse_positions(ppt, instrument_sase, laser_ppt_mask, ppl_offset)[0])
        has_ppt[index] = True

    psh.map(count_pulses, ppt_data)

    # Fill any missing values with the highest.
    pulse_counts[has_ppt == False] = pulse_counts.max()

    if trailing_trigger:
        # Add a single count to every train for trailing trigger.
        pulse_counts += 1

    # Compute offsets based on pulse counts.
    pulse_offsets = np.zeros_like(pulse_counts)
    np.cumsum(pulse_counts[:-1], out=pulse_offsets[1:])

    # Total number of pulses.
    num_pulses = int(pulse_counts.sum())
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, ax = plt.subplots(num=1, ncols=1, nrows=1, figsize=(9, 4), clear=True)

 ax.set_title('Pulse count')
 ax.plot(dc.train_ids, pulse_counts, lw=1)
 ax.set_xlabel('Train ID')
 ax.set_ylabel('Number of pulses')
 ax.set_ylim(0, max(300, pulse_counts.max() + 10))
 ax.ticklabel_format(style='plain')
 pass
 ```

 %% Cell type:markdown id: tags:

 ### Find triggers

 The trigger defines the boundary of a pulse on the digitizer trace, which is stored per train.

 %% Cell type:code id: tags:

 ``` python
 # A trigger defines the boundary of a pulse on the digitizer trace stored per train. This cell creates a
 # global variable:
 # * `triggers [(start: int32, stop: int32, offset: float64, fel: bool, ppl: bool): num_pulses]`
 #   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
 # 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
 # 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))

 triggers = psh.alloc(shape=(num_pulses,), dtype=trigger_dt, fill=(-1, -1, np.nan, -1, 0, 0))

 clock_factor = remi['digitizer']['clock_factor']

 min_trace_len = min([
    dc[src, key].entry_shape[0] for det_name in remi['detector'].keys()
    for src, key in remi.get_detector_sourcekeys(det_name)
 ])

 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)
    num_pulses = len(all_pos)

    if num_pulses > 0:
        if len(ppl_pos) == 0 and ppl_offset < 0:
            # 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
            # end, so we correct by adding the ppl_offset to relative positions
            # when computing trace positions.
            pos_corr = abs(ppl_offset)
        else:
            pos_corr = 0

        rel_pos = all_pos - all_pos[0]

        if num_pulses > 1:
            pulse_len = np.unique(rel_pos[1:] - rel_pos[:-1]).min()
        elif num_pulses == 1:
            pulse_len = single_pulse_length

        start_frac = first_pulse_offset + (rel_pos + pos_corr) * 2 * clock_factor
        start_int = start_frac.astype(int)

        train_triggers = triggers[pulse_offsets[index]:int(pulse_offsets[index]+num_pulses)]
        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['offset'] = start_frac - start_int
        train_triggers['pulse'] = all_pos.astype(np.int16)
        train_triggers['fel'] = [pos in fel_pos for pos in all_pos]
        train_triggers['ppl'] = [pos in ppl_pos for pos in all_pos]

        last_sample = train_triggers['stop'].max()

    else:
        last_sample = first_pulse_offset

    if trailing_trigger:
        # Add trailing trigger if required.
        trigger = triggers[int(pulse_offsets[index]+pulse_counts[index]-1)]
        trigger['start'] = last_sample
        trigger['stop'] = min_trace_len
        trigger['offset'] = 0.0
        trigger['pulse'] = -1
        trigger['fel'] = False
        trigger['ppl'] = False

 with timing('find_triggers'):
    psh.map(trigger_by_ppt, ppt_data)

 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
    # necessarily a problem, as the pattern could simply have changed in between trains
    # with each train being split properly.
    # 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.

    one = np.uint64(1)  # Because np.uint64 + int = np.float64
    pulse_deltas = set()

    for pulse_id, (offset, count) in enumerate(zip(
        pulse_offsets, pulse_counts - one if trailing_trigger else pulse_counts
    )):
        deltas = triggers['pulse'][offset+one:offset+count] - triggers['pulse'][offset:offset+count-one]

        if len(np.unique(deltas)) > 1:
            for delta in deltas:
                pulse_deltas.add(delta)

    if len(pulse_deltas) > 1:
        delta_str = ', '.join([str(x) for x in sorted(pulse_deltas)])
        warning(f'Different pulse lengths (PPT: {delta_str}) encountered within single trains, '
                f'separated pulse spectra may split up signals!')
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, (lx, rx) = plt.subplots(num=2, ncols=2, nrows=1, figsize=(9, 4), clear=True,
                             gridspec_kw=dict(top=0.75))

 # Display ~400 pulses or 10 trains, whatever is lower
 n_trains = max(abs(pulse_offsets - 200).argmin(), 5)

 visible_triggers = triggers[:pulse_offsets[n_trains]]

 pulse_index = np.arange(len(visible_triggers))
 pumped = visible_triggers['fel'] & visible_triggers['ppl']
 fel_only = visible_triggers['fel'] & ~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[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')

 max_start = visible_triggers['start'].max()

 lx.vlines(pulse_offsets[:n_trains], 0, max_start, color='grey', linewidth=1, alpha=0.2)
 lx.tick_params(right=True)

 lx.set_xlabel('Pulse index')
 lx.set_xlim(-15, pulse_offsets[n_trains]+15)

 lx.set_ylabel('Trigger position')
 lx.set_ylim(-max_start // 20, max_start + max_start // 20)

 lx.legend(fontsize='small', loc='lower right')

 train_lx = lx.twiny()
 train_lx.set_xlabel('Train ID', labelpad=8)
 train_lx.set_xlim(lx.get_xlim())
 train_lx.set_xticks(pulse_offsets[:n_trains])
 train_lx.set_xticklabels([str(int(x)) for x in dc.train_ids[:n_trains]],
                         rotation=-45, fontsize='x-small')

 rx.plot(triggers['start'], lw=0.2)

 rx.set_xlabel('Pulse index')
 rx.tick_params(left=False, labelleft=False, right=True, labelright=True)

 pass
 ```

 %% Cell type:markdown id: tags:

 # Analog signal to digital edges

 %% Cell type:markdown id: tags:

 ### Find edges in analog signal

 %% Cell type:code id: tags:

 ``` python
 psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_find_edges))
 threadpool_limits(limits=remi.get_num_workers(mt_avg_trace))

 det_data = {}

 for i, (det_name, det) in enumerate(remi['detector'].items()):
    det_sourcekeys = remi.get_detector_sourcekeys(det_name)
    det_get_traces = remi.get_traces_getter(det_name)
    trace_len = dc[next(iter(det_sourcekeys))].entry_shape[0]

    edges = psh.alloc(shape=(num_pulses, 7, det['max_hits']),
                      dtype=np.float64, fill=np.nan)
    amplitudes = psh.alloc(shape=(num_pulses, 7, det['max_hits']),
                           dtype=np.float64, fill=np.nan)
    avg_traces = psh.alloc_per_worker(shape=(7, trace_len), dtype=np.float64)

    def prepare_edge_worker(worker_id):
        correct_func = remi.get_baseline_corrector()
        discr_func, discr_params = remi.get_discriminator(det['channels'])

        source_name = remi['digitizer']['source']
        bl_start, bl_stop, _ = remi.get_baseline_limits(trace_len)
        bl_sym = remi['digitizer']['baseline_symmetry']
        time_cal = remi.get_time_calibration()

        traces_corr = np.empty((7, trace_len), dtype=np.float64)
        baselines = np.empty(bl_sym, dtype=np.float64)
        yield

    @psh.with_init(prepare_edge_worker)
    def find_edges(worker_id, index, train_id, data):
        try:
            data = det_get_traces(data[source_name])
        except KeyError:
            return

        for channel_idx in range(7):
            correct_func(data[channel_idx], traces_corr[channel_idx],
                         baselines, bl_start, bl_stop)

        avg_traces[worker_id] += traces_corr

        pulses_slice = np.s_[pulse_offsets[index]:pulse_offsets[index]+pulse_counts[index]]

        for trigger, pulse_edges, pulse_amplitudes in zip(
            triggers[pulses_slice], edges[pulses_slice], amplitudes[pulses_slice]
        ):
            trigger_slice = np.s_[trigger['start']:trigger['stop']]

            for trace, channel_params, channel_edges, channel_amplitudes in zip(
                traces_corr, discr_params, pulse_edges, pulse_amplitudes
            ):
                discr_func(trace[trigger_slice], edges=channel_edges,
                           amplitudes=channel_amplitudes, **channel_params)

            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"]}')

    with timing(f'find_edges, {det_name}'):
        psh.map(find_edges, dc.select(det_sourcekeys))

    if not np.isfinite(edges).any():
        warning(f'No edges found for {det_name}')

    fig, (ux, bx) = plt.subplots(num=110+i, 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))

    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']):
        ux.hist(finite_flattened_slice(amplitudes, np.s_[:, edge_idx, :]),
                bins=1000, range=(0, 2048), histtype='step', lw=1,
                color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)

        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',
                lw=1, color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)

    ux.legend()
    ux.set_title('Pulse height distributions')
    ux.set_xlabel('Pulse height')
    ux.set_yscale('log')
    ux.set_xlim(0, 2048)
    ux.set_ylim(10, 1.5*ux.get_ylim()[1])

    bx.set_title('Fractional edge distributions')
    bx.set_xlabel('Edge positions - ⌊edge positions⌋')
    bx.set_yscale('log')
    bx.set_xlim(-0.05, 1.2)
    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:

 ### Global average of analog signals

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    fig, axs = plt.subplots(num=10+i, nrows=7, figsize=(9.5, 8), clear=True,
                            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')

    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].tick_params(labelbottom=False)
        axs[edge_idx].set_ylabel(edge_name)

    axs[-1].tick_params(labelbottom=True)
 pass
 ```

 %% Cell type:markdown id: tags:

 ### Sample for found edges

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    edges = det_data[det_name]['edges']

    fig = plt.figure(num=100+i, figsize=(9.5, 8))
    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')

    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
        col = (signal_idx % 2) if signal_idx < 6 else np.s_[:]
        ax = fig.add_subplot(grid[row, col])

        finite_edges = np.isfinite(edges[:, signal_idx, 0])
        if not finite_edges.any():
            warning(f'No edges found for {det_name}/{signal_name}')
            continue

        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]
        trigger = triggers[pulse_idx]

        sourcekey = remi.get_channel_sourcekey(
            remi['detector'][det_name]['channels'][signal_idx])
        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)

        remi.get_baseline_corrector()(
            train_trace, corr_trace,
            np.empty(remi['digitizer']['baseline_symmetry'], dtype=np.float64),
            *remi.get_baseline_limits(len(train_trace))[:2])

        pulse_trace = corr_trace[np.s_[trigger['start']:trigger['stop']]]

        x_time = remi.get_time_calibration() * (np.arange(len(pulse_trace)) + trigger['offset'])

        ax.plot(x_time, pulse_trace, lw=1)
        ax.set_xlim(x_time[0], x_time[-1])
        ax.set_ylim(-200, pulse_trace.max()*1.1)
        ax.text(x_time[-1], pulse_trace.max(),
                f'T{train_idx} P{pulse_idx - pulse_offsets[train_idx]} ',
                va='top', ha='right')
        ax.tick_params(labelbottom=False)
        ax.set_ylabel(signal_name)

        ax.vlines(edges[pulse_idx, signal_idx, :], *ax.get_ylim(), color='red', linewidth=1)
 pass
 ```

 %% Cell type:markdown id: tags:

 ### Digitized channel spectra

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    fig = plt.figure(num=20+i, figsize=(9.5, 6))

    edges = det_data[det_name]['edges']

    min_edge = np.nanmin(edges)
    max_edge = np.nanmax(edges)

    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')

    numx = fig.add_subplot(grid[0, 0])
    numx.set_title('Edges per pulse')

    agg_window = num_pulses // 60
    max_num_edges = 0.0
    max_spectral_intensity = 0
    hist_axs = []

    for edge_idx, edge_name in enumerate(['u1', 'u2', 'v1', 'v2', 'w1', 'w2', 'mcp']):
        if edge_idx < 6:
            row = 1 + edge_idx % 2
            col = edge_idx // 2
        else:
            row = 0
            col = np.s_[1:3]

        ax = fig.add_subplot(grid[row, col])
        ax.set_title(f'TOF spectrum: {edge_name}')

        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.reshape(-1, agg_window).mean(axis=1)

        if (num_edges == 0).all():
            warning(f'No edges found for {det_name}/{edge_name}')
            continue

        if edge_idx < 6:
            plot_kwargs = dict(c=f'C{edge_idx}', ls='solid', lw=1.0)
        else:
            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)
        max_num_edges = max(max_num_edges, num_edges.max())

        y, _, _ = ax.hist(finite_flattened_slice(edges, np.s_[:, edge_idx, :]),
                          bins=int((max_edge - min_edge) // 5), range=(min_edge, max_edge),
                          color=plot_kwargs['c'], histtype='step', linewidth=1)
        hist_axs.append(ax)

        max_spectral_intensity = max(max_spectral_intensity, y.max())

    numx.tick_params(labelbottom=False)
    numx.set_ylim(0, 1.2*max_num_edges)

    for ax in hist_axs:
        ax.set_ylim(0, max_spectral_intensity*1.1)
        ax.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
 pass
 ```

 %% Cell type:markdown id: tags:

 # Detector diagnostics

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    edges = det_data[det_name]['edges']

    sort = remi.get_dld_sorter(det_name)

    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,
                                     sum_limit=max(sort.uncorrected_time_sum_half_widths),
                                     sum_shifts=sum_shifts)

    if not is_valid.any():
        warning(f'No valid preliminary edge combinations found for {det_name}')

    signals, sums = remi.get_signals_and_sums(edges, indices=sort.channel_indices, sum_shifts=sum_shifts,
                                              mask=is_valid)
    fig = plot_detector_diagnostics(signals=signals, sums=sums, fig_num=30+i, im_scale=1.5,
                                    sum_range=max(sort.uncorrected_time_sum_half_widths),
                                    sorter=sort)
    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']:
        n_masked = is_valid.sum()
        signals = 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)
        fig = plot_detector_diagnostics(signals=signals, sums=sums, fig_num=40+i, im_scale=1.5,
                                        sum_range=max(sort.uncorrected_time_sum_half_widths),
                                        sorter=sort)
        fig.text(0.02, 0.98, det_name.upper() + ' after corrections', rotation=90, ha='left', va='top', size='x-large')
 pass
 ```

 %% Cell type:markdown id: tags:

 # Hit reconstruction

 %% Cell type:code id: tags:

 ``` python
 psh.set_default_context('processes', num_workers=remi.get_num_workers(mp_rec_hits))

 for det_name, det in remi['detector'].items():
    edges = det_data[det_name]['edges']

    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))
    hit_counts = psh.alloc(shape=len(dc.train_ids), dtype=np.uint32)

    def prepare_hit_worker(worker_id):
        sort = remi.get_dld_sorter(det_name)
        yield

    @psh.with_init(prepare_hit_worker)
    def reconstruct_hits(worker_id, index, train_id):
        hit_counts[index] += sort.run_on_train(
            edges, signals, hits, pulse_offsets[index], pulse_offsets[index] + pulse_counts[index])

    with timing(f'rec_hits, {det_name}'):
        psh.map(reconstruct_hits, dc.train_ids)

    det_data[det_name].update(signals=signals, hits=hits, hit_counts=hit_counts)
 ```

 %% Cell type:code id: tags:

 ``` python
 fig, ax = plt.subplots(num=50+i, figsize=(9.5, 4), ncols=1, clear=True,
                       gridspec_kw=dict(top=0.92, right=0.98, left=0.05, bottom=0.12))

 max_num_hits = 0.0

 for det_name in remi['detector'].keys():
    agg_window = num_pulses // min(1000, num_pulses)

    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.reshape(-1, agg_window).mean(axis=1)
    max_num_hits = max(max_num_hits, num_hits.max())

    ax.plot(np.arange(0, (num_pulses // agg_window) * agg_window, agg_window), num_hits,
            lw=1, label=det_name.upper())

 ax.set_title('Hits per pulse')
 ax.set_xlabel('Pulse index')
 ax.set_ylim(0, max_num_hits*1.1)
 ax.legend()
 pass
 ```

 %% Cell type:markdown id: tags:

 ### 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.:

 * `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.
 * `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.

 | Method | `u+v+w +mcp` |
 | - | - |
 | 0 | `2+2+2 +1` |
 | 1 | `0+2+2 +1` |
 | 2 | `2+0+2 +1` |
 | 3 | `2+2+0 +1` |
 | 4 | `1+2+2 +1` (2 permutations) |
 | 5 | `2+1+2 +1` (2 permutations) |
 | 6 | `2+2+1 +1` (2 permutations) |
 | 7 | `2+2+2 +0` |
 | 8 | `0+2+2 +0` |
 | 9 | `2+0+2 +0` |
 | 10 | `2+2+0 +0` |
 | 11 | `1+2+2 +0` (2 permutations) |
 | 12 | `2+1+2 +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) |
 | 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) |
 | 17 | `2+1+1 +0` `1+2+1 +0` `1+1+2 +0` (12 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) |

 * 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.

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    hits = det_data[det_name]['hits']

    fig, ax = plt.subplots(num=60+i, figsize=(9.5, 5), ncols=1, clear=True,
                           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')

    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)
    ax.bar(np.arange(20), method_bins, width=0.5)

    ax.set_xlabel('Reconstruction method')
    ax.set_xlim(-0.5, 19.5)
    ax.set_xticks(np.arange(20))

    ax.set_ylabel('Number of hits')
    ax.set_ylim(0, method_bins.max()*1.05)
    ylims = ax.get_ylim()

    ax.tick_params(which='both', right=True, labelright=True)

    num_risky = method_bins[15:].sum()
    num_total = method_bins.sum()

    ax.text(14.2, method_bins.max(), f'{(100*(num_total-num_risky)/num_total):.2g}%',
            va='top', ha='right', color='black')
    ax.text(14.8, method_bins.max(), f'{(100*num_risky/num_total):.2g}%',
            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)

    labelx = ax.twiny()
    labelx.set_xlim(*ax.get_xlim())
    labelx.set_xticks(ax.get_xticks())
    labelx.set_xticklabels([
        '2+2+2 +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',
        '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',
        '2+1+1 +1',
        '2+1+0 +1',
        '1+1+1 +1',
        '2+1+1 +0',
        '1+1+1 +0',
        '2+1+0 +0',
    ], rotation=90)

    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

    rely = ax.twinx()
    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, 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_ylabel('Percentage of total hits')
 pass
 ```

 %% Cell type:markdown id: tags:

 ### Detector image and fishes

 %% Cell type:code id: tags:

 ``` python
 for i, det_name in enumerate(remi['detector'].keys()):
    flat_hits = det_data[det_name]['hits'].reshape(-1)
    flat_hits = flat_hits[np.isfinite(flat_hits[:]['x'])]
    flat_hits = flat_hits[flat_hits['m'] <= 10]

    fig = plt.figure(num=70+i, figsize=(9, 13.5))

    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')

    imp = fig.add_axes([0.1 + 0.25/2, 0.56, 0.6, 0.4])
    txp = fig.add_axes([0.1, 0.28, 0.85, 0.22])
    typ = fig.add_axes([0.1, 0.04, 0.85, 0.22])

    if flat_hits.size == 0:
        warning(f'No hits found for {det_name}')
        continue

    im_radius = remi['detector'][det_name]['mcp_radius']*1.1

    imp.hist2d(flat_hits['x'], flat_hits['y'], bins=(256, 256),
               range=[[-im_radius, im_radius], [-im_radius, im_radius]], norm=LogNorm())
    imp.xaxis.set_label_position('top')
    imp.set_xlabel('X / mm')
    imp.set_ylabel('Y / mm')
    imp.tick_params(right=True, labelright=True, top=True, labeltop=True)
    imp.grid()

    min_tof = flat_hits['t'].min()
    max_tof = flat_hits['t'].max()

-    num_tof_bins = int((max_tof - min_tof) // 5)
+    num_tof_bins = min(int((max_tof - min_tof) // 10), 500)

    if num_tof_bins == 0:
        warning(f'All TOFs limited to single bin for {det_name}')
        continue

    for ax, dim_label in zip([txp, typ], ['x', 'y']):
        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())
        ax.set_ylabel(f'{dim_label.upper()} / mm')

    typ.set_xlabel('Time-of-flight / ns')
    txp.tick_params(bottom=True, labelbottom=False, top=True, labeltop=True, right=True, labelright=True)
    typ.tick_params(right=True, labelright=True, top=True)
 pass
 ```

 %% Cell type:markdown id: tags:

 # Transformed data files

 %% Cell type:code id: tags:

 ``` python
 # 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])
 if not proposal and m is not None:
    proposal = int(m[1])

 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')}

 control_sources = [det_device_id.format(karabo_id=karabo_id, det_name=det_name.upper())
                   for det_name in remi['detector']]

 channels = []
 if save_raw_triggers or save_raw_edges:
    channels.append('raw')
 if save_rec_signals or save_rec_hits:
    channels.append('rec')

 instrument_channels = [
    f'{device_id}:{det_output_key}/{channel}'
    for device_id in control_sources
    for channel in channels
 ]
 ```

 %% Cell type:code id: tags:

 ``` python
 Path(out_folder).mkdir(parents=True, exist_ok=True)
 print('Writing sequence files', flush=True, end='')

 t_write = timing('write_files')
 t_write.__enter__()

 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]

    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,
                             control_sources=control_sources, instrument_channels=instrument_channels)
        outp.create_index(
            seq_train_ids,
            timestamps=dc.select_trains(by_id[seq_train_ids]).train_timestamps().astype(np.uint64)
        )

        for det_name in remi['detector']:
            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_control_data = outp.create_control_source(cur_device_id)
            # Manually manipulate the file here, still creates the index properly.
            remi.attach_detector_config(det_name, cur_control_data.get_run_group())
            cur_control_data.create_index(len(seq_train_ids))

            cur_fast_data = outp.create_instrument_source(f'{cur_device_id}:{det_output_key}')

            cur_data = det_data[det_name]

            if save_raw_triggers:
                cur_fast_data.create_key('raw.triggers', triggers[pulse_mask],
                                         maxshape=(None,) + triggers.shape[1:],
                                         chunks=tuple(chunks_triggers), **dataset_kwargs)

            if save_raw_edges:
                cur_fast_data.create_key('raw.edges', cur_data['edges'][pulse_mask],
                                         maxshape=(None,) + cur_data['edges'].shape[1:],
                                         chunks=tuple(chunks_edges if chunks_edges[-1] <= cur_max_hits
                                                      else chunks_edges[:-1] + [cur_max_hits]),
                                         **dataset_kwargs)

            if save_raw_amplitudes:
                cur_fast_data.create_key('raw.amplitudes', cur_data['amplitudes'][pulse_mask],
                                         maxshape=(None,) + cur_data['amplitudes'].shape[1:],
                                         chunks=tuple(chunks_amplitudes if chunks_amplitudes[-1] <= cur_max_hits
                                                      else chunks_amplitudes[:-1] + [cur_max_hits]),
                                         **dataset_kwargs)

            if save_rec_signals:
                cur_fast_data.create_key('rec.signals', cur_data['signals'][pulse_mask],
                                         maxshape=(None,) + cur_data['signals'].shape[1:],
                                         chunks=tuple(chunks_signals if chunks_signals[-1] <= cur_max_hits
                                                      else chunks_signals[:-1] + [cur_max_hits]),
                                         **dataset_kwargs)

            if save_rec_hits:
                cur_fast_data.create_key('rec.hits', cur_data['hits'][pulse_mask],
                                         maxshape=(None,) + hits.shape[1:],
                                         chunks=tuple(chunks_hits if chunks_hits[-1] <= cur_max_hits
                                                      else chunks_hits[:-1] + [cur_max_hits]),
                                         **dataset_kwargs)

            cur_fast_data.create_index(raw=pulse_counts[train_mask], rec=pulse_counts[train_mask])

    print('.', flush=True, end='')

 print('')
 t_write.__exit__()
 ```