Adds channel_peak_params to read digitizer parameters for peak integration

src/toolbox_scs/detectors/__init__.py

@@ -4,7 +4,8 @@ from .xgm import (
 from .tim import (
     load_TIM,)
 from .digitizers import(
-    get_peaks, get_tim_peaks, get_fast_adc_peaks, find_integ_params)
+    get_peaks, get_tim_peaks, get_fast_adc_peaks, find_integ_params,
+    channel_peak_params)
 from .dssc_data import (
     save_xarray, load_xarray, get_data_formatted, save_attributes_h5)
 from .dssc_misc import (
@@ -26,6 +27,7 @@ __all__ = (
     "get_tim_peaks",
     "get_fast_adc_peaks",
     "find_integ_params",
+    "channel_peak_params",
     "load_dssc_info",
     "create_dssc_bins",
     "calc_xgm_frame_indices",

src/toolbox_scs/detectors/digitizers.py

@@ -121,10 +121,10 @@ def get_peaks(run,
         required if data is a DataArray or None.
     key: str
         Key for digitizer data, e.g. 'digitizers.channel_1_A.raw.samples'.
-        Only required if data a DataArray or is None.
+        Only required if data is a DataArray or is None.
     digitizer: string
         name of digitizer, e.g. 'FastADC' or 'ADQ412'. Used to determine
-        the sampling rate when useRaw is True.
+        the sampling rate.
     useRaw: bool
         If True, extract peaks from raw traces. If False, uses the APD (or
         peaks) data from the digitizer.
@@ -157,13 +157,16 @@ def get_peaks(run,
         raise ValueError('At least data or source + key arguments ' +
                          'are required.')
     # load data
+    arr = None
     if data is None:
         arr = run.get_array(source, key)
-    elif isinstance(data, str):
+    if isinstance(data, str):
         arr = run.get_array(*_mnemonics[data].values())
         source = _mnemonics[data]['source']
         key = _mnemonics[data]['key']
-    else:
+    if arr is None:
+        if source is None or key is None:
+            raise ValueError('source and/or key arguments missing.')
         arr = data
     dim = [d for d in arr.dims if d != 'trainId'][0]
     # check if bunch pattern is provided
@@ -192,7 +195,7 @@ def get_peaks(run,
                                drop=True).trainId
         mask_on = mask.sel(trainId=valid_tid)
         if (mask_on == mask_on[0]).all().values is False:
-            log.debug('Pattern changed during the run!')
+            log.info('Pattern changed during the run!')
         pid = np.unique(np.where(mask_on)[1])
         npulses = len(pid)
         extra_coords = pid
@@ -202,6 +205,13 @@ def get_peaks(run,
     if extra_dim is None:
         extra_dim = 'pulseId'
 
+    # minimum pulse period, according to digitizer type
+    min_distance = 1
+    if digitizer == 'FastADC':
+        min_distance = 24
+    if digitizer == 'ADQ412':
+        min_distance = 440
+
     # 1. Use peak-integrated data from digitizer
     if useRaw is False:
         # 1.1 No bunch pattern provided
@@ -217,41 +227,26 @@ def get_peaks(run,
             return arr.sel({dim: 0}).expand_dims(extra_dim).T.assign_coords(
                 {extra_dim: extra_coords})
         # verify period used by APD and match it to pid from bunchPattern
-        if digitizer == 'FastADC':
-            adc_source = source.split(':')[0]
-            enable_key = (source.split(':')[1].split('.')[0]
-                          + '.enablePeakComputation.value')
-            if run.get_array(adc_source, enable_key)[0] is False:
-                raise ValueError('The digitizer did not record ' +
-                                 'peak-integrated data.')
-            period_key = (source.split(':')[1].split('.')[0] +
-                          '.pulsePeriod.value')
-            period = run.get_array(adc_source, period_key)[0].values/24
-        if digitizer == 'ADQ412':
-            board_source = source.split(':')[0]
-            board = key[19]
-            channel = key[21]
-            channel_to_number = {'A': 0, 'B': 1, 'C': 2, 'D': 3}
-            channel = channel_to_number[channel]
-            in_del_key = (f'board{board}.apd.channel_{channel}' +
-                          '.initialDelay.value')
-            initialDelay = run.get_array(board_source, in_del_key)[0].values
-            up_lim_key = (f'board{board}.apd.channel_{channel}' +
-                          '.upperLimit.value')
-            upperLimit = run.get_array(board_source, up_lim_key)[0].values
-            period = (upperLimit - initialDelay)/440
-        stride = (pid[1] - pid[0]) / period
-        if period < 1:
-            log.warning('Warning: the pulse period in digitizer was ' +
-                        'smaller than the separation of two pulses ' +
-                        'at 4.5 MHz.')
-            stride = 1
+        peak_params = channel_peak_params(run, source, key)
+        log.debug(f'peak_params: {peak_params}')
+        if peak_params['enable'] == 0:
+            raise ValueError('The digitizer did not record ' +
+                             'peak-integrated data.')
+        period = peak_params['period']
+        if period < min_distance:
+            log.warning(f'Warning: the pulse period ({period} samples) of '
+                        'digitizer was smaller than the pulse separation at ' +
+                        f'4.5 MHz ({min_distance} samples).')
+        step = period / min_distance
+        pid_diff = np.min(np.diff(pid))
+        stride = pid_diff / step
+        log.debug(f'dig step: {step}, pid diff: {pid_diff}, ' +
+                  f'stride: {stride}')
         if stride < 1:
-            raise ValueError('Pulse period in digitizer was too large ' +
-                             'compared to actual pulse separation. Some ' +
-                             'pulses were not recorded.')
+            raise ValueError(f'Pulse pattern in digitizer (1 out of {step} ' +
+                             'pulses @ 4.5 MHz) does not match the actual ' +
+                             'bunch pattern. Some pulses were not recorded.')
         stride = int(stride)
-        log.debug(f'period {period}, stride {stride}')
         if npulses*stride > arr.sizes[dim]:
             raise ValueError('The number of pulses recorded by digitizer ' +
                              f'that correspond to actual {pattern} pulses ' +
@@ -261,12 +256,6 @@ def get_peaks(run,
         return peaks.assign_coords({extra_dim: extra_coords})
 
     # 2. Use raw data from digitizer
-    # minimum samples between two pulses, according to digitizer type
-    min_distance = 1
-    if digitizer == 'FastADC':
-        min_distance = 24
-    if digitizer == 'ADQ412':
-        min_distance = 440
     if autoFind:
         integParams = find_integ_params(arr.mean(dim='trainId'),
                                         min_distance=min_distance)
@@ -311,6 +300,94 @@ def get_peaks(run,
     return peaks.assign_coords({extra_dim: extra_coords})
 
 
+def channel_peak_params(run, source, key=None, digitizer=None,
+                        channel=None, board=None):
+    """
+    Extract peak-integration parameters used by a channel of the digitizer.
+
+    Parameters
+    ----------
+    run: extra_data.DataCollection
+        DataCollection containing the digitizer data
+    source: str
+        ToolBox mnemonic of a digitizer data, e.g. "MCP2apd" or
+        "FastADC4peaks", or name of digitizer source, e.g.
+        'SCS_UTC1_ADQ/ADC/1:network'.
+    key: str
+        optional, used in combination of source (if source is not a ToolBox
+        mnemonics) instead of digitizer, channel and board.
+    digitizer: {"FastADC", "ADQ412"} str
+        Type of digitizer. If None, inferred from the source mnemonic.
+    channel: int or str
+        The digitizer channel for which to retrieve the parameters. If None,
+        inferred from the source mnemonic.
+    board: int
+        Board of the ADQ412. If None, inferred from the source mnemonic.
+
+    Returns
+    -------
+    dict with peak integration parameters.
+    """
+    if source in _mnemonics:
+        m = _mnemonics[source]
+        source = m['source']
+        key = m['key']
+    if key is not None:
+        if 'network' in source:
+            digitizer = 'ADQ412'
+            ch_to_int = {'A': 0, 'B': 1, 'C': 2, 'D': 3}
+            k = key.split('.')[1].split('_')
+            channel = ch_to_int[k[2]]
+            board = k[1]
+        else:
+            digitizer = 'FastADC'
+            channel = int(source.split(':')[1].split('.')[0].split('_')[1])
+    if digitizer is None:
+        raise ValueError('digitizer argument is missing.')
+    if channel is None:
+        raise ValueError('channel argument is missing.')
+    if isinstance(channel, str):
+        ch_to_int = {'A': 0, 'B': 1, 'C': 2, 'D': 3}
+        channel = ch_to_int[channel]
+    if board is None and digitizer == 'ADQ412':
+        raise ValueError('board argument is missing.')
+    keys = None
+    if digitizer == 'ADQ412':
+        baseKey = f'board{board}.apd.channel_{channel}.'
+        keys = ['baseStart', 'baseStop', 'pulseStart',
+                'pulseStop', 'initialDelay', 'upperLimit',
+                'enable']
+        keys = {k: baseKey + k + '.value' for k in keys}
+        keys['numPulses'] = f'board{board}.apd.numberOfPulses.value'
+    if digitizer == 'FastADC':
+        if ":" in source:
+            baseKey = source.split(':')[1].split('.')[0]+'.'
+        else:
+            baseKey = f'channel_{channel}.'
+        keys = ['baseStart', 'baseStop', 'initialDelay',
+                'peakSamples', 'numPulses', 'pulsePeriod',
+                'enablePeakComputation']
+        keys = {k: baseKey + k + '.value' for k in keys}
+    if ":" in source:
+        source = source.split(':')[0]
+    log.debug(f'retrieving data from source {source}, keys = ' +
+              f'{list(keys.values())}')
+    tid, data = run.select(source).train_from_index(0)
+    result = [data[source][k] for k in keys.values()]
+    result = dict(zip(keys.keys(), result))
+    if digitizer == 'ADQ412':
+        result['period'] = result.pop('upperLimit') - \
+                           result.pop('initialDelay')
+    if digitizer == 'FastADC':
+        result['period'] = result.pop('pulsePeriod')
+        result['pulseStart'] = result['initialDelay']
+        result['pulseStop'] = result.pop('initialDelay') + \
+                              result.pop('peakSamples')
+        result['enable'] = result.pop('enablePeakComputation')
+
+    return result
+
+
 def find_integ_params(trace, min_distance=1, height=None, width=1):
     """
     find integration parameters necessary for peak integration of a raw
@@ -524,7 +601,7 @@ def get_fast_adc_peaks(data=None, run=None, bunchPattern='scs_ppl',
         peaks = get_peaks(run, d,
                           source=mnemonic['source'],
                           key=mnemonic['key'],
-                          digitizer='ADQ412',
+                          digitizer='FastADC',
                           useRaw=useRaw,
                           autoFind=autoFind,
                           integParams=integParams,