Merge branch 'digitizers_LM' into 'master'

Simplify digitizer functions and pulse ID coordinates assignment for XGM, digitizers

See merge request !100

doc/bunch_pattern_decoding.rst

+The bunch pattern table is an array of 2700 values per train (the maximum number of pulses at 4.5 MHz provided by the machine) and contains information on how the pulses are distributed among SASE 1, 2, 3, and the various lasers at European XFEL.
+The data stored in the bunch pattern table (mnemonic *bunchPatternTable*) can be extracted using the wrappers to the `euxfel_bunch_pattern <https://pypi.org/project/euxfel-bunch-pattern/>`_ package as follows:
+
 .. code:: ipython3
 
     import toolbox_scs as tb
@@ -6,29 +9,17 @@
     proposalNB = 2511
     runNB = 176
 
-**option 1**
-
-This method uses function we implemented.
-
-.. code:: ipython3
+    run, data = tb.load(proposalNB, runNB, "bunchPatternTable")
+    ppl_mask = tbm.is_ppl(data.bunchPatternTable)
 
-    fields = ["bunchPatternTable"]
-    run = tb.load(fields, runNB, proposalNB)
-    bpt = run['bunchPatternTable']
+ppl_mask is a boolean DataArray of dimensions trainId x 2700, where True values indicate where a laser pulse from the PP laser was triggered.
 
-    bpt_dec = tbm.extractBunchPattern(
-                        run['bunchPatternTable'],'scs_ppl')
+.. note::
+   The position of the PP laser pulses with respect to that of the SASE 3 pulses is arbitrary. The PP laser pattern always starts at pulse Id 0, while that of SASE 3 can vary, depending on the machine parameters.
 
-
-**option 2**
-
-This method uses function from the euxfel_bunch_pattern package.
+From this mask, one can obtain the number of pulses per train by summing along the 'pulse_slot' dimension:
 
 .. code:: ipython3
 
+   ppl_npulses = ppl_mask.sum(dim='pulse_slot')
 
-    run = tb.load_run(proposalNB, runNB)
-    mnemonic = tb.mnemonics["bunchPatternTable"]
-    bpt = run.get_array(*mnemonic.values())
-
-    bpt_is_laser = tbm.is_ppl(bpt)

doc/conf.py

@@ -43,12 +43,17 @@ extensions = [
     'sphinx.ext.viewcode',
     'sphinx.ext.coverage',
     'sphinx.ext.napoleon',
-    'autoapi.extension',
+    # Comment out autoapi to be able to run nbsphinx
+    # 'autoapi.extension',
     'sphinx_rtd_theme',
+    'nbsphinx'
 ]
 autoapi_dirs = ['../src/toolbox_scs']
 autoapi_ignore = ['*/deprecated/*']
 
+# Don't add .txt suffix to source files:
+html_sourcelink_suffix = ''
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
 

doc/howtos.rst

@@ -4,7 +4,8 @@
 top
 ---
 
-* :doc:`load run and data <load>`.
+.. * :doc:`load run and data <load>`.
+* :doc:`load data in memory <Loading_data_in_memory>`.
 
 misc
 ----
@@ -12,15 +13,23 @@ misc
 * :doc:`bunch pattern decoding <bunch_pattern_decoding>`.
 
 
-detectors (dssc)
-----------------
+detectors
+---------
 
+DSSC
+++++
 Most of the functions within toolbox_scs.detectors can be accessed directly. This is useful during development, or when working in a non-standardized way, which is often neccessary during data evaluation. For frequent routines there is the possibility to use dssc objects that guarantee consistent data structure, and reduce the amount of recurring code within the notebook.
 
 * bin data using toolbox_scs.tbdet -> *to be documented*.
 * :doc:`bin data using the DSSCBinner <dssc/DSSCBinner>`.
 * post processing, data analysis -> *to be documented*
 
+Point detectors
++++++++++++++++
+Detectors that produce one point per pulse, or 0D detectors, are all handled in a similar way. Such detectors are, for instance, the X-ray Gas Monitor (XGM), the Transmitted Intensity Monitor (TIM), the electron Bunch Arrival Monitor (BAM) or the photo diodes monitoring the PP laser.
+
+* :doc:`extract data from point detectors <point_detectors/point_detectors>`.
+
 routines
 --------
 

doc/load.rst

+Loading data in memory is performed as follows:
+
 **Option 1**:
 
 .. code:: python3
@@ -13,9 +15,9 @@
     proposalNr = 2565
     runNr = 19
 
-    run_data = tb.load(fields, runNr, proposalNr)
+    run, data = tb.load(proposalNr, runNr, fields)
 
-run_data is an xarray dataArray. It has an attribute called 'run' containing the underlying extra_data dataCollection.
+run is an extra_data dataCollection and data is an xarray Dataset containing all variables listed in fields. For convinience, data also contains the variable bunchPatternTable, which is used by other functions of the ToolBox. All variables are aligned by train Id.
 
 **Option 2**:
 
@@ -28,7 +30,7 @@ run_data is an xarray dataArray. It has an attribute called 'run' containing the
     proposalNr = 2565
     runNr = 19
 
-    run = tb.load_run(proposalNr, runNr)
-    run_data = run.get_array(*mnemonic.values())
+    run, _ = tb.load(proposalNr, runNr)
+    data = run.get_array(*mnemonic.values())
 
-run is an extra_data dataCollection and run_data an xarray dataArray for a single data source.
+run is an extra_data dataCollection and data an xarray dataArray for a single data source.

doc/point_detectors/point_detectors.rst

+blablabla
\ No newline at end of file

doc/requirements.txt

 sphinx
 sphinx_rtd_theme
 autoapi
-sphinx-autoapi
\ No newline at end of file
+sphinx-autoapi
+nbsphinx
\ No newline at end of file

notebook_examples/XAS and XMCD energy shift investigation.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook shows an example of X-ray absorption spectroscopy (XAS) and X-ray magnetic circular dichroism (XMCD) experiment."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Import toolbox_scs and subpackages"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "%matplotlib notebook\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.rcParams['figure.constrained_layout.use'] = True\n",
+    "\n",
+    "import toolbox_scs as tb\n",
+    "import toolbox_scs.detectors as tbdet\n",
+    "import toolbox_scs.routines as tbr\n",
+    "import extra_data as ed\n",
+    "import xarray as xr\n",
+    "\n",
+    "import logging\n",
+    "logging.basicConfig(level=logging.INFO)\n",
+    "log_root = logging.getLogger(__name__)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Load data, extract pulse-resolved signals from detectors\n",
+    "\n",
+    "Here, we use the ToolBox mnemonics to get the arrays: MCP1apd is the pulse-resolved peak-integrated TIM data from MCP 1, SCS_SA3 is the XGM data for SASE 3 pulses in SCS hutch, nrj is the monochromator energy in eV and magnet is the current applied to the magnet in A.\n",
+    "\n",
+    "The function get_all_detectors() takes care of aligning the pulse-resolved data (TIM and XGM) according to the sase 3 bunch pattern."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "runNB = 491\n",
+    "proposalNB = 900094\n",
+    "mcp = 'MCP1apd'\n",
+    "fields = [mcp, 'SCS_SA3', 'nrj', 'magnet']\n",
+    "run, ds = tb.load(proposalNB, runNB, fields)\n",
+    "mcp = mcp.replace('apd', 'peaks')\n",
+    "ds\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Split positive and negative magnetic fields, calculate XAS and XMCD spectra\n",
+    "\n",
+    "Here we first define the monochromator energy bins, and apply the xas() subroutine to subsets of the data. The subsets correspond to positive magnet current and negative magnet current. The xas() routine returns (among other quantities) the binned energy and the absorption, which are then plotted for both positive and negative magnet currents."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nrj_bins_edges = np.linspace(704, 712, 31)\n",
+    "pxas = tbr.xas(ds.where(ds['magnet'] > 2, drop=True), bins=nrj_bins_edges, plot=False, Itkey=mcp)\n",
+    "nxas = tbr.xas(ds.where(ds['magnet'] < -2, drop=True), bins=nrj_bins_edges, plot=False, Itkey=mcp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(pxas['nrj'], pxas['muA'])\n",
+    "plt.plot(nxas['nrj'], nxas['muA'])\n",
+    "plt.ylabel('$\\mu_A$')\n",
+    "plt.xlabel('Energy [eV]')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now calculate the XMCD using the positive and negative data returned by xas() routine."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = tbr.xasxmcd(pxas, nxas)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(res['nrj'], res['muXAS'])\n",
+    "plt.ylabel('$\\mu_A$')\n",
+    "plt.xlabel('Energy [eV]')\n",
+    "plt.twinx()\n",
+    "plt.plot(res['nrj'], -res['muXMCD'], c='C1')\n",
+    "plt.ylabel('$\\mu_{XMCD}$')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Only 2000 first train"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "subds = ds.isel({'trainId':slice(1600)})\n",
+    "\n",
+    "pxas2 = tbr.xas(subds.where(subds['magnet'] > 2, drop=True), bins=nrj_bins_edges, plot=False, Itkey=mcp)\n",
+    "nxas2 = tbr.xas(subds.where(subds['magnet'] < -2, drop=True), bins=nrj_bins_edges, plot=False, Itkey=mcp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(subds['nrj'])\n",
+    "plt.ylabel('energy [eV]')\n",
+    "plt.xlabel('train number')\n",
+    "plt.twinx()\n",
+    "plt.plot(subds['magnet'], c='C1')\n",
+    "plt.ylabel('magnet current [A]')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.errorbar(pxas['nrj'], pxas['muA'], pxas['sterrA'], label='pos all trains')\n",
+    "plt.errorbar(pxas2['nrj'], pxas2['muA'], pxas2['sterrA'], c='C0', ls='--', label='pos first 2000 trains')\n",
+    "\n",
+    "plt.errorbar(nxas['nrj'], nxas['muA'], nxas['sterrA'], label='neg all trains')\n",
+    "plt.errorbar(nxas2['nrj'], nxas2['muA'], nxas['sterrA'], c='C1', ls='--', label='neg first 2000 trains')\n",
+    "plt.legend()\n",
+    "plt.ylabel('$\\mu_A$')\n",
+    "plt.xlabel('Energy [eV]')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res2 = tbr.xasxmcd(pxas2, nxas2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.plot(res['nrj'], res['muXAS'])\n",
+    "plt.plot(res2['nrj'], res2['muXAS'], ls='--', c='C0')\n",
+    "plt.ylabel('$\\mu_A$')\n",
+    "plt.xlabel('Energy [eV]')\n",
+    "plt.twinx()\n",
+    "plt.plot(res['nrj'], -res['muXMCD'], c='C1')\n",
+    "plt.plot(res2['nrj'], -res2['muXMCD'], ls='--', c='C1')\n",
+    "plt.ylabel('$\\mu_{XMCD}$')\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Mono up or down"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.signal import savgol_filter\n",
+    "\n",
+    "dE = np.gradient(savgol_filter(ds['nrj'], 9, 1))\n",
+    "mono_turn = dE > 0\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.plot(1e3*dE)\n",
+    "plt.ylabel('Delta energy (meV)')\n",
+    "plt.twinx()\n",
+    "plt.plot(mono_turn, ls='', marker='.', c='C1')\n",
+    "plt.ylabel('mono going up')\n",
+    "plt.title(f'run {runNB} p{proposalNB}')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds['mono_turn'] = xr.DataArray(np.gradient(savgol_filter(ds['nrj'], 9, 1))>0, dims=['trainId'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "subds = ds#.isel({'trainId':slice(1600)})\n",
+    "\n",
+    "Eshift = 0.05\n",
+    "shift_subds = subds\n",
+    "shift_subds['nrj'] = subds['nrj'] + Eshift\n",
+    "pxas3a = tbr.xas(shift_subds.where((subds['magnet'] > 2)*subds['mono_turn'], drop=True),\n",
+    "                bins=nrj_bins_edges, plot=False, Itkey=mcp)\n",
+    "nxas3a = tbr.xas(shift_subds.where((subds['magnet'] < -2)*subds['mono_turn'], drop=True),\n",
+    "                bins=nrj_bins_edges, plot=False, Itkey=mcp)\n",
+    "\n",
+    "shift_subds['nrj'] = subds['nrj'] - Eshift\n",
+    "pxas3b = tbr.xas(shift_subds.where((subds['magnet'] > 2)*np.logical_not(subds['mono_turn']), drop=True),\n",
+    "                bins=nrj_bins_edges, plot=False, Itkey=mcp)\n",
+    "nxas3b = tbr.xas(shift_subds.where((subds['magnet'] < -2)*np.logical_not(subds['mono_turn']), drop=True),\n",
+    "                bins=nrj_bins_edges, plot=False, Itkey=mcp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "#plt.errorbar(pxas['nrj'], pxas['muA'], pxas['sterrA'])\n",
+    "plt.errorbar(pxas3a['nrj'], pxas3a['muA'], pxas3a['sterrA'], c='C0', ls='-', label='+H, mono up')\n",
+    "plt.errorbar(pxas3b['nrj'], pxas3b['muA'], pxas3b['sterrA'], c='C0', ls='--', \n",
+    "             label=f'+H, mono down, dE = {-1e3*Eshift} meV')\n",
+    "\n",
+    "#plt.errorbar(nxas['nrj'], nxas['muA'], nxas['sterrA'])\n",
+    "plt.errorbar(nxas3a['nrj'], nxas3a['muA'], nxas3a['sterrA'], c='C1', ls='-', label='-H, mono up')\n",
+    "plt.errorbar(nxas3b['nrj'], nxas3b['muA'], nxas3b['sterrA'], c='C1', ls='--', \n",
+    "             label=f'-H, mono down, dE = {-1e3*Eshift} meV')\n",
+    "\n",
+    "plt.xlabel('Energy (eV)')\n",
+    "plt.ylabel('XAS (arb. units)')\n",
+    "plt.title(f'run {runNB} p{proposalNB}')\n",
+    "plt.xlim([707, 709])\n",
+    "plt.ylim([3, 4.5])\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure()\n",
+    "plt.title(f'run {runNB} p{proposalNB}')\n",
+    "plt.xlabel('Energy (eV)')\n",
+    "plt.ylabel('XAS (mono up) - XAS (mono down) (arb. units)')\n",
+    "plt.plot(pxas3a['nrj'], pxas3a['muA'] - pxas3b['muA'])\n",
+    "plt.plot(nxas3a['nrj'], nxas3a['muA'] - nxas3b['muA'])\n",
+    "plt.ylim([-0.1, 0.18])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "xfel",
+   "language": "python",
+   "name": "xfel"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

src/toolbox_scs/__init__.py

-from .load import (load, concatenateRuns, get_array,
-                   load_run, run_by_path)
+from .load import (load, concatenateRuns, get_array, run_by_path)
 
 from .constants import mnemonics
 
@@ -8,7 +7,6 @@ __all__ = (
     "load",
     "concatenateRuns",
     "get_array",
-    "load_run",
     "run_by_path",
     # Classes
     # Variables

src/toolbox_scs/detectors/__init__.py

 from .xgm import (
-    load_xgm, cleanXGMdata, matchXgmTimPulseId, calibrateXGMs,
-    autoFindFastAdcPeaks)
-from .tim import (
-    load_TIM,)
+    get_xgm, calibrate_xgm)
+from .digitizers import (
+    get_peaks, get_tim_peaks, get_laser_peaks, get_digitizer_peaks,
+    check_peak_params)
+from .bam_detectors import get_bam
 from .dssc_data import (
     save_xarray, load_xarray, get_data_formatted, save_attributes_h5)
 from .dssc_misc import (
@@ -12,25 +13,30 @@ from .dssc_processing import (
     process_dssc_data)
 from .dssc import (
     DSSCBinner, DSSCFormatter)
-from .azimuthal_integrator import AzimuthalIntegrator, AzimuthalIntegratorDSSC
+from .azimuthal_integrator import (
+    AzimuthalIntegrator, AzimuthalIntegratorDSSC)
 
 __all__ = (
     # Functions
-    "load_xgm",
-    "cleanXGMdata",
-    "load_TIM",
-    "matchXgmTimPulseId",
-    "load_dssc_info",
-    "create_dssc_bins",
-    "calc_xgm_frame_indices",
-    "get_xgm_formatted",
-    "process_dssc_data",
+    "get_xgm",
+    "calibrate_xgm",
+    "get_peaks",
+    "get_tim_peaks",
+    "get_laser_peaks",
+    "get_digitizer_peaks",
+    "check_peak_params",
+    "get_bam",
     "save_xarray",
     "load_xarray",
     "get_data_formatted",
     "save_attributes_h5",
+    "load_dssc_info",
+    "create_dssc_bins",
     "quickmask_DSSC_ASIC",
+    "get_xgm_formatted",
     "load_mask",
+    "calc_xgm_frame_indices",
+    "process_dssc_data",
     # Classes
     "DSSCBinner",
     "DSSCFormatter",
@@ -58,8 +64,9 @@ clean_ns = [
     'dssc_plot',
     'azimuthal_integrator',
     'FastCCD',
-    'tim',
     'xgm',
+    'digitizers',
+    'bam_detectors'
     ]
 
 

src/toolbox_scs/detectors/bam_detectors.py

+""" Beam Arrival Monitor related sub-routines
+
+    Copyright (2021) SCS Team.
+
+    (contributions preferrably comply with pep8 code structure
+    guidelines.)
+"""
+
+import logging
+
+import numpy as np
+import xarray as xr
+
+from ..constants import mnemonics as _mnemonics
+from ..misc.bunch_pattern_external import is_pulse_at
+import toolbox_scs.load as tbload
+
+
+log = logging.getLogger(__name__)
+
+
+def get_bam(run, mnemonics=None, merge_with=None, bunchPattern='sase3'):
+    """
+    Load beam arrival monitor (BAM) data and align their pulse ID
+    according to the bunch pattern. Sources can be loaded on the fly
+    via the mnemonics argument, or processed from an existing data set
+    (merge_with).
+
+    Parameters
+    ----------
+    run: extra_data.DataCollection
+        DataCollection containing the digitizer data.
+    mnemonics: str or list of str
+        mnemonics for BAM, e.g. "BAM1932M" or ["BAM414", "BAM1932M"].
+        If None, defaults to "BAM1932M" in case no merge_with dataset
+        is provided.
+    merge_with: xarray Dataset
+        If provided, the resulting Dataset will be merged with this
+        one. The BAM variables of merge_with (if any) will also be
+        selected, aligned and merged.
+    bunchPattern: str
+        'sase1' or 'sase3' or 'scs_ppl', bunch pattern
+        used to extract peaks. The pulse ID dimension will be named
+        'sa1_pId', 'sa3_pId' or 'ol_pId', respectively.
+    Returns
+    -------
+    xarray Dataset with pulse-resolved BAM variables aligned,
+        merged with Dataset *merge_with* if provided.
+
+    Example
+    -------
+    >>> import toolbox_scs as tb
+    >>> import toolbox_scs.detectors as tbdet
+    >>> run, _ = tb.load(2711, 303)
+    >>> bam = tbdet.get_bam(run)
+
+    """
+    # get the list of mnemonics to process
+    mnemonics = tbload.mnemonics_to_process(mnemonics, merge_with, 'BAM')
+
+    if len(mnemonics) == 0:
+        log.info('No array with unaligned BAM data to extract. Skipping.')
+        return merge_with
+    else:
+        log.info(f'Extracting BAM data from {mnemonics}.')
+
+    # Prepare the dataset of non-BAM data to merge with
+    if bool(merge_with):
+        mw_ds = merge_with.drop(mnemonics, errors='ignore')
+    else:
+        mw_ds = xr.Dataset()
+
+    # extract the XFEL pulses: slice(0,5400,2)
+    roi = np.s_[:5400:2]
+    ds = xr.Dataset()
+    for m in mnemonics:
+        if bool(merge_with) and m in merge_with:
+            val = merge_with[m].isel({_mnemonics[m]['dim'][0]: roi})
+            log.debug(f'Using {m} from merge_with dataset.')
+        else:
+            val = run.get_array(*_mnemonics[m].values(), roi=roi, name=m)
+            log.debug(f'Loading {m} from DataCollection.')
+        val[_mnemonics[m]['dim'][0]] = np.arange(2700)
+        ds = ds.merge(val, join='inner')
+
+    # check if bunch pattern table exists
+    if bool(merge_with) and 'bunchPatternTable' in merge_with:
+        bpt = merge_with['bunchPatternTable']
+        log.debug('Using bpt from merge_with dataset.')
+    elif _mnemonics['bunchPatternTable']['source'] in run.all_sources:
+        bpt = run.get_array(*_mnemonics['bunchPatternTable'].values())
+        log.debug('Loaded bpt from extra_data run.')
+    else:
+        bpt = None
+
+    # align the pulse Id
+    if bpt is not None and len(ds.variables) > 0:
+        dim_names = {'sase3': 'sa3_pId', 'sase1': 'sa1_pId',
+                     'scs_ppl': 'ol_pId'}
+        mask = is_pulse_at(bpt, bunchPattern)
+        mask = mask.rename({'pulse_slot': dim_names[bunchPattern]})
+        ds = ds.rename({_mnemonics['BAM1932M']['dim'][0]:
+                        dim_names[bunchPattern]})
+        ds = ds.where(mask, drop=True)
+
+    # merge with non-BAM dataset
+    ds = mw_ds.merge(ds, join='inner')
+
+    return ds

src/toolbox_scs/detectors/dssc.py

@@ -19,7 +19,7 @@ import joblib
 import numpy as np
 import xarray as xr
 
-from ..load import load_run
+import toolbox_scs as tb
 from ..util.exceptions import ToolBoxValueError, ToolBoxFileError
 from .dssc_data import (
     save_xarray, load_xarray, save_attributes_h5,
@@ -84,7 +84,7 @@ class DSSCBinner:
         self.proposal = proposal_nr
         self.runnr = run_nr
         self.info = load_dssc_info(proposal_nr, run_nr)
-        self.run = load_run(proposal_nr, run_nr)
+        self.run, _ = tb.load(proposal_nr, run_nr)
         self.binners = {}
         for b in binners:
             self.add_binner(b, binners[b])

src/toolbox_scs/detectors/dssc_misc.py

@@ -11,9 +11,8 @@ from imageio import imread
 
 import extra_data as ed
 
-from ..constants import mnemonics as _mnemonics
-from ..misc.bunch_pattern_external import is_sase_3
-
+from .xgm import get_xgm
+from .digitizers import get_tim_peaks
 
 log = logging.getLogger(__name__)
 
@@ -95,15 +94,6 @@ def create_dssc_bins(name, coordinates, bins):
                      'trainId, x, or y')
 
 
-def _bunch_pattern_sase3(run_obj):
-    log.debug('load sase3 bunch pattern info')
-    bpt = run_obj.get_array(*_mnemonics["bunchPatternTable"].values())
-    is_sase3 = is_sase_3(bpt.values[0])
-    ind_sase3 = np.where(is_sase3 == 1)[0]
-    n_pulses = len(ind_sase3)
-    return is_sase3, ind_sase3, n_pulses
-
-
 def get_xgm_formatted(run_obj, xgm_name, dssc_frame_coords):
     """
     Load the xgm data and define coordinates along the pulse dimension.
@@ -123,27 +113,16 @@ def get_xgm_formatted(run_obj, xgm_name, dssc_frame_coords):
         xgm data with coordinate 'pulse'.
     """
     log.debug('load raw xgm data')
-    xgm_raw = run_obj.get_array(*_mnemonics[xgm_name].values())
-
-    log.debug('create formatted xarray.dataArray')
-    _, ind_sase3, n_pulses = _bunch_pattern_sase3(run_obj)
-    coords = {'trainId': xgm_raw.trainId.values}
-    shape = (len(xgm_raw.trainId.values), n_pulses)
-    xgm_formatted = xr.DataArray(np.zeros(shape, dtype=np.float64),
-                                 coords=coords,
-                                 dims=['trainId', 'pulse'])
-    for i in range(len(xgm_formatted.trainId)):
-        xgm_formatted[i, :] = xgm_raw[i, 0:n_pulses].values
-
+    xgm = get_xgm(run_obj, xgm_name)[xgm_name]
+    pulse_dim = [d for d in xgm.dims if d != 'trainId'][0]
+    xgm = xgm.rename({pulse_dim: 'pulse'})
     if type(dssc_frame_coords) == int:
-        dssc_frame_coords = np.linspace(0,
-                                        (n_pulses-1)*dssc_frame_coords,
-                                        n_pulses,
-                                        dtype=np.uint64)
+        dssc_frame_coords = np.arange(xgm.sizes['pulse']*dssc_frame_coords,
+                                      step=dssc_frame_coords, dtype=np.uint64)
 
-    xgm_formatted['pulse'] = dssc_frame_coords
+    xgm['pulse'] = dssc_frame_coords
     log.info('loaded formatted xgm data.')
-    return xgm_formatted
+    return xgm
 
 
 def get_tim_formatted(run_obj, tim_names, dssc_frame_coords):
@@ -161,40 +140,22 @@ def get_tim_formatted(run_obj, tim_names, dssc_frame_coords):
 
     Returns
     -------
-    xgm: xarray.DataArray
-        xgm data with coordinate 'pulse'.
+    tim: xarray.DataArray
+        tim data with coordinate 'pulse'.
     """
     log.debug('load tim data')
-    data = []
-    for name in tim_names:
-        data.append(run_obj.get_array(*_mnemonics[name].values()))
-    tim_raw = xr.zeros_like(data[0])
-    for d in data:
-        tim_raw -= d
-    tim_raw /= len(data)
-
-    log.debug('calculate pulse indices in raw data')
-    _, ind_sase3, n_pulses = _bunch_pattern_sase3(run_obj)
-    ind_pulse_tim = ((np.array(ind_sase3)-ind_sase3[0])/4).astype(np.int64)
-
-    log.debug('create formatted xarray.dataArray')
-    coords = {'trainId': tim_raw.trainId.values}
-    shape = (len(tim_raw.trainId.values), n_pulses)
-    tim_formatted = xr.DataArray(np.zeros(shape, dtype=np.float64),
-                                 coords=coords,
-                                 dims=['trainId', 'pulse'])
-    for i in range(len(tim_formatted.trainId)):
-        tim_formatted[i, :] = np.take(tim_raw[i, :].values, ind_pulse_tim)
-
+    tim = get_tim_peaks(run_obj, tim_names)
+    # average over all tim sources
+    tim = -tim.to_array().mean(dim='variable')
+    pulse_dim = [d for d in tim.dims if d != 'trainId'][0]
+    tim = tim.rename({pulse_dim: 'pulse'})
     if type(dssc_frame_coords) == int:
-        dssc_frame_coords = np.linspace(0,
-                                        (n_pulses-1)*dssc_frame_coords,
-                                        n_pulses,
-                                        dtype=np.uint64)
+        dssc_frame_coords = np.arange(tim.sizes['pulse']*dssc_frame_coords,
+                                      step=dssc_frame_coords, dtype=np.uint64)
 
-    tim_formatted['pulse'] = dssc_frame_coords
+    tim['pulse'] = dssc_frame_coords
     log.info('loaded formatted tim data.')
-    return tim_formatted
+    return tim
 
 
 def quickmask_DSSC_ASIC(poslist):

src/toolbox_scs/detectors/tim.py

-""" Tim related sub-routines
-
-    Copyright (2019) SCS Team.
-    contributions preferrably comply with pep8 code structure 
-    guidelines.
-"""
-
-import logging
-
-import xarray as xr
-
-from .xgm import matchXgmTimPulseId
-from ..constants import mnemonics as _mnemonics
-
-
-
-log = logging.getLogger(__name__)
-
-
-def load_TIM(run, apd='MCP2apd'):
-    """
-    Load TIM traces and match them to SASE3 pulses.
-
-    Parameters
-    ----------
-    run: extra_data.DataCollection, extra_data.RunDirectory
-
-    Returns
-    -------
-    timdata : xarray.DataArray
-        xarray DataArray containing the tim data
-
-    Example
-    -------
-    >>> import toolbox_scs as tb
-    >>> import toolbox_scs.detectors as tbdet
-    >>> run = tb.run_by_proposal(2212, 235)
-    >>> data = tbdet.load_TIM(run)
-    """
-    
-    fields = ["sase1", "sase3", "npulses_sase3", 
-              "npulses_sase1", apd, "SCS_SA3", "nrj"]
-    timdata = xr.Dataset()
-    
-    for f in fields:
-        m = _mnemonics[f]
-        timdata[f] = run.get_array(m['source'],
-                                   m['key'], 
-                                   extra_dims=m['dim'])
-
-    timdata.attrs['run'] = run
-    timdata = matchXgmTimPulseId(timdata)
-    
-    return timdata.rename({'sa3_pId': 'pulse'})[apd]
\ No newline at end of file

src/toolbox_scs/load.py

@@ -12,24 +12,30 @@ import logging
 
 import numpy as np
 import xarray as xr
-
 import extra_data as ed
 
-from .misc.bunch_pattern import extractBunchPattern
-from .constants import mnemonics as _mnemonics_ld
+from .constants import mnemonics as _mnemonics
 from .util.exceptions import ToolBoxValueError, ToolBoxPathError
 from .util.data_access import find_run_dir
+import toolbox_scs.detectors as tbdet
 
 log = logging.getLogger(__name__)
 
 
-def load(fields, runNB, proposalNB,
+def load(proposalNB=None, runNB=None,
+         fields=None,
          subFolder='raw',
          display=False,
          validate=False,
          subset=ed.by_index[:],
          rois={},
-         useBPTable=True):
+         extract_tim=True,
+         extract_laser=True,
+         extract_xgm=True,
+         extract_bam=True,
+         tim_bp='sase3',
+         laser_bp='scs_ppl',
+         ):
     """
     Load a run and extract the data. Output is an xarray with aligned
     trainIds
@@ -37,16 +43,16 @@ def load(fields, runNB, proposalNB,
     Parameters
     ----------
 
-    fields: list of strings, list of dictionaries
-        list of mnemonic strings to load specific data such as "fastccd",
-        "SCS_XGM", or dictionnaries defining a custom mnemonic such as
-            {"extra":
-                {'SCS_CDIFFT_MAG/SUPPLY/CURRENT',
-                 'actual_current.value', None}}
+    proposalNB: (str, int)
+        proposal number e.g. 'p002252' or 2252
     runNB: (str, int)
         run number as integer
-    proposalNB: (str, int)
-        of the proposal number e.g. 'p002252' or 2252
+    fields: str, list of str, list of dict
+        list of mnemonics to load specific data such as "fastccd",
+        "SCS_XGM", or dictionnaries defining a custom mnemonic such as
+        {"extra": {'source: 'SCS_CDIFFT_MAG/SUPPLY/CURRENT',
+                   'key': 'actual_current.value',
+                   'dim': None}}
     subFolder: (str)
         sub-folder from which to load the data. Use 'raw' for raw data
         or 'proc' for processed data.
@@ -60,31 +66,49 @@ def load(fields, runNB, proposalNB,
     rois: dictionary
         a dictionnary of mnemonics with a list of rois definition and
         the desired names, for example:
-            {'fastccd':
-                {'ref':
-                    {'roi': by_index[730:890, 535:720],
-                     'dim': ['ref_x', 'ref_y']},
-                 'sam':
-                    {'roi':by_index[1050:1210, 535:720],
-                     'dim': ['sam_x', 'sam_y']}}}
-    useBPTable: boolean
-        If True, uses the raw bunch pattern table to extract sase pulse number
-        and indices in the trains. If false, load the data from BUNCH_DECODER
-        middle layer device.
+        {'fastccd': {'ref': {'roi': by_index[730:890, 535:720],
+                             'dim': ['ref_x', 'ref_y']},
+                     'sam': {'roi':by_index[1050:1210, 535:720],
+                             'dim': ['sam_x', 'sam_y']}}}
+    extract_tim: bool
+        If True, extracts the peaks from TIM variables (e.g. 'MCP2raw',
+        'MCP3apd') and aligns the pulse Id with the sase3 bunch pattern.
+    extract_laser: bool
+        If True, extracts the peaks from FastADC variables (e.g. 'FastADC5raw',
+        'FastADC3peaks') and aligns the pulse Id with the PP laser bunch pattern.
+    extract_xgm: bool
+        If True, extracts the values from XGM variables (e.g. 'SCS_SA3',
+        'XTD10_XGM') and aligns the pulse Id with the sase1 / sase3 bunch pattern.
+    extract_bam: bool
+        If True, extracts the values from BAM variables (e.g. 'BAM1932M')
+        and aligns the pulse Id with the sase3 bunch pattern.
+    tim_bp: str
+        bunch pattern used to extract the TIM pulses. Ignored if extract_tim=False.
+    laser_bp: str
+        bunch pattern used to extract the TIM pulses. Ignored if extract_tim=False.
 
     Returns
     -------
-    res: xarray.DataArray
-        an xarray DataSet with aligned trainIds
+    run, data: DataCollection, xarray.DataArray
+        extra_data DataCollection of the proposal and run number and an
+        xarray Dataset with aligned trainIds and pulseIds
+
+    Example
+    -------
+    >>> import toolbox_scs as tb
+    >>> run, data = tb.load(2212, 208, ['SCS_SA3', 'MCP2apd', 'nrj'])
+
     """
     try:
         runFolder = find_run_dir(proposalNB, runNB)
     except ToolBoxPathError as err:
         log.error(f"{err.message}")
         raise
-
     run = ed.RunDirectory(runFolder).select_trains(subset)
-
+    if fields is None:
+        return run, xr.Dataset()
+    if isinstance(fields, str):
+        fields = [fields]
     if validate:
         # get_ipython().system('extra-data-validate ' + runFolder)
         pass
@@ -96,25 +120,20 @@ def load(fields, runNB, proposalNB,
     vals = []
 
     # load pulse pattern info
-    if useBPTable:
-        bp_mnemo = _mnemonics_ld['bunchPatternTable']
-        if bp_mnemo['source'] not in run.all_sources:
-            print('Source {} not found in run. Skipping!'.format(
-                                _mnemonics_ld['bunchPatternTable']['source']))
-        else:
-            bp_table = run.get_array(bp_mnemo['source'], bp_mnemo['key'],
-                                     extra_dims=bp_mnemo['dim'])
-            sase1, npulses_sase1, dummy = extractBunchPattern(
-                                              bp_table, 'sase1')
-            sase3, npulses_sase3, dummy = extractBunchPattern(
-                                              bp_table, 'sase3')
-            keys += ["sase1", "npulses_sase1", "sase3", "npulses_sase3"]
-            vals += [sase1, npulses_sase1, sase3, npulses_sase3]
+    if _mnemonics['bunchPatternTable']['source'] in run.all_sources:
+        bpt = run.get_array(*_mnemonics['bunchPatternTable'].values())
+        keys.append("bunchPatternTable")
+        vals.append(bpt)
+    elif _mnemonics['bunchPatternTable_SA3']['source'] in run.all_sources:
+        log.info('Did not find SCS bunch pattern table but found the SA3 one.')
+        bpt = run.get_array(*_mnemonics['bunchPatternTable_SA3'].values())
+        keys.append("bunchPatternTable")
+        vals.append(bpt)
     else:
-        fields += ["sase1", "sase3", "npulses_sase3", "npulses_sase1"]
-
+        log.warning('Source {} and {} not found in run. Skipping!'.format(
+                            _mnemonics['bunchPatternTable']['source'],
+                            _mnemonics['bunchPatternTable_SA3']['source'],))
     for f in fields:
-
         if type(f) == dict:
             # extracting mnemomic defined on the spot
             if len(f.keys()) > 1:
@@ -124,73 +143,73 @@ def load(fields, runNB, proposalNB,
             v = f[k]
         else:
             # extracting mnemomic from the table
-            if f in _mnemonics_ld:
-                v = _mnemonics_ld[f]
+            if f in _mnemonics:
+                v = _mnemonics[f]
                 k = f
             else:
-                print('Unknow mnemonic "{}". Skipping!'.format(f))
+                print(f'Unknow mnemonic "{f}". Skipping!')
                 continue
-
         if k in keys:
             continue  # already loaded, skip
-
         if display:
-            print('Loading {}'.format(k))
-
+            print(f'Loading {k}')
         if v['source'] not in run.all_sources:
-            print('Source {} not found in run. Skipping!'.format(v['source']))
+            log.warning(f'Source {v["source"]} not found in run. Skipping!')
+            print(f'Source {v["source"]} not found in run. Skipping!')
             continue
-
         if k not in rois:
             # no ROIs selection, we read everything
-            vals.append(run.get_array(v['source'], v['key'],
-                                      extra_dims=v['dim']))
+            arr = run.get_array(*v.values())
+            if len(arr) == 0:
+                log.warning(f'Empty array for {f}: {v["source"]}, {v["key"]}. '
+                            'Skipping!')
+                print(f'Empty array for {f}: {v["source"]}, {v["key"]}. '
+                      'Skipping!')
+                continue
+            vals.append(arr)
             keys.append(k)
         else:
             # ROIs selection, for each ROI we select a region of the data and
             # save it with new name and dimensions
             for nk, nv in rois[k].items():
-                vals.append(run.get_array(v['source'], v['key'],
-                                          extra_dims=nv['dim'],
-                                          roi=nv['roi']))
+                arr = run.get_array(v['source'], v['key'],
+                                    extra_dims=nv['dim'],
+                                    roi=nv['roi'])
+                if len(arr) == 0:
+                    log.warning(f'Empty array for {f}: {v["source"]}, {v["key"]}. '
+                                'Skipping!')
+                    print(f'Empty array for {f}: {v["source"]}, {v["key"]}. '
+                          'Skipping!')
+                    continue
+                vals.append(arr)
                 keys.append(nk)
-
     aligned_vals = xr.align(*vals, join='inner')
-    result = dict(zip(keys, aligned_vals))
-    result = xr.Dataset(result)
-    result.attrs['run'] = run
-    result.attrs['runFolder'] = runFolder
-    return result
+    data = dict(zip(keys, aligned_vals))
+    data = xr.Dataset(data)
+    data.attrs['runFolder'] = runFolder
 
+    tim = [k for k in _mnemonics if 'MCP' in k and k in data]
+    if extract_tim and len(tim) > 0:
+        data = tbdet.get_tim_peaks(run, mnemonics=tim, merge_with=data,
+                                   bunchPattern=tim_bp)
 
-def load_run(proposal, run, **kwargs):
-    """
-    Get run in given proposal
+    laser = [k for k in _mnemonics if 'FastADC' in k and k in data]
+    if extract_laser and len(laser) > 0:
+        data = tbdet.get_laser_peaks(run, mnemonics=laser, merge_with=data,
+                                     bunchPattern=laser_bp)
 
-    Wraps the extra_data open_run routine, out of convenience for the toolbox
-    user. More information can be found in the karabo_data documentation.
+    xgm = ['XTD10_XGM', 'XTD10_XGM_sigma', 'XTD10_SA3', 'XTD10_SA3_sigma', 
+           'XTD10_SA1', 'XTD10_SA1_sigma', 'SCS_XGM', 'SCS_XGM_sigma', 
+           'SCS_SA1', 'SCS_SA1_sigma', 'SCS_SA3', 'SCS_SA3_sigma']
+    xgm = [k for k in xgm if k in data]
+    if extract_xgm and len(xgm) > 0:
+        data = tbdet.get_xgm(run, mnemonics=xgm, merge_with=data)
 
-    Parameters
-    ----------
-    proposal: str, int
-        Proposal number
-    run: str, int
-        Run number
-
-    **kwargs
-    --------
-    data: str
-        default -> 'raw'
-    include: str
-        default -> '*'
+    bam = [k for k in _mnemonics if 'BAM' in k and k in data]
+    if extract_bam and len(bam) > 0:
+        data = tbdet.get_bam(run, mnemonics=bam, merge_with=data)
 
-    Returns
-    -------
-    run : extra_data.DataCollection
-        DataCollection object containing information about the specified
-        run. Data can be loaded using built-in class methods.
-    """
-    return ed.open_run(proposal, run, **kwargs)
+    return run, data
 
 
 def run_by_path(path):
@@ -280,8 +299,8 @@ def get_array(run, mnemonic_key=None, stepsize=None):
             data = xr.DataArray(
                         np.ones(len(run.train_ids), dtype=np.int16),
                         dims=['trainId'], coords={'trainId': run.train_ids})
-        elif mnemonic_key in _mnemonics_ld:
-            mnem = _mnemonics_ld[mnemonic_key]
+        elif mnemonic_key in _mnemonics:
+            mnem = _mnemonics[mnemonic_key]
             data = run.get_array(*mnem.values())
         else:
             raise ToolBoxValueError("Invalid mnemonic", mnemonic_key)
@@ -295,3 +314,79 @@ def get_array(run, mnemonic_key=None, stepsize=None):
         raise
 
     return data
+
+
+def mnemonics_to_process(mnemo_list, merge_with, detector, func=None):
+    """
+    Finds the list of mnemonics, within mnemo_list and merge_with, that
+    correspond to arrays that are not yet loaded and/or processed by a
+    detector function. Removes the mnemonics of the already processed
+    arrays from the list.
+
+    Parameters
+    ----------
+    mnemo_list: str or list of str
+        ToolBox mnemonics of pulse-resolved detector arrays
+    merge_with: xarray Dataset
+        Dataset that may contain non-processed arrays
+    detector: str
+        One in {'ADQ412', 'FastADC', 'XGM', 'BAM'}
+    func: function
+        function that takes one argument, an unprocessed mnemonic string,
+        and converts it into a processed one, i.e. from 'MCP2apd' to
+        'MCP2peaks'. If None, the function returns the input mnemonic.
+
+    Returns
+    -------
+    mnemonics: list of str
+        the mnemonics to process
+    """
+    if func is None:
+        def func(x):
+            return x
+
+    if detector == 'BAM':
+        det_mnemos = [m for m in _mnemonics if 'BAM' in m]
+        default_mnemo = 'BAM1932M'
+        default_processed = 'BAM1932M'
+    if detector == 'XGM':
+        det_mnemos = ['XTD10_XGM', 'XTD10_XGM_sigma', 'XTD10_SA3',
+                      'XTD10_SA3_sigma', 'XTD10_SA1', 'XTD10_SA1_sigma',
+                      'SCS_XGM', 'SCS_XGM_sigma', 'SCS_SA1', 'SCS_SA1_sigma',
+                      'SCS_SA3', 'SCS_SA3_sigma']
+        default_mnemo = 'SCS_SA3'
+        default_processed = 'SCS_SA3'
+    if detector == 'ADQ412':
+        det_mnemos = [m for m in _mnemonics if 'MCP' in m]
+        default_mnemo = 'MCP2apd'
+        default_processed = 'MCP2peaks'
+    if detector == 'FastADC':
+        det_mnemos = [m for m in _mnemonics if 'FastADC' in m]
+        default_mnemo = 'FastADC5raw'
+        default_processed = 'FastADC5peaks'
+
+    dig_dims = list(set([_mnemonics[m]['dim'][0] for m in det_mnemos]))
+    processed_mnemos = list(set([func(m) for m in det_mnemos]))
+
+    # create a list of mnemonics to process from the provided mnemonics and
+    # merge_with Dataset
+    mw_mnemos = []
+    mw_processed = []
+    if bool(merge_with):
+        mw_mnemos = [m for m in merge_with if m in det_mnemos and
+                     any(dim in merge_with[m].dims for dim in dig_dims)]
+        mw_processed = [m for m in merge_with if m in processed_mnemos and
+                        any(dim in merge_with[m].dims for dim in dig_dims)
+                        is False]
+    if mnemo_list is None:
+        mnemonics = []
+        if len(mw_mnemos) == 0 and default_processed not in mw_processed:
+            mnemonics = [default_mnemo]
+    else:
+        mnemonics = [mnemo_list] if isinstance(mnemo_list, str) else mnemo_list
+    mnemonics = list(set(mnemonics + mw_mnemos))
+    for m in mnemonics[:]:
+        if func(m) in mw_processed:
+            mnemonics.remove(m)
+
+    return mnemonics