diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f60120ca5..480f2dad4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -24,12 +24,6 @@ repos: - id: requirements-txt-fixer - id: trailing-whitespace -- repo: https://github.com/Carreau/velin.git - rev: "0.0.12" - hooks: - - id: velin - args: ["--check", "--write"] - - repo: https://github.com/asottile/setup-cfg-fmt rev: "v1.20.1" hooks: @@ -41,7 +35,7 @@ repos: - id: isort - repo: https://github.com/asottile/pyupgrade - rev: "v2.34.0" + rev: "v2.37.1" hooks: - id: pyupgrade args: ["--py36-plus"] @@ -87,3 +81,9 @@ repos: - id: rst-backticks - id: rst-directive-colons - id: rst-inline-touching-normal + +- repo: https://github.com/pre-commit/mirrors-prettier + rev: "v2.7.1" + hooks: + - id: prettier + types_or: [json] diff --git a/docs/source/conf.py b/docs/source/conf.py index c115832a9..a83ae5ca0 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -64,7 +64,9 @@ # enforce consistent usage of NumPy-style docstrings napoleon_numpy_docstring = True napoleon_google_docstring = False -napoleon_use_rtype = False +napoleon_use_ivar = True +napoleon_custom_sections = ["JSON Configuration Example"] + # intersphinx intersphinx_mapping = { diff --git a/docs/source/manuals/build_dsp.rst b/docs/source/manuals/build_dsp.rst index 7ba6b8b46..8830594b7 100644 --- a/docs/source/manuals/build_dsp.rst +++ b/docs/source/manuals/build_dsp.rst @@ -39,6 +39,8 @@ Writing custom processors # 1) Import Python modules # + from __future__ import annotations + import numpy as np from numba import guvectorize @@ -67,14 +69,46 @@ Writing custom processors # - Use "w_" for waveforms, "t_" for indexes, "a_" for amplitudes # - Use underscore casing for the name of the processor and variables, e.g., # "a_trap_energy_in" or "t_trigger_in" + # - use type annotations # - def the_processor_template(w_in, t_in, a_in, w_out, t_out): + def the_processor_template(w_in: np.ndarray, + t_in: float, + a_in: float, + w_out: np.ndarray, + t_out: float) -> None: # 4) Document the algorithm # - """ - Add here a complete description of what the processor does, including the - meaning of input and output variables. + """One-liner description of the processor. + + Add here a more detailed description of what the processor does. + Document input parameters in the "Parameters" section. Add a JSON + example for ProcessingChain configuration in the last section. + + Parameters + ---------- + w_in + the input waveform. + t_in + a scalar parameter in the time domain + a_in + a scalar parameter in the amplitude domain + w_out + the output waveform + t_out + an output scalar value in the time domain + + JSON Configuration Example + -------------------------- + + .. code-block :: json + + "wf_bl": { + "function": "the_processor_template", + "module": "pygama.dsp.processors", + "args": ["waveform", "t_a", "a_b", "wf_filtered", "t_result"], + "unit": "ADC" + } """ # 5) Initialize output parameters diff --git a/src/pygama/dsp/__init__.py b/src/pygama/dsp/__init__.py index 8988108e8..540702519 100644 --- a/src/pygama/dsp/__init__.py +++ b/src/pygama/dsp/__init__.py @@ -1,14 +1,13 @@ -""" +r""" The pygama signal processing framework is contained in the :mod:`.dsp` submodule. This code is responsible for running a variety of discrete signal processors on data, producing tier 2 (dsp) files from tier 1 (raw) files. The main contents of this submodule are: * :mod:`.processors`: A collection of Numba functions that perform individual - DSP transforms and reductions on our data. Individual processors are held in - the ``_processors`` directory, which is where contributors will write new - processors for inclusion in the analysis. Available processors include all - numpy ufuncs as well. + DSP transforms and reductions on our data. Here is where contributors will + write new processors for inclusion in the analysis. Available processors + include all :class:`numpy.ufunc`\ s as well. * :class:`.ProcessingChain`: A class that manages and efficiently runs a list of DSP processors * :func:`.build_processing_chain`: A function that builds a :class:`.ProcessingChain` diff --git a/src/pygama/dsp/_processors/__init__.py b/src/pygama/dsp/_processors/__init__.py deleted file mode 100644 index 37e898bfb..000000000 --- a/src/pygama/dsp/_processors/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -""" -Package documentation -""" diff --git a/src/pygama/dsp/_processors/fftw.py b/src/pygama/dsp/_processors/fftw.py deleted file mode 100644 index c851bb335..000000000 --- a/src/pygama/dsp/_processors/fftw.py +++ /dev/null @@ -1,123 +0,0 @@ -import numpy as np -from numba import guvectorize -from pyfftw import FFTW - - -def dft(buf_in, buf_out): - """ - Perform discrete Fourier transforms using the FFTW library. FFTW optimizes - the FFT algorithm based on the size of the arrays, with SIMD parallelized - commands. This optimization requires initialization, so this is a factory - function that returns a Numba gufunc that performs the FFT. FFTW works on - fixed memory buffers, so you must tell it what memory to use ahead of time. - When using this with ProcessingChain, to ensure the correct buffers are used, - call ProcessingChain.get_variable('var_name') to give it the internal memory - buffer directly---with build_dsp, you can just give it the name, and it will - automatically happen. The possible dtypes for the input/outputs are: - - float32/float (size n) -> complex64 (size n/2+1) - - float64/double (size n) -> complex128 (size n/2+1) - - float128/longdouble (size n) -> complex256/clongdouble (size n/2+1) - - complex64 (size n) -> complex64 (size n ) - - complex128 (size n) -> complex128 (size n ) - - complex256/clongdouble (size n) -> complex256/clongdouble (size n ) - """ - try: - dft_fun = FFTW(buf_in, buf_out, axes=(-1,), direction='FFTW_FORWARD') - except ValueError: - raise ValueError("""Incompatible array types/shapes. Allowed: - - float32/float (size n) -> complex64 (size n/2+1) - - float64/double (size n) -> complex128 (size n/2+1) - - float128/longdouble (size n) -> complex256/clongdouble (size n/2+1) - - complex64 (size n) -> complex64 (size n) - - complex128 (size n) -> complex128 (size n) - - complex256/clongdouble (size n) -> complex256/clongdouble (size n)""") - - typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' - sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' - - @guvectorize([typesig], sizesig, forceobj=True) - def dft(wf_in, dft_out): - dft_fun(wf_in, dft_out) - - return dft - -def inv_dft(buf_in, buf_out): - """ - Perform inverse discrete Fourier transforms using the FFTW library. FFTW - optimizes the FFT algorithm based on the size of the arrays, with SIMD parallelized - commands. This optimization requires initialization, so this is a factory - function that returns a Numba gufunc that performs the FFT. FFTW works on - fixed memory buffers, so you must tell it what memory to use ahead of time. - When using this with ProcessingChain, to ensure the correct buffers are used, - call ProcessingChain.get_variable('var_name') to give it the internal memory - buffer directly---with build_dsp, you can just give it the name, and it will - automatically happen. The possible dtypes for the input/outputs are: - - complex64 (size n/2+1) -> float32/float (size n) - - complex128 (size n/2+1) -> float64/double (size n) - - complex256/clongdouble (size n/2+1) -> float128/longdouble (size n) - - complex64 (size n ) -> complex64 (size n) - - complex128 (size n ) -> complex128 (size n) - - complex256/clongdouble (size n ) -> complex256/clongdouble (size n) - """ - try: - idft_fun = FFTW(buf_in, buf_out, axes=(-1,), direction='FFTW_BACKWARD') - except ValueError: - raise ValueError("""Incompatible array types/shapes. Allowed: - - complex64 (size n/2+1) -> float32/float (size n) - - complex128 (size n/2+1) -> float64/double (size n) - - complex256/clongdouble (size n/2+1) -> float128/longdouble (size n) - - complex64 (size n) -> complex64 (size n) - - complex128 (size n) -> complex128 (size n) - - complex256/clongdouble (size n) -> complex256/clongdouble (size n)""") - - typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' - sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' - - @guvectorize([typesig], sizesig, forceobj=True) - def inv_dft(wf_in, dft_out): - idft_fun(wf_in, dft_out) - - return inv_dft - -def psd(buf_in, buf_out): - """ - Perform discrete Fourier transforms using the FFTW library, and use it to get - the power spectral density. FFTW optimizes the FFT algorithm based on the - size of the arrays, with SIMD parallelized commands. This optimization - requires initialization, so this is a factory function that returns a Numba gufunc - that performs the FFT. FFTW works on fixed memory buffers, so you must - tell it what memory to use ahead of time. When using this with ProcessingChain, - to ensure the correct buffers are used, call ProcessingChain.get_variable('var_name') - to give it the internal memory buffer directly---with build_dsp, you can just - give it the name, and it will automatically happen. The possible dtypes for the - input/outputs are: - - complex64 (size n) -> float32/float (size n ) - - complex128 (size n) -> float64/double (size n ) - - complex256/clongdouble (size n) -> float128/longdouble (size n ) - - float32/float (size n) -> float32/float (size n/2+1) - - float64/double (size n) -> float64/double (size n/2+1) - - float128/longdouble (size n) -> float128/longdouble (size n/2+1) - """ - - # build intermediate array for the dft, which will be abs'd to get the PSD - buf_dft = np.ndarray(buf_out.shape, np.dtype('complex' + str(buf_out.dtype.itemsize * 16))) - try: - dft_fun = FFTW(buf_in, buf_dft, axes=(-1,), direction='FFTW_FORWARD') - except ValueError: - raise ValueError("""Incompatible array types/shapes. Allowed: - - complex64 (size n) -> float32/float (size n) - - complex128 (size n) -> float64/double (size n) - - complex256/clongdouble (size n) -> float128/longdouble (size n) - - float32/float (size n) -> float32/float (size n/2+1) - - float64/double (size n) -> float64/double (size n/2+1) - - float128/longdouble (size n) -> float128/longdouble (size n/2+1)""") - - typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' - sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' - - @guvectorize([typesig], sizesig, forceobj=True) - def psd(wf_in, psd_out): - dft_fun(wf_in, buf_dft) - np.abs(buf_dft, psd_out) - - return psd diff --git a/src/pygama/dsp/_processors/param_lookup.py b/src/pygama/dsp/_processors/param_lookup.py deleted file mode 100644 index 47e767edc..000000000 --- a/src/pygama/dsp/_processors/param_lookup.py +++ /dev/null @@ -1,26 +0,0 @@ -import numpy as np -from numba import from_dtype, guvectorize, types - - -def param_lookup(param_dict, default_val, dtype): - """ - Generate the ufunc lookup(channel, val), which returns a NumPy array of - values corresponding to various channels that are looked up in the provided - param_dict. If there is no key, use default_val instead. - """ - out_type = from_dtype(np.dtype(dtype)) - - # convert types to avoid any necessity of casting - param_dict = {types.uint32(k): out_type(v) for k, v in param_dict.items()} - default_val = out_type(default_val) - - @guvectorize(["void(uint32, " + out_type.name + "[:])"], - "()->()", forceobj=True) - def lookup(channel, val): - """ - Look up a value for the provided channel from a dictionary provided at - compile time. - """ - val[0] = param_dict.get(channel, default_val) - - return lookup diff --git a/src/pygama/dsp/_processors/pulse_injector.py b/src/pygama/dsp/_processors/pulse_injector.py deleted file mode 100644 index 5959d9e54..000000000 --- a/src/pygama/dsp/_processors/pulse_injector.py +++ /dev/null @@ -1,47 +0,0 @@ -from math import exp, log - -import numpy as np -from numba import guvectorize - - -@guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])", - "void(float64[:], float64, float64,float64,float64,float64[:])"], - "(n),(),(),(), ()->(n)", nopython=True, cache=True) - -def inject_sig_pulse(wf_in,t0,rt,A,decay, wf_out): - """ - Inject sigmoid pulse into existing waveform to simulate pileup - """ - - - wf_out[:] = np.nan - - if np.isnan(wf_in).any() or np.isnan(rt) or np.isnan(t0) or np.isnan(A) or np.isnan(decay): - return - - rise = 4*log(99)/rt - - wf_out[:] = wf_in[:] - for T in range(len(wf_out)): - wf_out[T] = wf_out[T]+ (A/(1+exp((-rise)*(T-(t0+rt/2))))*exp(-(1/decay)*(T-t0))) - -@guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])", - "void(float64[:], float64, float64,float64,float64,float64[:])"], - "(n),(),(),(), ()->(n)", nopython=True, cache=True) - -def inject_exp_pulse(wf_in,t0,rt,A,decay, wf_out): - """ - Inject exponential pulse into existing waveform to simulate pileup - """ - - wf_out[:] = np.nan - - if np.isnan(wf_in).any() or np.isnan(rt) or np.isnan(t0) or np.isnan(A) or np.isnan(decay): - return - - wf_out[:] = wf_in[:] - for T in range(len(wf_out)): - if (t0<= T)& (T<= t0+rt): - wf_out[T] += (A*exp((T-t0-rt)/(rt))*exp(-(1/decay)*(T-t0))) - elif (T>t0+rt): - wf_out[T] += (A*exp(-(1/decay)*(T-t0))) diff --git a/src/pygama/dsp/build_dsp.py b/src/pygama/dsp/build_dsp.py index 85a8a708b..0640930a4 100644 --- a/src/pygama/dsp/build_dsp.py +++ b/src/pygama/dsp/build_dsp.py @@ -23,47 +23,53 @@ log = logging.getLogger(__name__) -def build_dsp(f_raw: str, f_dsp: str, dsp_config: str | dict = None, - lh5_tables: list[str] = None, database: str = None, - outputs: list[str] = None, n_max: int = np.inf, - write_mode: str = 'r', buffer_len: int = 3200, - block_width: int = 16, chan_config: dict = None) -> None: - """ - Convert raw-tier LH5 data into dsp-tier LH5 data by running a sequence of - processors via the :class:`~.processing_chain.ProcessingChain`. +def build_dsp(f_raw: str, + f_dsp: str, + dsp_config: str | dict = None, + lh5_tables: list[str] = None, + database: str | dict = None, + outputs: list[str] = None, + n_max: int = np.inf, + write_mode: str = None, + buffer_len: int = 3200, + block_width: int = 16, + chan_config: dict[str, str] = None) -> None: + """Convert raw-tier LH5 data into dsp-tier LH5 data by running a sequence + of processors via the :class:`~.processing_chain.ProcessingChain`. Parameters ---------- - f_raw : str - name of raw LH5 file to read from - f_dsp : str - name of dsp LH5 file to write to - dsp_config : str or dict - dict or name of JSON file containing + f_raw + name of raw-tier LH5 file to read from. + f_dsp + name of dsp-tier LH5 file to write to. + dsp_config + :class:`dict` or name of JSON file containing :class:`~.processing_chain.ProcessingChain` config. See - :func:`~.processing_chain.build_processing_chain` for details - lh5_tables : list of str, optional - list of HDF5 groups to consider in the input file. If None, process all - valid groups - database : str, optional - name of JSON file containing a parameter database. See - :func:`~.processing_chain.build_processing_chain` for details - outputs : list of str, optional + :func:`~.processing_chain.build_processing_chain` for details. + lh5_tables + list of HDF5 groups to consider in the input file. If ``None``, process + all valid groups. + database + dictionary or name of JSON file containing a parameter database. See + :func:`~.processing_chain.build_processing_chain` for details. + outputs list of parameter names to write to the output file. If not provided, - use list provided under ``"outputs"`` in the DSP configuration file - n_max : int, optional - Number of waveforms to process. Default all - write_mode : {'r', 'a', 'u'}, optional + use list provided under ``"outputs"`` in the DSP configuration file. + n_max + number of waveforms to process. + write_mode + - ``None`` -- create new output file if it does not exist - `'r'` -- delete existing output file with same name before writing - `'a'` -- append to end of existing output file - `'u'` -- update values in existing output file - buffer_len : int, optional, default=3200 - number of waveforms to read/write from disk at a time - block_width : int, optional, default=16 - number of waveforms to process at a time - chan_config : dict, optional + buffer_len + number of waveforms to read/write from/to disk at a time. + block_width + number of waveforms to process at a time. + chan_config contains JSON DSP configuration file names for every table in - `lh5_tables` + `lh5_tables`. """ if chan_config is not None: @@ -124,6 +130,9 @@ def build_dsp(f_raw: str, f_dsp: str, dsp_config: str | dict = None, database = None raise ValueError('input database is not a valid JSON file or dict') + if write_mode is None and os.path.isfile(f_dsp): + raise FileExistsError(f"output file {f_dsp} exists. Set the 'write_mode' keyword") + # clear existing output files if write_mode == 'r': if os.path.isfile(f_dsp): diff --git a/src/pygama/dsp/errors.py b/src/pygama/dsp/errors.py index 457c78287..4d18c24f1 100644 --- a/src/pygama/dsp/errors.py +++ b/src/pygama/dsp/errors.py @@ -1,3 +1,6 @@ +from __future__ import annotations + + class DSPError(Exception): """Base class for signal processors.""" pass diff --git a/src/pygama/dsp/processing_chain.py b/src/pygama/dsp/processing_chain.py index ed69117a2..86ee4b7d0 100644 --- a/src/pygama/dsp/processing_chain.py +++ b/src/pygama/dsp/processing_chain.py @@ -13,7 +13,7 @@ from abc import ABCMeta, abstractmethod from copy import deepcopy from dataclasses import dataclass -from typing import Any +from typing import Any, Union import numpy as np from numba import vectorize @@ -25,6 +25,8 @@ log = logging.getLogger(__name__) +LGDO = Union[lgdo.Scalar, lgdo.Array, lgdo.VectorOfVectors, lgdo.Struct] + # Filler value for variables to be automatically deduced later auto = 'auto' @@ -99,25 +101,26 @@ class ProcChainVar: Members can be set to ``auto`` to attempt to deduce these when adding this variable to a processor for the first time. """ - def __init__(self, proc_chain: ProcessingChain, name: str, shape: tuple = auto, + def __init__(self, proc_chain: ProcessingChain, name: str, + shape: int | tuple[int, ...] = auto, dtype: np.dtype = auto, grid: CoordinateGrid = auto, unit: str | Unit = auto, is_coord: bool = auto) -> None: """ Parameters ---------- - proc_chain : ProcessingChain + proc_chain DOCME - name : str + name DOCME - shape : tuple, optional, default='auto' + shape DOCME - dtype : numpy.dtype or str, optional, default='auto' + dtype DOCME - grid : CoordinateGrid + grid DOCME - unit : str or pint.Unit, optional, default='auto' + unit DOCME - is_coord : bool, optional, default='auto' + is_coord DOCME """ assert isinstance(proc_chain, ProcessingChain) and isinstance(name, str) @@ -169,7 +172,7 @@ def __setattr__(self, name: str, value: Any) -> None: super().__setattr__(name, value) - def get_buffer(self, unit: str | Unit = None) -> np.array: + def get_buffer(self, unit: str | Unit = None) -> np.ndarray: # If buffer needs to be created, do so now if self._buffer is None: if self.shape is auto: @@ -225,18 +228,18 @@ def period(self): def offset(self): return self.grid.offset if self.grid else None - def description(self): + def description(self) -> str: return (f"{self.name}(shape: {self.shape}, " f"dtype: {self.dtype}, grid: {self.grid}, " f"unit: {self.unit}, is_coord: {self.is_coord})") - def update_auto(self, shape: tuple = auto, dtype: np.dtype = auto, + def update_auto(self, shape: int | tuple[int, ...] = auto, + dtype: np.dtype = auto, grid: CoordinateGrid = auto, unit: str | Unit = auto, is_coord: bool = auto, period: period = None, offset: offset = 0) -> None: - """ - Update any variables set to auto; leave the others alone. Emit a message - only if anything was updated. + """Update any variables set to ``auto``; leave the others alone. Emit a + message only if anything was updated. """ updated = False @@ -299,9 +302,9 @@ def __init__(self, block_width: int = 8, buffer_len: int = None) -> None: """ Parameters ---------- - block_width : int + block_width number of entries to simultaneously process. - buffer_len : int + buffer_len length of input and output buffers. Should be a multiple of `block_width` """ @@ -316,9 +319,15 @@ def __init__(self, block_width: int = 8, buffer_len: int = None) -> None: self._block_width = block_width self._buffer_len = buffer_len - def add_variable(self, name: str, dtype: np.dtype | str = auto, shape: int | tuple = auto, - grid: CoordinateGrid = auto, unit: str | Unit = auto, is_coord: bool = auto, - period: CoordinateGrid.period = None, offset: CoordinateGrid.offset = 0) -> ProcChainVar: + def add_variable(self, + name: str, + dtype: np.dtype | str = auto, + shape: int | tuple[int, ...] = auto, + grid: CoordinateGrid = auto, + unit: str | Unit = auto, + is_coord: bool = auto, + period: CoordinateGrid.period = None, + offset: CoordinateGrid.offset = 0) -> ProcChainVar: """Add a named variable containing a block of values or arrays. Parameters @@ -359,24 +368,26 @@ def add_variable(self, name: str, dtype: np.dtype | str = auto, shape: int | tup self._vars_dict[name] = var return var - # TODO: type: numpy array or lgdo. need an lgdo type - def link_input_buffer(self, varname: str, buff=None): - """Link an input buffer to a variable + def link_input_buffer(self, + varname: str, + buff: np.ndarray | LGDO = None + ) -> np.ndarray | LGDO: + """Link an input buffer to a variable. Parameters ---------- - varname : str + varname name of internal variable to copy into buffer at the end of processor execution. If variable does not yet exist, it will - be created with a similar shape to the provided buffer - buff : numpy.array or lgdo class, optional, default=None + be created with a similar shape to the provided buffer. + buff object to use as input buffer. If ``None``, create a new buffer - with a similar shape to the variable + with a similar shape to the variable. Returns ------- - buffer: numpy.array or lgdo class - buff or newly allocated input buffer + buffer + `buff` or newly allocated input buffer. """ self._validate_name(varname, raise_exception=True) var = self.get_variable(varname, expr_only=True) @@ -418,23 +429,26 @@ def link_input_buffer(self, varname: str, buff=None): return buff - def link_output_buffer(self, varname: str, buff=None): + def link_output_buffer(self, + varname: str, + buff: np.ndarray | LGDO = None + ) -> np.ndarray | LGDO: """Link an output buffer to a variable. Parameters ---------- varname - name of internal variable to copy into buffer at the end - of processor execution. If variable does not yet exist, it will - be created with a similar shape to the provided buffer - buff : numpy.array or lgdo class - object to use as output buffer. If None, create a new buffer with a - similar shape to the variable + name of internal variable to copy into buffer at the end of + processor execution. If variable does not yet exist, it will be + created with a similar shape to the provided buffer. + buff + object to use as output buffer. If ``None``, create a new buffer + with a similar shape to the variable Returns ------- - numpy.array or lgdo object - buff or newly allocated output buffer + buffer + `buff` or newly allocated output buffer. """ self._validate_name(varname, raise_exception=True) var = self.get_variable(varname, expr_only=True) @@ -532,7 +546,7 @@ def get_variable(self, expr: str, get_names_only: bool = False, - ``keyword = expr``: return a ``dict`` with a single element pointing from keyword to the parsed `expr`. This is used for `kwargs`. If `expr_only` is ``True``, raise an exception if we see - this + this. If `get_names_only` is set to ``True``, do not fetch or allocate new arrays, instead return a list of variable names found in the expression. @@ -761,7 +775,7 @@ def _execute_procs(self, begin: int, end: int) -> str: for out_man in self._output_managers: out_man.write(begin, end) - def __str__(self): + def __str__(self) -> str: return 'Input variables:\n ' \ + '\n '.join([str(in_man) for in_man in self._input_managers]) \ + '\nProcessors:\n ' \ @@ -1054,7 +1068,7 @@ def __str__(self) -> str: # Ok, this one's not LGDO class NumpyIOManager(IOManager): """IO Manager for buffers that are numpy arrays""" - def __init__(self, io_buf: np.array, var: ProcChainVar): + def __init__(self, io_buf: np.ndarray, var: ProcChainVar) -> None: assert isinstance(io_buf, np.ndarray) \ and isinstance(var, ProcChainVar) @@ -1180,7 +1194,7 @@ def __str__(self) -> str: class LGDOWaveformIOManager(IOManager): - def __init__(self, wf_table: lgdo.WaveformTable, variable: ProcChainVar): + def __init__(self, wf_table: lgdo.WaveformTable, variable: ProcChainVar) -> None: assert isinstance(wf_table, lgdo.WaveformTable) \ and isinstance(variable, ProcChainVar) @@ -1249,20 +1263,24 @@ def __str__(self) -> str: f"t0(shape={self.wf_table.t0.nda.shape}, dtype={self.wf_table.t0.nda.dtype}, attrs={self.wf_table.t0.attrs}))") -def build_processing_chain(lh5_in: lgdo.Table, dsp_config: dict | str, db_dict: dict = None, - outputs: list[str] = None, block_width: int = 16) -> tuple[ProcessingChain, list[str], lgdo.Table]: +def build_processing_chain(lh5_in: lgdo.Table, + dsp_config: dict | str, + db_dict: dict = None, + outputs: list[str] = None, + block_width: int = 16 + ) -> tuple[ProcessingChain, list[str], lgdo.Table]: """Produces a :class:`ProcessingChain` object and an LH5 :class:`~.lgdo.table.Table` for output parameters from an input LH5 :class:`~.lgdo.table.Table` and a JSON recipe. Parameters ---------- - lh5_in : lgdo.Table + lh5_in HDF5 table from which raw data is read. At least one row of entries should be read in prior to calling this! - dsp_config : dict or str - A dict or json filename containing the recipes for computing DSP + dsp_config + A dictionary or JSON filename containing the recipes for computing DSP parameter from raw parameters. The format is as follows: .. code-block:: json @@ -1290,42 +1308,43 @@ def build_processing_chain(lh5_in: lgdo.Table, dsp_config: dict | str, db_dict: names of parameters computed - ``function`` -- string, name of function to call. Function - should implement the gufunc interface, a factory function - returning a gufunc, or an arbitrary function that can be - mapped onto a gufunc + should implement the :class:`numpy.gufunc` interface, a factory + function returning a ``gufunc``, or an arbitrary function that + can be mapped onto a ``gufunc`` - ``module`` -- string, name of module containing function - ``args``-- list of strings or numerical values. Contains list of names of computed and input parameters or constant values used as inputs to function. Note that outputs should be fed by reference as args! Arguments read - from the database are prepended with db. + from the database are prepended with ``db``. - ``kwargs`` -- dictionary. Keyword arguments for :meth:`ProcesssingChain.add_processor`. - ``init_args`` -- list of strings or numerical values. List of names of computed and input parameters or constant values - used to initialize a gufunc via a factory function + used to initialize a :class:`numpy.gufunc` via a factory + function - ``unit`` -- list of strings. Units for parameters - ``defaults`` -- dictionary. Default value to be used for arguments read from the database - db_dict : dict, optional, default=None - A nested dict pointing to values for db args. e.g. if a processor - uses arg db.trap.risetime, it will look up db_dict['trap']['risetime'] - and use the found value. If no value is found, use the default - defined in the config file. + db_dict + A nested :class:`dict` pointing to values for database arguments. As + instance, if a processor uses the argument ``db.trap.risetime``, it + will look up ``db_dict['trap']['risetime']`` and use the found value. + If no value is found, use the default defined in `dsp_config`. - outputs : list of str, optional, default=None - List of parameters to put in the output lh5 table. If None, - use the parameters in the 'outputs' list from config + outputs + List of parameters to put in the output LH5 table. If ``None``, + use the parameters in the ``"outputs"`` list from `dsp_config`. - block_width : int, optional, default=16 + block_width number of entries to process at once. To optimize performance, a multiple of 16 is preferred, but if performance is not an issue any value can be used. Returns ------- - (proc_chain, field_mask, lh5_out) : tuple + (proc_chain, field_mask, lh5_out) - `proc_chain` -- :class:`ProcessingChain` object that is executed - `field_mask` -- list of input fields that are used - `lh5_out` -- output :class:`~.lgdo.table.Table` containing processed diff --git a/src/pygama/dsp/processors.py b/src/pygama/dsp/processors/__init__.py similarity index 70% rename from src/pygama/dsp/processors.py rename to src/pygama/dsp/processors/__init__.py index 09ecda024..db56a5d87 100644 --- a/src/pygama/dsp/processors.py +++ b/src/pygama/dsp/processors/__init__.py @@ -41,7 +41,7 @@ vectorization on a vector-CPU. Modern CPUs typically have 256- or 512-bit wide processing units, which can accommodate multiple 32- or 64-bit numbers. Programming with these, however, is quite difficult and requires specialized - commands to be called. Luckily for us, many NumPy :class`~numpy.ufunc`\ s + commands to be called. Luckily for us, many NumPy :class:`~numpy.ufunc`\ s will automatically use these for us, speeding up our code! 3. :class:`~numpy.ufunc`\ s are capable of calculating their output in place, @@ -59,37 +59,32 @@ to use functions that operate in place as much as possible! """ -from ._processors.bl_subtract import bl_subtract -from ._processors.convolutions import cusp_filter, t0_filter, zac_filter -from ._processors.fftw import dft, inv_dft, psd -from ._processors.fixed_time_pickoff import fixed_time_pickoff -from ._processors.gaussian_filter1d import gaussian_filter1d -from ._processors.get_multi_local_extrema import get_multi_local_extrema -from ._processors.linear_slope_fit import linear_slope_fit -from ._processors.log_check import log_check -from ._processors.min_max import min_max -from ._processors.moving_windows import ( +from .bl_subtract import bl_subtract +from .convolutions import cusp_filter, t0_filter, zac_filter +from .fftw import dft, inv_dft, psd +from .fixed_time_pickoff import fixed_time_pickoff +from .gaussian_filter1d import gaussian_filter1d +from .get_multi_local_extrema import get_multi_local_extrema +from .linear_slope_fit import linear_slope_fit +from .log_check import log_check +from .min_max import min_max +from .moving_windows import ( avg_current, moving_window_left, moving_window_multi, moving_window_right, ) -from ._processors.multi_a_filter import multi_a_filter -from ._processors.multi_t_filter import multi_t_filter, remove_duplicates -from ._processors.optimize import optimize_1pz, optimize_2pz -from ._processors.param_lookup import param_lookup -from ._processors.pole_zero import double_pole_zero, pole_zero -from ._processors.presum import presum -from ._processors.pulse_injector import inject_exp_pulse, inject_sig_pulse -from ._processors.saturation import saturation -from ._processors.soft_pileup_corr import soft_pileup_corr, soft_pileup_corr_bl -from ._processors.time_point_thresh import time_point_thresh -from ._processors.trap_filters import ( - asym_trap_filter, - trap_filter, - trap_norm, - trap_pickoff, -) -from ._processors.upsampler import upsampler -from ._processors.wiener_filter import wiener_filter -from ._processors.windower import windower +from .multi_a_filter import multi_a_filter +from .multi_t_filter import multi_t_filter, remove_duplicates +from .optimize import optimize_1pz, optimize_2pz +from .param_lookup import param_lookup +from .pole_zero import double_pole_zero, pole_zero +from .presum import presum +from .pulse_injector import inject_exp_pulse, inject_sig_pulse +from .saturation import saturation +from .soft_pileup_corr import soft_pileup_corr, soft_pileup_corr_bl +from .time_point_thresh import time_point_thresh +from .trap_filters import asym_trap_filter, trap_filter, trap_norm, trap_pickoff +from .upsampler import upsampler +from .wiener_filter import wiener_filter +from .windower import windower diff --git a/src/pygama/dsp/_processors/bl_subtract.py b/src/pygama/dsp/processors/bl_subtract.py similarity index 56% rename from src/pygama/dsp/_processors/bl_subtract.py rename to src/pygama/dsp/processors/bl_subtract.py index 08e982b50..6988375a2 100644 --- a/src/pygama/dsp/_processors/bl_subtract.py +++ b/src/pygama/dsp/processors/bl_subtract.py @@ -1,35 +1,34 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp.errors import DSPFatal - @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),()->(n)", nopython=True, cache=True) -def bl_subtract(w_in, a_baseline, w_out): - """ - Subtract the constant baseline from the entire waveform. +def bl_subtract(w_in: np.ndarray, a_baseline: float, w_out: np.ndarray) -> None: + """Subtract the constant baseline from the entire waveform. Parameters ---------- - w_in : array-like - The input waveform - a_baseline : float - The baseline value to subtract - w_out : array-like - The output waveform with the baseline subtracted - - Examples - -------- + w_in + the input waveform. + a_baseline + the baseline value to subtract. + w_out + the output waveform with the baseline subtracted. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_bl": { "function": "bl_subtract", "module": "pygama.dsp.processors", "args": ["waveform", "baseline", "wf_bl"], - "unit": "ADC", - "prereqs": ["waveform", "baseline"] + "unit": "ADC" } """ w_out[:] = np.nan diff --git a/src/pygama/dsp/_processors/convolutions.py b/src/pygama/dsp/processors/convolutions.py similarity index 69% rename from src/pygama/dsp/_processors/convolutions.py rename to src/pygama/dsp/processors/convolutions.py index a60b31612..dd4740dc5 100644 --- a/src/pygama/dsp/_processors/convolutions.py +++ b/src/pygama/dsp/processors/convolutions.py @@ -1,28 +1,36 @@ +from __future__ import annotations + +from typing import Callable + import numpy as np from numba import guvectorize from pygama.dsp.errors import DSPFatal -def cusp_filter(length, sigma, flat, decay): - """ - Apply a CUSP filter to the waveform. Note that it is composed of a - factory function that is called using the init_args argument and that - the function the waveforms are passed to using args. +def cusp_filter(length: int, sigma: float, flat: int, decay: int) -> Callable: + """Apply a CUSP filter to the waveform. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. Parameters ---------- - length : int - The length of the filter to be convolved - sigma : float - The curvature of the rising and falling part of the kernel - flat : int - The length of the flat section - decay : int - The decay constant of the exponential to be convolved - - Examples - -------- + length + the length of the filter to be convolved. + sigma + the curvature of the rising and falling part of the kernel. + flat + the length of the flat section. + decay + the decay constant of the exponential to be convolved. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_cusp": { @@ -30,7 +38,6 @@ def cusp_filter(length, sigma, flat, decay): "module": "pygama.dsp.processors", "args": ["wf_bl", "wf_cusp(101,f)"], "unit": "ADC", - "prereqs": ["wf_bl"], "init_args": ["len(wf_bl)-100", "40*us", "3*us", "45*us"] } """ @@ -62,14 +69,14 @@ def cusp_filter(length, sigma, flat, decay): @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])"], "(n),(m)", forceobj=True) - def cusp_out(w_in, w_out): + def cusp_out(w_in: np.ndarray, w_out: np.ndarray) -> None: """ Parameters ---------- - w_in : array-like - The input waveform - w_out : array-like - The filtered waveform + w_in + the input waveform. + w_out + the filtered waveform. """ w_out[:] = np.nan @@ -83,25 +90,30 @@ def cusp_out(w_in, w_out): return cusp_out -def zac_filter(length, sigma, flat, decay): - """ - Apply a ZAC (Zero Area CUSP) filter to the waveform. Note that it is - composed of a factory function that is called using the init_args - argument and that the function the waveforms are passed to using args. + +def zac_filter(length: int, sigma: float, flat: int, decay: int) -> Callable: + """Apply a ZAC (Zero Area CUSP) filter to the waveform. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. Parameters ---------- - length : int - The length of the filter to be convolved - sigma : float - The curvature of the rising and falling part of the kernel - flat : int - The length of the flat section - decay : int - The decay constant of the exponential to be convolved - - Examples - -------- + length + the length of the filter to be convolved. + sigma + the curvature of the rising and falling part of the kernel. + flat + the length of the flat section. + decay + the decay constant of the exponential to be convolved. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_zac": { @@ -109,7 +121,6 @@ def zac_filter(length, sigma, flat, decay): "module": "pygama.dsp.processors", "args": ["wf_bl", "wf_zac(101,f)"], "unit": "ADC", - "prereqs": ["wf_bl"], "init_args": ["len(wf_bl)-100", "40*us", "3*us", "45*us"], } """ @@ -159,14 +170,14 @@ def zac_filter(length, sigma, flat, decay): @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])"], "(n),(m)", forceobj=True) - def zac_out(w_in, w_out): + def zac_out(w_in: np.ndarray, w_out: np.ndarray) -> None: """ Parameters ---------- - w_in : array-like - The input waveform - w_out : array-like - The filtered waveform + w_in + the input waveform. + w_out + the filtered waveform. """ w_out[:] = np.nan @@ -180,24 +191,28 @@ def zac_out(w_in, w_out): return zac_out -def t0_filter(rise, fall): - """ - Apply a modified, asymmetric trapezoidal filter to the waveform. Note - that it is composed of a factory function that is called using the init_args - argument and that the function the waveforms are passed to using args. +def t0_filter(rise: int, fall: int) -> Callable: + """Apply a modified, asymmetric trapezoidal filter to the waveform. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. Parameters ---------- - rise : int - The length of the rise section. This is the linearly increasing + rise + the length of the rise section. This is the linearly increasing section of the filter that performs a weighted average. - fall : int - The length of the fall section. This is the simple averaging part + fall + the length of the fall section. This is the simple averaging part of the filter. - Examples - -------- + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_t0_filter": { @@ -205,7 +220,6 @@ def t0_filter(rise, fall): "module": "pygama.dsp.processors", "args": ["wf_pz", "wf_t0_filter(3748,f)"], "unit": "ADC", - "prereqs": ["wf_pz"], "init_args": ["128*ns", "2*us"] } """ @@ -221,13 +235,13 @@ def t0_filter(rise, fall): @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])"], "(n),(m)", forceobj=True) - def t0_filter_out(w_in, w_out): + def t0_filter_out(w_in: np.ndarray, w_out: np.ndarray) -> None: """ Parameters ---------- - w_in : array-like + w_in The input waveform - w_out : array-like + w_out The filtered waveform """ w_out[:] = np.nan diff --git a/src/pygama/dsp/processors/fftw.py b/src/pygama/dsp/processors/fftw.py new file mode 100644 index 000000000..40e0b441a --- /dev/null +++ b/src/pygama/dsp/processors/fftw.py @@ -0,0 +1,159 @@ +from __future__ import annotations + +from typing import Callable + +import numpy as np +from numba import guvectorize +from pyfftw import FFTW + + +def dft(buf_in: np.ndarray, buf_out: np.ndarray) -> Callable: + """Perform discrete Fourier transforms using the FFTW library. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. + + Note + ---- + FFTW optimizes the FFT algorithm based on the size of the arrays, with SIMD + parallelized commands. This optimization requires initialization, so this + is a factory function that returns a Numba gufunc that performs the FFT. + FFTW works on fixed memory buffers, so you must tell it what memory to use + ahead of time. When using this with + :class:`~.dsp.processing_chain.ProcessingChain`, to ensure the correct + buffers are used, call + :meth:`~.dsp.processing_chain.ProcessingChain.get_variable` to give it the + internal memory buffer directly. With :func:`~.dsp.build_dsp.build_dsp`, + you can just give it the name, and it will automatically happen. The + possible `dtypes` for the input/outputs are: + + =============================== ========= =============================== ============= + :class:`numpy.dtype` Size :class:`numpy.dtype` Size + ------------------------------- --------- ------------------------------- ------------- + ``float32``/``float`` :math:`n` ``complex64`` :math:`n/2+1` + ``float64``/``double`` :math:`n` ``complex128`` :math:`n/2+1` + ``float128``/``longdouble`` :math:`n` ``complex256``/``clongdouble`` :math:`n/2+1` + ``complex64`` :math:`n` ``complex64`` :math:`n` + ``complex128`` :math:`n` ``complex128`` :math:`n` + ``complex256``/``clongdouble`` :math:`n` ``complex256``/``clongdouble`` :math:`n` + =============================== ========= =============================== ============= + """ + try: + dft_fun = FFTW(buf_in, buf_out, axes=(-1,), direction='FFTW_FORWARD') + except ValueError: + raise ValueError("incompatible array types/shapes. See function documentation for allowed values") + + typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' + sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' + + @guvectorize([typesig], sizesig, forceobj=True) + def dft(wf_in: np.ndarray, dft_out: np.ndarray) -> None: + dft_fun(wf_in, dft_out) + + return dft + + +def inv_dft(buf_in: np.ndarray, buf_out: np.ndarray) -> Callable: + """Perform inverse discrete Fourier transforms using the FFTW library. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. + + Note + ---- + FFTW optimizes the FFT algorithm based on the size of the arrays, with SIMD + parallelized commands. This optimization requires initialization, so this + is a factory function that returns a Numba gufunc that performs the FFT. + FFTW works on fixed memory buffers, so you must tell it what memory to use + ahead of time. When using this with + :class:`~.dsp.processing_chain.ProcessingChain`, to ensure the correct + buffers are used, call + :meth:`~.dsp.processing_chain.ProcessingChain.get_variable` to give it the + internal memory buffer directly. With :func:`~.dsp.build_dsp.build_dsp`, + you can just give it the name, and it will automatically happen. The + possible `dtypes` for the input/outputs are: + + =============================== ============= =============================== ========= + :class:`numpy.dtype` Size :class:`numpy.dtype` Size + ------------------------------- ------------- ------------------------------- --------- + ``complex64`` :math:`n/2+1` ``float32``/``float`` :math:`n` + ``complex128`` :math:`n/2+1` ``float64``/``double`` :math:`n` + ``complex256``/``clongdouble`` :math:`n/2+1` ``float128``/``longdouble`` :math:`n` + ``complex64`` :math:`n` ``complex64`` :math:`n` + ``complex128`` :math:`n` ``complex128`` :math:`n` + ``complex256``/``clongdouble`` :math:`n` ``complex256``/``clongdouble`` :math:`n` + =============================== ============= =============================== ========= + """ + try: + idft_fun = FFTW(buf_in, buf_out, axes=(-1,), direction='FFTW_BACKWARD') + except ValueError: + raise ValueError("incompatible array types/shapes. See function documentation for allowed values") + + typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' + sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' + + @guvectorize([typesig], sizesig, forceobj=True) + def inv_dft(wf_in: np.ndarray, dft_out: np.ndarray) -> None: + idft_fun(wf_in, dft_out) + + return inv_dft + + +def psd(buf_in: np.ndarray, buf_out: np.ndarray) -> Callable: + """Perform discrete Fourier transforms using the FFTW library, and use it to get + the power spectral density. + + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. + + Note + ---- + FFTW optimizes the FFT algorithm based on the size of the arrays, with SIMD + parallelized commands. This optimization requires initialization, so this + is a factory function that returns a Numba gufunc that performs the FFT. + FFTW works on fixed memory buffers, so you must tell it what memory to use + ahead of time. When using this with + :class:`~.dsp.processing_chain.ProcessingChain`, to ensure the correct + buffers are used, call + :meth:`~.dsp.processing_chain.ProcessingChain.get_variable` to give it the + internal memory buffer directly. With :func:`~.dsp.build_dsp.build_dsp`, + you can just give it the name, and it will automatically happen. The + possible `dtypes` for the input/outputs are: + + =============================== ========= ============================ ============= + :class:`numpy.dtype` Size :class:`numpy.dtype` Size + ------------------------------- --------- ---------------------------- ------------- + ``complex64`` :math:`n` ``float32``/``float`` :math:`n` + ``complex128`` :math:`n` ``float64``/``double`` :math:`n` + ``complex256``/``clongdouble`` :math:`n` ``float128``/``longdouble`` :math:`n` + ``float32``/``float`` :math:`n` ``float32``/``float`` :math:`n/2+1` + ``float64``/``double`` :math:`n` ``float64``/``double`` :math:`n/2+1` + ``float128``/``longdouble`` :math:`n` ``float128``/``longdouble`` :math:`n/2+1` + =============================== ========= ============================ ============= + """ + + # build intermediate array for the dft, which will be abs'd to get the PSD + buf_dft = np.ndarray(buf_out.shape, np.dtype('complex' + str(buf_out.dtype.itemsize * 16))) + try: + dft_fun = FFTW(buf_in, buf_dft, axes=(-1,), direction='FFTW_FORWARD') + except ValueError: + raise ValueError("incompatible array types/shapes. See function documentation for allowed values") + + typesig = 'void(' + str(buf_in.dtype) + '[:, :], ' + str(buf_out.dtype) + '[:, :])' + sizesig = '(m, n)->(m, n)' if buf_in.shape == buf_out.shape else '(m, n),(m, l)' + + @guvectorize([typesig], sizesig, forceobj=True) + def psd(wf_in: np.ndarray, psd_out: np.ndarray) -> None: + dft_fun(wf_in, buf_dft) + np.abs(buf_dft, psd_out) + + return psd diff --git a/src/pygama/dsp/_processors/fixed_time_pickoff.py b/src/pygama/dsp/processors/fixed_time_pickoff.py similarity index 62% rename from src/pygama/dsp/_processors/fixed_time_pickoff.py rename to src/pygama/dsp/processors/fixed_time_pickoff.py index c40036085..b9fad4137 100644 --- a/src/pygama/dsp/_processors/fixed_time_pickoff.py +++ b/src/pygama/dsp/processors/fixed_time_pickoff.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,30 +9,30 @@ @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),()->()", nopython=True, cache=True) -def fixed_time_pickoff(w_in, t_in, a_out): - """ - Pick off the waveform value at the provided index. If the - provided index is out of range, return NaN. +def fixed_time_pickoff(w_in: np.ndarray, t_in: int, a_out: np.ndarray) -> None: + """Pick off the waveform value at the provided index. + + If the provided index `t_in` is out of range, return :any:`numpy.nan`. Parameters ---------- - w_in : array-like - The input waveform - t_in : int - The waveform index to pick off - a_out : float - The output pick-off value - - Examples - -------- + w_in + the input waveform + t_in + the waveform index to pick off + a_out + the output pick-off value + + JSON Configuration Example + -------------------------- + .. code-block :: json "trapEftp": { "function": "fixed_time_pickoff", "module": "pygama.dsp.processors", "args": ["wf_trap", "tp_0+10*us", "trapEftp"], - "unit": "ADC", - "prereqs": ["wf_trap", "tp_0"] + "unit": "ADC" } """ a_out[0] = np.nan diff --git a/src/pygama/dsp/_processors/gaussian_filter1d.py b/src/pygama/dsp/processors/gaussian_filter1d.py similarity index 89% rename from src/pygama/dsp/_processors/gaussian_filter1d.py rename to src/pygama/dsp/processors/gaussian_filter1d.py index cb0f90111..3eb885276 100644 --- a/src/pygama/dsp/_processors/gaussian_filter1d.py +++ b/src/pygama/dsp/processors/gaussian_filter1d.py @@ -33,26 +33,27 @@ # The only thing changed was the calculation of the convulution, which # originally called a function from a C library. In this code, the convolution is # performed with NumPy's built in convolution function. +from __future__ import annotations import numpy from numba import guvectorize -def gaussian_filter1d(sigma, truncate): +def gaussian_filter1d(sigma: int, truncate: float = 4.0) -> numpy.ndarray: """1-D Gaussian filter. + Note + ---- + This processor is composed of a factory function that is called using the + `init_args` argument. The input and output waveforms are passed using + `args`. + Parameters ---------- - %(input)s - sigma : scalar + sigma standard deviation for Gaussian kernel - truncate : float, optional - Truncate the filter at this many standard deviations. - Default is 4.0. - - Returns - ------- - gaussian_filter1d : ndarray + truncate + truncate the filter at this many standard deviations. """ def _gaussian_kernel1d(sigma, radius): """ @@ -79,17 +80,6 @@ def _gaussian_kernel1d(sigma, radius): extension_length = int(len(weights) / 2) + 1 - """Calculate a 1-D correlation along the given axis. - The lines of the array along the given axis are correlated with the - given weights. - - Parameters - ---------- - %(input)s - weights : array - 1-D sequence of numbers. - %(mode_reflect)s - """ @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])", "void(int32[:], int32[:])", diff --git a/src/pygama/dsp/_processors/get_multi_local_extrema.py b/src/pygama/dsp/processors/get_multi_local_extrema.py similarity index 65% rename from src/pygama/dsp/_processors/get_multi_local_extrema.py rename to src/pygama/dsp/processors/get_multi_local_extrema.py index 939cd5b3a..91c187e06 100644 --- a/src/pygama/dsp/_processors/get_multi_local_extrema.py +++ b/src/pygama/dsp/processors/get_multi_local_extrema.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,32 +9,33 @@ @guvectorize(["void(float32[:], float32, float32[:], float32[:], float32[:], float32[:], float32[:])", "void(float64[:], float64, float64[:], float64[:],float64[:], float64[:], float64[:])"], "(n),(),(m),(m),(),(),()", nopython=True, cache=True) -def get_multi_local_extrema(w_in, a_delta_in, vt_max_out, vt_min_out, n_max_out, n_min_out, flag_out): - """ - Get lists of indices of the local maxima and minima of data - The "local" extrema are those maxima / minima that have heights / depths of - at least a_delta_in. - Converted from MATLAB script at: http://billauer.co.il/peakdet.html +def get_multi_local_extrema(w_in: np.ndarray, a_delta_in: float, + vt_max_out: np.ndarray, vt_min_out: np.ndarray, + n_max_out: int, n_min_out: int, flag_out: int) -> None: + """Get lists of indices of the local maxima and minima of data. + + The "local" extrema are those maxima (minima) that have heights (depths) of + at least `a_delta_in`. Converted from a `MATLAB script + `_ by E. Billauer. Parameters ---------- - w_in : array-like - The array of data within which extrema will be found - a_delta_in : scalar - The absolute level by which data must vary (in one direction) about an - extremum in order for it to be tagged - vt_max_out, vt_min_out : array-like, array-like - Arrays of fixed length (padded with nans) that hold the indices of - the identified local maxima and minima - n_max_out, n_min_out : scalar, scalar - The number of maxima and minima found in a waveform - flag_out : scalar - Returns 1 if there is only one maximum and it is a simple waveform, - Returns 0 if there are no peaks, or multiple peaks in a waveform + w_in + the array of data within which extrema will be found. + a_delta_in + the absolute level by which data must vary (in one direction) about an + extremum in order for it to be tagged. + vt_max_out, vt_min_out + arrays of fixed length (padded with :any:`numpy.nan`) that hold the + indices of the identified local maxima and minima. + n_max_out, n_min_out + the number of maxima and minima found in a waveform. + flag_out + returns ``1`` if there is only one maximum and it is a simple waveform, + returns ``0`` if there are no peaks, or multiple peaks in a waveform. """ # prepare output - vt_max_out[:]= np.nan vt_min_out[:] = np.nan n_max_out[0] = np.nan @@ -40,13 +43,10 @@ def get_multi_local_extrema(w_in, a_delta_in, vt_max_out, vt_min_out, n_max_out, flag_out[0] = np.nan # initialize internal counters - n_max_counter = 0 n_min_counter = 0 - # Checks - if (np.isnan(w_in).any() or np.isnan(a_delta_in)): return diff --git a/src/pygama/dsp/_processors/linear_slope_fit.py b/src/pygama/dsp/processors/linear_slope_fit.py similarity index 76% rename from src/pygama/dsp/_processors/linear_slope_fit.py rename to src/pygama/dsp/processors/linear_slope_fit.py index 5bde7cd00..e9e0e1f5d 100644 --- a/src/pygama/dsp/_processors/linear_slope_fit.py +++ b/src/pygama/dsp/processors/linear_slope_fit.py @@ -1,13 +1,14 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp.errors import DSPFatal - @guvectorize(["void(float32[:], float32[:], float32[:], float32[:], float32[:])", "void(float64[:], float64[:], float64[:], float64[:], float64[:])"], "(n)->(),(),(),()", nopython=True, cache=True) -def linear_slope_fit(w_in, mean, stdev, slope, intercept): +def linear_slope_fit(w_in: np.ndarray, mean: float, stdev: float, + slope: float, intercept: float) -> None: """ Calculate the mean and standard deviation of the waveform using Welford's method as well as the slope an intercept of the waveform @@ -15,19 +16,20 @@ def linear_slope_fit(w_in, mean, stdev, slope, intercept): Parameters ---------- - w_in : array-like - The input waveform - mean : float - The mean of the waveform - stdev : float - The standard deviation of the waveform - slope : float - The slope of the linear fit - intercept : float - The intercept of the linear fit + w_in + the input waveform. + mean + the mean of the waveform. + stdev + the standard deviation of the waveform. + slope + the slope of the linear fit. + intercept + the intercept of the linear fit. + + JSON Configuration Example + -------------------------- - Examples - -------- .. code-block :: json "bl_mean, bl_std, bl_slope, bl_intercept": { @@ -35,7 +37,6 @@ def linear_slope_fit(w_in, mean, stdev, slope, intercept): "module": "pygama.dsp.processors", "args": ["wf_blsub[0:1650]", "bl_mean", "bl_std", "bl_slope", "bl_intercept"], "unit": ["ADC", "ADC", "ADC", "ADC"], - "prereqs": ["wf_blsub"] } """ mean [0] = np.nan diff --git a/src/pygama/dsp/_processors/log_check.py b/src/pygama/dsp/processors/log_check.py similarity index 71% rename from src/pygama/dsp/_processors/log_check.py rename to src/pygama/dsp/processors/log_check.py index e2682b235..fbb7ab813 100644 --- a/src/pygama/dsp/_processors/log_check.py +++ b/src/pygama/dsp/processors/log_check.py @@ -1,34 +1,34 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp.errors import DSPFatal - @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])"], "(n)->(n)", nopython=True, cache=True) -def log_check(w_in, w_log): +def log_check(w_in: np.ndarray, w_log: np.ndarray) -> None: """ Calculate the logarithm of the waveform if all its values are positive; otherwise, return NaN. Parameters ---------- - w_in : array-like - The input waveform - w_log : array-like - The output waveform with logged values + w_in + the input waveform. + w_log + the output waveform with logged values. + + JSON Configuration Example + -------------------------- - Examples - -------- .. code-block :: json "wf_logged": { "function": "log_check", "module": "pygama.dsp.processors", "args": ["wf_blsub[2100:]", "wf_logged"], - "unit": "ADC", - "prereqs": ["wf_blsub"] + "unit": "ADC" } """ w_log[:] = np.nan diff --git a/src/pygama/dsp/_processors/min_max.py b/src/pygama/dsp/processors/min_max.py similarity index 61% rename from src/pygama/dsp/_processors/min_max.py rename to src/pygama/dsp/processors/min_max.py index 5fae1461a..511b67123 100644 --- a/src/pygama/dsp/_processors/min_max.py +++ b/src/pygama/dsp/processors/min_max.py @@ -1,41 +1,44 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp.errors import DSPFatal - @guvectorize(["void(float32[:], float32[:], float32[:], float32[:], float32[:])", "void(float64[:], float64[:], float64[:], float64[:], float64[:])"], "(n)->(),(),(),()", nopython=True, cache=True) -def min_max(w_in, t_min, t_max, a_min, a_max): - """ - Find the value and index of the minimum and maximum values - in the waveform. Note that the first instance of each extremum - in the waveform will be returned. +def min_max(w_in: np.ndarray, t_min: int, t_max: int, + a_min: float, a_max: float) -> None: + """Find the value and index of the minimum and maximum values + in the waveform. + + Note + ---- + The first found instance of each extremum in the waveform will be returned. Parameters ---------- - w_in : array-like - The input waveform - t_min : int - The index of the minimum value - t_max : int - The index of the maximum value - a_min : float - The minimum value - a_max : float - The maximum value - - Examples - -------- + w_in + the input waveform + t_min + the index of the minimum value + t_max + the index of the maximum value + a_min + the minimum value + a_max + the maximum value + + JSON Configuration Example + -------------------------- + .. code-block :: json "tp_min, tp_max, wf_min, wf_max": { "function": "min_max", "module": "pygama.dsp.processors", "args": ["waveform", "tp_min", "tp_max", "wf_min", "wf_max"], - "unit": ["ns", "ns", "ADC", "ADC"], - "prereqs": ["waveform"] + "unit": ["ns", "ns", "ADC", "ADC"] } """ a_min[0] = np.nan diff --git a/src/pygama/dsp/_processors/moving_windows.py b/src/pygama/dsp/processors/moving_windows.py similarity index 65% rename from src/pygama/dsp/_processors/moving_windows.py rename to src/pygama/dsp/processors/moving_windows.py index a594b2d27..4dae202fa 100644 --- a/src/pygama/dsp/_processors/moving_windows.py +++ b/src/pygama/dsp/processors/moving_windows.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,38 +9,41 @@ @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),()->(n)", nopython=True, cache=True) -def moving_window_left(w_in, length, w_out): - ''' - Applies a moving average window to the waveform from the left, assumes that the baseline is at 0. +def moving_window_left(w_in: np.ndarray, length: float, w_out: np.ndarray) -> None: + """Applies a moving average window to the waveform. + + Note + ---- + Starts from the left and assumes that the baseline is at zero. Parameters ---------- - w_in : array-like - The input waveform - length : float - Length of the moving window to be applied - w_out : array-like - Output waveform after moving window applied - - Examples - -------- + w_in + the input waveform. + length + length of the moving window to be applied. + w_out + output waveform after moving window applied. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_mw": { "function": "moving_window_left", "module": "pygama.dsp.processors", "args": ["wf_pz", "96*ns", "wf_mw"], - "unit": "ADC", - "prereqs": ["wf_pz"] + "unit": "ADC" } - ''' + """ w_out[:] = np.nan if (np.isnan(w_in).any()): return - if (not length >= 0 or not length< len(w_in)): + if (not length >= 0 or not length < len(w_in)): raise DSPFatal('length is out of range, must be between 0 and the length of the waveform') w_out[0] = w_in[0] @@ -51,42 +56,41 @@ def moving_window_left(w_in, length, w_out): @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),()->(n)", nopython=True, cache=True) -def moving_window_right(w_in, length, w_out): - ''' - Applies a moving average window to the waveform from the right. +def moving_window_right(w_in: np.ndarray, length: float, w_out: np.ndarray) -> None: + """Applies a moving average window to the waveform from the right. Parameters ---------- - w_in : array-like - The input waveform - length : float - Length of the moving window to be applied - w_out : array-like - Output waveform after moving window applied - - Examples - -------- + w_in + the input waveform. + length + length of the moving window to be applied. + w_out + output waveform after moving window applied. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_mw": { "function": "moving_window_right", "module": "pygama.dsp.processors", "args": ["wf_pz", "96*ns", "wf_mw"], - "unit": "ADC", - "prereqs": ["wf_pz"] + "unit": "ADC" } - ''' + """ w_out[:] = np.nan if (np.isnan(w_in).any()): return - if (not length >= 0 or not length< len(w_in)): + if (not length >= 0 or not length < len(w_in)): raise DSPFatal('length is out of range, must be between 0 and the length of the waveform') - w_out[-1]= w_in[-1] - for i in range(1, int(length),1): + w_out[-1] = w_in[-1] + for i in range(1, int(length), 1): w_out[len(w_in)-1-i] = w_out[len(w_in)-i] + (w_in[len(w_in)-1-i]-w_out[-1])/length for i in range(int(length), len(w_in), 1): w_out[len(w_in)-1-i] = w_out[len(w_in)-i] + (w_in[len(w_in)-1-i] - w_in[len(w_in)-1-i+int(length)])/length @@ -95,36 +99,36 @@ def moving_window_right(w_in, length, w_out): @guvectorize(["void(float32[:], float32, float32, int32, float32[:])", "void(float64[:], float64, float64, int32, float64[:])"], "(n),(),(),()->(n)", nopython=True, cache=True) -def moving_window_multi(w_in, length, num_mw, mw_type,w_out): - """ - Apply a series of moving-average windows to the waveform, alternating +def moving_window_multi(w_in: np.ndarray, length: float, + num_mw: int, mw_type: int, w_out: np.ndarray) -> None: + """Apply a series of moving-average windows to the waveform, alternating its application between the left and the right. Parameters ---------- - w_in : array-like - The input waveform - length : float - Length of the moving window to be applied - num_mw : int - The number of moving windows - mw_type : int - 0 = alternate moving windows right and left - 1 = only left - 2 = only right - w_out : array-like - The windowed waveform - - Examples - -------- + w_in + the input waveform. + length + length of the moving window to be applied. + num_mw + the number of moving windows. + mw_type + - ``0`` -- alternate moving windows right and left + - ``1`` -- only left + - ``2`` -- only right + w_out + the windowed waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "curr_av": { "function": "moving_window_multi", "module": "pygama.dsp.processors", "args": ["curr", "96*ns", "3", "0", "curr_av"], - "unit": "ADC/sample", - "prereqs": ["curr"] + "unit": "ADC/sample" } """ w_out[:] = np.nan @@ -161,35 +165,31 @@ def moving_window_multi(w_in, length, num_mw, mw_type,w_out): w_buf = w_out.copy() - - @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),(),(m)", nopython=True, cache=True) -def avg_current(w_in, length, w_out): - """ - Calculate the derivative of a WF, averaged across n samples. Dimension of - deriv should be len(wf) - n +def avg_current(w_in: np.ndarray, length: float, w_out: np.ndarray) -> None: + """Calculate the derivative of a waveform, averaged across `length` samples. Parameters ---------- - w_in : array-like - The input waveform - length : float - Length of the moving window to be applied - w_out : array-like - Output waveform after derivation - - Examples - -------- + w_in + the input waveform. + length + length of the moving window to be applied. + w_out + output waveform after derivation. + + JSON Configuration Example + -------------------------- + .. code-block :: json "curr": { "function": "avg_current", "module": "pygama.dsp.processors", "args": ["wf_pz", 1, "curr(len(wf_pz)-1, f)"], - "unit": "ADC/sample", - "prereqs": ["wf_pz"] + "unit": "ADC/sample" } """ @@ -198,9 +198,8 @@ def avg_current(w_in, length, w_out): if (np.isnan(w_in).any()): return - if (not length >= 0 or not length< len(w_in)): + if (not length >= 0 or not length < len(w_in)): raise DSPFatal('length is out of range, must be between 0 and the length of the waveform') - w_out[:] = w_in[int(length):] - w_in[:-int(length)] - w_out/=length + w_out /= length diff --git a/src/pygama/dsp/_processors/multi_a_filter.py b/src/pygama/dsp/processors/multi_a_filter.py similarity index 93% rename from src/pygama/dsp/_processors/multi_a_filter.py rename to src/pygama/dsp/processors/multi_a_filter.py index c1eaa2991..d11219e5b 100644 --- a/src/pygama/dsp/_processors/multi_a_filter.py +++ b/src/pygama/dsp/processors/multi_a_filter.py @@ -1,9 +1,10 @@ import numpy as np from numba import guvectorize -from pygama.dsp._processors.fixed_time_pickoff import fixed_time_pickoff from pygama.dsp.errors import DSPFatal +from .fixed_time_pickoff import fixed_time_pickoff + @guvectorize(["void(float32[:], float32[:], float32[:])", "void(float64[:], float64[:], float64[:])"], diff --git a/src/pygama/dsp/_processors/multi_t_filter.py b/src/pygama/dsp/processors/multi_t_filter.py similarity index 52% rename from src/pygama/dsp/_processors/multi_t_filter.py rename to src/pygama/dsp/processors/multi_t_filter.py index 4ebd2b81d..8fef16117 100644 --- a/src/pygama/dsp/_processors/multi_t_filter.py +++ b/src/pygama/dsp/processors/multi_t_filter.py @@ -1,33 +1,42 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp._processors.time_point_thresh import time_point_thresh from pygama.dsp.errors import DSPFatal +from .time_point_thresh import time_point_thresh + @guvectorize(["void(float32[:],float32[:],float32[:])", "void(float64[:],float64[:],float64[:])"], "(n),(n) -> (n)", nopython=True, cache=True) -def remove_duplicates(t_in, vt_min_in, t_out): - """ - time_point_thresh has issues with afterpulsing in waveforms that causes an - aferpulse peak's tp0 to be sent to 0 or the same index as the tp0 for the - first pulse. This only happens when the relative minimum between the first - pulse and the afterpulse is greater than the threshold. So, we sweep +def remove_duplicates(t_in: np.ndarray, vt_min_in: np.ndarray, t_out: np.ndarray) -> None: + """Helper function to remove duplicate peak positions. + + :func:`.time_point_thresh` has issues with afterpulsing in waveforms that + causes an aferpulse peak position to be sent to zero or the same index as + the first pulse. This only happens when the relative minimum between the + first pulse and the afterpulse is greater than the threshold. So, we sweep through the array again to ensure there are no duplicate indices. If there - are duplicate indices caused by a misidentified tp0 of an afterpulse, we - replace its index by that of the corresponding minimum found using the - get_multi_local_extrema function. It also checks to make sure that the + are duplicate indices caused by a misidentified position of an afterpulse, + we replace its index by that of the corresponding minimum found using + :func:`.get_multi_local_extrema`. It also checks to make sure that the maximum of a waveform isn't right at index 0. Parameters ---------- - t_in : array-like - The array of indices that we want to remove duplicates from - vt_min_in : array-like - List of indices of minima that we want to replace duplicates in t_out with - t_out : array-like - The array we want to return that will have no duplicate indices in it + t_in + the array of indices that we want to remove duplicates from. + vt_min_in + list of indices of minima that we want to replace duplicates in `t_out` + with. + t_out + the array we want to return that will have no duplicate indices in it. + + See Also + -------- + .multi_t_filter """ # initialize arrays t_out[:] = np.nan @@ -55,32 +64,39 @@ def remove_duplicates(t_in, vt_min_in, t_out): t_out[:] = np.append(t_out[1:],np.nan) - @guvectorize(["void(float32[:], float32[:], float32[:], float32[:], float32[:])", "void(float64[:], float64[:], float64[:], float64[:], float64[:])"], "(n),(),(m),(m),(m)", forceobj=True, cache=True) -def multi_t_filter(w_in, a_threshold_in, vt_max_in, vt_min_in, t_out): - """ - Gets list of indices of the start of leading edges of multiple peaks within a waveform. - Is built to handle afterpulses/delayed cross talk and trains of pulses. - The multi_t_filter works by calling the vectorized functions - "get_multi_local_extrema" which returns a list of the maxima and minima in a waveform, - and then the list of maxima is fed into "time_point_thresh" which returns - the final times that waveform is less than a specified threshold. +def multi_t_filter(w_in: np.ndarray, a_threshold_in: float, + vt_max_in: np.ndarray, vt_min_in: np.ndarray, t_out: np.ndarray) -> None: + """Gets list of indices of the start of leading edges of multiple peaks + within a waveform. + + Is built to handle afterpulses/delayed cross-talk and trains of pulses. + Works by calling the vectorized functions :func:`.get_multi_local_extrema` + which returns a list of the maxima and minima in a waveform, and then the + list of maxima is fed into :func:`.time_point_thresh` which returns the + final times that waveform is less than a specified threshold. Parameters ---------- - w_in : array-like - The array of data within which the list of tp0s will be found - a_threshold_in : scalar - Threshold to search for using time_point_thresh - vt_max_in : array-like - The array of max positions for each wf - vt_min_in : array-like - The array of min positions for each wf - t_out : array-like - Output array of fixed length (padded with nans) that hold the indices of - the identified initial rise times of peaks in the signal + w_in + the array of data within which the list of leading edge times will be + found. + a_threshold_in + threshold to search for using :func:`.time_point_thresh`. + vt_max_in + the array of maximum positions for each waveform. + vt_min_in + the array of minimum positions for each waveform. + t_out + output array of fixed length (padded with :any:`numpy.nan`) that hold + the indices of the identified initial rise times of peaks in the + signal. + + See Also + -------- + ~.time_point_thresh.time_point_thresh """ # initialize arrays, padded with the elements we want @@ -91,7 +107,7 @@ def multi_t_filter(w_in, a_threshold_in, vt_max_in, vt_min_in, t_out): return if (np.isnan(vt_max_in).all() and np.isnan(vt_min_in).all()): return - if (not len(t_out)<=len(w_in)): + if (not len(t_out) <= len(w_in)): raise DSPFatal('The length of your return array must be smaller than the length of your waveform') # Initialize an intermediate array to hold the tp0 values before we remove duplicates from it diff --git a/src/pygama/dsp/_processors/optimize.py b/src/pygama/dsp/processors/optimize.py similarity index 57% rename from src/pygama/dsp/_processors/optimize.py rename to src/pygama/dsp/processors/optimize.py index cfe8082ef..70bda1af2 100644 --- a/src/pygama/dsp/_processors/optimize.py +++ b/src/pygama/dsp/processors/optimize.py @@ -1,3 +1,7 @@ +from __future__ import annotations + +from typing import Callable + import numpy as np from iminuit import Minuit from numba import guvectorize @@ -8,13 +12,13 @@ class Model: - """ - The model class containing the function to minimize. + """A class representing the function to minimize. """ errordef = Minuit.LEAST_SQUARES # the constructor - def __init__(self, func, w_in, baseline, beg, end): + def __init__(self, func: Callable, w_in: np.ndarray, + baseline: float, beg: int, end: int) -> None: self.func = func self.x = np.arange(beg, end) self.y = np.asarray(w_in, dtype=np.float64) - baseline @@ -22,46 +26,51 @@ def __init__(self, func, w_in, baseline, beg, end): self.end = end # the function to minimize - def __call__(self, args): + def __call__(self, args) -> np.ndarray: y_pz = self.func(self.y, *args)[self.beg:self.end] return np.abs(np.sum(self.x) * np.sum(y_pz) - len(self.x) * np.sum(self.x * y_pz)) + @guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])", "void(float64[:], float64, float64, float64, float64, float64[:])"], "(n),(),(),(),()->()", forceobj=True) -def optimize_1pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, val0_out): - """ - Find the optimal, single pole-zero cancellation's parameter +def optimize_1pz(w_in: np.ndarray, a_baseline_in: float, t_beg_in: int, + t_end_in: int, p0_in: float, val0_out: float) -> None: + """Find the optimal, single pole-zero cancellation's parameter by minimizing the slope in the waveform's specified time range. Parameters ---------- - w_in : array-like - The input waveform - a_baseline_in : float - The resting baseline - t_beg_in : int - The lower bound's index for the time range over - which to optimize the pole-zero cancellation - t_end_in : int - The upper bound's index for the time range over - which to optimize the pole-zero cancellation - p0_in : float - The initial guess of the optimal time constant - val0_out : float - The output value of the best-fit time constant - - Examples - -------- + w_in + the input waveform. + a_baseline_in + the resting baseline. + t_beg_in + the lower bound's index for the time range over + which to optimize the pole-zero cancellation. + t_end_in + the upper bound's index for the time range over + which to optimize the pole-zero cancellation. + p0_in + the initial guess of the optimal time constant. + val0_out + the output value of the best-fit time constant. + + JSON Configuration Example + -------------------------- + .. code-block :: json "tau0": { "function": "optimize_1pz", "module": "pygama.dsp.processors", "args": ["waveform", "baseline", "0", "20*us", "500*us", "tau0"], - "prereqs": ["waveform", "baseline"], "unit": "us" } + + See Also + -------- + ~.pole_zero.pole_zero, .double_pole_zero """ val0_out[0] = np.nan @@ -84,48 +93,53 @@ def optimize_1pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, val0_out): m.migrad() val0_out[0] = m.values[0] + @guvectorize(["void(float32[:], float32, float32, float32, float32, float32, float32, float32[:], float32[:], float32[:])", "void(float64[:], float64, float64, float64, float64, float64, float64, float64[:], float64[:], float64[:])"], "(n),(),(),(),(),(),()->(),(),()", forceobj=True) -def optimize_2pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, p1_in, p2_in, val0_out, val1_out, val2_out): - """ - Find the optimal, double pole-zero cancellation's parameters - by minimizing the slope in the waveform's specified time range. +def optimize_2pz(w_in: np.ndarray, a_baseline_in: float, t_beg_in: int, + t_end_in: int, p0_in: float, p1_in: float, p2_in: float, + val0_out: float, val1_out: float, val2_out: float) -> None: + """Find the optimal, double pole-zero cancellation's parameters by + minimizing the slope in the waveform's specified time range. Parameters ---------- - w_in : array-like - The input waveform - a_baseline_in : float - The resting baseline - t_beg_in : int - The lower bound's index for the time range over - which to optimize the pole-zero cancellation - t_end_in : int - The upper bound's index for the time range over - which to optimize the pole-zero cancellation - p0_in : float - The initial guess of the optimal, longer time constant - p1_in : float - The initial guess of the optimal, shorter time constant - p2_in : float - The initial guess of the optimal fraction - val0_out : float - The output value of the best-fit, longer time constant - val1_out : float - The output value of the best-fit, shorter time constant - val2_out : float - The output value of the best-fit fraction - - Examples - -------- + w_in + the input waveform. + a_baseline_in + the resting baseline. + t_beg_in + the lower bound's index for the time range over + which to optimize the pole-zero cancellation. + t_end_in + the upper bound's index for the time range over + which to optimize the pole-zero cancellation. + p0_in + the initial guess of the optimal, longer time constant. + p1_in + the initial guess of the optimal, shorter time constant. + p2_in + the initial guess of the optimal fraction. + val0_out + the output value of the best-fit, longer time constant. + val1_out + the output value of the best-fit, shorter time constant. + val2_out + the output value of the best-fit fraction. + + JSON Configuration Example + -------------------------- + .. code-block :: json "tau1, tau2, frac": { "function": "optimize_2pz", "module": "pygama.dsp.processors", - "args": ["waveform", "baseline", "0", "20*us", "500*us", "20*us", "0.02", "tau1", "tau2", "frac"], - "prereqs": ["waveform", "baseline"], + "args": [ + "waveform", "baseline", "0", "20*us", "500*us", + "20*us", "0.02", "tau1", "tau2", "frac" + ], "unit": "us" } """ diff --git a/src/pygama/dsp/processors/param_lookup.py b/src/pygama/dsp/processors/param_lookup.py new file mode 100644 index 000000000..94ee1a40d --- /dev/null +++ b/src/pygama/dsp/processors/param_lookup.py @@ -0,0 +1,28 @@ +from __future__ import annotations + +import numpy as np +from numba import from_dtype, guvectorize, types + + +def param_lookup(param_dict: dict[int, float], default_val: float, + dtype: str | np.dtype) -> np.ufunc: + """Generate the :class:`numpy.ufunc` ``lookup(channel, val)``, which + returns a NumPy array of values corresponding to various channels that are + looked up in the provided `param_dict`. If there is no key, use + `default_val` instead. + """ + out_type = from_dtype(np.dtype(dtype)) + + # convert types to avoid any necessity of casting + param_dict = {types.uint32(k): out_type(v) for k, v in param_dict.items()} + default_val = out_type(default_val) + + @guvectorize(["void(uint32, " + out_type.name + "[:])"], + "()->()", forceobj=True) + def lookup(channel: int, val: np.ndarray) -> None: + """Look up a value for the provided channel from a dictionary provided + at compile time. + """ + val[0] = param_dict.get(channel, default_val) + + return lookup diff --git a/src/pygama/dsp/_processors/pole_zero.py b/src/pygama/dsp/processors/pole_zero.py similarity index 50% rename from src/pygama/dsp/_processors/pole_zero.py rename to src/pygama/dsp/processors/pole_zero.py index 3b34212d7..c733371ee 100644 --- a/src/pygama/dsp/_processors/pole_zero.py +++ b/src/pygama/dsp/processors/pole_zero.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,30 +9,29 @@ @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),()->(n)", nopython=True, cache=True) -def pole_zero(w_in, t_tau, w_out): - """ - Apply a pole-zero cancellation using the provided time +def pole_zero(w_in: np.ndarray, t_tau: float, w_out: np.ndarray) -> None: + """Apply a pole-zero cancellation using the provided time constant to the waveform. Parameters ---------- - w_in : array-like - The input waveform - t_tau : float - The time constant of the exponential to be deconvolved - w_out : array-like - The pole-zero cancelled waveform - - Examples - -------- + w_in + the input waveform. + t_tau + the time constant of the exponential to be deconvolved. + w_out + the pole-zero cancelled waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_pz": { "function": "pole_zero", "module": "pygama.dsp.processors", "args": ["wf_bl", "400*us", "wf_pz"], - "unit": "ADC", - "prereqs": ["wf_bl"] + "unit": "ADC" } """ w_out[:] = np.nan @@ -43,64 +44,75 @@ def pole_zero(w_in, t_tau, w_out): for i in range(1, len(w_in), 1): w_out[i] = w_out[i-1] + w_in[i] - w_in[i-1] * const + @guvectorize(["void(float32[:], float32, float32, float32, float32[:])", "void(float64[:], float64, float64, float64, float64[:])"], "(n),(),(),()->(n)", nopython=True, cache=True) -def double_pole_zero(w_in, t_tau1, t_tau2, frac, w_out): - """ +def double_pole_zero(w_in: np.ndarray, t_tau1: float, t_tau2: float, + frac: float, w_out: np.ndarray) -> np.ndarray: + r""" Apply a double pole-zero cancellation using the provided time constants to the waveform. Parameters ---------- - w_in : array-like - The input waveform - t_tau1 : float - The time constant of the first exponential to be deconvolved - t_tau2 : float - The time constant of the second exponential to be deconvolved - frac : float - The fraction which the second exponential contributes - w_out : array-like - The pole-zero cancelled waveform - - Examples - -------- + w_in + the input waveform. + t_tau1 + the time constant of the first exponential to be deconvolved. + t_tau2 + the time constant of the second exponential to be deconvolved. + frac + the fraction which the second exponential contributes. + w_out + the pole-zero cancelled waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_pz": { "function": "double_pole_zero", "module": "pygama.dsp.processors", "args": ["wf_bl", "400*us", "20*us", "0.02", "wf_pz"], - "unit": "ADC", - "prereqs": ["wf_bl"] + "unit": "ADC" } Notes ----- This algorithm is an IIR filter to deconvolve the function - .. math:: s(t) = A[frac*\\exp\\left(-\\frac{t}{t_{tau2}}\\right) + (1-frac)*\\exp\\left(-\\frac{t}{t_{tau1}}\\right)] + .. math:: + s(t) = A \left[ f \cdot \exp\left(-\frac{t}{\tau_2} \right) + + (1-f) \cdot \exp\left(-\frac{t}{\tau_1}\right) \right] - into a single step function of amplitude A. - This filter is derived by z-transforming the input (:math:`s(t)`) and output (step function) - signals, and then determining the transfer function. For shorthand, define :math:`a=\\exp\\left(-\\frac{1}{t_{tau1}}\\right)` - and :math:`b=\\exp\\left(-\\frac{1}{t_{tau2}}\\right)`, the transfer function is then: + (:math:`f` = `frac`) into a single step function of amplitude :math:`A`. + This filter is derived by :math:`z`-transforming the input (:math:`s(t)`) + and output (step function) signals, and then determining the transfer + function. For shorthand, define :math:`a=\exp(-1/\tau_1)` and + :math:`b=\exp(-1/\tau_2)`, the transfer function is then: - .. math:: H(z) = \\frac{1-(a+b)z^{-1}+abz^{-2}}{1+(fb-fa-b-1)z^{-1}-(fb-fa-b)z^{-2}} + .. math:: + H(z) = \frac{1 - (a+b)z^{-1} + abz^{-2}} + {1 + (fb - fa - b - 1)z^{-1}-(fb - fa - b)z^{-2}} - By equating the transfer function to the ratio of output to input waveforms :math:`H(z) = \\frac{w_{out}(z)}{w_{in}(z)}` - and then taking the inverse z-transform yields the filter as run in the function, where f is the frac parameter: + By equating the transfer function to the ratio of output to input waveforms + :math:`H(z) = w_\text{out}(z) / w_\text{in}(z)` and then taking the + inverse :math:`z`-transform yields the filter as run in the function, where + :math:`f` is the `frac` parameter: .. math:: - w_{out}[n] = &w_{in}[n]-(a+b)w_{in}[n-1]+(ab)w_{in}[n-2]\\\\ - &-(fb-fa-b-1)w_{out}[n-1]+(fb-fa-b)w_{out}[n-2] + w_\text{out}[n] =& w_\text{in}[n] - (a+b)w_\text{in}[n-1] + + abw_\text{in}[n-2] \\ + & -(fb - fa - b - 1)w_\text{out}[n-1] + + (fb - fa - b)w_\text{out}[n-2] """ w_out[:] = np.nan if np.isnan(w_in).any() or np.isnan(t_tau1) or np.isnan(t_tau2) or np.isnan(frac): return - if (len(w_in)<=3): + if (len(w_in) <= 3): raise DSPFatal('The length of the waveform must be larger than 3 for the filter to work safely') a = np.exp(-1 / t_tau1) diff --git a/src/pygama/dsp/_processors/presum.py b/src/pygama/dsp/processors/presum.py similarity index 60% rename from src/pygama/dsp/_processors/presum.py rename to src/pygama/dsp/processors/presum.py index c1e6430ff..489fb665e 100644 --- a/src/pygama/dsp/_processors/presum.py +++ b/src/pygama/dsp/processors/presum.py @@ -1,24 +1,24 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize -from pygama.dsp.errors import DSPFatal - @guvectorize(["void(float32[:], float32[:])", "void(float64[:], float64[:])"], "(n),(m)", nopython=True, cache=True) -def presum(w_in, w_out): - """ - Presum the waveform. Combine bins in chunks of len(w_in) / len(w_out), - which is hopefully an integer. If it isn't, then some samples at the end - will be omitted. +def presum(w_in: np.ndarray, w_out: np.ndarray) -> None: + """Presum the waveform. + + Combine bins in chunks of ``len(w_in) / len(w_out)``, which is hopefully an + integer. If it isn't, then some samples at the end will be omitted. Parameters ---------- - w_in : array-like - The input waveform - w_out : array-like - The output waveform + w_in + the input waveform. + w_out + the output waveform. """ w_out[:] = np.nan diff --git a/src/pygama/dsp/processors/pulse_injector.py b/src/pygama/dsp/processors/pulse_injector.py new file mode 100644 index 000000000..8cfeef4b3 --- /dev/null +++ b/src/pygama/dsp/processors/pulse_injector.py @@ -0,0 +1,81 @@ +from __future__ import annotations + +from math import exp, log + +import numpy as np +from numba import guvectorize + + +@guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])", + "void(float64[:], float64, float64,float64,float64,float64[:])"], + "(n),(),(),(), ()->(n)", nopython=True, cache=True) +def inject_sig_pulse(wf_in: np.ndarray, t0: int, rt: float, + A: float, decay: float, wf_out: np.ndarray) -> None: + r"""Inject sigmoid pulse into existing waveform to simulate pileup. + + .. math:: + s(t) = \frac{A}{1 + \exp[-4 \log(99) (t - t_0 - t_r/2) / t_r]} + e^{-(t-t_0)/\tau} + + Parameters + ---------- + wf_in + the input waveform. + t0 + the position :math:`t_0` of the injected waveform. + rt + the rise time :math:`t_r` of the injected waveform. + A + the amplitude :math:`A` of the injected waveform. + decay + the decay parameter :math:`\tau` of the injected waveform. + wf_out + the output waveform. + """ + + wf_out[:] = np.nan + + if np.isnan(wf_in).any() or np.isnan(rt) or np.isnan(t0) or np.isnan(A) or np.isnan(decay): + return + + rise = 4*log(99)/rt + + wf_out[:] = wf_in[:] + for T in range(len(wf_out)): + wf_out[T] = wf_out[T] + A/(1+exp(-rise*(T-(t0+rt/2))))*exp(-(1/decay)*(T-t0)) + + +@guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])", + "void(float64[:], float64, float64,float64,float64,float64[:])"], + "(n),(),(),(), ()->(n)", nopython=True, cache=True) +def inject_exp_pulse(wf_in: np.ndarray, t0: int, rt: float, A: float, + decay: float, wf_out: np.ndarray) -> None: + """Inject exponential pulse into existing waveform to simulate pileup. + + Parameters + ---------- + wf_in + the input waveform. + t0 + the position of the injected waveform. + rt + the rise time of the injected waveform. + A + the amplitude of the injected waveform. + decay + the exponential decay constant of the injected waveform. + wf_out + the output waveform. + """ + + wf_out[:] = np.nan + + if np.isnan(wf_in).any() or np.isnan(rt) or np.isnan(t0) or np.isnan(A) or np.isnan(decay): + return + + wf_out[:] = wf_in[:] + for T in range(len(wf_out)): + if T <= t0 and T <= t0+rt: + wf_out[T] += (A*exp((T-t0-rt)/(rt))*exp(-(1/decay)*(T-t0))) + elif (T > t0 + rt): + wf_out[T] += (A*exp(-(1/decay)*(T-t0))) diff --git a/src/pygama/dsp/_processors/saturation.py b/src/pygama/dsp/processors/saturation.py similarity index 64% rename from src/pygama/dsp/_processors/saturation.py rename to src/pygama/dsp/processors/saturation.py index 9f800ca93..700c84a84 100644 --- a/src/pygama/dsp/_processors/saturation.py +++ b/src/pygama/dsp/processors/saturation.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,33 +9,31 @@ @guvectorize(["void(float32[:], float32, float32[:], float32[:])", "void(float64[:], float64, float64[:], float64[:])"], "(n),()->(),()", nopython=True, cache=True) -def saturation(w_in, bit_depth_in, n_lo_out, n_hi_out): - """ - Count the number of samples in the waveform that are - saturated at the minimum and maximum possible values based - on the bit depth. +def saturation(w_in: np.ndarray, bit_depth_in: int, n_lo_out: int, n_hi_out: int) -> None: + """Count the number of samples in the waveform that are saturated at the + minimum and maximum possible values based on the bit depth. Parameters ---------- - w_in : array-like - The input waveform - bit_depth_in : int - The bit depth of the analog-to-digital converter - n_lo_out : int - The output number of samples at the minimum - n_hi_out : int - The output number of samples at the maximum - - Examples - -------- + w_in + the input waveform + bit_depth_in + the bit depth of the analog-to-digital converter + n_lo_out + the output number of samples at the minimum + n_hi_out + the output number of samples at the maximum + + JSON Configuration Example + -------------------------- + .. code-block :: json "sat_lo, sat_hi": { "function": "saturation", "module": "pygama.dsp.processors", "args": ["waveform", "16", "sat_lo", "sat_hi"], - "unit": "ADC", - "prereqs": ["waveform"] + "unit": "ADC" } """ n_lo_out[0] = np.nan diff --git a/src/pygama/dsp/_processors/soft_pileup_corr.py b/src/pygama/dsp/processors/soft_pileup_corr.py similarity index 69% rename from src/pygama/dsp/_processors/soft_pileup_corr.py rename to src/pygama/dsp/processors/soft_pileup_corr.py index f67599933..a847b4d80 100644 --- a/src/pygama/dsp/_processors/soft_pileup_corr.py +++ b/src/pygama/dsp/processors/soft_pileup_corr.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,34 +9,33 @@ @guvectorize(["void(float32[:], float32, float32, float32[:])", "void(float64[:], float64, float64, float64[:])"], "(n),(),()->(n)", nopython=True, cache=True) -def soft_pileup_corr(w_in, n_in, tau_in, w_out): - """ - Fit the baseline to an exponential with the provided time - constant and then subtract the best-fit function from the - entire waveform. +def soft_pileup_corr(w_in: np.ndarray, n_in: int, + tau_in: float, w_out: np.ndarray) -> None: + """Fit the baseline to an exponential with the provided time constant and + then subtract the best-fit function from the entire waveform. Parameters ---------- - w_in : array-like - The input waveform - n_in : int + w_in + The input waveform. + n_in The number of samples at the beginning of the waveform - to fit to an exponential - tau_in : float - The fixed, exponential time constant - w_out : array-like - The output waveform with the exponential subtracted - - Examples - -------- + to fit to an exponential. + tau_in + The fixed, exponential time constant. + w_out + The output waveform with the exponential subtracted. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_bl": { "function": "soft_pileup_corr", "module": "pygama.dsp.processors", "args": ["waveform", "1000", "500*us", "wf_bl"], - "unit": "ADC", - "prereqs": ["waveform"] + "unit": "ADC" } """ w_out[:] = np.nan @@ -67,39 +68,40 @@ def soft_pileup_corr(w_in, n_in, tau_in, w_out): for i in range(0, len(w_in), 1): w_out[i] = w_in[i] - (A * np.exp(-1.0 * i / tau_in) + B) + @guvectorize(["void(float32[:], float32, float32, float32, float32[:])", "void(float64[:], float64, float64, float64, float64[:])"], "(n),(),(),()->(n)", nopython=True, cache=True) -def soft_pileup_corr_bl(w_in, n_in, tau_in, B_in, w_out): - """ - Fit the baseline to an exponential with the provided time +def soft_pileup_corr_bl(w_in: np.ndarray, n_in: int, tau_in: float, + B_in: float, w_out: np.ndarray) -> None: + """Fit the baseline to an exponential with the provided time constant and then subtract the best-fit function from the entire waveform. Parameters ---------- - w_in : array-like - The input waveform - n_in : int - The number of samples at the beginning of the waveform - to fit to an exponential - tau_in : float - The fixed, exponential time constant - B_in : float - The fixed, exponential constant - w_out : array-like - The output waveform with the exponential subtracted - - Examples - -------- + w_in + the input waveform. + n_in + the number of samples at the beginning of the waveform + to fit to an exponential. + tau_in + the fixed, exponential time constant. + B_in + the fixed, exponential constant. + w_out + the output waveform with the exponential subtracted. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_bl": { "function": "soft_pileup_corr_bl", "module": "pygama.dsp.processors", "args": ["waveform", "1000", "500*us", "baseline", "wf_bl"], - "unit": "ADC", - "prereqs": ["waveform", "baseline"] + "unit": "ADC" } """ w_out[:] = np.nan diff --git a/src/pygama/dsp/_processors/time_point_thresh.py b/src/pygama/dsp/processors/time_point_thresh.py similarity index 67% rename from src/pygama/dsp/_processors/time_point_thresh.py rename to src/pygama/dsp/processors/time_point_thresh.py index 93b65bb48..5fdbb6be9 100644 --- a/src/pygama/dsp/_processors/time_point_thresh.py +++ b/src/pygama/dsp/processors/time_point_thresh.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,34 +9,34 @@ @guvectorize(["void(float32[:], float32, float32, float32, float32[:])", "void(float64[:], float64, float64, float64, float64[:])"], "(n),(),(),()->()", nopython=True, cache=True) -def time_point_thresh(w_in, a_threshold, t_start, walk_forward, t_out): - """ - Find the index where the waveform value crosses the threshold, - walking either forward or backward from the starting index. +def time_point_thresh(w_in: np.ndarray, a_threshold: float, t_start: int, + walk_forward: int, t_out: float) -> None: + """Find the index where the waveform value crosses the threshold, walking + either forward or backward from the starting index. Parameters ---------- - w_in : array-like - The input waveform - a_threshold : float - The threshold value - t_start : int - The starting index - walk_forward : int - The backward (0) or forward (1) search direction - t_out : float - The index where the waveform value crosses the threshold - - Examples - -------- + w_in + the input waveform. + a_threshold + the threshold value. + t_start + the starting index. + walk_forward + the backward (``0``) or forward (``1``) search direction. + t_out + the index where the waveform value crosses the threshold. + + JSON Configuration Example + -------------------------- + .. code-block :: json "tp_0": { "function": "time_point_thresh", "module": "pygama.dsp.processors", "args": ["wf_atrap", "bl_std", "tp_start", 0, "tp_0"], - "unit": "ns", - "prereqs": ["wf_atrap", "bl_std", "tp_start"] + "unit": "ns" } """ t_out[0] = np.nan diff --git a/src/pygama/dsp/_processors/trap_filters.py b/src/pygama/dsp/processors/trap_filters.py similarity index 71% rename from src/pygama/dsp/_processors/trap_filters.py rename to src/pygama/dsp/processors/trap_filters.py index d2d810f11..086d4a9a5 100644 --- a/src/pygama/dsp/_processors/trap_filters.py +++ b/src/pygama/dsp/processors/trap_filters.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,31 +9,30 @@ @guvectorize(["void(float32[:], int32, int32, float32[:])", "void(float64[:], int32, int32, float64[:])"], "(n),(),()->(n)", nopython=True, cache=True) -def trap_filter(w_in, rise, flat, w_out): - """ - Apply a symmetric trapezoidal filter to the waveform. +def trap_filter(w_in: np.ndarray, rise: int, flat: int, w_out: np.ndarray) -> None: + """Apply a symmetric trapezoidal filter to the waveform. Parameters ---------- - w_in : array-like - The input waveform - rise : int - The number of samples averaged in the rise and fall sections - flat : int - The delay between the rise and fall sections - w_out : array-like - The filtered waveform - - Examples - -------- + w_in + the input waveform. + rise + the number of samples averaged in the rise and fall sections. + flat + the delay between the rise and fall sections. + w_out + the filtered waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_tf": { "function": "trap_filter", "module": "pygama.dsp.processors", "args": ["wf_pz", "10*us", "3*us", "wf_tf"], - "unit": "ADC", - "prereqs": ["wf_pz"] + "unit": "ADC" } """ w_out[:] = np.nan @@ -58,35 +59,35 @@ def trap_filter(w_in, rise, flat, w_out): for i in range(2 * rise + flat, len(w_in), 1): w_out[i] = w_out[i-1] + w_in[i] - w_in[i-rise] - w_in[i-rise-flat] + w_in[i-2*rise-flat] + @guvectorize(["void(float32[:], int32, int32, float32[:])", "void(float64[:], int32, int32, float64[:])"], "(n),(),()->(n)", nopython=True, cache=True) -def trap_norm(w_in, rise, flat, w_out): - """ - Apply a symmetric trapezoidal filter to the waveform, normalized - by the number of samples averaged in the rise and fall sections. +def trap_norm(w_in: np.ndarray, rise: int, flat: int, w_out: np.ndarray) -> None: + """Apply a symmetric trapezoidal filter to the waveform, normalized by the + number of samples averaged in the rise and fall sections. Parameters ---------- - w_in : array-like - The input waveform - rise : int - The number of samples averaged in the rise and fall sections - flat : int - The delay between the rise and fall sections - w_out : array-like - The normalized, filtered waveform - - Examples - -------- + w_in + the input waveform. + rise + the number of samples averaged in the rise and fall sections. + flat + the delay between the rise and fall sections. + w_out + the normalized, filtered waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_tf": { "function": "trap_norm", "module": "pygama.dsp.processors", "args": ["wf_pz", "10*us", "3*us", "wf_tf"], - "unit": "ADC", - "prereqs": ["wf_pz"] + "unit": "ADC" } """ w_out[:] = np.nan @@ -113,37 +114,38 @@ def trap_norm(w_in, rise, flat, w_out): for i in range(2 * rise + flat, len(w_in), 1): w_out[i] = w_out[i-1] + (w_in[i] - w_in[i-rise] - w_in[i-rise-flat] + w_in[i-2*rise-flat]) / rise + @guvectorize(["void(float32[:], int32, int32, int32, float32[:])", "void(float64[:], int32, int32, int32, float64[:])"], "(n),(),(),()->(n)", nopython=True, cache=True) -def asym_trap_filter(w_in, rise, flat, fall, w_out): - """ - Apply an asymmetric trapezoidal filter to the waveform, normalized - by the number of samples averaged in the rise and fall sections. +def asym_trap_filter(w_in: np.ndarray, rise: int, flat: int, + fall: int, w_out: np.ndarray) -> None: + """Apply an asymmetric trapezoidal filter to the waveform, normalized by + the number of samples averaged in the rise and fall sections. Parameters ---------- - w_in : array-like - The input waveform - rise : int - The number of samples averaged in the rise section - flat : int - The delay between the rise and fall sections - fall : int - The number of samples averaged in the fall section - w_out : array-like - The normalized, filtered waveform - - Examples - -------- + w_in + the input waveform. + rise + the number of samples averaged in the rise section. + flat + the delay between the rise and fall sections. + fall + the number of samples averaged in the fall section. + w_out + the normalized, filtered waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_af": { "function": "asym_trap_filter", "module": "pygama.dsp.processors", "args": ["wf_pz", "128*ns", "64*ns", "2*us", "wf_af"], - "unit": "ADC", - "prereqs": ["wf_pz"] + "unit": "ADC" } """ w_out[:] = np.nan @@ -173,38 +175,39 @@ def asym_trap_filter(w_in, rise, flat, fall, w_out): for i in range(rise + flat + fall, len(w_in), 1): w_out[i] = w_out[i-1] + (w_in[i] - w_in[i-rise]) / rise - (w_in[i-rise-flat] - w_in[i-rise-flat-fall]) / fall + @guvectorize(["void(float32[:], int32, int32, float32, float32[:])", "void(float64[:], int32, int32, float64, float64[:])"], "(n),(),(),()->()", nopython=True, cache=True) -def trap_pickoff(w_in, rise, flat, t_pickoff, a_out): - """ - Pick off the value at the provided index of a symmetric trapezoidal +def trap_pickoff(w_in: np.ndarray, rise: int, flat: int, + t_pickoff: float, a_out: float) -> None: + """Pick off the value at the provided index of a symmetric trapezoidal filter to the input waveform, normalized by the number of samples averaged in the rise and fall sections. Parameters ---------- - w_in : array-like - The input waveform - rise : int - The number of samples averaged in the rise and fall sections - flat : int - The delay between the rise and fall sections - t_pickoff : float - The waveform index to pick off - a_out : float - The output pick-off value of the filtered waveform - - Examples - -------- + w_in + the input waveform. + rise + the number of samples averaged in the rise and fall sections. + flat + the delay between the rise and fall sections. + t_pickoff + the waveform index to pick off. + a_out + the output pick-off value of the filtered waveform. + + JSON Configuration Example + -------------------------- + .. code-block :: json "ct_corr": { "function": "trap_pickoff", "module": "pygama.dsp.processors", "args": ["wf_pz", "1.5*us", 0, "tp_0", "ct_corr"], - "unit": "ADC", - "prereqs": ["wf_pz", "tp_0"] + "unit": "ADC" } """ a_out[0] = np.nan diff --git a/src/pygama/dsp/_processors/upsampler.py b/src/pygama/dsp/processors/upsampler.py similarity index 66% rename from src/pygama/dsp/_processors/upsampler.py rename to src/pygama/dsp/processors/upsampler.py index a9dfea0e0..1d122c30e 100644 --- a/src/pygama/dsp/_processors/upsampler.py +++ b/src/pygama/dsp/processors/upsampler.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,20 +9,18 @@ @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),(),(m)", nopython=True, cache=True) - -def upsampler(w_in, upsample, w_out): - """ - Upsamples the waveform by the number specified, - a series of moving windows should be applied afterwards for smoothing +def upsampler(w_in: np.ndarray, upsample: float, w_out: np.ndarray) -> None: + """Upsamples the waveform by the number specified, a series of moving + windows should be applied afterwards for smoothing. Parameters ---------- - w_in : array-like + w_in waveform to upsample. - upsample : float - Number of samples to increase each sample to. - w_out : array-like - Output array for upsampled waveform + upsample + number of samples to increase each sample to. + w_out + output array for upsampled waveform. """ w_out[:] = np.nan diff --git a/src/pygama/dsp/_processors/wiener_filter.py b/src/pygama/dsp/processors/wiener_filter.py similarity index 74% rename from src/pygama/dsp/_processors/wiener_filter.py rename to src/pygama/dsp/processors/wiener_filter.py index 9bec0a18a..49eea8270 100644 --- a/src/pygama/dsp/_processors/wiener_filter.py +++ b/src/pygama/dsp/processors/wiener_filter.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -5,21 +7,28 @@ from pygama.dsp.errors import DSPFatal -def wiener_filter(file_name_array): - """Apply a wiener filter to the waveform. +def wiener_filter(file_name_array: list[str]) -> np.ndarray: + """Apply a Wiener filter to the waveform. - Note that this convolution is performed in the frequency domain + Note + ---- + The convolution is performed in the frequency domain. This processor is + composed of a factory function that is called using the `init_args` + argument. The input and output waveforms are passed using `args`. The input + must be the Fourier transform of the waveform. The output is the filtered + waveform in the frequency domain. Parameters ---------- - file_name_array : string - Array with path to an lh5 file containing the time domain version - of the superpulse in one column and noise waveform in another, - the superpulse group must be titled 'spms/processed/superpulse' and - the noise waveform must be called 'spms/processed/noise_wf' - - Examples - -------- + file_name_array + Array with path to an HDF5 file containing the time domain version of + the superpulse in one column and noise waveform in another, the + superpulse HDF5 group must be titled ``spms/processed/superpulse`` and + the noise waveform must be called ``spms/processed/noise_wf``. + + JSON Configuration Example + -------------------------- + .. code-block :: json "wf_wiener": { @@ -27,7 +36,6 @@ def wiener_filter(file_name_array): "module": "pygama.dsp.processors", "args": ["wf_bl_fft", "wf_wiener(2000,f)"], "unit": "dB", - "prereqs": ["wf_bl_fft"], "init_args": ["/path/to/file/wiener.lh5"] } """ @@ -97,14 +105,14 @@ def PSF(superpulse, fft_superpulse): @guvectorize(["void(complex64[:], complex64[:])", "void(complex128[:], complex128[:])"], "(n)->(n)", forceobj=True) - def wiener_out(fft_w_in, fft_w_out): + def wiener_out(fft_w_in: np.ndarray, fft_w_out: np.ndarray) -> None: """ Parameters ---------- - fft_w_in : array-like - The fourier transform input waveform - fft_w_out : array-like - The filtered waveform, in the frequency domain + fft_w_in + the Fourier transformed input waveform. + fft_w_out + the filtered waveform, in the frequency domain. """ fft_w_out[:] = np.nan diff --git a/src/pygama/dsp/_processors/windower.py b/src/pygama/dsp/processors/windower.py similarity index 57% rename from src/pygama/dsp/_processors/windower.py rename to src/pygama/dsp/processors/windower.py index feae15521..3a9b272a8 100644 --- a/src/pygama/dsp/_processors/windower.py +++ b/src/pygama/dsp/processors/windower.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np from numba import guvectorize @@ -7,23 +9,25 @@ @guvectorize(["void(float32[:], float32, float32[:])", "void(float64[:], float64, float64[:])"], "(n),(),(m)", nopython=True, cache=True) -def windower(w_in, t0_in, w_out): - """ - Return a shorter sample of the waveform, starting at the - specified index. Note that the length of the output waveform - is determined by the length of "w_out" rather than an input - parameter. If the the length of "w_out" plus "t0_in" extends - past the end of "w_in" or if "t0_in" is negative, remaining - values are padded with NaN. +def windower(w_in: np.ndarray, t0_in: int, w_out: np.ndarray) -> None: + """Return a shorter sample of the waveform, starting at the + specified index. + + Note + ---- + The length of the output waveform is determined by the length of `w_out` + rather than an input parameter. If the the length of `w_out` plus `t0_in` + extends past the end of `w_in` or if `t0_in` is negative, remaining values + are padded with :any:`numpy.nan`. Parameters ---------- - w_in : array-like - The input waveform - t0_in : int - The starting index of the window - w_out : array-like - The windowed waveform + w_in + the input waveform. + t0_in + the starting index of the window. + w_out + the windowed waveform. """ w_out[:] = np.nan diff --git a/src/pygama/math/peak_fitting.py b/src/pygama/math/peak_fitting.py index eb524a388..bf2fb44ce 100644 --- a/src/pygama/math/peak_fitting.py +++ b/src/pygama/math/peak_fitting.py @@ -303,8 +303,7 @@ def gauss_mode_width_max(hist, bins, var=None, mode_guess=None, n_bins=5, Returns ------- - pars : ndarray containing the parameters (mode, sigma, maximum) of the - gaussian fit + pars : ndarray containing the parameters (mode, sigma, maximum) of the gaussian fit - mode : the estimated x-position of the maximum - sigma : the estimated width of the peak. Equivalent to a gaussian width (sigma), but based only on the curvature within n_bins of diff --git a/src/pygama/raw/build_raw.py b/src/pygama/raw/build_raw.py index e6fad99cd..5befc955a 100644 --- a/src/pygama/raw/build_raw.py +++ b/src/pygama/raw/build_raw.py @@ -22,7 +22,7 @@ def build_raw(in_stream: int, in_stream_type: str = None, out_spec: str | dict | RawBufferLibrary = None, buffer_size: int = 8192, - n_max: int = np.inf, overwrite: bool = True, **kwargs) -> None: + n_max: int = np.inf, overwrite: bool = False, **kwargs) -> None: """Convert data into LEGEND HDF5 raw-tier format. Takes an input stream of a given type and writes to output file(s) @@ -171,7 +171,7 @@ def build_raw(in_stream: int, in_stream_type: str = None, if len(out_file_glob) > 1: raise RuntimeError(f'got multiple matches for out_file {out_file}: {out_file_glob}') if not overwrite: - raise RuntimeError(f'file {out_file_glob[0]} exists. Use option overwrite to proceed.') + raise FileExistsError(f'file {out_file_glob[0]} exists. Use option overwrite to proceed.') os.remove(out_file_glob[0]) diff --git a/tests/dsp/configs/icpc-dsp-config.json b/tests/dsp/configs/icpc-dsp-config.json new file mode 100644 index 000000000..0168da525 --- /dev/null +++ b/tests/dsp/configs/icpc-dsp-config.json @@ -0,0 +1,327 @@ +{ + "outputs": [ + "tp_min", + "tp_max", + "wf_min", + "wf_max", + "bl_mean", + "bl_std", + "bl_slope", + "bl_intercept", + "pz_slope", + "pz_std", + "pz_mean", + "trapTmax", + "tp_0_est", + "tp_0_atrap", + "tp_10", + "tp_20", + "tp_50", + "tp_80", + "tp_90", + "tp_99", + "tp_100", + "tp_01", + "tp_95", + "A_max", + "QDrift", + "dt_eff", + "tp_aoe_max", + "tp_aoe_samp", + "trapEmax", + "trapEftp", + "cuspEmax", + "zacEmax", + "zacEftp", + "cuspEftp" + ], + "processors": { + "tp_min, tp_max, wf_min, wf_max": { + "function": "min_max", + "module": "pygama.dsp.processors", + "args": ["waveform", "tp_min", "tp_max", "wf_min", "wf_max"], + "unit": ["ns", "ns", "ADC", "ADC"] + }, + "wf_blsub": { + "function": "bl_subtract", + "module": "pygama.dsp.processors", + "args": ["waveform", "baseline", "wf_blsub"], + "unit": "ADC" + }, + "bl_mean , bl_std, bl_slope, bl_intercept": { + "function": "linear_slope_fit", + "module": "pygama.dsp.processors", + "args": [ + "wf_blsub[0:750]", + "bl_mean", + "bl_std", + "bl_slope", + "bl_intercept" + ], + "unit": ["ADC", "ADC", "ADC", "ADC"] + }, + "wf_pz": { + "function": "pole_zero", + "module": "pygama.dsp.processors", + "args": ["wf_blsub", "db.pz.tau", "wf_pz"], + "unit": "ADC", + "defaults": { "db.pz.tau": "27460.5" } + }, + "pz_mean , pz_std, pz_slope, pz_intercept": { + "function": "linear_slope_fit", + "module": "pygama.dsp.processors", + "args": ["wf_pz[1500:]", "pz_mean", "pz_std", "pz_slope", "pz_intercept"], + "unit": ["ADC", "ADC", "ADC", "ADC"] + }, + "wf_t0_filter": { + "function": "t0_filter", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "wf_t0_filter(len(wf_pz), 'f', grid=wf_pz.grid)"], + "init_args": ["128*ns/wf_pz.period", "2*us/wf_pz.period"], + "unit": "ADC" + }, + "wf_atrap": { + "function": "asym_trap_filter", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "128*ns", "4", "2*us", "wf_atrap"], + "unit": "ADC" + }, + "conv_tmin ,tp_start, conv_min, conv_max": { + "function": "min_max", + "module": "pygama.dsp.processors", + "args": ["wf_t0_filter", "conv_tmin", "tp_start", "conv_min", "conv_max"], + "unit": ["ns", "ns", "ADC", "ADC"] + }, + "tp_0_atrap": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_atrap", "bl_std", "tp_start", 0, "tp_0_atrap"], + "unit": "ns" + }, + "tp_0_est": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_t0_filter", "bl_std", "tp_start", 0, "tp_0_est(unit=ns)"], + "unit": "ns" + }, + "wf_trap": { + "function": "trap_norm", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "db.ttrap.rise", "db.ttrap.flat", "wf_trap"], + "unit": "ADC", + "defaults": { "db.ttrap.rise": "10*us", "db.ttrap.flat": "3.008*us" } + }, + "trapTmax": { + "function": "amax", + "module": "numpy", + "args": ["wf_trap", 1, "trapTmax"], + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, + "unit": "ADC" + }, + "wf_etrap": { + "function": "trap_norm", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "db.etrap.rise", "db.etrap.flat", "wf_etrap"], + "unit": "ADC", + "defaults": { "db.etrap.rise": "10*us", "db.etrap.flat": "3.008*us" } + }, + "trapEmax": { + "function": "amax", + "module": "numpy", + "args": ["wf_etrap", 1, "trapEmax"], + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, + "unit": "ADC" + }, + "trapEftp": { + "function": "fixed_time_pickoff", + "module": "pygama.dsp.processors", + "args": [ + "wf_etrap", + "tp_0_est+db.etrap.rise+db.etrap.flat*db.etrap.sample", + "trapEftp" + ], + "unit": "ADC", + "defaults": { + "db.etrap.rise": "10*us", + "db.etrap.flat": "3*us", + "db.etrap.sample": "0.8" + } + }, + "wf_cusp": { + "function": "cusp_filter", + "module": "pygama.dsp.processors", + "args": ["wf_blsub", "wf_cusp(101, 'f')"], + "init_args": [ + "len(wf_blsub)-100", + "db.cusp.sigma/wf_blsub.period", + "db.cusp.flat/wf_blsub.period", + "db.pz.tau" + ], + "defaults": { + "db.cusp.sigma": "20*us", + "db.cusp.flat": "3*us", + "db.pz.tau": "27460.5" + }, + "unit": "ADC" + }, + "cuspEmax": { + "function": "amax", + "module": "numpy", + "args": ["wf_cusp", 1, "cuspEmax"], + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, + "unit": "ADC" + }, + "cuspEftp": { + "function": "fixed_time_pickoff", + "module": "pygama.dsp.processors", + "args": ["wf_cusp", "db.cusp.sample", "cuspEftp"], + "unit": "ADC", + "defaults": { "db.cusp.sample": "50" } + }, + "wf_zac": { + "function": "zac_filter", + "module": "pygama.dsp.processors", + "args": ["wf_blsub", "wf_zac(101, 'f')"], + "init_args": [ + "len(wf_blsub)-100", + "db.zac.sigma/wf_blsub.period", + "db.zac.flat/wf_blsub.period", + "db.pz.tau" + ], + "defaults": { + "db.zac.sigma": "20*us", + "db.zac.flat": "3*us", + "db.pz.tau": "27460.5" + }, + "unit": "ADC" + }, + "zacEmax": { + "function": "amax", + "module": "numpy", + "args": ["wf_zac", 1, "zacEmax"], + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, + "unit": "ADC" + }, + "zacEftp": { + "function": "fixed_time_pickoff", + "module": "pygama.dsp.processors", + "args": ["wf_zac", "db.zac.sample", "zacEftp"], + "defaults": { "db.zac.sample": "50" }, + "unit": "ADC" + }, + "tp_100": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax", "tp_0_est", 1, "tp_100"], + "unit": "ns" + }, + "tp_99": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "0.99*trapTmax", "tp_0_est", 1, "tp_99"], + "unit": "ns" + }, + "tp_95": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.95", "tp_99", 0, "tp_95"], + "unit": "ns" + }, + "tp_90": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.9", "tp_95", 0, "tp_90"], + "unit": "ns" + }, + "tp_80": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.8", "tp_90", 0, "tp_80"], + "unit": "ns" + }, + "tp_50": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.5", "tp_80", 0, "tp_50"], + "unit": "ns" + }, + "tp_20": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.2", "tp_50", 0, "tp_20"], + "unit": "ns" + }, + "tp_10": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.1", "tp_20", 0, "tp_10"], + "unit": "ns" + }, + "tp_01": { + "function": "time_point_thresh", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "trapTmax*0.01", "tp_10", 0, "tp_01"], + "unit": "ns" + }, + "wf_trap2": { + "function": "trap_norm", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "4*us", "96*ns", "wf_trap2"], + "unit": "ADC" + }, + "trapQftp": { + "function": "fixed_time_pickoff", + "module": "pygama.dsp.processors", + "args": ["wf_trap2", "tp_0_est + 8.096*us", "trapQftp"], + "unit": "ADC" + }, + "QDrift": { + "function": "multiply", + "module": "numpy", + "args": ["trapQftp", "(4*us)/waveform.period", "QDrift"], + "unit": "ADC" + }, + "dt_eff": { + "function": "divide", + "module": "numpy", + "args": ["QDrift", "trapTmax", "dt_eff"], + "unit": "ns" + }, + "wf_le": { + "function": "windower", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "tp_0_est", "wf_le(301, 'f')"], + "unit": "ADC" + }, + "curr": { + "function": "avg_current", + "module": "pygama.dsp.processors", + "args": ["wf_le", 1, "curr(len(wf_le)-1, 'f')"], + "unit": "ADC/sample" + }, + "curr_up": { + "function": "upsampler", + "module": "pygama.dsp.processors", + "args": ["curr", "16", "curr_up(4784, 'f')"], + "unit": "ADC/sample" + }, + "curr_av": { + "function": "moving_window_multi", + "module": "pygama.dsp.processors", + "args": ["curr_up", "48", 3, 0, "curr_av"], + "unit": "ADC/sample" + }, + "aoe_t_min, tp_aoe_max, A_min, A_max": { + "function": "min_max", + "module": "pygama.dsp.processors", + "args": ["curr_av", "aoe_t_min", "tp_aoe_max", "A_min", "A_max"], + "unit": ["ns", "ns", "ADC/sample", "ADC/sample"] + }, + "tp_aoe_samp": { + "function": "add", + "module": "numpy", + "args": ["tp_0_est", "tp_aoe_max/16", "tp_aoe_samp"], + "unit": "ns" + } + } +} diff --git a/tests/dsp/configs/sipm-dsp-config.json b/tests/dsp/configs/sipm-dsp-config.json new file mode 100644 index 000000000..3b015024c --- /dev/null +++ b/tests/dsp/configs/sipm-dsp-config.json @@ -0,0 +1,23 @@ +{ + "outputs": ["bl_mean", "bl_std"], + "processors": { + "wf_blsub": { + "function": "bl_subtract", + "module": "pygama.dsp.processors", + "args": ["waveform", "baseline", "wf_blsub"], + "unit": "ADC" + }, + "bl_mean , bl_std, bl_slope, bl_intercept": { + "function": "linear_slope_fit", + "module": "pygama.dsp.processors", + "args": [ + "wf_blsub[0:750]", + "bl_mean", + "bl_std", + "bl_slope", + "bl_intercept" + ], + "unit": ["ADC", "ADC", "ADC", "ADC"] + } + } +} diff --git a/tests/dsp/conftest.py b/tests/dsp/conftest.py new file mode 100644 index 000000000..93bb01bb3 --- /dev/null +++ b/tests/dsp/conftest.py @@ -0,0 +1,36 @@ +import pytest + +import pygama.dsp.processors # noqa: F401 +from pygama.lgdo import LH5Store +from pygama.raw.build_raw import build_raw + + +@pytest.fixture(scope='session') +def geds_raw_tbl(lgnd_test_data): + store = LH5Store() + obj, _ = store.read_object( + '/geds/raw', + lgnd_test_data.get_path('lh5/LDQTA_r117_20200110T105115Z_cal_geds_raw.lh5'), + n_rows=10) + return obj + + +@pytest.fixture(scope='session') +def spms_raw_tbl(lgnd_test_data): + + out_file = '/tmp/L200-comm-20211130-phy-spms.lh5' + out_spec = { + 'FCEventDecoder': { + 'raw': { + 'key_list': [[0, 6]], + 'out_stream': f"{out_file}:/spms" + } + } + } + + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=out_spec, n_max=10, overwrite=True) + + store = LH5Store() + obj, _ = store.read_object('/spms/raw', out_file, n_rows=10) + return obj diff --git a/tests/dsp/test_build_dsp.py b/tests/dsp/test_build_dsp.py new file mode 100644 index 000000000..322a785cd --- /dev/null +++ b/tests/dsp/test_build_dsp.py @@ -0,0 +1,77 @@ +import os +from pathlib import Path + +import pytest + +from pygama import lgdo +from pygama.dsp import build_dsp +from pygama.lgdo.lh5_store import LH5Store, ls +from pygama.raw import build_raw + +config_dir = Path(__file__).parent/'configs' + + +@pytest.fixture(scope="module") +def multich_raw_file(lgnd_test_data): + out_file = '/tmp/L200-comm-20211130-phy-spms.lh5' + out_spec = { + 'FCEventDecoder': { + 'ch{key}': { + 'key_list': [[0, 6]], + 'out_stream': out_file + ":{name}", + 'out_name': 'raw' + } + } + } + + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=out_spec, overwrite=True) + + return out_file + + +def test_build_dsp_basics(lgnd_test_data): + build_dsp(lgnd_test_data.get_path('lh5/LDQTA_r117_20200110T105115Z_cal_geds_raw.lh5'), + '/tmp/LDQTA_r117_20200110T105115Z_cal_geds_dsp.lh5', + dsp_config=f'{config_dir}/icpc-dsp-config.json', + database={'pz': {'tau': 27460.5}}, + write_mode='r') + + assert os.path.exists('/tmp/LDQTA_r117_20200110T105115Z_cal_geds_dsp.lh5') + + with pytest.raises(FileExistsError): + build_dsp(lgnd_test_data.get_path('lh5/LDQTA_r117_20200110T105115Z_cal_geds_raw.lh5'), + '/tmp/LDQTA_r117_20200110T105115Z_cal_geds_dsp.lh5', + dsp_config=f'{config_dir}/icpc-dsp-config.json') + + with pytest.raises(FileNotFoundError): + build_dsp('non-existent-file.lh5', + '/tmp/LDQTA_r117_20200110T105115Z_cal_geds_dsp.lh5', + dsp_config=f'{config_dir}/icpc-dsp-config.json', + write_mode='r') + + +def test_build_dsp_channelwise(multich_raw_file): + chan_config = { + 'ch0/raw': f'{config_dir}/sipm-dsp-config.json', + 'ch1/raw': f'{config_dir}/sipm-dsp-config.json', + 'ch2/raw': f'{config_dir}/sipm-dsp-config.json' + } + + out_file = '/tmp/L200-comm-20211130-phy-spms_dsp.lh5' + build_dsp(multich_raw_file, + out_file, + {}, + n_max=5, + lh5_tables=chan_config.keys(), + chan_config=chan_config, + write_mode='r') + + assert ls(out_file) == ['ch0', 'ch1', 'ch2', 'dsp_info'] + assert ls(out_file, 'ch0/') == ['ch0/dsp'] + assert ls(out_file, 'ch0/dsp/') == ['ch0/dsp/bl_mean', 'ch0/dsp/bl_std'] + + store = LH5Store() + lh5_obj, n_rows = store.read_object('/ch0/dsp/bl_mean', out_file) + assert isinstance(lh5_obj, lgdo.Array) + assert len(lh5_obj) == 5 diff --git a/tests/dsp/test_processing_chain.py b/tests/dsp/test_processing_chain.py new file mode 100644 index 000000000..09677f328 --- /dev/null +++ b/tests/dsp/test_processing_chain.py @@ -0,0 +1,121 @@ +import pytest + +from pygama import lgdo +from pygama.dsp.processing_chain import build_processing_chain + + +def test_waveform_slicing(geds_raw_tbl): + dsp_config = { + "outputs": ["wf_blsub"], + "processors": { + "wf_blsub": { + "function": "bl_subtract", + "module": "pygama.dsp.processors", + "args": ["waveform[0:100]", "baseline", "wf_blsub"], + "unit": "ADC" + }, + } + } + proc_chain, _, tbl_out = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + assert list(tbl_out.keys()) == ["wf_blsub"] + assert isinstance(tbl_out["wf_blsub"], lgdo.WaveformTable) + assert tbl_out["wf_blsub"].wf_len == 100 + + +def test_processor_none_arg(geds_raw_tbl): + dsp_config = { + "outputs": ["wf_cum"], + "processors": { + "wf_cum": { + "function": "cumsum", + "module": "numpy", + "args": ["waveform", 1, None, "wf_cum"], + "kwargs": {"signature": "(n),(),()->(n)", "types": ["fii->f"]}, + "unit": "ADC" + } + } + } + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + dsp_config["processors"]["wf_cum"]["args"][2] = "None" + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + +def test_processor_kwarg_assignment(geds_raw_tbl): + dsp_config = { + "outputs": ["wf_cum"], + "processors": { + "wf_cum": { + "function": "cumsum", + "module": "numpy", + "args": ["waveform", "axis=1", "out=wf_cum"], + "kwargs": {"signature": "(n),()->(n)", "types": ["fi->f"]}, + "unit": "ADC" + } + } + } + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + dsp_config["processors"]["wf_cum"]["args"][1] = "dtypo=None" + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + with pytest.raises(TypeError): + proc_chain.execute(0, 1) + + +def test_processor_dtype_arg(geds_raw_tbl): + dsp_config = { + "outputs": ["wf_cum"], + "processors": { + "wf_cum": { + "function": "cumsum", + "module": "numpy", + "args": ["waveform", "axis=0", "dtype='int32'", "out=wf_cum"], + "kwargs": {"signature": "(n),(),()->(n)", "types": ["fiU->i"]}, + "unit": "ADC" + } + } + } + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + +def test_scipy_gauss_filter(geds_raw_tbl): + dsp_config = { + "outputs": ["wf_gaus"], + "processors": { + "wf_gaus": { + "function": "gaussian_filter1d", + "module": "scipy.ndimage", + "args": [ + "waveform", "0.1*us", "mode='reflect'", "truncate=3", + "output=wf_gaus" + ], + "kwargs": {"signature": "(n),(),(),()->(n)", "types": ["ffUf->f"]}, + "unit": "ADC" + } + } + } + proc_chain, _, _ = build_processing_chain(geds_raw_tbl, dsp_config) + proc_chain.execute(0, 1) + + +def test_processor_variable_array_output(spms_raw_tbl): + dsp_config = { + "outputs": ["wf_cum"], + "processors": { + "vt_max_out, vt_min_out, n_max_out, n_min_out, flag_out": { + "function": "get_multi_local_extrema", + "module": "pygama.dsp.processors", + "args": ["curr", "3*bl_std", "vt_max_out(10)", "vt_min_out(10)", "n_max_out", "n_min_out", "flag_out"], + "unit": "ADC" + } + } + } + + proc_chain, _, _ = build_processing_chain(spms_raw_tbl, dsp_config) + proc_chain.execute(0, 1) diff --git a/tests/raw/configs/fc-out-spec.json b/tests/raw/configs/fc-out-spec.json index 260c8930f..ebfbb4343 100644 --- a/tests/raw/configs/fc-out-spec.json +++ b/tests/raw/configs/fc-out-spec.json @@ -1,8 +1,8 @@ { - "FCEventDecoder": { - "spms": { - "key_list": [[2, 4]], - "out_stream": "/tmp/L200-comm-20211130-phy-spms.lh5" - } + "FCEventDecoder": { + "spms": { + "key_list": [[2, 4]], + "out_stream": "/tmp/L200-comm-20211130-phy-spms.lh5" } + } } diff --git a/tests/raw/configs/orca-out-spec-cli.json b/tests/raw/configs/orca-out-spec-cli.json index 141f32698..3e536c71b 100644 --- a/tests/raw/configs/orca-out-spec-cli.json +++ b/tests/raw/configs/orca-out-spec-cli.json @@ -1,8 +1,8 @@ { - "ORFlashCamADCWaveformDecoder": { - "g{key}": { - "key_list": [[0, 3]], - "out_stream": "/tmp/{orig_basename}.lh5" - } + "ORFlashCamADCWaveformDecoder": { + "g{key}": { + "key_list": [[0, 3]], + "out_stream": "/tmp/{orig_basename}.lh5" } + } } diff --git a/tests/raw/configs/orca-out-spec.json b/tests/raw/configs/orca-out-spec.json index da2be4985..208a506c7 100644 --- a/tests/raw/configs/orca-out-spec.json +++ b/tests/raw/configs/orca-out-spec.json @@ -1,8 +1,8 @@ { - "ORFlashCamADCWaveformDecoder": { - "geds": { - "key_list": [[0, 3]], - "out_stream": "/tmp/L200-comm-20220519-phy-geds.lh5" - } + "ORFlashCamADCWaveformDecoder": { + "geds": { + "key_list": [[0, 3]], + "out_stream": "/tmp/L200-comm-20220519-phy-geds.lh5" } + } } diff --git a/tests/raw/test_build_raw.py b/tests/raw/test_build_raw.py index 592889b92..6ce83402a 100644 --- a/tests/raw/test_build_raw.py +++ b/tests/raw/test_build_raw.py @@ -1,32 +1,33 @@ -import os.path +import os from pathlib import Path import pytest -import pygama.raw as raw from pygama.lgdo.lh5_store import LH5Store, ls +from pygama.raw import build_raw config_dir = Path(__file__).parent/'configs' def test_build_raw_basics(lgnd_test_data): with pytest.raises(FileNotFoundError): - raw.build_raw(in_stream='non-existent-file') + build_raw(in_stream='non-existent-file') with pytest.raises(FileNotFoundError): - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - out_spec='non-existent-file.json') + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec='non-existent-file.json') def test_build_raw_fc(lgnd_test_data): - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio')) + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + overwrite=True) assert lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.lh5') != '' out_file = '/tmp/L200-comm-20211130-phy-spms.lh5' - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - out_spec=out_file) + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=out_file, overwrite=True) assert os.path.exists('/tmp/L200-comm-20211130-phy-spms.lh5') @@ -42,44 +43,48 @@ def test_build_raw_fc_out_spec(lgnd_test_data): } } - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - out_spec=out_spec, n_max=10) + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=out_spec, n_max=10, overwrite=True) store = LH5Store() lh5_obj, n_rows = store.read_object('/spms', out_file) assert n_rows == 10 assert (lh5_obj['channel'].nda == [2, 3, 4, 2, 3, 4, 2, 3, 4, 2]).all() - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - out_spec=f'{config_dir}/fc-out-spec.json', n_max=10) + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=f'{config_dir}/fc-out-spec.json', n_max=10, overwrite=True) def test_build_raw_fc_channelwise_out_spec(lgnd_test_data): out_file = '/tmp/L200-comm-20211130-phy-spms.lh5' out_spec = { 'FCEventDecoder': { - 's{key}': { + 'ch{key}': { 'key_list': [[0, 6]], - 'out_stream': out_file + 'out_stream': out_file + ":{name}", + 'out_name': 'raw' } } } - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - out_spec=out_spec) + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), + out_spec=out_spec, overwrite=True) - assert ls(out_file) == ['s0', 's1', 's2', 's3', 's4', 's5'] + assert ls(out_file) == ['ch0', 'ch1', 'ch2', 'ch3', 'ch4', 'ch5'] + assert ls(out_file, 'ch0/') == ['ch0/raw'] + assert ls(out_file, 'ch0/raw/waveform') == ['ch0/raw/waveform'] def test_build_raw_orca(lgnd_test_data): - raw.build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca')) + build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), + overwrite=True) assert lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.lh5') != '' out_file = '/tmp/L200-comm-20220519-phy-geds.lh5' - raw.build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), - out_spec=out_file) + build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), + out_spec=out_file, overwrite=True) assert os.path.exists('/tmp/L200-comm-20220519-phy-geds.lh5') @@ -95,19 +100,18 @@ def test_build_raw_orca_out_spec(lgnd_test_data): } } - raw.build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), - out_spec=out_spec, n_max=10) + build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), + out_spec=out_spec, n_max=10, overwrite=True) store = LH5Store() lh5_obj, n_rows = store.read_object('/geds', out_file) assert n_rows == 10 assert (lh5_obj['channel'].nda == [2, 3, 4, 2, 3, 4, 2, 3, 4, 2]).all() - raw.build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), - out_spec=f'{config_dir}/orca-out-spec.json', n_max=10) + build_raw(in_stream=lgnd_test_data.get_path('orca/fc/L200-comm-20220519-phy-geds.orca'), + out_spec=f'{config_dir}/orca-out-spec.json', n_max=10, overwrite=True) def test_build_raw_overwrite(lgnd_test_data): - with pytest.raises(RuntimeError): - raw.build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio'), - overwrite=False) + with pytest.raises(FileExistsError): + build_raw(in_stream=lgnd_test_data.get_path('fcio/L200-comm-20211130-phy-spms.fcio')) diff --git a/tests/vis/configs/hpge-dsp-config.json b/tests/vis/configs/hpge-dsp-config.json index 74912a5a7..797044bd4 100644 --- a/tests/vis/configs/hpge-dsp-config.json +++ b/tests/vis/configs/hpge-dsp-config.json @@ -1,29 +1,29 @@ { - "processors": { - "wf_blsub": { - "function": "bl_subtract", - "module": "pygama.dsp.processors", - "args": ["waveform", "baseline", "wf_blsub"], - "unit": "ADC" - }, - "wf_pz": { - "function": "pole_zero", - "module": "pygama.dsp.processors", - "args": ["wf_blsub", "150*us", "wf_pz"], - "unit": "ADC" - }, - "wf_trap": { - "function": "trap_norm", - "module": "pygama.dsp.processors", - "args": ["wf_pz", "10*us", "3.008*us", "wf_trap"], - "unit": "ADC" - }, - "trapEmax": { - "function": "amax", - "module": "numpy", - "args": ["wf_trap", 1, "trapEmax"], - "kwargs": {"signature": "(n),()->()", "types":["fi->f"]}, - "unit": "ADC" - } + "processors": { + "wf_blsub": { + "function": "bl_subtract", + "module": "pygama.dsp.processors", + "args": ["waveform", "baseline", "wf_blsub"], + "unit": "ADC" + }, + "wf_pz": { + "function": "pole_zero", + "module": "pygama.dsp.processors", + "args": ["wf_blsub", "150*us", "wf_pz"], + "unit": "ADC" + }, + "wf_trap": { + "function": "trap_norm", + "module": "pygama.dsp.processors", + "args": ["wf_pz", "10*us", "3.008*us", "wf_trap"], + "unit": "ADC" + }, + "trapEmax": { + "function": "amax", + "module": "numpy", + "args": ["wf_trap", 1, "trapEmax"], + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, + "unit": "ADC" } + } } diff --git a/tutorials/metadata/LPGTA_dsp.json b/tutorials/metadata/LPGTA_dsp.json index 9db719341..e243f3531 100644 --- a/tutorials/metadata/LPGTA_dsp.json +++ b/tutorials/metadata/LPGTA_dsp.json @@ -1,15 +1,28 @@ { - "outputs": [ "trapE", "trapEftp", "ct_corr", "bl", "bl_sig", "A_10", "AoE", - "tp_max", "tp_0", "tp_10", "tp_50", "tp_80", "tp_90", - "dcr_raw"], - "processors":{ - "bl , bl_sig, slope, intercept":{ + "outputs": [ + "trapE", + "trapEftp", + "ct_corr", + "bl", + "bl_sig", + "A_10", + "AoE", + "tp_max", + "tp_0", + "tp_10", + "tp_50", + "tp_80", + "tp_90", + "dcr_raw" + ], + "processors": { + "bl , bl_sig, slope, intercept": { "function": "linear_slope_fit", "module": "pygama.dsp.processors", - "args" : ["waveform[:1650]", "bl","bl_sig", "slope","intercept"], - "unit": ["ADC","ADC","ADC","ADC"] + "args": ["waveform[:1650]", "bl", "bl_sig", "slope", "intercept"], + "unit": ["ADC", "ADC", "ADC", "ADC"] }, - "wf_blsub":{ + "wf_blsub": { "function": "subtract", "module": "numpy", "args": ["waveform", "bl", "wf_blsub"], @@ -22,7 +35,7 @@ "args": ["wf_blsub", "db.pz_const", "wf_pz"], "prereqs": ["wf_blsub"], "unit": "ADC", - "defaults": { "db.pz_const":"400*us" } + "defaults": { "db.pz_const": "400*us" } }, "wf_trap": { "function": "trap_norm", @@ -35,15 +48,25 @@ "function": "asym_trap_filter", "module": "pygama.dsp.processors", "prereqs": ["wf_pz"], - "args": ["wf_pz", "db.atrap.rise", "db.atrap.flat", "db.atrap.fall", "wf_atrap"], - "defaults" : {"db.atrap.rise":"round(20*ns)", "db.atrap.flat":"round(1*us)","db.atrap.fall":"round(4*us)"}, + "args": [ + "wf_pz", + "db.atrap.rise", + "db.atrap.flat", + "db.atrap.fall", + "wf_atrap" + ], + "defaults": { + "db.atrap.rise": "round(20*ns)", + "db.atrap.flat": "round(1*us)", + "db.atrap.fall": "round(4*us)" + }, "unit": "ADC" }, "trapE": { "function": "amax", "module": "numpy", "args": ["wf_trap", 1, "trapE"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC", "prereqs": ["wf_trap"] }, @@ -51,7 +74,7 @@ "function": "fixed_time_pickoff", "module": "pygama.dsp.processors", "args": ["wf_trap", "tp_0+13.5*us", "trapEftp"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC", "prereqs": ["wf_trap", "tp_0"] }, @@ -67,14 +90,14 @@ "function": "amax", "module": "numpy", "args": ["curr10", 1, "A_10"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC", "prereqs": ["curr10"] }, "AoE": { "function": "divide", "module": "numpy", - "args":["A_10", "trapE", "AoE"], + "args": ["A_10", "trapE", "AoE"], "unit": "1/sample", "prereqs": ["A_10", "trapE"] }, @@ -82,7 +105,7 @@ "function": "argmax", "module": "numpy", "args": ["wf_blsub", 1, "tp_max"], - "kwargs": {"signature":"(n),()->()", "types":["fi->i"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->i"] }, "unit": "ns", "prereqs": ["wf_blsub"] }, @@ -130,7 +153,7 @@ "ct_corr": { "function": "trap_pickoff", "module": "pygama.dsp.processors", - "args":["wf_pz", "4*us", 0, "tp_0", "ct_corr"], + "args": ["wf_pz", "4*us", 0, "tp_0", "ct_corr"], "unit": "ADC", "prereqs": ["wf_pz", "tp_0"] } diff --git a/tutorials/metadata/cage.json b/tutorials/metadata/cage.json index 8dfdd2258..50e8304ab 100644 --- a/tutorials/metadata/cage.json +++ b/tutorials/metadata/cage.json @@ -1,25 +1,25 @@ { - "experiment" : "cage", + "experiment": "cage", - "runDB" : "/global/project/projectdirs/legend/software/CAGE/processing/runDB.json", - "fileDB" : "/global/project/projectdirs/legend/software/CAGE/processing/fileDB.h5", - "runSelectionDB" : "/global/project/projectdirs/legend/software/CAGE/processing/metadata/run_selection.json", + "runDB": "/global/project/projectdirs/legend/software/CAGE/processing/runDB.json", + "fileDB": "/global/project/projectdirs/legend/software/CAGE/processing/fileDB.h5", + "runSelectionDB": "/global/project/projectdirs/legend/software/CAGE/processing/metadata/run_selection.json", - "daq" : "ORCA", - "daq_dir" : "/global/cfs/cdirs/m2676/data/cage/Data", - "daq_ignore" : [".log",".DS_Store", "openFiles"], - "lh5_dir" : "/global/cfs/cdirs/m2676/data/cage/LH5", + "daq": "ORCA", + "daq_dir": "/global/cfs/cdirs/m2676/data/cage/Data", + "daq_ignore": [".log", ".DS_Store", "openFiles"], + "lh5_dir": "/global/cfs/cdirs/m2676/data/cage/LH5", - "nersc_dir":"/global/project/projectdirs/m2676/data/cage", + "nersc_dir": "/global/project/projectdirs/m2676/data/cage", - "tier_dirs" : ["raw","dsp","hit"], - "subsystems" : [""], - "run_types" : ["",""], - "evt_dirs" : [""], + "tier_dirs": ["raw", "dsp", "hit"], + "subsystems": [""], + "run_types": ["", ""], + "evt_dirs": [""], - "unique_key" : "cage-cyc{cycle}-{YYYY}{mm}{dd}", - "daq_template" : "{YYYY}-{mm}-{dd}-CAGERun{cycle}", - "lh5_template" : "cage_run{run}_cyc{cycle}_{tier}.lh5", + "unique_key": "cage-cyc{cycle}-{YYYY}{mm}{dd}", + "daq_template": "{YYYY}-{mm}-{dd}-CAGERun{cycle}", + "lh5_template": "cage_run{run}_cyc{cycle}_{tier}.lh5", - "ecal_default" : "/global/project/projectdirs/legend/software/CAGE/processing/metadata/config_ecal.json" + "ecal_default": "/global/project/projectdirs/legend/software/CAGE/processing/metadata/config_ecal.json" } diff --git a/tutorials/metadata/dsp_config.json b/tutorials/metadata/dsp_config.json index 509dab4a2..b2201c668 100644 --- a/tutorials/metadata/dsp_config.json +++ b/tutorials/metadata/dsp_config.json @@ -1,15 +1,32 @@ { "outputs": [ - "timestamp", "channel", "trapEmax", "trapEftp", "triE", "bl_mean", "bl_sig", "A_10", "AoE", "dcr", "zacE", "cuspE" + "timestamp", + "channel", + "trapEmax", + "trapEftp", + "triE", + "bl_mean", + "bl_sig", + "A_10", + "AoE", + "dcr", + "zacE", + "cuspE" ], - "processors":{ - "bl_mean , bl_sig, bl_slope, bl_intercept":{ + "processors": { + "bl_mean , bl_sig, bl_slope, bl_intercept": { "function": "linear_slope_fit", "module": "pygama.dsp.processors", - "args" : ["waveform[0: 1000]", "bl_mean","bl_sig", "bl_slope","bl_intercept"], - "unit": ["ADC","ADC","ADC","ADC"] + "args": [ + "waveform[0: 1000]", + "bl_mean", + "bl_sig", + "bl_slope", + "bl_intercept" + ], + "unit": ["ADC", "ADC", "ADC", "ADC"] }, - "wf_blsub":{ + "wf_blsub": { "function": "subtract", "module": "numpy", "args": ["waveform", "bl_mean", "wf_blsub"], @@ -20,7 +37,7 @@ "module": "pygama.dsp.processors", "args": ["wf_blsub", "db.pz_const", "wf_pz"], "unit": "ADC", - "defaults": { "db.pz_const":"72*us" } + "defaults": { "db.pz_const": "72*us" } }, "wf_trap": { "function": "trap_norm", @@ -40,78 +57,78 @@ "args": ["wf_pz", "0.05*us", "2*us", "3*us", "wf_atrap"], "unit": "ADC" }, - "tp_min, tp_max, wf_min, wf_max":{ + "tp_min, tp_max, wf_min, wf_max": { "function": "min_max", "module": "pygama.dsp.processors", "args": ["waveform", "tp_min", "tp_max", "wf_min", "wf_max"], - "unit": ["ns","ns","ADC", "ADC"] + "unit": ["ns", "ns", "ADC", "ADC"] }, "tp_0": { "function": "time_point_thresh", "module": "pygama.dsp.processors", "args": ["wf_atrap", 0, "tp_max", 0, "tp_0"], "unit": "ns" - }, + }, "tp_95": { "function": "time_point_thresh", "module": "pygama.dsp.processors", - "args": ["wf_blsub", "0.95*tp_max", "tp_max", 0,"tp_95"], + "args": ["wf_blsub", "0.95*tp_max", "tp_max", 0, "tp_95"], "unit": "ns" - }, + }, "tp_80": { "function": "time_point_thresh", "module": "pygama.dsp.processors", - "args": ["wf_blsub", "0.8*tp_max", "tp_max", 0,"tp_80"], + "args": ["wf_blsub", "0.8*tp_max", "tp_max", 0, "tp_80"], "unit": "ns" - }, + }, "tp_50": { "function": "time_point_thresh", "module": "pygama.dsp.processors", - "args": ["wf_blsub", "0.50*tp_max", "tp_max", 0,"tp_50"], + "args": ["wf_blsub", "0.50*tp_max", "tp_max", 0, "tp_50"], "unit": "ns" - }, + }, "tp_20": { "function": "time_point_thresh", "module": "pygama.dsp.processors", - "args": ["wf_blsub", "0.2*tp_max", "tp_max", 0,"tp_20"], + "args": ["wf_blsub", "0.2*tp_max", "tp_max", 0, "tp_20"], "unit": "ns" - }, + }, "tp_05": { "function": "time_point_thresh", "module": "pygama.dsp.processors", - "args": ["wf_blsub", "0.05*tp_max", "tp_max", 0,"tp_05"], + "args": ["wf_blsub", "0.05*tp_max", "tp_max", 0, "tp_05"], "unit": "ns" - }, - "trapEftp": { - "function": "fixed_time_pickoff", - "module": "pygama.dsp.processors", - "args": ["wf_trap", "tp_0+(5*us+9*us)", "trapEftp"], - "unit": "ADC" - }, - "dcr_raw": { - "function": "trap_pickoff", - "module": "pygama.dsp.processors", - "args": ["wf_pz", 200, 1000, "tp_0+1.504*us", "dcr_raw"], - "unit": "ADC" - }, - "dcr": { - "function": "divide", - "module": "numpy", - "args": ["dcr_raw", "trapEftp", "dcr"], - "unit": "ADC" - }, + }, + "trapEftp": { + "function": "fixed_time_pickoff", + "module": "pygama.dsp.processors", + "args": ["wf_trap", "tp_0+(5*us+9*us)", "trapEftp"], + "unit": "ADC" + }, + "dcr_raw": { + "function": "trap_pickoff", + "module": "pygama.dsp.processors", + "args": ["wf_pz", 200, 1000, "tp_0+1.504*us", "dcr_raw"], + "unit": "ADC" + }, + "dcr": { + "function": "divide", + "module": "numpy", + "args": ["dcr_raw", "trapEftp", "dcr"], + "unit": "ADC" + }, "trapEmax": { "function": "amax", "module": "numpy", "args": ["wf_trap", 1, "trapEmax"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC" }, "triE": { "function": "amax", "module": "numpy", "args": ["wf_triangle", 1, "triE"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC" }, "curr10": { @@ -125,41 +142,51 @@ "function": "amax", "module": "numpy", "args": ["curr10", 1, "A_10"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC" }, "AoE": { "function": "divide", "module": "numpy", - "args":["A_10", "trapEmax", "AoE"], + "args": ["A_10", "trapEmax", "AoE"], "unit": "1/sample" }, "zac_wf": { "function": "zac_filter", "module": "pygama.dsp.processors", "args": ["waveform", "zac_wf(101, 'f')"], - "init_args": ["len(waveform)", "round(3.5*us/waveform.period)", "round(2.5*us/waveform.period)", "round(72*us/waveform.period)"], + "init_args": [ + "len(waveform)", + "round(3.5*us/waveform.period)", + "round(2.5*us/waveform.period)", + "round(72*us/waveform.period)" + ], "unit": "ADC" }, "zacE": { "function": "amax", "module": "numpy", "args": ["zac_wf", 1, "zacE"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC" }, "cusp_wf": { "function": "cusp_filter", "module": "pygama.dsp.processors", "args": ["waveform", "cusp_wf(101, 'f')"], - "init_args": ["len(waveform)", "round(10*us/waveform.period)", "round(1*us/waveform.period)", "round(200*us/waveform.period)"], + "init_args": [ + "len(waveform)", + "round(10*us/waveform.period)", + "round(1*us/waveform.period)", + "round(200*us/waveform.period)" + ], "unit": "ADC" }, "cuspE": { "function": "amax", "module": "numpy", "args": ["cusp_wf", 1, "cuspE"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "ADC" }, "curr": { @@ -179,7 +206,7 @@ "function": "amax", "module": "numpy", "args": ["curr_psd[10*mhz*len(curr):]", 1, "hf_max"], - "kwargs": {"signature":"(n),()->()", "types":["fi->f"]}, + "kwargs": { "signature": "(n),()->()", "types": ["fi->f"] }, "unit": "" }, "power_spect": { @@ -189,6 +216,5 @@ "init_args": ["wf_blsub", "power_spect(len(wf_blsub)//2+1, 'f')"], "unit": "MHz" } - } } diff --git a/tutorials/metadata/optimizer_results.json b/tutorials/metadata/optimizer_results.json index fdcbc39d1..dca6f86a2 100644 --- a/tutorials/metadata/optimizer_results.json +++ b/tutorials/metadata/optimizer_results.json @@ -1 +1,9 @@ -{"pz2": {"tau1": "52.00*us", "tau2": "6.40*us", "frac": "0.043"}, "etrap": {"rise": "11.00*us", "flat": "3.00*us"}, "atrap": {"rise": "20*ns", "flat": "1*us", "fall": "4*us"}, "tp_ftp": {"ftp": "tp_0 + 12.5*us"}, "dcr_pz": {"tau1": "52.00*us", "tau2": "6.40*us", "frac": "0.043"}, "dcr_trap": {"flat": "20.33*us", "rise": "7.00*us"}, "dcr": {"ftp": "80*us"}} +{ + "pz2": { "tau1": "52.00*us", "tau2": "6.40*us", "frac": "0.043" }, + "etrap": { "rise": "11.00*us", "flat": "3.00*us" }, + "atrap": { "rise": "20*ns", "flat": "1*us", "fall": "4*us" }, + "tp_ftp": { "ftp": "tp_0 + 12.5*us" }, + "dcr_pz": { "tau1": "52.00*us", "tau2": "6.40*us", "frac": "0.043" }, + "dcr_trap": { "flat": "20.33*us", "rise": "7.00*us" }, + "dcr": { "ftp": "80*us" } +}