[REMI] Save pulse amplitudes during discrimination

Description

The first step in actual signal processing for the Quad/Hex DLD detectors is discrimination of the analog signal into digital edges. In a recent REMI beamtime, one of the detectors seem to require constant tuning of their discrimination thresholds indicating unusually variation in the pulse heights. So far it was never necessary to look at the pulse heights, so this MR now adds tracking, plotting (see below) and subsequent storage of the pulse amplitudes corresponding to each edge found. The tracking part lives in the actual Cython implementations currently still in the EXtra-remi package.

How Has This Been Tested?

run = 50  # Run ID.
in_folder = '/gpfs/exfel/exp/SQS/202301/p003090/raw'  # Partial input path appended with run ID.
calib_config_path = '/gpfs/exfel/exp/SQS/202301/p003090/usr/remi_proc_config.yaml'  # Path to correction and transform configuration
ppl_offset = -96  # In units of the PPT.
first_pulse_offset = 18000  # Sample position where the first pulse begins, ignored when PPT is reconstructed.

Relevant Documents (optional)

Top plot shows pulse height distribution, bottom plot shows fractional component of edge positions after interpolation. The latter plot is useful to judge the quality of interpolation used (currently linear, but could be extended to slower sinc interpolation in the future if necessary). Given enough samples, this distribution is expected to be constant (i.e. any fractional position occurs with equal probability). Linear interpolation for very steep slopes can e.g. cause this distribution to be non-constant, indicating poor interpolation quality.

Types of changes

• New feature (non-breaking change which adds functionality)

Reviewers

@ahmedk

notebooks/REMI/REMI_Digitize_and_Transform.ipynb

@@ -28,10 +28,12 @@
     "det_output_key = 'output'  # Pipeline name for fast data output.\n",
     "save_raw_triggers = True  # Whether to save trigger position in files.\n",
     "save_raw_edges = True  # Whether to save digitized edge positions in files.\n",
+    "save_raw_amplitudes = True  # Whether to save analog pulse amplitudes in files.\n",
     "save_rec_signals = True  # Whether to save reconstructed signals (u1-w2, mcp) in files.\n",
     "save_rec_hits = True  # Whether to save reoncstructed hits (x,y,t,m) in files.\n",
     "chunks_triggers = [500]  # HDF chunk size for triggers.\n",
     "chunks_edges = [500, 7, 50]  # HDF chunk size for edges.\n",
+    "chunks_amplitudes = [500, 7, 50]  # HDF chunk size for amplitudes.\n",
     "chunks_hits = [50, 50]  # HDF chunk size for hits.\n",
     "chunks_signals = [50, 50]  # HDF chunk size for signals.\n",
     "dataset_compression = 'gzip'  # HDF compression method.\n",
@@ -188,7 +190,7 @@
     "    discr_func, discr_params = remi.get_discriminator([trigger_edge_channel])\n",
     "\n",
     "    edges = np.zeros(1000, dtype=np.float64)\n",
-    "    num_pulses = discr_func(trigger_trace, edges, **discr_params[0])\n",
+    "    num_pulses = discr_func(trigger_trace, edges=edges, **discr_params[0])\n",
     "    edges = edges[:num_pulses]\n",
     "\n",
     "    first_edge = edges[0]\n",
@@ -514,13 +516,15 @@
     "\n",
     "det_data = {}\n",
     "\n",
-    "for det_name, det in remi['detector'].items():\n",
+    "for i, (det_name, det) in enumerate(remi['detector'].items()):\n",
     "    det_sourcekeys = remi.get_detector_sourcekeys(det_name)\n",
     "    det_get_traces = remi.get_traces_getter(det_name)\n",
     "    trace_len = dc[next(iter(det_sourcekeys))].entry_shape[0]\n",
     "    \n",
     "    edges = psh.alloc(shape=(num_pulses, 7, det['max_hits']),\n",
     "                      dtype=np.float64, fill=np.nan)\n",
+    "    amplitudes = psh.alloc(shape=(num_pulses, 7, det['max_hits']),\n",
+    "                           dtype=np.float64, fill=np.nan)\n",
     "    avg_traces = psh.alloc_per_worker(shape=(7, trace_len), dtype=np.float64)\n",
     "    \n",
     "    def prepare_edge_worker(worker_id):\n",
@@ -551,24 +555,55 @@
     "\n",
     "        pulses_slice = np.s_[pulse_offsets[index]:pulse_offsets[index]+pulse_counts[index]]\n",
     "\n",
-    "        for trigger, pulse_edges in zip(triggers[pulses_slice], edges[pulses_slice]):\n",
+    "        for trigger, pulse_edges, pulse_amplitudes in zip(\n",
+    "            triggers[pulses_slice], edges[pulses_slice], amplitudes[pulses_slice]\n",
+    "        ):\n",
     "            trigger_slice = np.s_[trigger['start']:trigger['stop']]\n",
     "                                                 \n",
-    "            for trace, channel_params, channel_edges in zip(traces_corr, discr_params, pulse_edges):\n",
-    "                discr_func(trace[trigger_slice], channel_edges, **channel_params)\n",
-    "\n",
-    "            pulse_edges += trigger['offset']\n",
-    "            pulse_edges *= time_cal\n",
+    "            for trace, channel_params, channel_edges, channel_amplitudes in zip(\n",
+    "                traces_corr, discr_params, pulse_edges, pulse_amplitudes\n",
+    "            ):\n",
+    "                discr_func(trace[trigger_slice], edges=channel_edges,\n",
+    "                           amplitudes=channel_amplitudes, **channel_params)\n",
     "            \n",
     "    with timing(f'find_edges, {det_name}'):\n",
     "        psh.map(find_edges, dc.select(det_sourcekeys))\n",
     "    \n",
+    "    fig, (ux, bx) = plt.subplots(num=110+i, ncols=1, nrows=2, figsize=(9.5, 8), clear=True,\n",
+    "                                 gridspec_kw=dict(left=0.1, right=0.98, top=0.98, bottom=0.1, hspace=0.25))\n",
+    "    \n",
+    "    fig.text(0.02, 0.98, det_name.upper(), rotation=90, ha='left', va='top', size='x-large')\n",
+    "\n",
+    "    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",
+    "                bins=1000, range=(0, 2048), histtype='step', lw=1,\n",
+    "                color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)\n",
+    "        \n",
+    "        cur_edges = finite_flattened_slice(edges, np.s_[:, edge_idx, :])\n",
+    "        bx.hist(cur_edges - np.floor(cur_edges), bins=500, range=(0, 1), histtype='step',\n",
+    "                lw=1, color=f'C{edge_idx}' if edge_idx < 6 else 'k', label=edge_name)\n",
+    "    \n",
+    "    ux.legend()\n",
+    "    ux.set_title('Pulse height distributions')\n",
+    "    ux.set_xlabel('Pulse height')\n",
+    "    ux.set_yscale('log')\n",
+    "    ux.set_xlim(0, 2048)\n",
+    "    ux.set_ylim(10, 1.5*ux.get_xlim()[1])\n",
     "    \n",
-    "    with np.printoptions(precision=2, suppress=True):\n",
-    "        print(edges[:5, :, :8])"
+    "    bx.set_title('Fractional edge distributions')\n",
+    "    bx.set_xlabel('Edge positions - ⌊edge positions⌋')\n",
+    "    bx.set_yscale('log')\n",
+    "    bx.set_xlim(-0.05, 1.2)\n",
+    "    bx.legend()\n",
+    "    \n",
+    "    # Properly offset edges to their trigger offset and convert to time.\n",
+    "    # This is not done earlier to preserve the information\n",
+    "    edges += triggers['offset'][:, None, None]\n",
+    "    edges *= remi.get_time_calibration()\n",
     "    \n",
     "    det_data[det_name] = {\n",
     "        'edges': edges,\n",
+    "        'amplitudes': amplitudes,\n",
     "        'avg_trace': avg_traces.sum(axis=0) / len(dc.train_ids)\n",
     "    }"
    ]
@@ -1081,6 +1116,10 @@
     "                cur_fast_data.create_key('raw.edges', cur_data['edges'][pulse_mask],\n",
     "                                         chunks=tuple(chunks_edges), **dataset_kwargs)\n",
     "                \n",
+    "            if save_raw_amplitudes:\n",
+    "                cur_fast_data.create_key('raw.amplitudes', cur_data['amplitudes'][pulse_mask],\n",
+    "                                         chunks=tuple(chunks_amplitudes), **dataset_kwargs)\n",
+    "                \n",
     "            if save_rec_signals:\n",
     "                cur_fast_data.create_key('rec.signals', cur_data['signals'][pulse_mask],\n",
     "                                         chunks=tuple(chunks_signals), **dataset_kwargs)\n",