Repository: legend-exp/pygama (language: python)

Issue:
Unit testing
Here's a to-do list, I'll kep it up-to-date. Good reference: https://scikit-hep.org/developer. Most of this will be implemented on the `refactor` branch only.

- [ ] Add micro-tests for each module/function separately. Read the [contribution guide (`refactor` branch)](https://legend-exp.github.io/pygama/refactor/contributing.html). Add @gipert as a reviewer for your pull requests. Assignments (put your name!):
  - `dsp` @iguinn 
  - `flow`
  - `hit` @wisecg 
  - `lgdo` @jasondet 
  - `math`
  - `pargen` @ggmarshall 
  - `raw` @jasondet 
  - `vis` @iguinn 
- [ ] Develop/import synthetic pulse injector
- [ ] Add high-level tests of the full production chain by using real data from https://github.com/legend-exp/legend-testdata @gipert, ...
  - Include all necessary data in legend-testdata
- [x] Let a Python linter run as part of tests to enforce uniform coding conventions (see also #242)
- [x] Add test coverage reporting through https://codecov.io and GitHub bot
Project packaging
Here's a to-do list, I'll keep it up-to-date. Good reference: https://scikit-hep.org/developer

- [x] Switch to more modern packaging practices: `pyproject.toml` (and `setup.cfg` + `setup.py`, if needed)
- [x] Do we want to publish the package on PyPI? This could be done automatically with GitHub actions
- [ ] Do we need to be more careful about licensing?
- [x] Add Dependabot (https://github.com/dependabot) to keep GitHub actions up to date
- [x] Enforce uniform coding conventions with a Python linter? Do it with https://pre-commit.ci bot
- [ ] Create separate repository for code that uses pygama, to keep code base clean? Apps, sandbox...
- [ ] Find right place for third-party code. e.g. Eigen (header-only C++)
- [x] How to handle optional dependencies (e.g. pyfcutils)?


Discussion / hints:

Fix patch:
diff --git a/.gitattributes b/.gitattributes
deleted file mode 100644
index 81d924c45..000000000
--- a/.gitattributes
+++ /dev/null
@@ -1,1 +0,0 @@
-*.ipynb    filter=jupyter_clear_output
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 000000000..f08ed35a5
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,9 @@
+version: 2
+updates:
+  # Maintain dependencies for GitHub Actions
+  - package-ecosystem: "github-actions"
+    directory: "/"
+    schedule:
+      interval: "daily"
+    ignore:
+      - dependency-name: "actions/*"
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index d320da5e5..b414b6a9d 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -32,7 +32,7 @@ jobs:
         python -m pip install .
     - name: Run unit tests
       run: |
-        pytest tests
+        pytest
 
   # test-in-legend-container:
   #   name: Test pygama in LEGEND container
@@ -48,7 +48,28 @@ jobs:
   #   - name: Install and run unit tests
   #     run: |
   #       python -m pip install .
-  #       pytest tests
+  #       pytest
+
+  test-coverage:
+    name: Calculate and upload test coverage
+    needs: build-and-test
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+        with:
+          fetch-depth: 2
+      - uses: actions/setup-python@v2
+        with:
+          python-version: '3.10'
+
+      - name: Generate Report
+        run: |
+          .github/workflows/get-dependencies.sh
+          python -m pip install .
+          python -m pip install coverage
+          coverage run -m pytest
+      - name: Upload Coverage to codecov.io
+        uses: codecov/codecov-action@v2
 
   deploy-docs:
     name: Deploy documentation on legend-exp.github.io/pygama
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 000000000..bad199496
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,82 @@
+ci:
+  autoupdate_commit_msg: "chore: update pre-commit hooks"
+  autofix_commit_msg: "style: pre-commit fixes"
+
+exclude: ^attic
+repos:
+- repo: https://github.com/pre-commit/pre-commit-hooks
+  rev: "v4.2.0"
+  hooks:
+  - id: check-added-large-files
+  - id: check-case-conflict
+  - id: check-executables-have-shebangs
+  - id: check-merge-conflict
+  - id: check-symlinks
+  - id: check-yaml
+  - id: check-json
+  - id: check-toml
+  - id: check-docstring-first
+  - id: debug-statements
+  - id: end-of-file-fixer
+  - id: forbid-new-submodules
+  - id: mixed-line-ending
+  - id: requirements-txt-fixer
+  - id: trailing-whitespace
+
+- repo: https://github.com/asottile/setup-cfg-fmt
+  rev: "v1.20.1"
+  hooks:
+  - id: setup-cfg-fmt
+
+- repo: https://github.com/PyCQA/isort
+  rev: "5.10.1"
+  hooks:
+  - id: isort
+
+- repo: https://github.com/asottile/pyupgrade
+  rev: "v2.32.0"
+  hooks:
+  - id: pyupgrade
+    args: ["--py36-plus"]
+
+# - repo: https://github.com/psf/black
+#   rev: "22.3.0"
+#   hooks:
+#   - id: black
+
+# - repo: https://github.com/PyCQA/flake8
+#   rev: "4.0.1"
+#   hooks:
+#   - id: flake8
+#     additional_dependencies: [flake8-bugbear, flake8-print]
+
+- repo: https://github.com/kynan/nbstripout
+  rev: "0.5.0"
+  hooks:
+    - id: nbstripout
+      args: ["--strip-empty-cells",
+             "--extra-keys", "metadata.kernelspec metadata.language_info"]
+
+- repo: https://github.com/mgedmin/check-manifest
+  rev: "0.48"
+  hooks:
+  - id: check-manifest
+    stages: [manual]
+
+- repo: https://github.com/codespell-project/codespell
+  rev: "v2.1.0"
+  hooks:
+  - id: codespell
+    args: ["-L", "hist,gaus,nd,ans,crate,nin,puls,spms"]
+
+- repo: https://github.com/shellcheck-py/shellcheck-py
+  rev: "v0.8.0.4"
+  hooks:
+  - id: shellcheck
+
+- repo: https://github.com/pre-commit/pygrep-hooks
+  rev: "v1.9.0"
+  hooks:
+  - id: rst-backticks
+  - id: rst-directive-colons
+  - id: rst-inline-touching-normal
diff --git a/LICENSE b/LICENSE
index e4f139cb4..cd482d897 100644
--- a/LICENSE
+++ b/LICENSE
@@ -198,4 +198,4 @@ Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
-limitations under the License.
\ No newline at end of file
+limitations under the License.
diff --git a/README.md b/README.md
index 8e8b63647..204cf72f8 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,9 @@
 # pygama
 
 ![GitHub tag (latest by date)](https://img.shields.io/github/v/tag/legend-exp/pygama?logo=git)
-![Documentation (main)](https://img.shields.io/badge/documentation-online-purple?logo=readthedocs&link=https%3A%2F%2Flegend-exp.github.io%2Fpygama)
-![GitHub Workflow Status (main)](https://img.shields.io/github/workflow/status/legend-exp/pygama/pygama/main?label=main%20branch&logo=github)
-![GitHub Workflow Status (dev)](https://img.shields.io/github/workflow/status/legend-exp/pygama/pygama/dev?label=dev%20branch&logo=github)
+![Documentation](https://img.shields.io/badge/documentation-online-purple?logo=readthedocs&link=https%3A%2F%2Flegend-exp.github.io%2Fpygama)
+![GitHub Workflow Status](https://img.shields.io/github/workflow/status/legend-exp/pygama/pygama/main?label=main%20branch&logo=github)
+![Codecov](https://img.shields.io/codecov/c/github/legend-exp/pygama?logo=codecov)
 ![GitHub issues](https://img.shields.io/github/issues/legend-exp/pygama?logo=github)
 ![GitHub pull requests](https://img.shields.io/github/issues-pr/legend-exp/pygama?logo=github)
 ![License](https://img.shields.io/github/license/legend-exp/pygama)
@@ -15,3 +15,4 @@
 * optimizing DSP routines and tuning associated analysis parameters
 * generating and selecting high-level event data for further analysis
 
+Check out the [online documentation](https://legend-exp.github.io/pygama).
diff --git a/codecov.yml b/codecov.yml
new file mode 100644
index 000000000..3010946f7
--- /dev/null
+++ b/codecov.yml
@@ -0,0 +1,16 @@
+codecov:
+  require_ci_to_pass: true
+
+coverage:
+  status:
+    project:
+      default:
+        target: 20%
+    patch: false
+
+github_checks:
+    annotations: false
+
+ignore:
+  - tests
+  - docs
diff --git a/docs/source/contributing.rst b/docs/source/contributing.rst
index 1f4e15de9..005a2bf71 100644
--- a/docs/source/contributing.rst
+++ b/docs/source/contributing.rst
@@ -8,7 +8,23 @@ that do not comply to the following directives will be rejected.
 Code style
 ----------
 
-*In progress...*
+A set of `pre-commit <https://pre-commit.com>`_ hooks is configured to make
+sure that *pygama* coherently follows standard coding style conventions. The
+pre-commit tool is able to identify common style problems and automatically fix
+them, wherever possible. Configured hooks are listed in the
+``.pre-commit-config.yaml`` file at the project root folder. They are run
+remotely on the GitHub repository through the `pre-commit bot
+<https://pre-commit.ci>`_, but can also be run locally before submitting a
+pull request (recommended):
+
+.. code-block:: console
+
+   $ pip install pre-commit      # required
+   $ pre-commit run --all-files  # analyse the source code and fix it wherever possible
+   $ pre-commit install          # install a Git pre-commit hook (optional)
+
+For a more comprehensive guide, check out the `Scikit-HEP documentation about
+code style <https://scikit-hep.org/developer/style>`_.
 
 Testing
 -------
@@ -22,11 +38,11 @@ Testing
   overview.
 * *pygama* tests belong to three categories:
 
-  :micro-tests: Should ensure the correct behaviour of each function
+  :unit tests: Should ensure the correct behaviour of each function
       independently, possibly without relying on other *pygama* methods. The
-      existence of these microscopic tests makes it possible to promptly
-      identify and fix the source of a bug. An example of this are tests for
-      each single DSP processor
+      existence of these micro-tests makes it possible to promptly identify and
+      fix the source of a bug. An example of this are tests for each single DSP
+      processor
 
   :integration tests: Should ensure that independent parts of the code base
       work well together and are integrated in a cohesive framework. An example
@@ -40,10 +56,19 @@ Testing
 * Unit tests are automatically run for every push event and pull request to the
   remote Git repository on a remote server (currently handled by GitHub
   actions). Every pull request must pass all tests before being approved for
-  merging
+  merging.
 * Additionally, pull request authors are required to provide tests with
   sufficient code coverage for every proposed change or addition. If necessary,
-  high-level functional tests should be updated.
+  high-level functional tests should be updated. We currently rely on
+  `codecov.io <https://app.codecov.io/gh/legend-exp/pygama>`_ to keep track of
+  the test coverage. To generate a local coverage report (recommended before
+  submitting pull requests), the ``coverage`` Python package is needed:
+
+  .. code-block:: console
+
+    $ pip install coverage
+    $ coverage run -m pytest
+    $ coverage report
 
 Documentation
 -------------
diff --git a/pygama/dsp/_processors/Wiener_filter.py b/pygama/dsp/_processors/Wiener_filter.py
index b69544762..8d8c6f54d 100644
--- a/pygama/dsp/_processors/Wiener_filter.py
+++ b/pygama/dsp/_processors/Wiener_filter.py
@@ -1,13 +1,14 @@
 import numpy as np
-from pygama import lh5
 from numba import guvectorize
+
+from pygama import lh5
 from pygama.dsp.errors import DSPFatal
 
 
 def Wiener_filter(file_name_array):
     """
     Apply a Wiener filter to the waveform.  Note that this convolution is performed in the frequency domain
-    
+
     Initialization Parameters
     -------------------------
     file_name_array : string
@@ -19,8 +20,8 @@ def Wiener_filter(file_name_array):
     fft_w_in : array-like
            The fourier transform input waveform
     fft_w_out: array-like
-           The filtered waveform, in the frequency domain     
-           
+           The filtered waveform, in the frequency domain
+
     Processing Chain Example
     ------------------------
     "wf_wiener": {
@@ -32,69 +33,69 @@ def Wiener_filter(file_name_array):
         "init_args": ["/path/to/file/wiener.lh5"]
     }
     """
-    
+
     sto = lh5.Store()
 
     # Check that the file is valid and the data is in the correct format
 
     try: file_name_array[0]
     except: raise DSPFatal('init_args must be an array with the filename')
-    
+
     file_name = file_name_array[0]
-     
+
     try: f = sto.gimme_file(file_name, 'r')
-    except: raise DSPFatal('File must be a valid lh5 file') 
+    except: raise DSPFatal('File must be a valid lh5 file')
 
     if 'spms/processed/superpulse' not in f:
         raise DSPFatal('lh5 file must have \'spms/processed/superpulse\' as a group')
-        
+
     if 'spms/processed/noise_wf' not in f:
         raise DSPFatal('lh5 file must have \'spms/processed/noise_wf\' as a group')
-        
-    # Read in the data 
-    
+
+    # Read in the data
+
     superpulse, _ = sto.read_object('spms/processed/superpulse', file_name)
     superpulse = superpulse.nda
-    
+
     noise_wf, _ = sto.read_object('spms/processed/noise_wf', file_name)
     noise_wf = noise_wf.nda
- 
+
     # Now check that the data are valid
-    
+
     if len(superpulse) <= 0:
         raise DSPFatal('The length of the filter must be positive')
-        
-    if len(superpulse) != len(noise_wf): 
+
+    if len(superpulse) != len(noise_wf):
         raise DSPFatal('The length of the superpulse must be equal to the length of the noise waveform')
 
     if np.argmax(superpulse) <= 0 or np.argmax(superpulse) > len(superpulse):
         raise DSPFatal('The index of the maximum of the superpulse must occur within the waveform')
-    
-    # Transform these to the frequency domain to eventually create the Wiener filter 
-    
+
+    # Transform these to the frequency domain to eventually create the Wiener filter
+
     fft_superpulse = np.fft.fft(superpulse)
     fft_noise_wf = np.fft.fft(noise_wf)
-    
-    # Create the point spread function for the detector's response 
-    
-    def PSF(superpulse, fft_superpulse):  
-        
+
+    # Create the point spread function for the detector's response
+
+    def PSF(superpulse, fft_superpulse):
+
         delta = np.zeros_like(superpulse)
         arg_max = np.argmax(superpulse)
         delta[arg_max] = np.amax(superpulse)
 
         return fft_superpulse/np.fft.fft(delta)
-    
-    # Now create the Wiener filter in the frequency domain 
-    
+
+    # Now create the Wiener filter in the frequency domain
+
     fft_PSF = PSF(superpulse, fft_superpulse)
     PSD_noise_wf = fft_noise_wf*np.conj(fft_noise_wf)
     PSD_superpulse = fft_superpulse*np.conj(fft_superpulse)
 
-    w_filter = (np.conj(fft_PSF))/((fft_PSF*np.conj(fft_PSF))+(PSD_noise_wf/PSD_superpulse)) 
-    
-    # Create a factory function that performs the convolution with the Wiener filter, the output is still in the frequency domain 
-    
+    w_filter = (np.conj(fft_PSF))/((fft_PSF*np.conj(fft_PSF))+(PSD_noise_wf/PSD_superpulse))
+
+    # Create a factory function that performs the convolution with the Wiener filter, the output is still in the frequency domain
+
     @guvectorize(["void(complex64[:], complex64[:])",
                   "void(complex128[:], complex128[:])"],
                  "(n)->(n)", forceobj=True)
@@ -109,4 +110,4 @@ def Wiener_out(fft_w_in, fft_w_out):
 
         fft_w_out[:] = fft_w_in * w_filter
 
-    return Wiener_out
\ No newline at end of file
+    return Wiener_out
diff --git a/pygama/dsp/_processors/bl_subtract.py b/pygama/dsp/_processors/bl_subtract.py
index 7c56ba79b..5302f7485 100644
--- a/pygama/dsp/_processors/bl_subtract.py
+++ b/pygama/dsp/_processors/bl_subtract.py
@@ -1,7 +1,9 @@
 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)
diff --git a/pygama/dsp/_processors/convolutions.py b/pygama/dsp/_processors/convolutions.py
index 4d6f96823..a1af0ca38 100644
--- a/pygama/dsp/_processors/convolutions.py
+++ b/pygama/dsp/_processors/convolutions.py
@@ -1,7 +1,9 @@
 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
@@ -61,7 +63,7 @@ def cusp_filter(length, sigma, flat, decay):
 
     den   = [1, -np.exp(-1 / decay)]
     cuspd = np.convolve(cusp, den, 'same')
-    
+
     @guvectorize(["void(float32[:], float32[:])",
                   "void(float64[:], float64[:])"],
                  "(n),(m)", forceobj=True)
diff --git a/pygama/dsp/_processors/fftw.py b/pygama/dsp/_processors/fftw.py
index 13959da06..c851bb335 100644
--- a/pygama/dsp/_processors/fftw.py
+++ b/pygama/dsp/_processors/fftw.py
@@ -2,6 +2,7 @@
 from numba import guvectorize
 from pyfftw import FFTW
 
+
 def dft(buf_in, buf_out):
     """
     Perform discrete Fourier transforms using the FFTW library.  FFTW optimizes
@@ -30,10 +31,10 @@ def dft(buf_in, buf_out):
     - 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)
@@ -62,16 +63,16 @@ def inv_dft(buf_in, buf_out):
         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) 
+        - 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)
@@ -104,19 +105,19 @@ def psd(buf_in, buf_out):
         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) 
+        - 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/pygama/dsp/_processors/fixed_time_pickoff.py b/pygama/dsp/_processors/fixed_time_pickoff.py
index dec3e4e3e..8dfc2d1ba 100644
--- a/pygama/dsp/_processors/fixed_time_pickoff.py
+++ b/pygama/dsp/_processors/fixed_time_pickoff.py
@@ -1,7 +1,9 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 @guvectorize(["void(float32[:], float32, float32[:])",
               "void(float64[:], float64, float64[:])"],
              "(n),()->()", nopython=True, cache=True)
@@ -39,5 +41,5 @@ def fixed_time_pickoff(w_in, t_in, a_out):
 
     if int(t_in) < 0 or int(t_in) >= len(w_in):
         return
-  
+
     a_out[0] = w_in[int(t_in)]
diff --git a/pygama/dsp/_processors/gaussian_filter1d.py b/pygama/dsp/_processors/gaussian_filter1d.py
index c5ec43936..5b02695b3 100644
--- a/pygama/dsp/_processors/gaussian_filter1d.py
+++ b/pygama/dsp/_processors/gaussian_filter1d.py
@@ -34,9 +34,10 @@
 # originally called a function from a C library.  In this code, the convolution is
 # performed with NumPy's built in convolution function.
 
-import numpy 
+import numpy
 from numba import guvectorize
 
+
 def gaussian_filter1d(sigma, truncate):
     """1-D Gaussian filter.
     Parameters
@@ -47,7 +48,7 @@ def gaussian_filter1d(sigma, truncate):
     truncate : float, optional
         Truncate the filter at this many standard deviations.
         Default is 4.0.
-    Returns 
+    Returns
     -------
     gaussian_filter1d : ndarray
     """
@@ -62,18 +63,18 @@ def _gaussian_kernel1d(sigma, radius):
         return phi_x
 
     sd = float(sigma)
-    
+
     # Make the radius of the filter equal to truncate standard deviations
-    
+
     lw = int(truncate * sd + 0.5)
-    
+
     # Since we are calling correlate, not convolve, revert the kernel
-    
+
     weights = _gaussian_kernel1d(sigma, lw)[::-1]
     weights = numpy.asarray(weights, dtype=numpy.float64)
 
     # Find the length of the kernel so we can reflect the signal an appropriate amount
-    
+
     extension_length = int(len(weights) / 2) + 1
 
     """Calculate a 1-D correlation along the given axis.
@@ -97,30 +98,30 @@ def gaussian_filter1d_out(wf_in, wf_out):
         # Extend the signal on the left and right by at least half of the length of the kernel
 
         wf_in = numpy.asarray(wf_in)
-        
+
         # Short warning message if kernel is larger than signal, in which case signal can't be convolved
-        
+
         if len(wf_in) < extension_length:
             print('Kernel calculated was larger than signal, try again with smaller parameters')
             return 0
-        
+
         # This mode extends as a reflection
         # ‘reflect’ (d c b a | a b c d | d c b a)
         # The input is extended by reflecting about the edge of the last pixel. This mode is also
         # sometimes referred to as half-sample symmetric.
-       
+
         reflected_front = numpy.flip(wf_in[0:extension_length])
         reflected_end   = numpy.flip(wf_in[-extension_length:])
-        
+
         # Extend the signal
-        
+
         extended_signal = wf_in
         extended_signal = numpy.concatenate((extended_signal, reflected_end  ), axis=None)
         extended_signal = numpy.concatenate((reflected_front, extended_signal), axis=None)
         output = numpy.correlate(extended_signal, weights, mode='same')
-        
+
         # Now extract the original signal length
-        
+
         wf_out[:] = output[extension_length:-extension_length]
-        
+
     return gaussian_filter1d_out
diff --git a/pygama/dsp/_processors/get_multi_local_extrema.py b/pygama/dsp/_processors/get_multi_local_extrema.py
index c6ac057d6..aedec2798 100644
--- a/pygama/dsp/_processors/get_multi_local_extrema.py
+++ b/pygama/dsp/_processors/get_multi_local_extrema.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numba import guvectorize
-from pygama.dsp.errors import DSPFatal
 
+from pygama.dsp.errors import DSPFatal
 
 
 @guvectorize(["void(float32[:], float32, float32[:], float32[:], float32[:], float32[:], float32[:])",
@@ -25,40 +25,40 @@ def get_multi_local_extrema(w_in, a_delta_in, vt_max_out, vt_min_out, n_max_out,
     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 
+    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 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 
-    
+    # prepare output
+
     vt_max_out[:]= np.nan
-    vt_min_out[:] = np.nan    
+    vt_min_out[:] = np.nan
     n_max_out[0] = np.nan
     n_min_out[0] = np.nan
     flag_out[0] = np.nan
-    
+
     # initialize internal counters
-    
-    n_max_counter = 0 
-    n_min_counter = 0 
-    
-    
-    # Checks 
-    
+
+    n_max_counter = 0
+    n_min_counter = 0
+
+
+    # Checks
+
     if (np.isnan(w_in).any() or np.isnan(a_delta_in)):
         return
 
     if (not len(vt_max_out)<len(w_in) or not len(vt_min_out)<len(w_in)):
         raise DSPFatal('The length of your return array must be smaller than the length of your waveform')
-    if (not a_delta_in >= 0): 
+    if (not a_delta_in >= 0):
         raise DSPFatal('a_delta_in must be positive')
-        
+
     # now loop over data
-    
+
     imax, imin = 0, 0
     find_max = True
     for i in range(len(w_in)):
@@ -82,11 +82,11 @@ def get_multi_local_extrema(w_in, a_delta_in, vt_max_out, vt_min_out, n_max_out,
                 n_min_counter += 1
                 imax = i
                 find_max = True
-                
+
     # set output
     n_max_out[0] = n_max_counter
     n_min_out[0] = n_min_counter
-    
+
     if n_max_out[0] == 1:
         flag_out[0] = 1
     else:
diff --git a/pygama/dsp/_processors/linear_slope_fit.py b/pygama/dsp/_processors/linear_slope_fit.py
index e39ead2c5..60aabd053 100644
--- a/pygama/dsp/_processors/linear_slope_fit.py
+++ b/pygama/dsp/_processors/linear_slope_fit.py
@@ -1,16 +1,18 @@
 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, mean, stdev, slope, intercept):
     """
     Calculate the mean and standard deviation of the waveform using
     Welford's method as well as the slope an intercept of the waveform
     using linear regression.
-    
+
     Parameters
     ----------
     w_in     : array-like
diff --git a/pygama/dsp/_processors/log_check.py b/pygama/dsp/_processors/log_check.py
index 24453e000..20932a268 100644
--- a/pygama/dsp/_processors/log_check.py
+++ b/pygama/dsp/_processors/log_check.py
@@ -1,7 +1,9 @@
 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)
@@ -28,11 +30,11 @@ def log_check(w_in, w_log):
     }
     """
     w_log[:] = np.nan
-    
+
     if np.isnan(w_in).any():
         return
 
     if np.any(w_in <= 0):
         return
-        
+
     w_log[:] = np.log(w_in[:])
diff --git a/pygama/dsp/_processors/min_max.py b/pygama/dsp/_processors/min_max.py
index b83e56c4f..8275da9b4 100644
--- a/pygama/dsp/_processors/min_max.py
+++ b/pygama/dsp/_processors/min_max.py
@@ -1,7 +1,9 @@
 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)
diff --git a/pygama/dsp/_processors/moving_windows.py b/pygama/dsp/_processors/moving_windows.py
index 658060d3c..02237c61f 100644
--- a/pygama/dsp/_processors/moving_windows.py
+++ b/pygama/dsp/_processors/moving_windows.py
@@ -1,7 +1,9 @@
 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)
diff --git a/pygama/dsp/_processors/multi_a_filter.py b/pygama/dsp/_processors/multi_a_filter.py
index ee487a0c7..adb619959 100644
--- a/pygama/dsp/_processors/multi_a_filter.py
+++ b/pygama/dsp/_processors/multi_a_filter.py
@@ -1,19 +1,21 @@
 import numpy as np
 from numba import guvectorize
-from pygama.dsp.errors import DSPFatal
+
 from pygama.dsp._processors.fixed_time_pickoff import fixed_time_pickoff
+from pygama.dsp.errors import DSPFatal
+
 
 @guvectorize(["void(float32[:], float32[:], float32[:])",
               "void(float64[:], float64[:], float64[:])"],
              "(n),(m),(m)", forceobj=True, cache=True)
-def multi_a_filter(w_in, vt_maxs_in, va_max_out): 
+def multi_a_filter(w_in, vt_maxs_in, va_max_out):
     """
     Finds the maximums in a waveform and returns the amplitude of the wave at those points
-    
+
     Parameters
     ----------
     w_in : array-like
-        The array of data within which amplitudes of extrema will be found 
+        The array of data within which amplitudes of extrema will be found
     vt_maxs_in : array-like
         The array of max positions for each wf
     Returns
@@ -21,17 +23,17 @@ def multi_a_filter(w_in, vt_maxs_in, va_max_out):
     va_max_out: array-like
         An array of the amplitudes of the maximums of the waveform
     """
-    
-    # Initialize output parameters 
-    
+
+    # Initialize output parameters
+
     va_max_out[:] = np.nan
-    
-    # Check inputs 
-    
+
+    # Check inputs
+
     if (np.isnan(w_in).any()):
         return
-        
+
     if (not len(vt_maxs_in)<len(w_in)):
         raise DSPFatal('The length of your return array must be smaller than the length of your waveform')
-    
-    fixed_time_pickoff(w_in, vt_maxs_in, va_max_out)
\ No newline at end of file
+
+    fixed_time_pickoff(w_in, vt_maxs_in, va_max_out)
diff --git a/pygama/dsp/_processors/multi_t_filter.py b/pygama/dsp/_processors/multi_t_filter.py
index 03d852f26..4851cfe33 100644
--- a/pygama/dsp/_processors/multi_t_filter.py
+++ b/pygama/dsp/_processors/multi_t_filter.py
@@ -1,52 +1,53 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp._processors.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  
+    """
+    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 through the array again 
-    to ensure there are no duplicate indices. If there are duplicate indicies caused by a
+    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 maximum of a waveform
     isn't right at index 0.
     ----------
     t_in : array-like
-        The array of indices that we want to remove duplicates from 
+        The array of indices that we want to remove duplicates from
     vt_min_in : array-like
-        List of indicies of minima that we want to replace duplicates in t_out with
+        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
     """
     # initialize arrays
     t_out[:] = np.nan
-    
+
     # checks
     if (np.isnan(t_in).all() and np.isnan(vt_min_in).all()): # we pad these with NaNs, so only return if there is nothing to analyze
         return
-    
+
     # check if any later indexed values are equal to the earliest instance
     k=0
     for index1 in range(len(t_in)):
         for index2 in range(len(t_in[index1+1:])):
-            if t_in[index1] == t_in[index2+index1+1]: 
+            if t_in[index1] == t_in[index2+index1+1]:
                 t_out[index2+index1+1] = vt_min_in[k]
         k+=1 # this makes sure that the index of the misidentified afterpulse tp0 is replaced with the correct corresponding minimum
-        
-    # Fill up the output with the rest of the values from the input that weren't repeats    
-    for index in range(len(t_in)): 
+
+    # Fill up the output with the rest of the values from the input that weren't repeats
+    for index in range(len(t_in)):
         if (np.isnan(t_out[index]) and not np.isnan(t_in[index])):
             t_out[index] = t_in[index]
-            
+
     # makes sure that the first maximum found isn't the start of the waveform
     if not np.isnan(t_out[0]):
-        if int(t_out[0]) == 0: 
+        if int(t_out[0]) == 0:
             t_out[:] = np.append(t_out[1:],np.nan)
 
 
@@ -58,15 +59,15 @@ 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 
+    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. 
+    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.
     Parameters
     ----------
     w_in : array-like
         The array of data within which the list of tp0s will be found
-    a_threshold_in: scalar 
+    a_threshold_in: scalar
         Threshold to search for using time_point_thresh
     vt_maxs_in : array-like
         The array of max positions for each wf
@@ -77,25 +78,25 @@ def multi_t_filter(w_in, a_threshold_in, vt_max_in, vt_min_in, t_out):
     t_out : array-like
         Array of fixed length (padded with nans) that hold the indices of
         the identified initial rise times of peaks in the signal
-        
+
     """
-     
+
     # initialize arrays, padded with the elements we want
-    t_out[:] = np.nan 
-    
-    # checks 
+    t_out[:] = np.nan
+
+    # checks
     if (np.isnan(w_in).any() or np.isnan(a_threshold_in)):
         return
     if (np.isnan(vt_max_in).all() and np.isnan(vt_min_in).all()):
-        return 
+        return
     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
     intermediate_t_out = np.full_like(t_out, np.nan, dtype=np.float32)
-    
+
     # Go through the list of maxima, calling time_point_thresh (the refactored version ignores the nan padding)
     time_point_thresh(w_in, a_threshold_in, vt_max_in, 0, intermediate_t_out)
 
     # Remove duplicates from the t_out list
-    remove_duplicates(intermediate_t_out, vt_min_in, t_out)
\ No newline at end of file
+    remove_duplicates(intermediate_t_out, vt_min_in, t_out)
diff --git a/pygama/dsp/_processors/optimize.py b/pygama/dsp/_processors/optimize.py
index 62cf67e0e..b8ae79ee3 100644
--- a/pygama/dsp/_processors/optimize.py
+++ b/pygama/dsp/_processors/optimize.py
@@ -1,8 +1,11 @@
 import numpy as np
+from iminuit import Minuit
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
-from .pole_zero import pole_zero, double_pole_zero
-from iminuit import Minuit
+
+from .pole_zero import double_pole_zero, pole_zero
+
 
 class Model:
     """
@@ -30,7 +33,7 @@ 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
     by minimizing the slope in the waveform's specified time range.
-    
+
     Parameters
     ----------
     w_in         : array-like
@@ -47,7 +50,7 @@ def optimize_1pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, val0_out):
                    The initial guess of the optimal time constant
     val0_out     : float
                    The output value of the best-fit time constant
-              
+
     Processing Chain Example
     ------------------------
     "tau0": {
@@ -78,7 +81,7 @@ def optimize_1pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, val0_out):
     m.errordef = Minuit.LEAST_SQUARES
     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)
@@ -86,7 +89,7 @@ def optimize_2pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, p1_in, p2_in, v
     """
     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
@@ -111,7 +114,7 @@ def optimize_2pz(w_in, a_baseline_in, t_beg_in, t_end_in, p0_in, p1_in, p2_in, v
                    The output value of the best-fit, shorter time constant
     val2_out     : float
                    The output value of the best-fit fraction
-              
+
     Processing Chain Example
     ------------------------
     "tau1, tau2, frac": {
diff --git a/pygama/dsp/_processors/param_lookup.py b/pygama/dsp/_processors/param_lookup.py
index a5e6bc397..47e767edc 100644
--- a/pygama/dsp/_processors/param_lookup.py
+++ b/pygama/dsp/_processors/param_lookup.py
@@ -1,5 +1,6 @@
 import numpy as np
-from numba import guvectorize, types, from_dtype
+from numba import from_dtype, guvectorize, types
+
 
 def param_lookup(param_dict, default_val, dtype):
     """
@@ -12,7 +13,7 @@ def param_lookup(param_dict, default_val, 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):
diff --git a/pygama/dsp/_processors/pole_zero.py b/pygama/dsp/_processors/pole_zero.py
index 4d6e5b617..41fc75cfb 100644
--- a/pygama/dsp/_processors/pole_zero.py
+++ b/pygama/dsp/_processors/pole_zero.py
@@ -1,7 +1,9 @@
 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)
@@ -29,7 +31,7 @@ def pole_zero(w_in, t_tau, w_out):
         "prereqs": ["wf_bl"]
     }
     """
-    w_out[:] = np.nan 
+    w_out[:] = np.nan
 
     if np.isnan(w_in).any() or np.isnan(t_tau):
         return
@@ -70,7 +72,7 @@ def double_pole_zero(w_in, t_tau1, t_tau2, frac, w_out):
         "prereqs": ["wf_bl"]
     }
     """
-    w_out[:] = np.nan 
+    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
diff --git a/pygama/dsp/_processors/presum.py b/pygama/dsp/_processors/presum.py
index f67b2faa9..d776463e3 100644
--- a/pygama/dsp/_processors/presum.py
+++ b/pygama/dsp/_processors/presum.py
@@ -1,7 +1,9 @@
 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)
diff --git a/pygama/dsp/_processors/pulse_injector.py b/pygama/dsp/_processors/pulse_injector.py
index b72078817..9853469e5 100644
--- a/pygama/dsp/_processors/pulse_injector.py
+++ b/pygama/dsp/_processors/pulse_injector.py
@@ -1,6 +1,8 @@
+from math import exp, log
+
 import numpy as np
 from numba import guvectorize
-from math import log,exp
+
 
 @guvectorize(["void(float32[:], float32, float32, float32, float32, float32[:])",
               "void(float64[:], float64, float64,float64,float64,float64[:])"],
@@ -40,6 +42,6 @@ def inject_exp_pulse(wf_in,t0,rt,A,decay, wf_out):
     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)/(rise))*exp(-(1/decay)*(T-t0))))
+            wf_out[T] += (A*exp((T-t0-rt)/(rise))*exp(-(1/decay)*(T-t0)))
         elif (T>t0+rt):
-            wf_out[T] += (A*exp(-(1/decay)*(T-t0)))
\ No newline at end of file
+            wf_out[T] += (A*exp(-(1/decay)*(T-t0)))
diff --git a/pygama/dsp/_processors/saturation.py b/pygama/dsp/_processors/saturation.py
index 5251c0dd5..da134e3e4 100644
--- a/pygama/dsp/_processors/saturation.py
+++ b/pygama/dsp/_processors/saturation.py
@@ -1,7 +1,9 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 @guvectorize(["void(float32[:], float32, float32[:], float32[:])",
               "void(float64[:], float64, float64[:], float64[:])"],
              "(n),()->(),()", nopython=True, cache=True)
@@ -10,7 +12,7 @@ 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.
-    
+
     Parameters
     ----------
     w_in        : array-like
@@ -21,7 +23,7 @@ def saturation(w_in, bit_depth_in, n_lo_out, n_hi_out):
                   The output number of samples at the minimum
     n_hi_out    : int
                   The output number of samples at the maximum
-    
+
     Processing Chain Example
     ------------------------
     "sat_lo, sat_hi": {
diff --git a/pygama/dsp/_processors/soft_pileup_corr.py b/pygama/dsp/_processors/soft_pileup_corr.py
index 099e74d15..1a801d533 100644
--- a/pygama/dsp/_processors/soft_pileup_corr.py
+++ b/pygama/dsp/_processors/soft_pileup_corr.py
@@ -1,7 +1,9 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 @guvectorize(["void(float32[:], float32, float32, float32[:])",
               "void(float64[:], float64, float64, float64[:])"],
              "(n),(),()->(n)", nopython=True, cache=True)
@@ -10,7 +12,7 @@ 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.
-    
+
     Parameters
     ----------
     w_in  : array-like
@@ -22,7 +24,7 @@ def soft_pileup_corr(w_in, n_in, tau_in, w_out):
             The fixed, exponential time constant
     w_out : array-like
             The output waveform with the exponential subtracted
-    
+
     Processing Chain Example
     ------------------------
     "wf_bl": {
@@ -71,7 +73,7 @@ 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
     constant and then subtract the best-fit function from the
     entire waveform.
-    
+
     Parameters
     ----------
     w_in  : array-like
@@ -85,7 +87,7 @@ def soft_pileup_corr_bl(w_in, n_in, tau_in, B_in, w_out):
             The fixed, exponential constant
     w_out : array-like
             The output waveform with the exponential subtracted
-    
+
     Processing Chain Example
     ------------------------
     "wf_bl": {
diff --git a/pygama/dsp/_processors/template.py b/pygama/dsp/_processors/template.py
index fe7030500..334daaa1c 100644
--- a/pygama/dsp/_processors/template.py
+++ b/pygama/dsp/_processors/template.py
@@ -2,10 +2,12 @@
 #
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 # 2) Provide instructions to Numba
-# 
+#
 # Documentation about Numba guvectorize decorator:
 # https://numba.pydata.org/numba-doc/latest/user/vectorize.html#the-guvectorize-decorator
 #
@@ -65,7 +67,7 @@ def the_processor_template(w_in, t_in, a_in, w_out, t_out):
         raise DSPFatal('The error message goes here')
 
     # 7) Algorithm
-    # 
+    #
     # Loop over waveforms by using a "for i in range(.., .., ..)" instruction.
     # Avoid loops based on a while condition which might lead to segfault.  Avoid
     # also enumerate/ndenumerate to keep code as similar as possible among all
diff --git a/pygama/dsp/_processors/time_point_thresh.py b/pygama/dsp/_processors/time_point_thresh.py
index 9a24406c7..5c44b99a8 100644
--- a/pygama/dsp/_processors/time_point_thresh.py
+++ b/pygama/dsp/_processors/time_point_thresh.py
@@ -1,7 +1,9 @@
 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)
@@ -37,7 +39,7 @@ def time_point_thresh(w_in, a_threshold, t_start, walk_forward, t_out):
 
     if np.isnan(w_in).any() or np.isnan(a_threshold) or np.isnan(t_start) or np.isnan(walk_forward):
         return
-    
+
     if np.floor(t_start) != t_start:
         raise DSPFatal('The starting index must be an integer')
 
diff --git a/pygama/dsp/_processors/trap_filters.py b/pygama/dsp/_processors/trap_filters.py
index aa6567ab6..8ab2a62a2 100644
--- a/pygama/dsp/_processors/trap_filters.py
+++ b/pygama/dsp/_processors/trap_filters.py
@@ -1,7 +1,9 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 @guvectorize(["void(float32[:], int32, int32, float32[:])",
               "void(float64[:], int32, int32, float64[:])"],
              "(n),(),()->(n)", nopython=True, cache=True)
@@ -37,13 +39,13 @@ def trap_filter(w_in, rise, flat, w_out):
 
     if int(rise) < 0:
         raise DSPFatal('The number of samples in the rise section must be positive')
-    
+
     if int(flat) < 0:
         raise DSPFatal('The number of samples in the flat section must be positive')
-    
+
     if 2 * int(rise) + int(flat) > len(w_in):
         raise DSPFatal('The trapezoid width is wider than the waveform')
-    
+
     w_out[0] = w_in[0]
     for i in range(1, rise, 1):
         w_out[i] = w_out[i-1] + w_in[i]
@@ -96,7 +98,7 @@ def trap_norm(w_in, rise, flat, w_out):
 
     if 2 * int(rise) + int(flat) > len(w_in):
         raise DSPFatal('The trapezoid width is wider than the waveform')
-    
+
     w_out[0] = w_in[0] / rise
     for i in range(1, rise, 1):
         w_out[i] = w_out[i-1] + w_in[i] / rise
@@ -173,7 +175,7 @@ def trap_pickoff(w_in, rise, flat, t_pickoff, a_out):
     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
@@ -186,7 +188,7 @@ def trap_pickoff(w_in, rise, flat, t_pickoff, a_out):
                The waveform index to pick off
     a_out    : float
                The output pick-off value of the filtered waveform
-    
+
     Processing Chain Example
     ------------------------
     "ct_corr": {
@@ -220,7 +222,7 @@ def trap_pickoff(w_in, rise, flat, t_pickoff, a_out):
 
     if not len(w_in) >= start_time >= 2*rise + flat:
         return
-    
+
     for i in range(start_time - rise, start_time, 1):
         I_1 += w_in[i]
     for i in range(start_time - 2 * rise - flat, start_time - rise - flat, 1):
diff --git a/pygama/dsp/_processors/windower.py b/pygama/dsp/_processors/windower.py
index 6d7a96627..2f34d02cd 100644
--- a/pygama/dsp/_processors/windower.py
+++ b/pygama/dsp/_processors/windower.py
@@ -1,7 +1,9 @@
 import numpy as np
 from numba import guvectorize
+
 from pygama.dsp.errors import DSPFatal
 
+
 @guvectorize(["void(float32[:], float32, float32[:])",
               "void(float64[:], float64, float64[:])"],
              "(n),(),(m)", nopython=True, cache=True)
@@ -13,7 +15,7 @@ def windower(w_in, t0_in, w_out):
     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.
-    
+
     Parameters
     ----------
     w_in : array-like
@@ -30,7 +32,7 @@ def windower(w_in, t0_in, w_out):
 
     if len(w_out) >= len(w_in):
         raise DSPFatal('The windowed waveform must be smaller than the input waveform')
-    
+
     beg = min(int(t0_in), len(w_in))
     end = max(beg + len(w_out), 0)
     if beg < 0:
diff --git a/pygama/dsp/build_dsp.py b/pygama/dsp/build_dsp.py
index 9eea30f5e..4c8307798 100755
--- a/pygama/dsp/build_dsp.py
+++ b/pygama/dsp/build_dsp.py
@@ -1,26 +1,28 @@
 #!/usr/bin/env python3
+import argparse
+import json
 import os
 import sys
-import json
-import h5py
 import time
+
+import h5py
 import numpy as np
-import argparse
 
 import pygama
+import pygama.lgdo.lh5_store as lh5
 from pygama import git
-from pygama.dsp.processing_chain import build_processing_chain
 from pygama.dsp.errors import DSPFatal
-import pygama.lgdo.lh5_store as lh5
+from pygama.dsp.processing_chain import build_processing_chain
 from pygama.math.utils import tqdm_range
 
+
 def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
                outputs=None, n_max=np.inf, write_mode='r', buffer_len=3200,
                block_width=16, verbose=1, chan_config=None):
     """
     Convert raw-tier LH5 data into dsp-tier LH5 data by running a sequence of
     processors via the ProcessingChain.
-    
+
     Parameters
     ----------
     f_raw : str
@@ -54,7 +56,7 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
     t_start = time.time()
 
     if isinstance(dsp_config, str):
-        with open(dsp_config, 'r') as config_file:
+        with open(dsp_config) as config_file:
             dsp_config = json.load(config_file)
 
     if not isinstance(dsp_config, dict):
@@ -90,13 +92,13 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
 
     # load DSP config (default: one config file for all tables)
     if isinstance(dsp_config, str):
-        with open(dsp_config, 'r') as config_file:
+        with open(dsp_config) as config_file:
             dsp_config = json.load(config_file, object_pairs_hook=OrderedDict)
 
     # get the database parameters. For now, this will just be a dict in a json
     # file, but eventually we will want to interface with the metadata repo
     if isinstance(database, str):
-        with open(database, 'r') as db_file:
+        with open(database) as db_file:
             database = json.load(db_file)
 
     if database and not isinstance(database, dict):
@@ -131,7 +133,7 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
         # if we have separate DSP files for each table, read them in here
         if chan_config is not None:
             f_config = chan_config[tb]
-            with open(f_config, 'r') as config_file:
+            with open(f_config) as config_file:
                 dsp_config = json.load(config_file, object_pairs_hook=OrderedDict)
             print('Processing table:', tb, 'with DSP config file:\n  ', f_config)
 
@@ -156,15 +158,15 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
             if proc_chain is None:
                 proc_chain, lh5_it.field_mask, tb_out = build_processing_chain(lh5_in, dsp_config, db_dict, outputs, verbose, block_width)
                 tqdm_it = iter(tqdm_range(0, tot_n_rows, buffer_len, verbose))
-            
+
             n_rows = min(tot_n_rows-start_row, n_rows)
             try:
                 proc_chain.execute(0, n_rows)
             except DSPFatal as e:
                 # Update the wf_range to reflect the file position
-                e.wf_range = "{}-{}".format(e.wf_range[0]+start_row, e.wf_range[1]+start_row)
+                e.wf_range = f"{e.wf_range[0]+start_row}-{e.wf_range[1]+start_row}"
                 raise e
-            
+
             raw_store.write_object(obj=tb_out,
                                    name=tb_name,
                                    lh5_file=f_dsp,
@@ -173,10 +175,10 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
                                    write_start=write_offset+start_row,
             )
             next(tqdm_it)
-            
+
             if start_row+n_rows > tot_n_rows:
                 break
-            
+
         print(f'Done.  Writing to file {f_dsp}')
 
     raw_store.write_object(dsp_info, 'dsp_info', f_dsp)
@@ -187,9 +189,9 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
 
 
 if __name__=="__main__":
-    parser = argparse.ArgumentParser(description=
-"""Process a single tier 1 LH5 file and produce a tier 2 LH5 file using a
-json config file and build_dsp.""")
+    parser = argparse.ArgumentParser(description="""Process a single tier 1 LH5
+    file and produce a tier 2 LH5 file using a json config file and
+    build_dsp.""")
 
     arg = parser.add_argument
     arg('file', help="Input raw LH5 file.")
@@ -206,17 +208,17 @@ def build_dsp(f_raw, f_dsp, dsp_config, lh5_tables=None, database=None,
 
     arg('-n', '--nevents', default=None, type=int,
         help="Number of waveforms to process. By default do the whole file")
-    
+
     arg('-v', '--verbose', action='store_const', const=2, default=1,
         help="Verbose output, useful for debugging")
     arg('-q', '--quiet', action='store_const', const=0, dest='verbose',
         help="Silent output, print only exceptions thrown")
-    
+
     arg('-b', '--block', default=16, type=int,
         help="Number of waveforms to process simultaneously. Default is 8")
     arg('-c', '--chunk', default=3200, type=int,
         help="Number of waveforms to read from disk at a time. Default is 256. THIS IS NOT IMPLEMENTED YET!")
-    
+
     arg('-r', '--recreate', action='store_const', const='r', dest='writemode', default='r',
         help="Overwrite file if it already exists. Default option. Multually exclusive with --update and --append")
     arg('-u', '--update', action='store_const', const='u', dest='writemode',
diff --git a/pygama/dsp/errors.py b/pygama/dsp/errors.py
index d16f435bf..7fd925557 100644
--- a/pygama/dsp/errors.py
+++ b/pygama/dsp/errors.py
@@ -16,7 +16,7 @@ class DSPFatal(DSPError):
         is caught, and appended to the error message
     """
     def __init__(self, *args):
-        super(DSPFatal, self).__init__(*args)
+        super().__init__(*args)
         self.wf_range = None
         self.processor = None
 
@@ -26,7 +26,7 @@ def __str__(self):
             suffix += '\nThrown while processing entries ' + str(self.wf_range)
         if self.processor:
             suffix += '\nThrown by ' + self.processor
-        return super(DSPFatal, self).__str__() + suffix
+        return super().__str__() + suffix
 
 class ProcessingChainError(DSPError):
     """Error thrown when there is a problem setting up a processing chain"""
diff --git a/pygama/dsp/processing_chain.py b/pygama/dsp/processing_chain.py
index 802f48e3f..0b5db96f9 100644
--- a/pygama/dsp/processing_chain.py
+++ b/pygama/dsp/processing_chain.py
@@ -1,28 +1,28 @@
 from __future__ import annotations
 
-import numpy as np
-import json
-import re
 import ast
-import itertools as it
 import importlib
+import itertools as it
+import json
+import re
 from abc import ABCMeta, abstractmethod
-
+from copy import copy, deepcopy
 from dataclasses import dataclass, field
 from typing import Union
-from copy import deepcopy, copy
 
+import numpy as np
 from numba import vectorize
 
-from .errors import ProcessingChainError, DSPFatal
 import pygama.lgdo.lh5_store as lh5
+from pygama.math.units import Quantity, Unit
 from pygama.math.units import unit_registry as ureg
-from pygama.math.units import Unit, Quantity
+
+from .errors import DSPFatal, ProcessingChainError
 
 # Filler value for variables to be automatically deduced later
 auto = 'auto'
 
-# Map from ast interpretter operations to functions to call and format string
+# Map from ast interpreter operations to functions to call and format string
 ast_ops_dict = {ast.Add: (np.add, '{}+{}'),
                 ast.Sub: (np.subtract, '{}-{}'),
                 ast.Mult: (np.multiply, '{}*{}'),
@@ -36,8 +36,8 @@ class CoordinateGrid:
     and offset. Period is a unitted Quantity, offset is a scalar in units of
     period, a Unit or a ProcChainVar. In the last case, a ProcessingChain
     variable is used to store a different offset for each event"""
-    period: Union[Quantity, Unit, str]
-    offset: Union[Quantity, ProcChainVar, float, int]
+    period: Quantity | Unit | str
+    offset: Quantity | ProcChainVar | float | int
 
     def __post_init__(self):
         if isinstance(self.period, str):
@@ -66,7 +66,7 @@ def get_offset(self, unit=None):
         # Get the offset (convert)ed to unit. If unit is None use period
         if unit is None: unit = self.period
         elif isinstance(unit, str): unit = ureg.Quantity(unit)
-        
+
         if isinstance(self.offset, ProcChainVar):
             return self.offset.get_buffer(CoordinateGrid(unit, 0))
         else:
@@ -75,7 +75,7 @@ def get_offset(self, unit=None):
     def __str__(self):
         offset = self.offset.name if isinstance(self.offset, ProcChainVar) \
             else str(self.offset)
-        return "({},{})".format(str(self.period), offset)
+        return f"({str(self.period)},{offset})"
 
 class ProcChainVar:
     """Helper data class with buffer and information for internal variables in
@@ -92,11 +92,11 @@ def __init__(self,
         assert isinstance(proc_chain, ProcessingChain) and isinstance(name, str)
         self.proc_chain = proc_chain
         self.name = name
-        
+
         # ndarray containing data buffer of size block_width x shape
         # list of ndarrays in different coordinate systems if is_coord is true
-        self._buffer: Union[list, np.ndarray] = None
-        
+        self._buffer: list | np.ndarray = None
+
         self.shape = shape
         self.dtype = dtype
         self.grid = grid
@@ -109,7 +109,7 @@ def __init__(self,
     def __setattr__(self, name, value):
         if value is auto:
             pass
-        
+
         elif name=='shape':
             if hasattr(value, '__iter__'):
                 value = tuple(value)
@@ -127,7 +127,7 @@ def __setattr__(self, name, value):
 
         elif name=='unit' and not value is None:
             value = str(value)
-            
+
         elif name=='is_coord':
             value = bool(value)
             if value:
@@ -136,8 +136,8 @@ def __setattr__(self, name, value):
                 elif isinstance(self._buffer, np.ndarray):
                     self._buffer = [(self._buffer, CoordinateGrid(self.unit, 0))]
 
-        super(ProcChainVar, self).__setattr__(name, value)
-        
+        super().__setattr__(name, value)
+
     def get_buffer(self, unit=None):
         # If buffer needs to be created, do so now
         if self._buffer is None:
@@ -166,11 +166,11 @@ def get_buffer(self, unit=None):
                 raise ProcessingChainError("Cannot deduce shape of "+self.name)
             if self.dtype is auto:
                 raise ProcessingChainError("Cannot deduce shape of "+self.name)
-            
+
             buff = np.zeros(shape=(self.proc_chain._block_width,)+self.shape, dtype=self.dtype)
             self._buffer.append((buff, unit))
             return buff
-        
+
         # check if coordinate conversion has been done already
         for buff, grid in self._buffer:
             if grid == unit:
@@ -199,17 +199,17 @@ def description(self):
             str(self.shape), str(self.dtype), str(self.grid), \
             str(self.unit), " (coord)" if self.is_coord is True else '' )
 
-    # Update any variables set to auto; leave the others alone. Emit a messsage
+    # Update any variables set to auto; leave the others alone. Emit a message
     # only if anything was updated
     def update_auto(self, shape=auto, dtype=auto, grid=auto, unit=auto, is_coord=auto, period=None, offset=0):
         updated = False
-        
+
         # Construct coordinate grid from period/offset if given
         if grid is auto and period is not None:
             if isinstance(offset, str):
                 offset = get_variable(offset, expr_only=True)
             grid = CoordinateGrid(period, offset)
-        
+
         if self.shape is auto and shape is not auto:
             self.shape = shape
             updated = True
@@ -227,10 +227,10 @@ def update_auto(self, shape=auto, dtype=auto, grid=auto, unit=auto, is_coord=aut
             updated = True
         if updated:
             self.proc_chain._print(2, 'Updated variable:', self.description())
-        
+
     def __str__(self):
         return self.name
-        
+
 ###############################################################################
 ############################### The Main Class ################################
 ###############################################################################
@@ -305,13 +305,13 @@ def add_variable(self, name, dtype=auto, shape=auto, grid=auto, unit=auto, is_co
             if isinstance(offset, str):
                 offset = get_variable(offset, expr_only=True)
             grid = CoordinateGrid(period, offset)
-        
+
         var = ProcChainVar(self, name, shape=shape, dtype=dtype, grid=grid,
                            unit=unit, is_coord=is_coord)
         self._vars_dict[name] = var
         return var
-    
-    
+
+
     def link_input_buffer(self, varname, buff=None):
         """Link an input buffer to a variable
         Parameters
@@ -321,14 +321,14 @@ def link_input_buffer(self, varname, buff=None):
             be created with a similar shape to the provided buffer
         buff : numpy array or lgdo class to use as input buffer. If None,
             create a new buffer with a similar shape to the variable
-        
+
         Return : buff or newly allocated input buffer
         """
         self._validate_name(varname, raise_exception=True)
         var = self.get_variable(varname, expr_only=True)
         if var is None:
             var = self.add_variable(varname)
-        
+
         if not isinstance(var, ProcChainVar):
             raise ProcessingChainError("Must link an input buffer to a processing chain variable")
 
@@ -359,8 +359,8 @@ def link_input_buffer(self, varname, buff=None):
         self._input_managers.append(out_man)
 
         return buff
-    
-    
+
+
     def link_output_buffer(self, varname, buff=None):
         """Link an output buffer to a variable
         Parameters
@@ -370,14 +370,14 @@ def link_output_buffer(self, varname, buff=None):
             be created with a similar shape to the provided buffer
         buff : numpy array or lgdo class to use as output buffer. If None,
             create a new buffer with a similar shape to the variable
-        
+
         Return : buff or newly allocated output buffer
         """
         self._validate_name(varname, raise_exception=True)
         var = self.get_variable(varname, expr_only=True)
         if var is None:
             var = self.add_variable(varname)
-        
+
         if not isinstance(var, ProcChainVar):
             raise ProcessingChainError("Must link an output buffer to a processing chain variable")
 
@@ -430,7 +430,7 @@ def execute(self, start=0, stop=None):
         if stop is None: stop = self._buffer_len
         for i in range(start, stop, self._block_width):
             self._execute_procs(i, min(i+self._block_width, self._buffer_len))
-        
+
 
 
     def get_variable(self, expr, get_names_only=False, expr_only=False):
@@ -497,7 +497,7 @@ def _parse_expr(self, node, expr, dry_run, var_name_list):
 
         elif isinstance(node, ast.Str):
             return node.s
-        
+
         elif isinstance(node, ast.Constant):
             return node.value
 
@@ -510,7 +510,7 @@ def _parse_expr(self, node, expr, dry_run, var_name_list):
             #check if it is a variable
             var_name_list.append(node.id)
             if dry_run: return None
-            
+
             val = self._vars_dict.get(node.id, None)
             if val is None:
                 val = self.add_variable(node.id)
@@ -534,7 +534,7 @@ def _parse_expr(self, node, expr, dry_run, var_name_list):
                 out = ProcChainVar(self, name, lhs.shape, lhs.dtype, lhs.grid, lhs.unit, is_coord=lhs.is_coord)
             else:
                 out = ProcChainVar(self, name, rhs.shape, rhs.dtype, rhs.grid, rhs.unit, is_coord=rhs.is_coord)
-                
+
             self._proc_managers.append(ProcessorManager(self, op, [lhs, rhs, out]))
             return out
 
@@ -544,7 +544,7 @@ def _parse_expr(self, node, expr, dry_run, var_name_list):
             if operand is None: return None
             op, op_form = ast_ops_dict[type(node.op)]
             name = '('+op_form.format(str(operand))+')'
-            
+
             if isinstance(operand, ProcChainVar):
                 out = ProcChainVar(self, name, operand.shape, operand.dtype,
                               operand.grid, operand.unit, operand.is_coord)
@@ -578,9 +578,9 @@ def get_index(slice_value):
             if isinstance(node.slice, ast.Index):
                 index = get_index(node.slice.value)
                 out_buf = val[..., index]
-                out_name='{}[{}]'.format(str(val), index),
+                out_name=f'{str(val)}[{index}]',
                 out_grid = None
-            
+
             elif isinstance(node.slice, ast.Slice):
                 sl = slice(get_index(node.slice.lower),
                            get_index(node.slice.upper),
@@ -591,7 +591,7 @@ def get_index(slice_value):
                     '' if sl.start is None else str(sl.start),
                     '' if sl.stop is None else str(sl.stop),
                     '' if sl.step is None else ':'+str(sl.step) )
-                
+
                 if val.grid is None:
                     out_grid = None
                 else:
@@ -603,7 +603,7 @@ def get_index(slice_value):
                         start = sl.start*val.period
                         if isinstance(off, ProcChainVar):
                             new_off = ProcChainVar(self,
-                                                   name="({}+{})".format(str(off), str(start)),
+                                                   name=f"({str(off)}+{str(start)})",
                                                    is_coord=True)
                             self._proc_managers.append(ProcessorManager(self, np.add, [off, start, new_off]))
                             off = new_off
@@ -614,7 +614,7 @@ def get_index(slice_value):
             elif isinstance(node.slice, ast.ExtSlice):
                 # TODO: implement this...
                 raise ProcessingChainError("ExtSlice still isn't implemented...")
-            
+
             # Create our return variable and set the buffer to the slice
             out = ProcChainVar(self, out_name,
                                shape = out_buf.shape[1:],
@@ -659,7 +659,7 @@ def get_index(slice_value):
 
     def _validate_name(self, name, raise_exception=False):
         """Check that name is alphanumeric, and not an already used keyword"""
-        isgood = re.match("\A\w+$", name) and name not in self.func_list \
+        isgood = re.match(r"\A\w+$", name) and name not in self.func_list \
             and not name in ureg
         if raise_exception and not isgood:
             raise ProcessingChainError(name+" is not a valid variable name")
@@ -733,7 +733,7 @@ def _round(var):
                  'round':_round, }
 
 
-########################################################################### 
+###########################################################################
 # The class that calls processors and makes sure variables are compatible #
 ###########################################################################
 class ProcessorManager:
@@ -742,7 +742,7 @@ class ProcessorManager:
     class DimInfo:
         length : int # length of arrays in this dimension
         grid : CoordinateGrid # period and offset of arrays in this dimension
-    
+
     def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types=None):
         assert isinstance(proc_chain, ProcessingChain) and callable(func) \
             and isinstance(params, list)
@@ -759,13 +759,13 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
         self.args = []
         # dict of kws -> raw values and buffers from params; we will fill this soon
         self.kwargs = {}
-        
+
         # Get the signature and list of valid types for the function
         self.signature = func.signature if signature is None else signature
         if self.signature is None:
             self.signature = ','.join(['()']*func.nin) + '->' \
                            + ','.join(['()']*func.nout)
-        
+
         # Get list of allowed type signatures
         if types is None:
             types = func.types.copy()
@@ -774,14 +774,14 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
         if not isinstance(types, list):
             types = [types]
         found_types = [ typestr.replace('->', '') for typestr in types]
-        
+
         # Make sure arrays obey the broadcasting rules, and make a dictionary
         # of the correct dimensions and unit system
-        dims_list = re.findall("\((.*?)\)", self.signature)
+        dims_list = re.findall(r"\((.*?)\)", self.signature)
 
         if not len(dims_list)==len(params)+len(kw_params):
             raise ProcessingChainError("Expected {} arguments from signature {}; found {}: ({})".format(len(dims_list), self.signature, len(params), ', '.join([str(par) for par in params])))
-        
+
         dims_dict = {} # map from dim name -> DimInfo
         outerdims = [] # list of DimInfo
         grid = None # period/offset to use for unit and coordinate conversions
@@ -809,7 +809,7 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
             for i in range(max(len(fun_dims), len(arr_dims))):
                 fd = fun_dims[-i-1] if i<len(fun_dims) else None
                 ad = arr_dims[-i-1] if i<len(arr_dims) else None
-                
+
                 if(isinstance(fd, str)):
                     if fd in dims_dict:
                         this_dim = dims_dict[fd]
@@ -821,14 +821,14 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
                             self._print(0, "Arrays of dimension", fd, "for", func.__name__, "do not have consistent period and offset!")
                     else:
                         dims_dict[fd] = self.DimInfo(ad, arr_grid)
-                
+
                 elif not fd:
                     # if we ran out of function dimensions, add a new outer dim
                     outerdims.insert(0, self.DimInfo(ad, arr_grid))
-                    
+
                 elif not ad:
                     continue
-                
+
                 elif fd.length != ad:
                     # If dimensions disagree, either insert a broadcasted array dimension or raise an exception
                     if len(fun_dims)>len(arr_dims):
@@ -837,10 +837,10 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
                         outerdims.insert(len(fun_dims)-i, self.DimInfo(ad, arr_grid))
                         fun_dims.insert(len(fun_dims)-i, ad)
                     else:
-                        raise ProcessingChainError("Failed to broadcast array dimensions for "+func.__name__+". Input arrays do not have consistent outer dimensions. Require: "+str(tuple([dim.length for dim in outerdims+fun_dims]))+"; found "+str(tuple(arr_dims)) + " for "+str(param))
+                        raise ProcessingChainError("Failed to broadcast array dimensions for "+func.__name__+". Input arrays do not have consistent outer dimensions. Require: "+str(tuple(dim.length for dim in outerdims+fun_dims))+"; found "+str(tuple(arr_dims)) + " for "+str(param))
                 elif not fd.grid:
                     outerdims[len(fun_dims)-i].grid = arr_grid
-                
+
                 elif arr_grid and fd.grid!=arr_grid:
                     self._print(0, "Arrays of dimension", fd, "for", func.__name__, "do not have consistent period and offset!")
 
@@ -871,9 +871,9 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
                 if not d in dims_dict:
                     raise ProcessingChainError("Could not deduce dimension " + d + " for " + str(param))
                 dim_list.append(dims_dict[d])
-            shape = tuple([d.length for d in dim_list])
+            shape = tuple(d.length for d in dim_list)
             this_grid = dim_list[-1].grid if dim_list else None
-            
+
             if isinstance(param, ProcChainVar):
                 # Deduce any automated descriptions of parameter
                 unit = None
@@ -892,14 +892,14 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
                 if param.is_coord and not grid:
                     grid = param._buffer[0][1]
                 param = param.get_buffer(grid)
-                
+
                 # reshape just in case there are some missing dimensions
                 arshape = list(param.shape)
                 for idim in range(-1, -1-len(shape), -1):
                     if arshape[idim]!=shape[idim]:
                         arshape.insert(len(arshape)+idim+1, 1)
                 param = param.reshape(tuple(arshape))
-                
+
             elif param is not None:
                 # Convert scalar to right type, including units
                 if isinstance(param, (Quantity, Unit)):
@@ -910,10 +910,10 @@ def __init__(self, proc_chain, func, params, kw_params={}, signature=None, types
                     param = dtype.type(round(param))
                 else:
                     param = dtype.type(param)
-            
+
             if arg_name is None: self.args.append(param)
             else: self.kwargs[arg_name] = param
-            
+
         self.proc_chain._print(2, 'Added processor:', str(self))
 
     def execute(self):
@@ -921,8 +921,8 @@ def execute(self):
 
     def __str__(self):
         return self.processor.__name__ + '(' \
-            + ", ".join([str(par) for par in self.params] + ["{}={}".format(k, str(v)) for k, v in self.kw_params.items()]) + ')'
-        
+            + ", ".join([str(par) for par in self.params] + [f"{k}={str(v)}" for k, v in self.kw_params.items()]) + ')'
+
 
 
 # A special processor manager for handling converting variables between unit systems
@@ -930,7 +930,7 @@ class UnitConversionManager(ProcessorManager):
     @vectorize(nopython=True, cache=True)
     def convert(buf_in, offset_in, offset_out, period_ratio):
         return (buf_in - offset_in) * period_ratio + offset_out
-    
+
     def __init__(self, var, unit):
         # reference back to our processing chain
         self.proc_chain = var.proc_chain
@@ -946,9 +946,9 @@ def __init__(self, var, unit):
         self.args = [from_buffer, from_grid.get_offset(),
                            unit.get_offset(), period_ratio, self.out_buffer]
         self.kwargs = {}
-        
+
         self.proc_chain._print(2, 'Added conversion:', str(self))
-        
+
 ##########################################################################
 # Now, classes to manage I/O from buffers into ProcessingChain variables #
 ##########################################################################
@@ -976,12 +976,12 @@ class NumpyIOManager(IOManager):
     def __init__(self, io_buf, var):
         assert isinstance(io_buf, np.ndarray) \
             and isinstance(var, ProcChainVar)
-        
+
         var.update_auto(dtype = io_buf.dtype,
                         shape = io_buf.shape[1:])
-        
+
         if var.shape != io_buf.shape[1:] or var.dtype != io_buf.dtype:
-            raise ProcessingChainError("numpy.array<{}>\{{}\}@{}) is not compatible with variable {}".format(self.io_buf.shape, self.io_buf.dtype, self.io_buf.data, str(self.var)))
+            raise ProcessingChainError(r"numpy.array<{}>\{{}\}@{}) is not compatible with variable {}".format(self.io_buf.shape, self.io_buf.dtype, self.io_buf.data, str(self.var)))
 
         self.io_buf = io_buf
         self.var = var
@@ -996,21 +996,21 @@ def write(self, start, end):
                   self.raw_var[0:end-start, ...], 'unsafe')
 
     def __str__(self):
-        return '{} linked to numpy.array({}, {})@{})'.format(str(self.var), self.io_buf.shape, self.io_buf.dtype, self.io_buf.data)
+        return f'{str(self.var)} linked to numpy.array({self.io_buf.shape}, {self.io_buf.dtype})@{self.io_buf.data})'
 
 
 class LGDOArrayIOManager(IOManager):
     def __init__(self, io_array, var):
         assert isinstance(io_array, lh5.Array) \
             and isinstance(var, ProcChainVar)
-        
+
         unit = io_array.attrs.get('units', None)
         var.update_auto(dtype = io_array.dtype,
                         shape = io_array.nda.shape[1:],
                         unit = unit )
-        
+
         if var.shape != io_array.nda.shape[1:] or var.dtype != io_array.dtype:
-            raise ProcessingChainError('LGDO object {}@{} is incompatible with {}'.format(self.io_buf.form_datatype(), self.raw_buf.data, str(self.var)))
+            raise ProcessingChainError(f'LGDO object {self.io_buf.form_datatype()}@{self.raw_buf.data} is incompatible with {str(self.var)}')
 
         if isinstance(var.unit, CoordinateGrid):
             if unit is None:
@@ -1018,16 +1018,16 @@ def __init__(self, io_array, var):
             elif ureg.is_compatible_with(var.unit.period, unit):
                 unit = ureg.Quantity(unit).u
             else:
-                raise ProcessingChainError("LGDO array and variable {} have incompatible units ({} and {})".format(str(var), str(var.unit.period.u), str(unit)))
-        
+                raise ProcessingChainError(f"LGDO array and variable {str(var)} have incompatible units ({str(var.unit.period.u)} and {str(unit)})")
+
         if unit is None and not var.unit is None:
             io_array.attrs['units'] = str(var.unit)
-            
+
         self.io_array = io_array
         self.raw_buf = io_array.nda
         self.var = var
         self.raw_var = var.get_buffer(unit)
-        
+
 
     def read(self, start, end):
         np.copyto(self.raw_var[0:end-start, ...],
@@ -1038,7 +1038,7 @@ def write(self, start, end):
                   self.raw_var[0:end-start, ...], 'unsafe')
 
     def __str__(self):
-        return '{} linked to {}'.format(str(self.var), str(self.io_array))
+        return f'{str(self.var)} linked to {str(self.io_array)}'
 
 
 # Waveforms
@@ -1056,13 +1056,13 @@ def __init__(self, wf_table: lh5.WaveformTable, variable: ProcChainVar):
         t0_units = wf_table.t0_units
         if dt_units is None: dt_units = t0_units
         elif t0_units is None: t0_units = dt_units
-        
+
         period = wf_table.dt_units
         if isinstance(period, str) and period in ureg:
             grid = CoordinateGrid(ureg.Quantity(self.dt_buf[0], period), self.t0_buf[0])
         else:
             grid = None
-        
+
         self.var = variable
         self.var.update_auto(shape = self.wf_buf.shape[1:],
                              dtype = self.wf_buf.dtype,
@@ -1073,15 +1073,15 @@ def __init__(self, wf_table: lh5.WaveformTable, variable: ProcChainVar):
         if dt_units is None:
             dt_units = self.var.grid.unit_str()
             t0_units = self.var.grid.unit_str()
-        
+
         self.wf_var = self.var.buffer
-        
+
         self.t0_var = self.var.grid.get_offset(t0_units)
         self.variable_t0 = isinstance(self.t0_var, np.ndarray)
         if not self.variable_t0:
             self.t0_buf[:] = self.t0_var
         self.wf_table.t0_units = t0_units
-        
+
         self.dt_buf[:] = self.var.grid.get_period(dt_units)
         self.wf_table.dt_units = dt_units
 
@@ -1092,11 +1092,11 @@ def write(self, start, end):
         self.wf_buf[start:end, ...] = self.wf_var[0:end-start, ...]
         if self.variable_t0:
             self.t0_buf[start:end, ...] = self.t0_var[0:end-start, ...]
-        
+
 
     def __str__(self):
-        return  '{} linked to <pygama.lgdo.WaveformTable: values: {}, dt: {}, t0: {}>'.format(str(self.var), str(self.wf_table.values), str(self.wf_table.dt), str(self.wf_table.t0))
-    
+        return  f'{str(self.var)} linked to <pygama.lgdo.WaveformTable: values: {str(self.wf_table.values)}, dt: {str(self.wf_table.dt)}, t0: {str(self.wf_table.t0)}>'
+
 
 
 def build_processing_chain(lh5_in, dsp_config, db_dict = None,
@@ -1104,7 +1104,7 @@ def build_processing_chain(lh5_in, dsp_config, db_dict = None,
     """
     Produces a ProcessingChain object and an lh5 table for output parameters
     from an input lh5 table and a json recipe.
-    
+
     Parameters
     ----------
     lh5_in : lgdo.Table
@@ -1152,12 +1152,12 @@ def build_processing_chain(lh5_in, dsp_config, db_dict = None,
         0: Print nothing (except errors...)
         1: Print basic warnings (default)
         2: Print basic debug info
-        3: Print friggin' everything!    
+        3: Print friggin' everything!
     block_width : int (optional)
         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
@@ -1183,7 +1183,7 @@ def build_processing_chain(lh5_in, dsp_config, db_dict = None,
 
     # prepare the processor list
     multi_out_procs = {}
-    db_parser = re.compile("db.[\w_.]+")
+    db_parser = re.compile(r"db.[\w_.]+")
     for key, node in processors.items():
         # if we have multiple outputs, add each to the processesors list
         keys = [k for k in re.split(",| ", key) if k!='']
@@ -1212,7 +1212,7 @@ def build_processing_chain(lh5_in, dsp_config, db_dict = None,
                 if arg==db_var: arg = db_node
                 else: arg = arg.replace(db_var, str(db_node))
             args[i] = arg
-        
+
         # parse the arguments list for prereqs, if not included explicitly
         if not 'prereqs' in node:
             prereqs = []
@@ -1237,7 +1237,7 @@ def resolve_dependencies(par, resolved, leafs, unresolved=[]):
         if par in resolved:
             return
         elif par in unresolved:
-            raise ProcessingChainError('Circular references detected: %s -> %s' % (par, edge))
+            raise ProcessingChainError(f'Circular references detected: {par} -> {edge}')
 
         # if we don't find a node, this is a leaf
         node = processors.get(par)
@@ -1295,12 +1295,12 @@ def resolve_dependencies(par, resolved, leafs, unresolved=[]):
             for i, name in enumerate(new_vars):
                 unit = recipe.get('unit', auto)
                 if isinstance(unit, list): unit = unit[i]
-                
+
                 proc_chain.add_variable(name, unit=unit)
 
         # get this list of kwargs
         kwargs = recipe.get('kwargs', {}) # might also need db lookup here
-        
+
         # if init_args are defined, parse any strings and then call func
         # as a factory/constructor function
         try:
@@ -1322,7 +1322,7 @@ def resolve_dependencies(par, resolved, leafs, unresolved=[]):
                                 print("Database lookup: using default value of", db_node, "for", db_var)
                         except (KeyError, TypeError):
                             raise ProcessingChainError('Did not find', db_var, 'in database, and could not find default value.')
-                
+
                     if arg==db_var: arg = db_node
                     else: arg = arg.replace(db_var, str(db_node))
 
@@ -1335,11 +1335,11 @@ def resolve_dependencies(par, resolved, leafs, unresolved=[]):
                     init_args.append(arg)
 
             if(verbosity>1):
-                print("Building function from init_args: {}({})".format(func.__name__, ", ".join(["{}".format(a) for a in init_args] + ["{}={}".format(k, v) for k, v in init_kwargs.items()])))
+                print("Building function from init_args: {}({})".format(func.__name__, ", ".join([f"{a}" for a in init_args] + [f"{k}={v}" for k, v in init_kwargs.items()])))
             func = func(*init_args)
         except KeyError:
             pass
-        
+
         proc_chain.add_processor(func, *args, **kwargs)
 
 
@@ -1353,12 +1353,11 @@ def resolve_dependencies(par, resolved, leafs, unresolved=[]):
             print("I don't know what to do with " + copy_par + ". Building output without it!")
         else:
             lh5_out.add_field(copy_par, buf_in)
-    
+
     # finally, add the output buffers to lh5_out and the proc chain
     for out_par in out_par_list:
         buf_out = proc_chain.link_output_buffer(out_par)
         lh5_out.add_field(out_par, buf_out)
-    
+
     field_mask = input_par_list + copy_par_list
     return (proc_chain, field_mask, lh5_out)
-
diff --git a/pygama/dsp/processors.py b/pygama/dsp/processors.py
index 29bc40164..2145188aa 100644
--- a/pygama/dsp/processors.py
+++ b/pygama/dsp/processors.py
@@ -1,4 +1,4 @@
-"""
+r"""
 Contains a list of DSP processors, implemented using Numba's
 :func:`numba.guvectorize` to implement NumPy's :class:`numpy.ufunc` interface.
 In other words, all of the functions are void functions whose outputs are given
@@ -59,26 +59,18 @@
    operate in place as much as possible!
 """
 
-from ._processors.pole_zero import pole_zero, double_pole_zero
-from ._processors.fixed_time_pickoff import fixed_time_pickoff
-from ._processors.time_point_thresh import time_point_thresh
-from ._processors.param_lookup import param_lookup
+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.presum import presum
-from ._processors.windower import windower
-from ._processors.bl_subtract import bl_subtract
-from ._processors.convolutions import cusp_filter, zac_filter, t0_filter
-from ._processors.trap_filters import trap_filter, trap_norm, asym_trap_filter, trap_pickoff
-from ._processors.moving_windows import moving_window_left, moving_window_right, moving_window_multi, avg_current
-from ._processors.soft_pileup_corr import soft_pileup_corr, soft_pileup_corr_bl
-from ._processors.optimize import optimize_1pz, optimize_2pz
-from ._processors.saturation import saturation
-from ._processors.gaussian_filter1d import gaussian_filter1d
-from ._processors.get_multi_local_extrema import get_multi_local_extrema
-from ._processors.multi_t_filter import multi_t_filter, remove_duplicates
+from ._processors.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.pulse_injector import inject_exp_pulse, inject_sig_pulse
 from ._processors.Wiener_filter import Wiener_filter
-from ._processors.pulse_injector import inject_sig_pulse,inject_exp_pulse
\ No newline at end of file
diff --git a/pygama/flow/datagroup.py b/pygama/flow/datagroup.py
index 2dfe38459..a1bdab996 100644
--- a/pygama/flow/datagroup.py
+++ b/pygama/flow/datagroup.py
@@ -1,12 +1,14 @@
 #!/usr/bin/env python3
-import os, json
-import pandas as pd
 import itertools
+import json
+import os
 from collections import OrderedDict
-from string import Formatter
 from pathlib import Path
-from parse import parse
 from pprint import pprint
+from string import Formatter
+
+import pandas as pd
+from parse import parse
 
 import pygama.utils as pu
 
diff --git a/pygama/hit/hit_to_evt.py b/pygama/hit/hit_to_evt.py
index 01d3518e5..ae431a239 100755
--- a/pygama/hit/hit_to_evt.py
+++ b/pygama/hit/hit_to_evt.py
@@ -1,14 +1,16 @@
 #!/usr/bin/env python3
+import argparse
 import os
 import time
-import argparse
+from pprint import pprint
+
 import h5py
 import numpy as np
 import pandas as pd
-from pprint import pprint
 
 import pygama.lh5 as lh5
 
+
 def hit_to_evt(f_hit:str, f_evt:str=None, lh5_tables:list=None, copy_cols:list=None,
                overwrite:bool=True, evt_tb_name:str=None, builder_config:dict=None):
     """
@@ -202,7 +204,7 @@ def cluster_events(tb_list:list, ts_unit:float=1e-8, ch_col:str='channel',
     #     copy_cols.extend(data_cols)
     # copy_cols = sorted(list(set(copy_cols))) # drop duplicates
 
-    # create a new dataframe where we SORT ALL ROWS by a striclty ascending timestamp
+    # create a new dataframe where we SORT ALL ROWS by a strictly ascending timestamp
     dfs = tb_list
     tcm = pd.concat(dfs).sort_values('tcm_sec')
     tcm.reset_index(inplace=True, drop=True)
diff --git a/pygama/lgdo/__init__.py b/pygama/lgdo/__init__.py
index a38a2b28d..4f9b7358b 100644
--- a/pygama/lgdo/__init__.py
+++ b/pygama/lgdo/__init__.py
@@ -29,13 +29,13 @@
 `h5py <https://www.h5py.org>`_.
 """
 
-from .lgdo_utils import *
-from .scalar import Scalar
 from .array import Array
-from .fixedsizearray import FixedSizeArray
 from .arrayofequalsizedarrays import ArrayOfEqualSizedArrays
-from .vectorofvectors import VectorOfVectors
+from .fixedsizearray import FixedSizeArray
+from .lgdo_utils import *
+from .lh5_store import LH5Store, load_dfs, load_nda
+from .scalar import Scalar
 from .struct import Struct
 from .table import Table
+from .vectorofvectors import VectorOfVectors
 from .waveform_table import WaveformTable
-from .lh5_store import LH5Store, load_nda, load_dfs
diff --git a/pygama/lgdo/array.py b/pygama/lgdo/array.py
index 6d49cffe5..63cafad7b 100644
--- a/pygama/lgdo/array.py
+++ b/pygama/lgdo/array.py
@@ -1,6 +1,8 @@
 import numpy as np
+
 from .lgdo_utils import *
 
+
 class Array:
     """
     Holds an ndarray and attributes
@@ -26,7 +28,7 @@ def __init__(self, nda=None, shape=None, dtype=None, fill_val=None, attrs={}):
         nda : ndarray (optional)
             An ndarray to be used for this object's internal array. Note: the
             array is used directly, not copied. If not supplied, internal memory
-            is newly allocated based on the shape and dtype arguments. 
+            is newly allocated based on the shape and dtype arguments.
         shape : tuple of ints (optional)
             A numpy-format shape specification for shape of the internal
             ndarray. Required if nda is None, otherwise unused.
@@ -44,7 +46,7 @@ def __init__(self, nda=None, shape=None, dtype=None, fill_val=None, attrs={}):
             if fill_val is None: nda = np.empty(shape, dtype=dtype)
             elif fill_val == 0: nda = np.zeros(shape, dtype=dtype)
             else: nda = np.full(shape, fill_val, dtype=dtype)
-        self.nda = nda 
+        self.nda = nda
         self.dtype = self.nda.dtype
         self.attrs = dict(attrs)
         if 'datatype' in self.attrs:
@@ -78,4 +80,3 @@ def resize(self, new_size):
         """Resize the array to new_size (int)"""
         new_shape = (new_size,) + self.nda.shape[1:]
         self.nda.resize(new_shape)
-
diff --git a/pygama/lgdo/arrayofequalsizedarrays.py b/pygama/lgdo/arrayofequalsizedarrays.py
index 49a15f80e..95217bf8a 100644
--- a/pygama/lgdo/arrayofequalsizedarrays.py
+++ b/pygama/lgdo/arrayofequalsizedarrays.py
@@ -1,5 +1,6 @@
-from .lgdo_utils import *
 from .array import Array
+from .lgdo_utils import *
+
 
 class ArrayOfEqualSizedArrays(Array):
     """
@@ -14,14 +15,14 @@ class ArrayOfEqualSizedArrays(Array):
 
 
     def __init__(self, dims=None, nda=None, shape=None, dtype=None, fill_val=None, attrs={}):
-        """ 
+        """
         Parameters
         ----------
         dims : tuple of ints (optional)
             specifies the dimensions required for building the
             ArrayOfEqualSizedArrays' datatype attribute
 
-        See Array.__init__ for optional args 
+        See Array.__init__ for optional args
         """
         self.dims = dims
         super().__init__(nda=nda, shape=shape, dtype=dtype, fill_val=fill_val, attrs=attrs)
@@ -44,4 +45,3 @@ def form_datatype(self):
     def __len__(self):
         """Provides __len__ for this array-like class"""
         return len(self.nda)
-
diff --git a/pygama/lgdo/fixedsizearray.py b/pygama/lgdo/fixedsizearray.py
index 171fb79f6..12d36faea 100644
--- a/pygama/lgdo/fixedsizearray.py
+++ b/pygama/lgdo/fixedsizearray.py
@@ -1,5 +1,6 @@
 from .array import Array
 
+
 class FixedSizeArray(Array):
     """
     An array of fixed-size arrays
@@ -16,9 +17,8 @@ class FixedSizeArray(Array):
     def __init__(self, nda=None, shape=None, dtype=None, fill_val=None, attrs={}):
         """ See Array.__init__ for optional args """
         super().__init__(nda=nda, shape=shape, dtype=dtype, fill_val=fill_val, attrs=attrs)
-        
+
 
     def dataype_name(self):
         """The name for this object's lh5 datatype attribute"""
         return 'fixedsize_array'
-
diff --git a/pygama/lgdo/lgdo_utils.py b/pygama/lgdo/lgdo_utils.py
index 336c7b7aa..696a4fdd5 100644
--- a/pygama/lgdo/lgdo_utils.py
+++ b/pygama/lgdo/lgdo_utils.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def get_element_type(obj):
     """Get the lgdo element type of a scalar or array
 
@@ -63,10 +64,8 @@ def parse_datatype(datatype):
     # for other datatypes, need to parse the datatype string
     from parse import parse
     datatype, element_description = parse('{}{{{}}}', datatype)
-    if datatype.endswith('>'): 
+    if datatype.endswith('>'):
         datatype, dims = parse('{}<{}>', datatype)
         dims = [int(i) for i in dims.split(',')]
         return datatype, tuple(dims), element_description
     else: return datatype, None, element_description.split(',')
-
-
diff --git a/pygama/lgdo/lh5_store.py b/pygama/lgdo/lh5_store.py
index 383c92ad8..1ea568bd8 100644
--- a/pygama/lgdo/lh5_store.py
+++ b/pygama/lgdo/lh5_store.py
@@ -1,20 +1,21 @@
-import sys, os
-import numpy as np
-import h5py
 import fnmatch
-from collections import defaultdict
+import glob
+import os
+import sys
 from bisect import bisect_left
+from collections import defaultdict
 
+import h5py
+import numpy as np
 import pandas as pd
-import glob
 
+from .array import Array
+from .arrayofequalsizedarrays import ArrayOfEqualSizedArrays
+from .fixedsizearray import FixedSizeArray
 from .lgdo_utils import *
 from .scalar import Scalar
 from .struct import Struct
 from .table import Table
-from .array import Array
-from .fixedsizearray import FixedSizeArray
-from .arrayofequalsizedarrays import ArrayOfEqualSizedArrays
 from .vectorofvectors import VectorOfVectors
 from .waveform_table import WaveformTable
 
@@ -33,7 +34,7 @@ def gimme_file(self, lh5_file, mode='r', verbosity=0):
         else: full_path = lh5_file
         if mode != 'r':
             directory = os.path.dirname(full_path)
-            if directory != '' and not os.path.exists(directory): 
+            if directory != '' and not os.path.exists(directory):
                 if verbosity > 0: print(f'making path {directory}')
                 os.makedirs(directory)
         if mode == 'r' and not os.path.exists(full_path):
@@ -47,9 +48,9 @@ def gimme_file(self, lh5_file, mode='r', verbosity=0):
 
 
     def gimme_group(self, group, base_group, grp_attrs=None, overwrite=False, verbosity=0):
-        if not isinstance(group, h5py.Group): 
+        if not isinstance(group, h5py.Group):
             if group in base_group: group = base_group[group]
-            else: 
+            else:
                 group = base_group.create_group(group)
                 if grp_attrs is not None: group.attrs.update(grp_attrs)
                 return group
@@ -63,7 +64,7 @@ def gimme_group(self, group, base_group, grp_attrs=None, overwrite=False, verbos
                 group.attrs.update(grp_attrs)
         return group
 
-            
+
     def get_buffer(self, name, lh5_file, size=None):
         """
         Returns an lh5 object appropriate for use as a pre-allocated buffer
@@ -75,7 +76,7 @@ def get_buffer(self, name, lh5_file, size=None):
 
 
 
-    def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None, 
+    def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                     field_mask=None, obj_buf=None, obj_buf_start=0, verbosity=0):
         """ Read LH5 object data from a file
 
@@ -85,7 +86,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
             Name of the lh5 object to be read (including its group path)
         lh5_file : str or h5py File object, or a list of either
             The file(s) containing the object to be read oad out. If a list of
-            files, array-like object data will be contatenated into the output
+            files, array-like object data will be concatenated into the output
             object
         start_row : int (optional)
             Starting entry for the object read (for array-like objects). For a
@@ -141,7 +142,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 if isinstance(idx, list) and len(idx) > 0 and not np.isscalar(idx[0]):
                     # a list of lists: must be one per file
                     idx_i = idx[i]
-                elif idx is not None: 
+                elif idx is not None:
                     # make idx a proper tuple if it's not one already
                     if not (isinstance(idx, tuple) and len(idx) == 1): idx = (idx,)
                     # idx is a long continuous array
@@ -180,7 +181,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
             return None, 0
 
         # make idx a proper tuple if it's not one already
-        if not (isinstance(idx, tuple) and len(idx) == 1): 
+        if not (isinstance(idx, tuple) and len(idx) == 1):
             if idx is not None: idx = (idx,)
 
         # get the object's datatype
@@ -192,7 +193,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
 
         # Scalar
         # scalars are dim-0 datasets
-        if datatype == 'scalar': 
+        if datatype == 'scalar':
             value = h5f[name][()]
             if elements == 'bool': value = np.bool(value)
             if obj_buf is not None:
@@ -234,10 +235,10 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 # fields. If they all had shared indexing, they should be in a
                 # table... Maybe should emit a warning? Or allow them to be
                 # dicts keyed by field name?
-                obj_dict[field], _ = self.read_object(name+'/'+field, 
-                                                      h5f, 
-                                                      start_row=start_row, 
-                                                      n_rows=n_rows, 
+                obj_dict[field], _ = self.read_object(name+'/'+field,
+                                                      h5f,
+                                                      start_row=start_row,
+                                                      n_rows=n_rows,
                                                       idx=idx,
                                                       verbosity=verbosity)
             # modify datatype in attrs if a field_mask was used
@@ -281,13 +282,13 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 fld_buf = None
                 if obj_buf is not None:
                     if not isinstance(obj_buf, Table) or field not in obj_buf:
-                        print("obj_buf for Table", name, 
+                        print("obj_buf for Table", name,
                               "not formatted correctly. returning new object")
                         obj_buf = None
                     else: fld_buf = obj_buf[field]
-                col_dict[field], n_rows_read = self.read_object(name+'/'+field, 
-                                                                h5f, 
-                                                                start_row=start_row, 
+                col_dict[field], n_rows_read = self.read_object(name+'/'+field,
+                                                                h5f,
+                                                                start_row=start_row,
                                                                 n_rows=n_rows,
                                                                 idx=idx,
                                                                 obj_buf=fld_buf,
@@ -312,13 +313,13 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 attrs['datatype'] =  'table' + '{' + ','.join(selected_fields) + '}'
 
             # fields have been read out, now return a table
-            if obj_buf is None: 
+            if obj_buf is None:
                 # if col_dict contains just 3 objects called t0, dt, and values,
                 # return a WaveformTable
                 if len(col_dict) == 3:
                     if 't0' in col_dict and 'dt' in col_dict and 'values' in col_dict:
-                        table = WaveformTable(t0=col_dict['t0'], 
-                                              dt=col_dict['dt'], 
+                        table = WaveformTable(t0=col_dict['t0'],
+                                              dt=col_dict['dt'],
                                               values=col_dict['values'])
                 else: table = Table(col_dict=col_dict, attrs=attrs)
                 # set (write) loc to end of tree
@@ -353,9 +354,9 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
 
             # read out cumulative_length
             cumulen_buf = None if obj_buf is None else obj_buf.cumulative_length
-            cumulative_length, n_rows_read = self.read_object(name+'/cumulative_length', 
-                                                              h5f, 
-                                                              start_row=start_row, 
+            cumulative_length, n_rows_read = self.read_object(name+'/cumulative_length',
+                                                              h5f,
+                                                              start_row=start_row,
                                                               n_rows=n_rows,
                                                               obj_buf=cumulen_buf,
                                                               obj_buf_start=obj_buf_start,
@@ -365,7 +366,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
 
             # determine the start_row and n_rows for the flattened_data readout
             da_start = 0
-            if start_row > 0 and n_rows_read > 0: 
+            if start_row > 0 and n_rows_read > 0:
                 # need to read out the cumulen sample -before- the first sample
                 # read above in order to get the starting row of the first
                 # vector to read out in flattened_data
@@ -373,7 +374,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
 
                 # check limits for values that will be used subsequently
                 if this_cumulen_nda[-1] < da_start:
-                    print("warning: cumulative_length non-increasing between entries", 
+                    print("warning: cumulative_length non-increasing between entries",
                           start_row, "and", start_row+n_rows_read, "??")
                     print(this_cumulen_nda[-1], da_start, start_row, n_rows_read)
 
@@ -385,7 +386,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
             # operations on the view change the original array because they are
             # numpy arrays, not lists.
             #
-            # First we need to substract off the in-file offset for the start of
+            # First we need to subtract off the in-file offset for the start of
             # read for flattened_data
             this_cumulen_nda -= da_start
 
@@ -398,7 +399,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 this_cumulen_nda += da_buf_start
 
             # Now prepare the object buffer if necessary
-            da_buf = None 
+            da_buf = None
             if obj_buf is not None:
                 da_buf = obj_buf.flattened_data
                 # grow da_buf if necessary to hold the data
@@ -406,17 +407,17 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 if len(da_buf) < dab_size: da_buf.resize(dab_size)
 
             # now read
-            flattened_data, dummy_rows_read = self.read_object(name+'/flattened_data', 
-                                                               h5f, 
-                                                               start_row=da_start, 
+            flattened_data, dummy_rows_read = self.read_object(name+'/flattened_data',
+                                                               h5f,
+                                                               start_row=da_start,
                                                                n_rows=da_nrows,
                                                                idx=idx,
                                                                obj_buf=da_buf,
                                                                obj_buf_start=da_buf_start,
                                                                verbosity=verbosity)
             if obj_buf is not None: return obj_buf, n_rows_read
-            return VectorOfVectors(flattened_data=flattened_data, 
-                                   cumulative_length=cumulative_length, 
+            return VectorOfVectors(flattened_data=flattened_data,
+                                   cumulative_length=cumulative_length,
                                    attrs=h5f[name].attrs), n_rows_read
 
 
@@ -457,26 +458,26 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
                 # have to apply this patch to h5py (or update h5py, if it's
                 # fixed): https://github.com/h5py/h5py/issues/1792
                 h5f[name].read_direct(obj_buf.nda, source_sel, dest_sel)
-            else: 
-                if n_rows == 0: 
+            else:
+                if n_rows == 0:
                     tmp_shape = (0,) + h5f[name].shape[1:]
                     nda = np.empty(tmp_shape, h5f[name].dtype)
                 else: nda = h5f[name][source_sel]
 
-            # special handling for bools 
+            # special handling for bools
             # (c and Julia store as uint8 so cast to bool)
             if elements == 'bool': nda = nda.astype(np.bool)
 
             # Finally, set attributes and return objects
             attrs=h5f[name].attrs
             if obj_buf is None:
-                if datatype == 'array': 
+                if datatype == 'array':
                     return Array(nda=nda, attrs=attrs), n_rows_to_read
-                if datatype == 'fixedsize_array': 
+                if datatype == 'fixedsize_array':
                     return FixedSizeArray(nda=nda, attrs=attrs), n_rows_to_read
-                if datatype == 'array_of_equalsized_arrays': 
-                    return ArrayOfEqualSizedArrays(nda=nda, 
-                                                   dims=shape, 
+                if datatype == 'array_of_equalsized_arrays':
+                    return ArrayOfEqualSizedArrays(nda=nda,
+                                                   dims=shape,
                                                    attrs=attrs), n_rows_to_read
             else:
                 if set(obj_buf.attrs.keys()) != set(attrs.keys()):
@@ -490,7 +491,7 @@ def read_object(self, name, lh5_file, start_row=0, n_rows=sys.maxsize, idx=None,
         return None
 
 
-    def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None, 
+    def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
                      wo_mode='append', write_start=0, verbosity=0):
         """Write an object into an lh5_file
 
@@ -508,12 +509,12 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
                 array is the same as append
             'overwrite_file' or 'of': delete file if present prior to writing to
                 it. write_start should be 0 (it's ignored)
-            
+
         """
         if wo_mode == 'write_safe':  wo_mode = 'w'
         if wo_mode == 'append':  wo_mode = 'a'
         if wo_mode == 'overwrite': wo_mode = 'o'
-        if wo_mode == 'overwrite_file': 
+        if wo_mode == 'overwrite_file':
             wo_mode = 'of'
             write_start = 0
         if wo_mode != 'w' and wo_mode != 'a' and wo_mode != 'o' and wo_mode != 'of':
@@ -532,10 +533,10 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
             group = self.gimme_group(name, group, grp_attrs=obj.attrs, overwrite=(wo_mode=='o'), verbosity=verbosity)
             fields = obj.keys()
             for field in obj.keys():
-                self.write_object(obj[field], 
-                                  field, 
-                                  lh5_file, 
-                                  group=group, 
+                self.write_object(obj[field],
+                                  field,
+                                  lh5_file,
+                                  group=group,
                                   start_row=start_row,
                                   n_rows=n_rows,
                                   wo_mode=wo_mode,
@@ -561,16 +562,16 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
             # cumulative lengths as appropriate for the in-file object
             offset = 0 # declare here because we have to subtract it off at the end
             if (wo_mode == 'a' or wo_mode == 'o') and 'cumulative_length' in group:
-                len_cl = len(group['cumulative_length']) 
+                len_cl = len(group['cumulative_length'])
                 if wo_mode == 'a': write_start = len_cl
                 if len_cl > 0: offset = group['cumulative_length'][write_start-1]
-            # Add offset to obj.cumulative_length itself to avoid memory allocation. 
+            # Add offset to obj.cumulative_length itself to avoid memory allocation.
             # Then subtract it off after writing!
             obj.cumulative_length.nda += offset
             self.write_object(obj.cumulative_length,
-                              'cumulative_length', 
-                              lh5_file, 
-                              group=group, 
+                              'cumulative_length',
+                              lh5_file,
+                              group=group,
                               start_row=start_row,
                               n_rows=n_rows,
                               wo_mode=wo_mode,
@@ -582,9 +583,9 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
             da_start = 0 if start_row == 0 else obj.cumulative_length.nda[start_row-1]
             da_n_rows = obj.cumulative_length.nda[n_rows-1] - da_start
             self.write_object(obj.flattened_data,
-                              'flattened_data', 
-                              lh5_file, 
-                              group=group, 
+                              'flattened_data',
+                              lh5_file,
+                              group=group,
                               start_row=da_start,
                               n_rows=da_n_rows,
                               wo_mode=wo_mode,
@@ -593,7 +594,7 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
             return
 
         # if we get this far, must be one of the Array types
-        elif isinstance(obj, Array): 
+        elif isinstance(obj, Array):
             if n_rows is None or n_rows > obj.nda.shape[0] - start_row:
                 n_rows = obj.nda.shape[0] - start_row
             nda = obj.nda[start_row:start_row+n_rows]
@@ -610,7 +611,7 @@ def write_object(self, obj, name, lh5_file, group='/', start_row=0, n_rows=None,
                 ds = group.create_dataset(name, data=nda, maxshape=maxshape)
                 ds.attrs.update(obj.attrs)
                 return
-            
+
             # Now append or overwrite
             ds = group[name]
             old_len = ds.shape[0]
@@ -642,13 +643,13 @@ def read_n_rows(self, name, lh5_file):
         datatype, shape, elements = parse_datatype(datatype)
 
         # scalars are dim-0 datasets
-        if datatype == 'scalar': 
+        if datatype == 'scalar':
             return None
 
         # recursively build a struct, return as a dictionary
         if datatype == 'struct':
             return None
-        
+
         # read a table into a dataframe
         if datatype == 'table':
             # read out each of the fields
@@ -661,12 +662,12 @@ def read_n_rows(self, name, lh5_file):
                     print('table', name, 'got strange n_rows_read', rows_read)
                     print(n_rows_read, 'was expected')
             return rows_read
-        
+
         # read out vector of vectors of different size
         if elements.startswith('array'):
             cumulen_buf = None
             return self.read_n_rows(name+'/cumulative_length', h5f)
-        
+
         # read out all arrays by slicing
         if 'array' in datatype:
             # compute the number of rows to read
@@ -678,7 +679,7 @@ def read_n_rows(self, name, lh5_file):
 def ls(lh5_file, lh5_group=''):
     """Return a list of lh5 groups in the input file and group, similar
     to ls or h5ls. Supports wildcards in group names.
-    
+
     Parameters
     ----------
     lh5_file : str
@@ -687,7 +688,7 @@ def ls(lh5_file, lh5_group=''):
         group to search
     lh5_st : LH5Store
         lh5_st object to use for searching
-    
+
     Returns
     -------
     groups : list of strs
@@ -695,23 +696,23 @@ def ls(lh5_file, lh5_group=''):
     """
 
     lh5_st = LH5Store()
-    
+
     # To use recursively, make lh5_file a h5group instead of a string
     if isinstance(lh5_file, str):
         lh5_file = lh5_st.gimme_file(lh5_file, 'r')
-    
+
     if lh5_group=='':
         lh5_group='*'
-    
+
     splitpath = lh5_group.split('/', 1)
     matchingkeys = fnmatch.filter(lh5_file.keys(), splitpath[0])
-    
+
     # if we gave a group name, go one deeper
     if len(matchingkeys)==1 and matchingkeys[0] == splitpath[0] \
        and isinstance(lh5_file[matchingkeys[0]], h5py.Group):
         splitpath.append('')
     ret = []
-    
+
     if len(splitpath)==1:
         return matchingkeys
     else:
@@ -747,9 +748,9 @@ def load_nda(f_list, par_list, lh5_group='', idx_list=None, verbose=True):
         Each entry contains the data for the specified parameter concatenated
         over all files in f_list
     """
-    if isinstance(f_list, str): 
+    if isinstance(f_list, str):
         f_list = [f_list]
-        if idx_list is not None: 
+        if idx_list is not None:
             idx_list = [idx_list]
     if idx_list is not None and len(f_list) != len(idx_list):
         print(f"load_nda: f_list len ({len(f_list)}) != idx_list len ({len(idx_list)})!")
@@ -801,17 +802,17 @@ class LH5Iterator:
     A class for iterating through one or more LH5 files, one block of entries
     at a time. This also accepts an entry list/mask to enable event selection,
     and a field mask.
-    
+
     This class can be used either for random access:
         lh5_obj, n_rows = lh5_it.read(entry)
     to read the block of entries starting at entry. In case of multiple files
     or the use of an event selection, entry refers to a global event index
-    accros files and does not count events that are excluded by the selection.
-    
+    across files and does not count events that are excluded by the selection.
+
     This can also be used as an iterator:
         for lh5_obj, entry, n_rows in LH5Iterator(...):
             do the thing!
-    
+
     This is intended for if you are reading a large quantity of data but
     want to limit your memory usage (particularly when reading in waveforms!).
     The lh5_obj that is read by this class is reused in order to avoid
@@ -839,7 +840,7 @@ def __init__(self, lh5_files, group, base_path='', entry_list=None, entry_mask=N
             Number of entries to read at a time while iterating through files
         """
         self.lh5_st = LH5Store(base_path=base_path, keep_open=True)
-        
+
         # List of files, with wildcards and env vars expanded
         if isinstance(lh5_files, str): lh5_files = [lh5_files]
         self.lh5_files = [f for f_wc in lh5_files for f in sorted(glob.glob(os.path.expandvars(f_wc)))]
@@ -847,13 +848,13 @@ def __init__(self, lh5_files, group, base_path='', entry_list=None, entry_mask=N
         self.file_map = np.array([self.lh5_st.read_n_rows(group, f) for f in self.lh5_files], 'int64').cumsum()
         self.group = group
         self.buffer_len = buffer_len
-        
+
         self.lh5_buffer = self.lh5_st.get_buffer(self.group, self.lh5_files[0], self.buffer_len) if len(self.lh5_files)>0 else None
         self.n_rows = 0
         self.current_entry = 0
 
         self.field_mask = field_mask
-        
+
         # List of entry indices from each file
         self.entry_list = None
         if entry_list is not None:
@@ -871,7 +872,7 @@ def __init__(self, lh5_files, group, base_path='', entry_list=None, entry_mask=N
                 self.entry_list = [[]]*len(self.file_map)
                 for i_file, local_list in enumerate(entry_list):
                     self.entry_list[i_file] = list(local_list)
-                    
+
         elif entry_mask is not None:
             # Convert entry mask into an entry list
             if isinstance(entry_mask, pd.Series):
@@ -899,20 +900,20 @@ def read(self, entry):
         local_entry = entry
         if i_file>0: local_entry -= self.entry_map[i_file-1]
         self.n_rows = 0
-        
+
         while(self.n_rows < self.buffer_len and i_file < len(self.file_map)):
             # Loop through files
             local_idx = self.entry_list[i_file] if self.entry_list is not None else None
             i_local = local_idx[local_entry] if local_idx is not None else local_entry
             self.lh5_buffer, n_rows = self.lh5_st.read_object(self.group, self.lh5_files[i_file], start_row = i_local, n_rows = self.buffer_len - self.n_rows, idx = local_idx, field_mask = self.field_mask, obj_buf = self.lh5_buffer, obj_buf_start = self.n_rows )
-            
+
             self.n_rows += n_rows
             i_file += 1
             local_entry = 0
 
         self.current_entry = entry
         return (self.lh5_buffer, self.n_rows)
-    
+
     def __len__(self):
         """Total number of entries"""
         return self.entry_map[-1] if len(self.entry_map)>0 else 0
diff --git a/pygama/lgdo/scalar.py b/pygama/lgdo/scalar.py
index aba55c777..9229318d0 100644
--- a/pygama/lgdo/scalar.py
+++ b/pygama/lgdo/scalar.py
@@ -1,6 +1,8 @@
 import numpy as np
+
 from .lgdo_utils import *
 
+
 class Scalar:
     """
     Holds just a value and some attributes (datatype, units, ...)
diff --git a/pygama/lgdo/struct.py b/pygama/lgdo/struct.py
index cdf2064f2..19f55fb6d 100644
--- a/pygama/lgdo/struct.py
+++ b/pygama/lgdo/struct.py
@@ -49,6 +49,6 @@ def add_field(self, name, obj):
         self.update_datatype()
 
 
-    def __len__(self): 
+    def __len__(self):
         """Structs are considered length=1 """
         return 1
diff --git a/pygama/lgdo/table.py b/pygama/lgdo/table.py
index 3b4214716..66446562d 100644
--- a/pygama/lgdo/table.py
+++ b/pygama/lgdo/table.py
@@ -1,4 +1,5 @@
 import pandas as pd
+
 from .struct import Struct
 
 
@@ -10,7 +11,7 @@ class Table(Struct):
     table I/O using functions like push_row(), is_full(), and clear()
 
     Note: if you write to a table and don't fill it up to its total size, be
-    sure to resize it before passing to data processing fucntions, as they will
+    sure to resize it before passing to data processing functions, as they will
     call len(table) to access valid data, which returns table.size.
     """
     # TODO: overload getattr to allow access to fields as object attributes?
@@ -44,7 +45,7 @@ def __init__(self, size=None, col_dict={}, attrs={}):
         # if col_dict is not empty, set size according to it
         # if size is also supplied, resize all fields to match it
         # otherwise, warn if the supplied fields have varying size
-        if len(col_dict) > 0: 
+        if len(col_dict) > 0:
             do_warn = True if size is None else False
             self.resize(new_size=size, do_warn=do_warn)
 
@@ -55,7 +56,7 @@ def __init__(self, size=None, col_dict={}, attrs={}):
         self.loc = 0
 
 
-    def datatype_name(self): 
+    def datatype_name(self):
         """The name for this lgdo's datatype attribute"""
         return 'table'
 
@@ -71,7 +72,7 @@ def resize(self, new_size = None, do_warn = False):
             if new_size is None: new_size = len(obj)
             elif len(obj) != new_size:
                 if do_warn:
-                    print('warning: resizing field', field, 
+                    print('warning: resizing field', field,
                           'with size', len(obj), '!=', new_size)
                 if isinstance(obj, Table): obj.resize(new_size, do_warn)
                 else: obj.resize(new_size)
@@ -169,7 +170,7 @@ def get_dataframe(self, cols=None, copy=False):
         """
         df = pd.DataFrame(copy=copy)
         if cols is None: cols = self.keys()
-        for col in cols: 
+        for col in cols:
             if not hasattr(self[col],'nda'):
                 print(f'column {col} does not have an nda, skipping...')
             else: df[col] = self[col].nda
diff --git a/pygama/lgdo/vectorofvectors.py b/pygama/lgdo/vectorofvectors.py
index 628b254b5..036d26e76 100644
--- a/pygama/lgdo/vectorofvectors.py
+++ b/pygama/lgdo/vectorofvectors.py
@@ -1,6 +1,8 @@
 import numpy as np
-from .lgdo_utils import *
+
 from .array import Array
+from .lgdo_utils import *
+
 
 class VectorOfVectors:
     """
@@ -8,8 +10,8 @@ class VectorOfVectors:
 
     For now only a 1D vector of 1D vectors is supported. Internal representation
     is as two ndarrays, one to store the flattened data contiguosly and one to
-    store the cumulative sum of lengths of each vector. 
-    """ 
+    store the cumulative sum of lengths of each vector.
+    """
 
 
     def __init__(self, flattened_data=None, cumulative_length=None, shape_guess=None, dtype=None, attrs={}):
@@ -82,7 +84,7 @@ def set_vector(self, i_vec, nda):
         """
         if i_vec<0 or i_vec>len(self.cumulative_length.nda)-1:
             print('VectorOfVectors: Error: bad i_vec', i_vec)
-            return 
+            return
         if len(nda.shape) != 1:
             print('VectorOfVectors: Error: nda had bad shape', nda.shape)
             return
diff --git a/pygama/lgdo/waveform_table.py b/pygama/lgdo/waveform_table.py
index 288854946..4020bb4ad 100644
--- a/pygama/lgdo/waveform_table.py
+++ b/pygama/lgdo/waveform_table.py
@@ -1,7 +1,8 @@
 import numpy as np
-from .table import Table
+
 from .array import Array
 from .arrayofequalsizedarrays import ArrayOfEqualSizedArrays
+from .table import Table
 from .vectorofvectors import VectorOfVectors
 
 
@@ -119,7 +120,7 @@ def values_units(self, units):
     @property
     def wf_len(self):
         return self.values.nda.shape[1]
-    
+
     @property
     def t0(self):
         return self['t0']
diff --git a/pygama/math/histogram.py b/pygama/math/histogram.py
index 00a5bb441..c0bb56c43 100644
--- a/pygama/math/histogram.py
+++ b/pygama/math/histogram.py
@@ -14,11 +14,12 @@
 2D hists).
 """
 
-import numpy as np
 import matplotlib.pyplot as plt
-import pygama.math.utils as pgu
+import numpy as np
 from pylab import rcParams
 
+import pygama.math.utils as pgu
+
 
 def get_hist(data, bins=None, range=None, dx=None, wts=None):
     """ return hist, bins, var after binning data
@@ -44,7 +45,7 @@ def get_hist(data, bins=None, range=None, dx=None, wts=None):
         The array of data to be histogrammed
     bins: int, array, or str (optional)
         int: the number of bins to be used in the histogram
-        array: an array of bin edges to use 
+        array: an array of bin edges to use
         str: the name of the np.histogram automatic binning algorithm to use
         If not provided, np.histogram's default auto-binning routine is used
     range : tuple (float, float) (optional)
@@ -79,7 +80,7 @@ def get_hist(data, bins=None, range=None, dx=None, wts=None):
     if wts is not None and np.shape(wts) == (): wts = np.full_like(data, wts)
     hist, bins = np.histogram(data, bins=bins, range=range, weights=wts)
 
-    if wts is None: 
+    if wts is None:
         # no weights: var = hist, but return a copy so that mods to var don't
         # modify hist.
         # Note: If you don't want a var copy, just call np.histogram()
@@ -225,7 +226,7 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
     Typically used by sending slices around a peak. Searches about argmax(hist)
     for the peak to fall by [fraction] from mx to bl
 
-    Paramaters
+    Parameters
     ----------
     fraction : float
         The fractional amplitude at which to evaluate the full width
@@ -237,16 +238,16 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
         An array of histogram variances. Used with the 'fit_slopes' method
     mx : float or tuple(float, float) (optional)
         The value to use for the max of the peak. If None, np.amax(hist) is
-        used. 
+        used.
     dmx : float (optional)
         The uncertainty in mx
     bl : float or tuple (float, float) (optional)
-        Used to specify an offset from which to estimate the FWFM. 
+        Used to specify an offset from which to estimate the FWFM.
     dbl : float (optional)
         The uncertainty in the bl
     method : string
-        'bins_over_f' : the simplest method: just take the diffence in the bin
-            centers that are over [fraction] of max. Only works for high stats and 
+        'bins_over_f' : the simplest method: just take the difference in the bin
+            centers that are over [fraction] of max. Only works for high stats and
             FWFM/bin_width >> 1
         'interpolate' : interpolate between the bins that cross the [fration]
             line.  Works well for high stats and a reasonable number of bins.
@@ -279,13 +280,13 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
     """
 
     # find bins over [fraction]
-    if mx is None: 
+    if mx is None:
         mx = np.amax(hist)
-        if var is not None and dmx == 0: 
+        if var is not None and dmx == 0:
             dmx = np.sqrt(var[np.argmax(hist)])
     idxs_over_f = hist > (bl + fraction * (mx-bl))
 
-    # argmax will return the index of the first occurence of a maximum
+    # argmax will return the index of the first occurrence of a maximum
     # so we can use it to find the first and last time idxs_over_f is "True"
     bin_lo = np.argmax(idxs_over_f)
     bin_hi = len(idxs_over_f) - np.argmax(idxs_over_f[::-1])
@@ -295,7 +296,7 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
     dheight2 = (fraction*dmx)**2 + ((1-fraction)*dbl)**2
 
     if method == 'bins_over_f':
-        # the simplest method: just take the diffence in the bin centers
+        # the simplest method: just take the difference in the bin centers
         fwfm = bin_centers[bin_hi] - bin_centers[bin_lo]
 
         # compute rough uncertainty as [bin width] (+) [dheight / slope]
@@ -311,7 +312,7 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
 
     elif method == 'interpolate':
         # interpolate between the two bins that cross the [fraction] line
-        # works well for high stats 
+        # works well for high stats
         if bin_lo < 1 or bin_hi >= len(hist)-1:
             print(f"get_fwhm: can't interpolate ({bin_lo}, {bin_hi})")
             return 0, 0
@@ -325,7 +326,7 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
         x_lo = bin_centers[bin_lo-1] + dx * dhf/dh
         # uncertainty
         dx2_lo = 0
-        if var is not None: 
+        if var is not None:
             dx2_lo = (dhf/dh)**2 * var[bin_lo] + ((dh-dhf)/dh)**2 * var[bin_lo-1]
             dx2_lo *= (dx/dh)**2
         dDdh = -dx/dh
@@ -336,12 +337,12 @@ def get_fwfm(fraction, hist, bins, var=None, mx=None, dmx=0, bl=0, dbl=0, method
         dh = hist[bin_hi] - hist[bin_hi+1]
         if dh == 0:
             raise ValueError(f"get_fwhm: interpolation failed, dh == 0")
-        x_hi = bin_centers[bin_hi] + dx * dhf/dh 
+        x_hi = bin_centers[bin_hi] + dx * dhf/dh
         if x_hi < x_lo:
             raise ValueError(f"get_fwfm: interpolation produced negative fwfm")
         # uncertainty
         dx2_hi = 0
-        if var is not None: 
+        if var is not None:
             dx2_hi = (dhf/dh)**2 * var[bin_hi+1] + ((dh-dhf)/dh)**2 * var[bin_hi]
             dx2_hi *= (dx/dh)**2
         dDdh += dx/dh
@@ -406,7 +407,7 @@ def plot_hist(hist, bins, var=None, show_stats=False, stats_hloc=0.75, stats_vlo
     """
     plot a step histogram, with optional error bars
     """
-    if fill: 
+    if fill:
         # the concat calls get the steps to draw correctly at the range boundaries
         # where="post" tells plt to draw the step y[i] between x[i] and x[i+1]
         save_color = None
@@ -434,7 +435,7 @@ def plot_hist(hist, bins, var=None, show_stats=False, stats_hloc=0.75, stats_vlo
         dmean = stddev/np.sqrt(N)
 
         mean, dmean = pgu.get_formatted_stats(mean, dmean, 2)
-        stats = r'$\mu=%s \pm %s$\n$\sigma=%#.3g$' % (mean, dmean, stddev)
+        stats = fr'$\mu={mean} \pm {dmean}$\n$\sigma={stddev:#.3g}$'
         stats_fontsize = rcParams['legend.fontsize']
         plt.text(stats_hloc, stats_vloc, stats, transform=plt.gca().transAxes, fontsize = stats_fontsize)
 
diff --git a/pygama/math/peak_fitting.py b/pygama/math/peak_fitting.py
index b6208e229..10ca8c3df 100644
--- a/pygama/math/peak_fitting.py
+++ b/pygama/math/peak_fitting.py
@@ -1,12 +1,13 @@
-import numpy as np
-from scipy.optimize import minimize, curve_fit, minimize_scalar, brentq
-from scipy.special import erf, erfc, gammaln
-from scipy.stats import crystalball
-from scipy.integrate import simps
+import math
 import sys
+
 import numba as nb
-import math
+import numpy as np
 from iminuit import Minuit, cost, util
+from scipy.integrate import simps
+from scipy.optimize import brentq, curve_fit, minimize, minimize_scalar
+from scipy.special import erf, erfc, gammaln
+from scipy.stats import crystalball
 
 import pygama.math.histogram as pgh
 
@@ -19,7 +20,7 @@ def fit_hist(func, hist, bins, var=None, guess=None,
     DEPRECATED USE FIT_BINNED
 
     do a binned fit to a histogram (nonlinear least squares).
-    can either do a poisson log-likelihood fit (jason's fave) or
+    can either do a poisson log-likelihood fit (jason's favourite) or
     use curve_fit w/ an arbitrary function.
     - hist, bins, var : as in return value of pygama.histograms.get_hist()
     - guess : initial parameter guesses. Should be optional -- we can auto-guess
@@ -44,7 +45,7 @@ def fit_hist(func, hist, bins, var=None, guess=None,
 def fit_binned(func, hist, bins, var=None, guess=None,
              cost_func='LL', Extended=True,  simplex=False, bounds=None, fixed = None):
     """
-    Do a binned fit to a histogram. 
+    Do a binned fit to a histogram.
     Default is Extended Log Likelihood fit, with option for either Least Squares or other cost function.
 
     Inputs
@@ -56,7 +57,7 @@ def fit_binned(func, hist, bins, var=None, guess=None,
 
     guess : initial guess parameters
 
-    cost_func: cost function to use 
+    cost_func: cost function to use
 
     Extended: run extended or non extended fit
 
@@ -64,20 +65,20 @@ def fit_binned(func, hist, bins, var=None, guess=None,
 
     bounds : list of tuples with bounds can be None, e.g. [(0,None), (0,10)]
 
-    fixed: list of parameter indices to fix 
+    fixed: list of parameter indices to fix
 
     Returns
     ------
     coeff, error cov_matrix : tuple(array, array, matrix)
     """
-        
-    
+
+
     if guess is None:
         print("auto-guessing not yet implemented, you must supply a guess.")
         return None, None
 
-    if cost_func == 'LL': 
-        
+    if cost_func == 'LL':
+
         if var is not None:
             t_arr = np.zeros((len(hist),2))
             t_arr[:,0] = hist
@@ -86,7 +87,7 @@ def fit_binned(func, hist, bins, var=None, guess=None,
 
         if Extended ==True:
             cost_func = cost.ExtendedBinnedNLL(hist,bins,  func)
-            
+
         else:
             cost_func = cost.BinnedNLL( hist,bins, func)
 
@@ -94,11 +95,11 @@ def fit_binned(func, hist, bins, var=None, guess=None,
 
         if var is None:
             var = hist # assume Poisson stats if variances are not provided
-        
+
         if len(bins) == len(hist)+1:
             bin_centres = pgh.get_bin_centers(bins)
 
-        # skip "okay" bins with content 0 +/- 0 
+        # skip "okay" bins with content 0 +/- 0
         # if bin content is non-zero but var = 0 let the user see the warning
         zeros = (hist == 0)
         zero_errors = (var == 0)
@@ -107,7 +108,7 @@ def fit_binned(func, hist, bins, var=None, guess=None,
         var = np.sqrt(var[mask])
         xvals = bin_centres[mask]
         cost_func = cost.LeastSquares(xvals, hist,var, func)
-            
+
     m = Minuit(cost_func, *guess)
     if bounds is not None:
         m.limits = bounds
@@ -122,11 +123,11 @@ def fit_binned(func, hist, bins, var=None, guess=None,
     return m.values, m.errors, m.covariance
 
 def fit_unbinned(func, data, guess=None,
-             Extended=True, cost_func = 'LL',simplex=False, 
+             Extended=True, cost_func = 'LL',simplex=False,
              bounds=None, fixed=None):
     """
 
-    Do a unbinned fit to data. 
+    Do a unbinned fit to data.
     Default is Extended Log Likelihood fit, with option for other cost functions.
 
     Inputs
@@ -134,13 +135,13 @@ def fit_unbinned(func, data, guess=None,
 
     func : the function to fit
 
-    data: 
+    data:
 
     guess : initial guess parameters
 
     Extended: run extended or non extended fit
 
-    cost_func: cost function to use 
+    cost_func: cost function to use
 
     simplex: whether to include a round of simpson minimisation before main minimisation
 
@@ -159,10 +160,10 @@ def fit_unbinned(func, data, guess=None,
     if cost_func =='LL':
         if Extended ==True:
             cost_func = cost.ExtendedUnbinnedNLL(data, func)
-            
+
         else:
             cost_func = cost.UnbinnedNLL(data, func)
-            
+
     m = Minuit(cost_func, *guess)
     if bounds is not None:
         m.limits = bounds
@@ -209,16 +210,16 @@ def goodness_of_fit(hist, bins, var, func, pars, method='var'):
     if method == 'Neyman' and np.any(hist==0):
         print("goodness_of_fit: hist cannot contain zeros for Neyman method")
         return 0, 0
-    
-    
+
+
     # compute expected values
     yy = func(pgh.get_bin_centers(bins), *pars) * pgh.get_bin_widths(bins)
-    
+
     if method == 'LR':
         log_lr = 2*np.sum(np.where(hist>0 , yy-hist + hist*np.log((hist+1.e-99) / (yy+1.e-99)), yy-hist))
         dof = len(hist) - len(pars)
-        return log_lr, dof   
-        
+        return log_lr, dof
+
     else:
         # compute chi2 numerator and denominator
         numerator = (hist - yy)**2
@@ -232,7 +233,7 @@ def goodness_of_fit(hist, bins, var, func, pars, method='var'):
             print(f"goodness_of_fit: unknown method {method}")
             return 0, 0
 
-        # compute chi2 and dof 
+        # compute chi2 and dof
         chisq = np.sum(numerator/denominator)
         dof = len(hist) - len(pars)
         return chisq, dof
@@ -261,7 +262,7 @@ def poisson_gof(pars, func, hist, bins, integral=None, **kwargs):
     return 2.*np.sum(mu + hist*(np.log( (hist+1.e-99) / (mu+1.e-99) ) + 1))
 
 
-def gauss_mode_width_max(hist, bins, var=None, mode_guess=None, n_bins=5, 
+def gauss_mode_width_max(hist, bins, var=None, mode_guess=None, n_bins=5,
                          cost_func='Least Squares', inflate_errors=False, gof_method='var'):
     """
     Get the max, mode, and width of a peak based on gauss fit near the max
@@ -274,7 +275,7 @@ def gauss_mode_width_max(hist, bins, var=None, mode_guess=None, n_bins=5,
     The advantage of using a gaussian over a polynomial directly is that the
     gaussian parameters are the ones we care about most for a peak, whereas for
     a poly we would have to extract them after the fit, accounting for
-    covariances. The guassian also better approximates most peaks farther down
+    covariances. The gaussian also better approximates most peaks farther down
     the peak. However, the gauss fit is nonlinear and thus less stable.
     Parameters
     ----------
@@ -305,7 +306,7 @@ def gauss_mode_width_max(hist, bins, var=None, mode_guess=None, n_bins=5,
         pars : 3-tuple 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 guassian
+            sigma : the estimated width of the peak. Equivalent to a gaussian
                 width (sigma), but based only on the curvature within n_bins of
                 the peak.  Note that the Taylor-approxiamted curvature of the
                 underlying function in the vicinity of the max is given by max /
@@ -422,7 +423,7 @@ def taylor_mode_max(hist, bins, var=None, mode_guess=None, n_bins=5, poissonLL=F
 def nb_erf(x):
 
     """
-    Numba version of error function 
+    Numba version of error function
     """
 
     y = np.empty_like(x)
@@ -434,7 +435,7 @@ def nb_erf(x):
 def nb_erfc(x):
 
     """
-    Numba version of complementary error function 
+    Numba version of complementary error function
     """
 
     y = np.empty_like(x)
@@ -464,7 +465,7 @@ def gauss_norm(x, mu, sigma):
     if sigma ==0: invs=np.nan
     else: invs = 1.0 / sigma
     z = (x - mu) * invs
-    invnorm = invs/ np.sqrt(2 * np.pi) 
+    invnorm = invs/ np.sqrt(2 * np.pi)
     return np.exp(-0.5 * z ** 2) * invnorm
 
 @nb.njit(**kwd)
@@ -483,7 +484,7 @@ def gauss_amp(x, mu, sigma, a):
 
 @nb.njit(**kwd)
 def gauss_pdf(x, mu, sigma, n_sig):
-    
+
     """
     Basic Gaussian pdf args; mu, sigma, n_sig (number of signal events)
     """
@@ -492,10 +493,10 @@ def gauss_pdf(x, mu, sigma, n_sig):
 
 
 def gauss_uniform(x, n_sig, mu, sigma, n_bkg, components = False):
-    
+
     """
-    define a gaussian signal on a uniform background, 
-    args: n_sig mu, sigma for the signal and n_bkg for the backgorund
+    define a gaussian signal on a uniform background,
+    args: n_sig mu, sigma for the signal and n_bkg for the background
     """
 
     if components==False:
@@ -505,12 +506,12 @@ def gauss_uniform(x, n_sig, mu, sigma, n_bkg, components = False):
 
 
 def gauss_linear(x, n_sig, mu, sigma, n_bkg, b, m, components=False):
-    
+
     """
     gaussian signal + linear background function
-    args: n_sig mu, sigma for the signal and n_bkg,b,m for the backgorund
+    args: n_sig mu, sigma for the signal and n_bkg,b,m for the background
     """
-    
+
 
     norm = (m/2 *np.nanmax(x)**2 + b*np.nanmax(x)) - (m/2 *np.nanmin(x)**2 + b*np.nanmin(x))
 
@@ -532,7 +533,7 @@ def step_int(x,mu,sigma, hstep):
 
 @nb.njit(**kwd)
 def unnorm_step_pdf(x,  mu, sigma, hstep):
-    
+
     """
     Unnormalised step function for use in pdfs
     """
@@ -555,7 +556,7 @@ def step_pdf(x,  mu, sigma, hstep, lower_range=np.inf , upper_range=np.inf):
         integral = step_int(np.array([np.nanmin(x), np.nanmax(x)]), mu, sigma, hstep)
     else:
         integral = step_int(np.array([lower_range, upper_range]), mu, sigma, hstep)
-    
+
     norm = integral[1]-integral[0]
     return step_f/norm
 
@@ -579,8 +580,8 @@ def step_cdf(x,mu,sigma, hstep, lower_range=np.inf , upper_range=np.inf):
 def gauss_step_pdf(x,  n_sig, mu, sigma, n_bkg, hstep, lower_range=np.inf , upper_range=np.inf, components=False):
 
     """
-    Pdf for Gaussian on step background 
-    args: n_sig mu, sigma for the signal and n_bkg,hstep for the backgorund
+    Pdf for Gaussian on step background
+    args: n_sig mu, sigma for the signal and n_bkg,hstep for the background
     """
 
     try:
@@ -600,7 +601,7 @@ def extended_gauss_step_pdf(x,  n_sig, mu, sigma, n_bkg, hstep, lower_range=np.i
 
     """
     Pdf for Gaussian on step background for Compton spectrum, returns also the total number of events for extended unbinned fits
-    args: n_sig mu, sigma for the signal and n_bkg, hstep for the backgorund
+    args: n_sig mu, sigma for the signal and n_bkg, hstep for the background
     """
 
     if components ==False:
@@ -610,10 +611,10 @@ def extended_gauss_step_pdf(x,  n_sig, mu, sigma, n_bkg, hstep, lower_range=np.i
         return n_sig+n_bkg, sig, bkg
 
 def gauss_step_cdf(x,  n_sig, mu, sigma,n_bkg, hstep, lower_range=np.inf , upper_range=np.inf, components=False):
-    
+
     """
-    Cdf for Gaussian on step background 
-    args: n_sig mu, sigma for the signal and n_bkg,hstep for the backgorund
+    Cdf for Gaussian on step background
+    args: n_sig mu, sigma for the signal and n_bkg,hstep for the background
     """
     try:
         bkg = step_cdf(x, mu, sigma, hstep, lower_range, upper_range)
@@ -630,7 +631,7 @@ def gauss_step_cdf(x,  n_sig, mu, sigma,n_bkg, hstep, lower_range=np.inf , upper
 
 @nb.njit(**kwd)
 def gauss_tail_pdf(x, mu, sigma, tau):
-    
+
     """
     A gaussian tail function template
     Can be used as a component of other fit functions w/args mu,sigma,tau
@@ -638,8 +639,8 @@ def gauss_tail_pdf(x, mu, sigma, tau):
 
     x = np.asarray(x)
     tmp = ((x-mu)/tau) + ((sigma**2)/(2*tau**2))
-    tail_f = np.where(tmp < limit, 
-                      gauss_tail_exact(x, mu, sigma, tau), 
+    tail_f = np.where(tmp < limit,
+                      gauss_tail_exact(x, mu, sigma, tau),
                       gauss_tail_approx(x, mu, sigma, tau))
     return tail_f
 
@@ -662,7 +663,7 @@ def gauss_tail_approx(x, mu, sigma, tau):
 def gauss_tail_integral(x,mu,sigma,tau):
 
     """
-    Integral for gaussian tail 
+    Integral for gaussian tail
     """
 
     abstau = np.abs(tau)
@@ -689,7 +690,7 @@ def gauss_tail_norm(x,mu,sigma,tau, lower_range=np.inf , upper_range=np.inf):
 def gauss_tail_cdf(x,mu,sigma,tau, lower_range=np.inf , upper_range=np.inf):
 
     """
-    CDF for gaussian tail 
+    CDF for gaussian tail
     """
 
     cdf = gauss_tail_integral(x,mu,sigma,tau)
@@ -702,10 +703,10 @@ def gauss_tail_cdf(x,mu,sigma,tau, lower_range=np.inf , upper_range=np.inf):
     c = 1-cdf[-1]
     return cdf+c
 
-def gauss_with_tail_pdf(x, mu, sigma,  htail,tau, components=False): 
+def gauss_with_tail_pdf(x, mu, sigma,  htail,tau, components=False):
 
     """
-    Pdf for gaussian with tail 
+    Pdf for gaussian with tail
     """
 
     if htail < 0 or htail > 1:
@@ -715,19 +716,19 @@ def gauss_with_tail_pdf(x, mu, sigma,  htail,tau, components=False):
             return np.full_like(x, np.nan, dtype='float64'), np.full_like(x, np.nan, dtype='float64')
 
     peak = gauss_norm(x,mu,sigma)
-    try: 
+    try:
         tail = gauss_tail_pdf(x, mu, sigma, tau)
     except ZeroDivisionError:
         tail = np.zeros_like(x, dtype=np.float64)
     if components ==False:
         return (1-htail)*peak + htail*tail
-    else: 
+    else:
         return (1-htail)*peak, htail*tail
 
 def gauss_with_tail_cdf(x, mu, sigma, htail,  tau, components=False):
 
     """
-    Cdf for gaussian with tail 
+    Cdf for gaussian with tail
     """
 
     if htail < 0 or htail > 1:
@@ -737,7 +738,7 @@ def gauss_with_tail_cdf(x, mu, sigma, htail,  tau, components=False):
             return np.full_like(x, np.nan, dtype='float64'), np.full_like(x, np.nan, dtype='float64')
 
     peak = gauss_cdf(x,mu,sigma)
-    try: 
+    try:
         tail = gauss_tail_cdf(x, mu, sigma, tau)
     except  ZeroDivisionError:
         tail = np.zeros_like(x, dtype=np.float64)
@@ -746,11 +747,11 @@ def gauss_with_tail_cdf(x, mu, sigma, htail,  tau, components=False):
     else:
         return (1-htail)*peak, htail*tail
 
-def radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep, 
+def radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep,
                 lower_range=np.inf , upper_range=np.inf,  components=False):
 
     """
-    David Radford's HPGe peak shape PDF consists of a gaussian with tail signal on a step background 
+    David Radford's HPGe peak shape PDF consists of a gaussian with tail signal on a step background
     """
 
     try:
@@ -766,9 +767,9 @@ def radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep,
         return pdf
     else:
         peak, tail = gauss_with_tail_pdf(x, mu, sigma, htail,  tau, components=components)
-        return n_sig *peak, n_sig*tail, n_bkg * bkg 
+        return n_sig *peak, n_sig*tail, n_bkg * bkg
 
-def extended_radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep, 
+def extended_radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep,
                          lower_range=np.inf , upper_range=np.inf, components=False):
 
     """
@@ -778,14 +779,14 @@ def extended_radford_pdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep,
     if components ==False:
         return n_sig + n_bkg, radford_pdf(x, n_sig,  mu, sigma, htail, tau, n_bkg, hstep, lower_range, upper_range)
     else:
-        peak, tail, bkg = radford_pdf(x, n_sig,  mu, sigma, htail, tau, n_bkg, hstep, 
+        peak, tail, bkg = radford_pdf(x, n_sig,  mu, sigma, htail, tau, n_bkg, hstep,
                                       lower_range, upper_range,components=components)
         return n_sig + n_bkg, peak, tail, bkg
 
 def radford_cdf(x, n_sig, mu, sigma, htail, tau, n_bkg, hstep, lower_range=np.inf , upper_range=np.inf,  components=False):
 
     """
-    Cdf for gaussian with tail signal and step background 
+    Cdf for gaussian with tail signal and step background
     """
     try:
         bkg = step_cdf(x, mu, sigma, hstep, lower_range, upper_range)
@@ -814,7 +815,7 @@ def neg_radford_peak_bgfree(E, sigma, htail, tau):
     if htail<0 or htail>1:
         print("htail outside allowed limits of 0 and 1")
         raise ValueError
-    
+
     res = minimize_scalar( neg_radford_peak_bgfree,
                            args=(sigma, htail, tau),
                            bounds=(-sigma-htail, sigma+htail) )
@@ -824,7 +825,7 @@ def neg_radford_peak_bgfree(E, sigma, htail, tau):
     # root find this to find the half-max energies
     def radford_peak_bgfree_halfmax(E, sigma, htail, tau, half_max):
         return gauss_with_tail_pdf(np.array([E]), 0, sigma, htail, tau)[0] - half_max
-    
+
     try:
         lower_hm = brentq( radford_peak_bgfree_halfmax,
                        -(2.5*sigma/2 + htail*tau), Emax,
@@ -841,9 +842,9 @@ def radford_peak_bgfree_halfmax(E, sigma, htail, tau, half_max):
         upper_hm = brentq( radford_peak_bgfree_halfmax,
                    Emax, 5*sigma,
                    args = (sigma, htail, tau, half_max) )
-    
+
     if cov is None: return upper_hm - lower_hm
-    
+
     #calculate uncertainty
     #amp set to 1, mu to 0, hstep+bg set to 0
     pars = [1,0, sigma, htail, tau,0,0]
@@ -874,10 +875,10 @@ def radford_peakshape_derivative(E, pars, step_norm):
     return n_sig*ret - n_bkg*hstep*gaus/step_norm #need norm factor for bkg
 
 def radford_parameter_gradient(E, pars, step_norm):
-    n_sig, mu, sigma, htail, tau, n_bkg, hstep = pars 
+    n_sig, mu, sigma, htail, tau, n_bkg, hstep = pars
 
-    gaus = gauss_norm(np.array([E, E-1]), mu, sigma)[0] 
-    tailL = gauss_tail_pdf(np.array([E, E-1]), mu, sigma, tau)[0] 
+    gaus = gauss_norm(np.array([E, E-1]), mu, sigma)[0]
+    tailL = gauss_tail_pdf(np.array([E, E-1]), mu, sigma, tau)[0]
     if n_bkg ==0:
         step_f = 0
     else:
@@ -887,11 +888,11 @@ def radford_parameter_gradient(E, pars, step_norm):
     y = (E-mu)/sigma
     sigtauL = sigma/tau
 
-    g_n_sig = 0.5*(htail*tailL + (1-htail)*gaus) 
+    g_n_sig = 0.5*(htail*tailL + (1-htail)*gaus)
     g_n_bkg = step_f
-    
+
     g_hs = n_bkg*math.erfc(y/np.sqrt(2))/step_norm
-    
+
     g_ht = (n_sig/2)*(tailL-gaus)
 
     #gradient of gaussian part
@@ -1000,7 +1001,7 @@ def Am_double(x,  n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3,
      - two steps (for the two lines)
      - two tails (for the two lines)
     """
-    bkg1 = n_bkg1*step_pdf(x, mu1, sigma1, hstep1, lower_range, upper_range ) 
+    bkg1 = n_bkg1*step_pdf(x, mu1, sigma1, hstep1, lower_range, upper_range )
     bkg2 = n_bkg2*step_pdf(x, mu2, sigma2, hstep2, lower_range, upper_range)
     if np.any(bkg1<0) or np.any(bkg2<0):
         return 0, np.zeros_like(x)
@@ -1011,23 +1012,23 @@ def Am_double(x,  n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3,
         return sig1+sig2+sig3+bkg1+bkg2
     else:
         return sig1,sig2,sig3,bkg1,bkg2
-    
-def extended_Am_double(x,  n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3, 
+
+def extended_Am_double(x,  n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3,
                        n_bkg1, hstep1, n_bkg2, hstep2,
                      lower_range=np.inf , upper_range=np.inf, components=False):
     if components ==False:
-        return n_sig1+n_sig2+n_sig3 + n_bkg1+n_bkg2, Am_double(n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, 
-                                                               n_sig3, mu3,sigma3, 
+        return n_sig1+n_sig2+n_sig3 + n_bkg1+n_bkg2, Am_double(n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2,
+                                                               n_sig3, mu3,sigma3,
                                                                n_bkg1, hstep1, n_bkg2, hstep2,
                                                                  lower_range, upper_range)
     else:
-        sig1,sig2,sig3,bkg1,bkg2 = Am_double(n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3, 
+        sig1,sig2,sig3,bkg1,bkg2 = Am_double(n_sig1, mu1, sigma1,  n_sig2, mu2,sigma2, n_sig3, mu3,sigma3,
                                              n_bkg1, hstep1, n_bkg2, hstep2,
                                              lower_range , upper_range,components=components)
         return n_sig1+n_sig2+n_sig3 + n_bkg1+n_bkg2, sig1,sig2,sig3,bkg1,bkg2
 
 
-def double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep, 
+def double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep,
                      lower_range=np.inf, upper_range=np.inf, components=False):
     """
     A Fit function exclusevly for a 133Ba 81keV peak situation
@@ -1041,11 +1042,11 @@ def double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep,
     sig1 = n_sig1*gauss_norm(x,mu1,sigma1)
     sig2 = n_sig2* gauss_norm(x,mu2,sigma2)
     if components == False:
-        return sig1 + sig2 + bkg 
+        return sig1 + sig2 + bkg
     else:
         return sig1, sig2, bkg
 
-def extended_double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep, 
+def extended_double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep,
                      lower_range=np.inf , upper_range=np.inf, components=False):
     """
     A Fit function exclusevly for a 133Ba 81keV peak situation
@@ -1053,13 +1054,13 @@ def extended_double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg
      - two gaussian peaks (two lines)
      - one step
      """
-    
+
     if components == False:
-        pdf = double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep, 
+        pdf = double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep,
                      lower_range, upper_range)
         return n_sig1+n_sig2+n_bkg, pdf
     else:
-        sig1, sig2, bkg = double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep, 
+        sig1, sig2, bkg = double_gauss_pdf(x,  n_sig1,  mu1, sigma1, n_sig2, mu2,sigma2,n_bkg,hstep,
                      lower_range, upper_range,components=components)
         return n_sig1+n_sig2+n_bkg, sig1, sig2, bkg
 
diff --git a/pygama/math/units.py b/pygama/math/units.py
index 224de4e10..494bcc093 100644
--- a/pygama/math/units.py
+++ b/pygama/math/units.py
@@ -1,4 +1,5 @@
-from pint import UnitRegistry, Unit, Quantity, Measurement, set_application_registry
+from pint import (Measurement, Quantity, Unit, UnitRegistry,
+                  set_application_registry)
 
 unit_registry = UnitRegistry(auto_reduce_dimensions=True)
 set_application_registry(unit_registry)
diff --git a/pygama/math/utils.py b/pygama/math/utils.py
index 78a60c11e..bafc4ddbf 100644
--- a/pygama/math/utils.py
+++ b/pygama/math/utils.py
@@ -2,10 +2,10 @@
 pygama convenience functions.
 """
 import sys
-import tqdm
-import numpy as np
-import matplotlib.pyplot as plt
 
+import matplotlib.pyplot as plt
+import numpy as np
+import tqdm
 
 
 def get_dataset_from_cmdline(args, run_db, cal_db):
@@ -69,7 +69,7 @@ def tqdm_range(start, stop, step=1, verbose=False, text=None, bar_length=20, uni
     stop : int
         ending iteration value
     step : int
-        step size inbetween each iteration
+        step size in between each iteration
     verbose : int
         verbose = 0 hides progress bar verbose > 0 displays progress bar
     text : str
@@ -110,7 +110,7 @@ def sizeof_fmt(num, suffix='B'):
     """
     for unit in ['', 'K', 'M', 'G', 'T', 'P', 'E', 'Z']:
         if abs(num) < 1024.0:
-            return "{:.3f} {}{}".format(num, unit, suffix)
+            return f"{num:.3f} {unit}{suffix}"
         num /= 1024.0
     return "{:.1f} {} {}".format(num, 'Y', suffix)
 
@@ -200,7 +200,7 @@ def tree_draw(tree, vars, tcut):
 
         for j, var in enumerate(tmp_list):
             # print(i, j, var, "getting V", j+1, "writing to np_arrs", i + j)
-            tmp = getattr(tree, "GetV{}".format(j + 1))()
+            tmp = getattr(tree, f"GetV{j + 1}")()
             np_arrs[i + j] = np.array([tmp[k] for k in range(n)])
 
     return tuple(np_arrs)
@@ -313,4 +313,3 @@ def fit_simple_scaling(x, y, var=1):
     scale_var = 1/np.sum(x*x/var)
     scale = np.sum(y*x/var) * scale_var
     return scale, scale_var
-
diff --git a/pygama/pargen/data_cleaning.py b/pygama/pargen/data_cleaning.py
index cfb71a055..0bb0fb700 100644
--- a/pygama/pargen/data_cleaning.py
+++ b/pygama/pargen/data_cleaning.py
@@ -5,10 +5,11 @@
 - find_pulser_properties (find pulser by looking for which peak has a constant time between events)
 - tag_pulsers
 """
-import numpy as np
-import matplotlib.pyplot as plt
 import matplotlib.gridspec as gs
+import matplotlib.pyplot as plt
+import numpy as np
 from scipy import stats
+
 from .peak_fitting import *
 
 
@@ -47,7 +48,7 @@ def gaussian_cut(data, cut_sigma=3, plotAxis=None):
         plotAxis.plot(bin_centers, fit, label="gaussian fit")
         plotAxis.axvline(result[0], color="g", label="fit mean")
         plotAxis.axvline(
-            cut_lo, color="r", label="+/- {} sigma".format(cut_sigma))
+            cut_lo, color="r", label=f"+/- {cut_sigma} sigma")
         plotAxis.axvline(cut_hi, color="r")
         plotAxis.legend()
         # plt.xlabel(params[i])
@@ -92,7 +93,7 @@ def xtalball_cut(data, cut_sigma=3, plotFigure=None):
         fit = xtalball(bin_centers, *result)
         plt.plot(bin_centers, fit, label="xtalball fit")
         plt.axvline(result[0], color="g", label="fit mean")
-        plt.axvline(cut_lo, color="r", label="+/- {} sigma".format(cut_sigma))
+        plt.axvline(cut_lo, color="r", label=f"+/- {cut_sigma} sigma")
         plt.axvline(cut_hi, color="r")
         plt.legend()
         # plt.xlabel(params[i])
@@ -105,7 +106,6 @@ def find_pulser_properties(df, energy="trap_max"):
 
     # print (df[energy])
     # exit()
-
     #find pulser by looking for which peak has a constant time between events
     #df should already be grouped by channel
 
diff --git a/pygama/pargen/dsp_optimize.py b/pygama/pargen/dsp_optimize.py
index d995b8b2e..6bed71868 100644
--- a/pygama/pargen/dsp_optimize.py
+++ b/pygama/pargen/dsp_optimize.py
@@ -1,14 +1,16 @@
-import numpy as np
-from .build_processing_chain import build_processing_chain
-from collections import namedtuple
-from pprint import pprint
 import multiprocessing as mp
+from collections import namedtuple
 from multiprocessing import get_context
+from pprint import pprint
+
+import numpy as np
+
+from .build_processing_chain import build_processing_chain
 
 
 def run_one_dsp(tb_data, dsp_config, db_dict=None, fom_function=None, verbosity=0, fom_kwargs=None):
     """
-    run one iteration of DSP on tb_data 
+    run one iteration of DSP on tb_data
 
     Optionally returns a value for optimization
 
@@ -40,13 +42,13 @@ def run_one_dsp(tb_data, dsp_config, db_dict=None, fom_function=None, verbosity=
     tb_out : lh5 Table
         If fom_function is None, returns the output lh5 table for the DSP iteration
     """
-    
+
     pc, tb_out = build_processing_chain(tb_data, dsp_config, db_dict=db_dict, verbosity=verbosity)
     pc.execute()
-    if fom_function is not None: 
+    if fom_function is not None:
         if fom_kwargs is not None:
             return fom_function(tb_out, verbosity, fom_kwargs)
-        else: 
+        else:
             return fom_function(tb_out, verbosity)
     else: return tb_out
 
@@ -56,7 +58,7 @@ def run_one_dsp(tb_data, dsp_config, db_dict=None, fom_function=None, verbosity=
 class ParGrid():
     """ Parameter Grid class
     Each ParGrid entry corresponds to a dsp parameter to be varied.
-    The ntuples must follow the pattern: 
+    The ntuples must follow the pattern:
     ( name parameter value_strs) : ( str, str, list of str)
     where name and parameter are the same as 'db.name.parameter' in the processing chain,
     value_strs is the array of strings to set the argument to.
@@ -87,7 +89,7 @@ def get_par_meshgrid(self, copy=False, sparse=False):
         Always uses Matrix indexing (natural for par grid) so that
         mg[i1][i2][...] corresponds to index order in self.dims
         Note copy is False by default as opposed to numpy default of True
-        """     
+        """
         axes = []
         for i in range(self.get_n_dimensions()):
             axes.append(self.dims[i].values_strs)
@@ -127,20 +129,20 @@ def print_data(self, indices):
             name, parameter, value_str = self.get_data(i_dim, i_par)
             print(f"{name}.{parameter} = {value_str}")
 
-    def set_dsp_pars(self, db_dict, indices):        
+    def set_dsp_pars(self, db_dict, indices):
         if db_dict is None:
-            db_dict = {}          
+            db_dict = {}
         for i_dim, i_par in enumerate(indices):
             name, parameter, value_str= self.get_data(i_dim, i_par)
             if name not in db_dict.keys():
                 db_dict[name] = {parameter:value_str}
             else:
-                db_dict[name][parameter] = value_str        
+                db_dict[name][parameter] = value_str
         return db_dict
 
 
 def run_grid(tb_data, dsp_config, grid, fom_function, db_dict=None, verbosity=1, **fom_kwargs):
-    """Extract a table of optimization values for a grid of DSP parameters 
+    """Extract a table of optimization values for a grid of DSP parameters
     The grid argument defines a list of parameters and values over which to run
     the DSP defined in dsp_config on tb_data. At each point, a scalar
     figure-of-merit is extracted
@@ -157,7 +159,7 @@ def run_grid(tb_data, dsp_config, grid, fom_function, db_dict=None, verbosity=1,
         list of output variables to appear in the output table for each grid point
     grid : ParGrid
         See ParGrid class for format
-    fom_function : function 
+    fom_function : function
         When given the output lh5 table of this DSP iteration, the fom_function
         must return a scalar figure-of-merit. Should accept verbosity as a
         second keyword argument
@@ -166,7 +168,7 @@ def run_grid(tb_data, dsp_config, grid, fom_function, db_dict=None, verbosity=1,
     verbosity : int (optional)
         verbosity for the processing chain and fom_function calls
 
-    **fom_kwargs : 
+    **fom_kwargs :
         Any keyword arguments for fom_function
 
 
@@ -203,7 +205,7 @@ def run_grid_point(tb_data, dsp_config, grids, fom_function, iii, db_dict=None,
         db_dict = grid.set_dsp_pars(db_dict, index)
 
     if verbosity > 1: pprint(dsp_config)
-    if verbosity > 0: 
+    if verbosity > 0:
         [grid.print_data(iii[i]) for i,grid in enumerate(grids)]
         print(' ')
     tb_out = run_one_dsp(tb_data,
@@ -226,7 +228,7 @@ def run_grid_point(tb_data, dsp_config, grids, fom_function, iii, db_dict=None,
         print("value:",res)
         out = {'indexes': [tuple(ii) for ii in iii], 'results': res}
 
-    else: 
+    else:
         out = {'indexes': [tuple(ii) for ii in iii], 'results': tb_out}
     return out
 
@@ -239,7 +241,7 @@ def get_grid_points(grid):
     complete =np.full(len(grid), False)
     while True:
         out.append([tuple(ii) for ii in iii])
-        
+
         for i, gri in enumerate(zip(iii,grid)):
             if not gri[1].iterate_indices(gri[0]):
                 print(f"{i} grid end")
@@ -248,7 +250,7 @@ def get_grid_points(grid):
         if all(complete): break
     return out
 
-def run_grid_multiprocess_parallel(tb_data, dsp_config, grid, fom_function, db_dict=None, verbosity=1, 
+def run_grid_multiprocess_parallel(tb_data, dsp_config, grid, fom_function, db_dict=None, verbosity=1,
                           processes=5, fom_kwargs=None):
 
     """
@@ -273,13 +275,13 @@ def run_grid_multiprocess_parallel(tb_data, dsp_config, grid, fom_function, db_d
     fom_function : function or None (optional)
         When given the output lh5 table of this DSP iteration, the
         fom_function must return a scalar figure-of-merit value upon which the
-        optimization will be based. Should accept verbosity as a second argument. 
-        If multiple grids provided can either pass one fom to have it run for each grid 
+        optimization will be based. Should accept verbosity as a second argument.
+        If multiple grids provided can either pass one fom to have it run for each grid
         or a list of fom to run different fom on each grid.
     verbosity : int (optional)
         verbosity for the processing chain and fom_function calls
-    fom_kwargs: 
-        any keyword arguments to pass to the fom, 
+    fom_kwargs:
+        any keyword arguments to pass to the fom,
         if multiple grids given will need to be a list of the fom_kwargs for each grid
 
     Returns:
@@ -320,7 +322,7 @@ def run_grid_multiprocess_parallel(tb_data, dsp_config, grid, fom_function, db_d
         else:
             grid_values[0][indexes[0]] = {f'{indexes[0]}':res['results']}
 
-        
+
     pool.close()
     pool.join()
-    return grid_values
\ No newline at end of file
+    return grid_values
diff --git a/pygama/pargen/ecal_th.py b/pygama/pargen/ecal_th.py
index cb56b811b..7dace63a6 100644
--- a/pygama/pargen/ecal_th.py
+++ b/pygama/pargen/ecal_th.py
@@ -1,17 +1,20 @@
-import pygama.lh5 as lh5
+import json
+import math
+import os
+import pathlib
+
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
-import os,json
-import pathlib
+import scipy.stats
+from matplotlib.backends.backend_pdf import PdfPages
 from scipy.optimize import curve_fit
+
+import pygama.analysis.calibration as cal
 #import pygama.genpar_tmp.cuts as cut
 import pygama.analysis.histograms as pgh
-import pygama.analysis.calibration as cal
 import pygama.analysis.peak_fitting as pgf
-import scipy.stats
-import math
-import matplotlib as mpl
-from matplotlib.backends.backend_pdf import PdfPages
+
 
 def fwhm_slope(x, m0, m1):
     """
@@ -26,7 +29,7 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
     """
     This is an example script for calibrating Th data.
     """
-    
+
     if isinstance(energy_params, str): energy_params = [energy_params]
 
     mpl.use('pdf')
@@ -42,7 +45,7 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
 
     Nevents = len(uncal_pass[energy_params[0]])
     print(f'{Nevents} events pass')
-    
+
     glines    = [583.191, 727.330, 860.564,1592.53,1620.50,2103.53,2614.50] # gamma lines used for calibration
     range_keV = [(20,20),(30,30), (40,40),(40,25),(25,40),(40,40),(60,60)] # side bands width
     funcs = [pgf.extended_radford_pdf,pgf.extended_radford_pdf,pgf.extended_radford_pdf,pgf.extended_radford_pdf,
@@ -74,12 +77,12 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
             pk_pars      = results['pk_pars']
             found_peaks = results['got_peaks_locs']
             fitted_peaks = results['fitted_keV']
-            
+
             for i, peak in enumerate(glines):
-                if peak not in fitted_peaks: 
+                if peak not in fitted_peaks:
                     kev_ranges[i] = (kev_ranges[i][0]-5,  kev_ranges[i][1]-5)
             for i, peak in enumerate(glines):
-                if peak not in fitted_peaks: 
+                if peak not in fitted_peaks:
                     kev_ranges[i] = (kev_ranges[i][0]-5,  kev_ranges[i][1]-5)
             for i, peak in enumerate(fitted_peaks):
                 try:
@@ -109,10 +112,10 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
             fitted_funcs = []
             fitted_gof_funcs = []
             for i, peak in enumerate(glines):
-                if peak in fitted_peaks: 
+                if peak in fitted_peaks:
                     fitted_funcs.append(funcs[i])
                     fitted_gof_funcs.append(gof_funcs[i])
-                    
+
             ecal_pass = pgf.poly(uncal_pass[energy_param], pars)
             xpb = 1
             xlo = 0
@@ -123,10 +126,10 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
             fitted_peaks = results['fitted_keV']
             pk_pars      = results['pk_pars']
             pk_covs      = results['pk_covs']
-        
+
             plot_title = f'{detector}-{measurement}-{run}'
             peaks_kev = results['got_peaks_keV']
-            
+
             pk_ranges = results['pk_ranges']
             p_vals = results['pk_pvals']
             mus = [pgf.get_mu_func(func_i, pars_i) for func_i, pars_i in zip(fitted_funcs, pk_pars)]
@@ -146,7 +149,7 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
                 counts, bs, bars = plt.hist(energies, bins=binning, histtype='step')
                 fit_vals = fitted_gof_funcs[i](bin_cs, *pk_pars[i])*np.diff(bs)
                 plt.plot(bin_cs, fit_vals)
-                plt.step(bin_cs, [(fval-count)/count if count != 0 else  (fval-count) for count, fval in zip(counts, fit_vals)] ) 
+                plt.step(bin_cs, [(fval-count)/count if count != 0 else  (fval-count) for count, fval in zip(counts, fit_vals)] )
                 plt.plot([bin_cs[10]],[0],label=get_peak_label(fitted_peaks[i]), linestyle='None' )
                 plt.plot([bin_cs[10]],[0],label = f'{fitted_peaks[i]:.1f} keV', linestyle='None')
                 plt.plot([bin_cs[10]],[0],label = f'{fwhms[i]:.2f} +- {dfwhms[i]:.2f} keV', linestyle='None')
@@ -169,11 +172,11 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
             fwhm_peaks   = np.array([], dtype=np.float32)
             indexes=[]
             for i,peak in enumerate(fitted_peaks):
-                if(peak==2103.53): 
+                if(peak==2103.53):
                     print(f"Tl SEP found at index {i}")
                     indexes.append(i)
                     continue
-                elif(peak==1592.53): 
+                elif(peak==1592.53):
                     print(f"Tl DEP found at index {i}")
                     indexes.append(i)
                     continue
@@ -189,9 +192,9 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
                 print(f'FWHM of {peak} keV peak is: {fwhms[i]:1.2f} +- {dfwhms[i]:1.2f} keV')
             param_guess  = [0.2,0.001]
             param_bounds = (0, [10., 1.])
-            fit_pars, fit_covs = curve_fit(fwhm_slope, fwhm_peaks, fwhms, sigma=dfwhms, 
+            fit_pars, fit_covs = curve_fit(fwhm_slope, fwhm_peaks, fwhms, sigma=dfwhms,
                                 p0=param_guess, bounds=param_bounds, absolute_sigma=True)
-            
+
             rng = np.random.default_rng(1)
             pars_b = rng.multivariate_normal(fit_pars, fit_covs, size=1000)
             fits = np.array([fwhm_slope(fwhm_peaks, *pars) for pars in pars_b])
@@ -223,15 +226,15 @@ def energy_cal_th(files, energy_params,  save_path, lh5_path='raw',n_events=1500
             pdf.savefig()
             plt.close()
 
-        
+
 
         output_dict[energy_param] = {'Qbb_fwhm': round(fit_qbb,2), 'Qbb_fwhm_err': round(qbb_err,2),
-                                    '2.6_fwhm': round(fwhms[-1],2), '2.6_fwhm_err': round(dfwhms[-1],2), 
-                                    "m0":fit_pars[0], "m1":fit_pars[1], 
+                                    '2.6_fwhm': round(fwhms[-1],2), '2.6_fwhm_err': round(dfwhms[-1],2),
+                                    "m0":fit_pars[0], "m1":fit_pars[1],
                                     "Calibration_pars":pars.tolist(),
                                     "Number_events": Nevents}
 
-    
+
     dict_save_path = os.path.join(save_path, f'{detector}.json')
     with open(dict_save_path,'w') as fp:
         json.dump(output_dict,fp, indent=4)
@@ -248,7 +251,7 @@ def get_peak_labels(labels, pars):
     return out_labels, out
 
 def get_peak_label(peak):
-    if peak == 583.191: 
+    if peak == 583.191:
         return 'Tl 583'
     elif peak == 727.33:
         return 'Bi 727'
@@ -261,4 +264,4 @@ def get_peak_label(peak):
     elif peak == 2103.53:
         return 'Tl SEP'
     elif peak == 2614.5:
-        return 'Tl FEP'
\ No newline at end of file
+        return 'Tl FEP'
diff --git a/pygama/pargen/energy_cal.py b/pygama/pargen/energy_cal.py
index ee6d61839..b67534bf6 100644
--- a/pygama/pargen/energy_cal.py
+++ b/pygama/pargen/energy_cal.py
@@ -6,14 +6,17 @@
 - hpge_E_calibration (main routine -- finds and fits peaks specified)
 """
 import sys
-import numpy as np
-import pygama.math.utils as pgu
-import pygama.math.peak_fitting as pgf
-import pygama.math.histogram as pgh
-import matplotlib.pyplot as plt
+
 import matplotlib.gridspec as gs
-from scipy.signal import medfilt, find_peaks_cwt
+import matplotlib.pyplot as plt
+import numpy as np
 import scipy.stats
+from scipy.signal import find_peaks_cwt, medfilt
+
+import pygama.math.histogram as pgh
+import pygama.math.peak_fitting as pgf
+import pygama.math.utils as pgu
+
 
 def hpge_find_E_peaks(hist, bins, var, peaks_keV, n_sigma=5, deg=0, Etol_keV=None, var_zero=1, verbose=False):
     """ Find uncalibrated E peaks whose E spacing matches the pattern in peaks_keV
@@ -62,7 +65,7 @@ def hpge_find_E_peaks(hist, bins, var, peaks_keV, n_sigma=5, deg=0, Etol_keV=Non
     if Etol_keV is None:
         #estimate Etol_keV
         pt_pars, pt_covs = hpge_fit_E_peak_tops(hist, bins, var, detected_max_locs, n_to_fit=15)
-        if sum(sum([sum(c) if c is not None else 0 for c in pt_covs])) == np.inf or sum(sum([sum(c) if c is not None else 0 for c in pt_covs])) == 0: print('hpge_find_E_peaks: can safely ignore previous covariance warning, not used')
+        if sum(sum(sum(c) if c is not None else 0 for c in pt_covs)) == np.inf or sum(sum(sum(c) if c is not None else 0 for c in pt_covs)) == 0: print('hpge_find_E_peaks: can safely ignore previous covariance warning, not used')
         pt_pars = pt_pars[np.array([x is not None for x in pt_pars])]
         med_sigma_ratio = np.median(np.stack(pt_pars)[:,1]/np.stack(pt_pars)[:,0])
 
@@ -193,7 +196,7 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
         The function to be fit to the peak in the (windowed) hist
     """
     if  func == pgf.gauss_step_cdf or func == pgf.gauss_step_pdf or func == pgf.extended_gauss_step_pdf:
-        # get mu and hieght from a gauss fit, also sigma as fallback
+        # get mu and height from a gauss fit, also sigma as fallback
         pars, cov = pgf.gauss_mode_width_max(hist, bins, var)
         bin_centres = pgh.get_bin_centers(bins)
         if pars is None:
@@ -224,9 +227,9 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
                     return []
 
         # now compute amp and return
-        n_sig = np.sum(hist[(bin_centres>mu-3*sigma)&(bin_centres<mu+3*sigma)]) 
+        n_sig = np.sum(hist[(bin_centres>mu-3*sigma)&(bin_centres<mu+3*sigma)])
         n_bkg = np.sum(hist)-n_sig
-        
+
         hstep = step /(bg + np.mean(hist[:10]))
         return [n_sig, mu, sigma/2,n_bkg, hstep,bins[0], bins[-1],0]
 
@@ -249,7 +252,7 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
 
         # get bg and step from edges of hist
         bg0 = np.mean(hist[-10:])
-        step = bg0 -np.mean(hist[:10]) 
+        step = bg0 -np.mean(hist[:10])
 
         # get sigma from fwfm with f = 1/sqrt(e)
         try:
@@ -257,7 +260,7 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
             if sigma == 0: raise ValueError
         except:
             sigma = pgh.get_fwfm(0.6065, hist, bins, var, mx=height, bl=bg0+step/2, method='fit_slopes')[0]
-            if sigma == 0: 
+            if sigma == 0:
                 print("get_hpge_E_peak_par_guess: sigma estimation failed")
                 if sigma_guess is not None:
                     sigma = sigma_guess
@@ -268,10 +271,10 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
         # for now hard-coded
         htail = 1./5
         tau = 0.5*sigma
-        
+
         hstep = step /(bg0 + np.mean(hist[:10]))
-        
-        n_sig = np.sum(hist[(bin_centres>mu-3*sigma)&(bin_centres<mu+3*sigma)]) 
+
+        n_sig = np.sum(hist[(bin_centres>mu-3*sigma)&(bin_centres<mu+3*sigma)])
         n_bkg = np.sum(hist)-n_sig
 
         parguess = [n_sig,  mu, sigma, htail, tau, n_bkg,hstep,bins[0], bins[-1],0]
@@ -281,7 +284,7 @@ def get_hpge_E_peak_par_guess(hist, bins, var, func):
     else:
         print(f'get_hpge_E_peak_par_guess not implemented for {func.__name__}')
         return []
-        
+
 
 def get_hpge_E_fixed( func):
     """
@@ -305,16 +308,16 @@ def get_hpge_E_bounds(func):
     if  func == pgf.radford_cdf or func == pgf.radford_pdf or func == pgf.extended_radford_pdf:
         return [(0,None), (None,None), (None,None), (0,1),(None,None),(0,None), (None,None)
                 ,(None,None),(None,None),(None,None)]
-    
+
     elif  func ==  pgf.gauss_step_cdf or func ==  pgf.gauss_step_pdf or func == pgf.extended_gauss_step_pdf :
         return [(0,None), (None,None), (None,None),(0,None), (None,None),(None,None)
                 ,(None,None),(None,None)]
-    
+
     else:
         print(f'get_hpge_E_bounds not implemented for {func.__name__}')
         return []
 
-def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_step_cdf, 
+def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_step_cdf,
                      method = 'unbinned', gof_funcs=None, n_events=15000, allowed_p_val= 0.05,
                      uncal_is_int=False, simplex=False):
     """ Fit the Energy peaks specified using the function given
@@ -331,7 +334,7 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
         array of number of bins to use for the fit window histogramming
     funcs : function or array of functions
         funcs to be used to fit each region
-    method: str 
+    method: str
             default is unbinned fit can specify to use binned fit method instead
     gof_funcs : function or array of functions
             functions to use for calculation goodness of fit if unspecified will use same func as fit
@@ -340,7 +343,7 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
         E_uncal
     simplex : bool determining whether to do a round of simpson minimisation before gradient minimisation
     n_events : int number of events to use for unbinned fit
-    allowed_p_val: lower limit on p_val of fit 
+    allowed_p_val: lower limit on p_val of fit
     Returns
     -------
     pars : list of array
@@ -370,13 +373,13 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
             gof_func_i = gof_funcs[i_peak] if hasattr(gof_funcs, '__len__') else gof_funcs
         else:
             gof_func_i = func_i
-        
+
         try:
             # bin a histogram
             Euc_min = mode_guesses[i_peak] - wleft_i
             Euc_max = mode_guesses[i_peak] + wright_i
             Euc_min, Euc_max, n_bins_i = pgh.better_int_binning(x_lo=Euc_min, x_hi=Euc_max, n_bins=n_bins_i)
-            
+
             if method =='unbinned':
                 energies = E_uncal[(E_uncal>Euc_min)&(E_uncal<Euc_max)][:n_events]
                 hist, bins, var = pgh.get_hist(energies, bins=n_bins_i, range=(Euc_min,Euc_max))
@@ -391,25 +394,25 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
                 par_guesses = get_hpge_E_peak_par_guess(hist, bins, var, func_i)
                 bounds = get_hpge_E_bounds(func_i)
                 fixed, mask = get_hpge_E_fixed(func_i)
-                pars_i,errs_i, cov_i = pgf.fit_binned(func_i, hist,  bins, var=var, 
-                                        guess=par_guesses, cost_func=method, Extended=True, 
+                pars_i,errs_i, cov_i = pgf.fit_binned(func_i, hist,  bins, var=var,
+                                        guess=par_guesses, cost_func=method, Extended=True,
                                                     fixed=fixed, simplex=simplex, bounds = bounds)
-            
+
             csqr = pgf.goodness_of_fit(hist, bins, None, gof_func_i, pars_i , method='LR')
             p_val = scipy.stats.chi2.sf(csqr[0], csqr[1])
 
-            
+
             pars_i = np.array(pars_i)[mask]
             errs_i = np.array(errs_i)[mask]
             cov_i = np.array(cov_i)[mask,:][:,mask]
-            
+
             total_events = pgf.get_total_events_func(func_i, pars_i, errors=errs_i)
-            if (sum([sum(c) if c is not None else 0 for c in cov_i]) == np.inf or 
-                sum([sum(c) if c is not None else 0 for c in cov_i]) == 0 or 
-                np.isnan(sum([sum(c) if c is not None else 0 for c in cov_i]))) :
+            if (sum(sum(c) if c is not None else 0 for c in cov_i) == np.inf or
+                sum(sum(c) if c is not None else 0 for c in cov_i) == 0 or
+                np.isnan(sum(sum(c) if c is not None else 0 for c in cov_i))) :
                 print(f'hpge_fit_E_peaks: cov estimation failed for i_peak={i_peak} at loc {mode_guesses[i_peak]:g}')
                 pars_i, errs_i, cov_i, p_val = None, None, None, None
-            
+
             elif (np.abs(errs_i/pars_i) < 1e-7).any():
                 print(f'hpge_fit_E_peaks: cov estimation failed for i_peak={i_peak} at loc {mode_guesses[i_peak]:g}, parameter error too low')
                 pars_i, errs_i, cov_i, p_val = None, None, None, None
@@ -417,13 +420,13 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
             elif np.abs(total_events[0] -np.sum(hist))/np.sum(hist)>0.1:
                 print(f'hpge_fit_E_peaks: fit failed for i_peak={i_peak} at loc {mode_guesses[i_peak]:g}, total_events is outside limit')
                 pars_i, errs_i, cov_i, p_val = None, None, None, None
-                
+
             elif p_val<allowed_p_val:
                 print(f'hpge_fit_E_peaks: fit failed for i_peak={i_peak}, p-value too low: {p_val}')
                 pars_i, errs_i, cov_i, p_val = None, None, None, None
 
         except: pars_i, errs_i, cov_i, p_val = None, None, None, None
-        
+
         #get binning
         binw_1 = (bins[-1]-bins[0])/(len(bins)-1)
 
@@ -434,7 +437,7 @@ def hpge_fit_E_peaks(E_uncal, mode_guesses, wwidths, n_bins=50, funcs=pgf.gauss_
         ranges.append([Euc_min, Euc_max])
         p_vals.append(p_val)
 
-    return (np.array(pars, dtype = object), np.array(errors, dtype = object), 
+    return (np.array(pars, dtype = object), np.array(errors, dtype = object),
         np.array(covs, dtype = object), np.array(binws), np.array(ranges), np.array(p_vals))
 
 
@@ -509,7 +512,7 @@ def hpge_fit_E_cal_func(mus, mu_vars, Es_keV, E_scale_pars, deg=0):
     return pars, cov
 
 
-def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False, range_keV=None, 
+def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False, range_keV=None,
                         funcs=pgf.gauss_step_cdf, gof_funcs = None, method = 'unbinned', gof_func =None,
                         n_events=15000, simplex=False, allowed_p_val=0.05, verbose=True):
     """ Calibrate HPGe data to a set of known peaks
@@ -526,7 +529,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
         degree of the polynomial for the E_cal function E_keV = poly(E_uncal).
         deg = 0 corresponds to a simple scaling E_keV = scale * E_uncal.
         Otherwise follows the convention in np.polyfit
-    method: str 
+    method: str
             default is unbinned fit can specify to use binned fit method instead
     gof_funcs : function or array of functions
             functions to use for calculation goodness of fit if unspecified will use same func as fit
@@ -537,7 +540,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
         ranges around which the peak fitting is performed
         if tuple(s) are supplied, they provide the left and right ranges
     n_events : int number of events to use for unbinned fit
-    allowed_p_val: lower limit on p_val of fit 
+    allowed_p_val: lower limit on p_val of fit
     Returns
     -------
     pars, cov : array, 2D array
@@ -599,7 +602,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
             print(f'\t{i}'.ljust(4) + str(Ei).ljust(9) + f'| {Li:g}'.ljust(5))
 
     # re-bin the histogram in ~0.2 keV bins with updated E scale par for peak-top fits
-    Euc_min, Euc_max = [(np.poly1d(roughpars)-i).roots for i in (peaks_keV[0]*.9, peaks_keV[-1]*1.1)]
+    Euc_min, Euc_max = ((np.poly1d(roughpars)-i).roots for i in (peaks_keV[0]*.9, peaks_keV[-1]*1.1))
     Euc_min = Euc_min[np.logical_and(Euc_min >= 0, Euc_min <= max(Euc_max))][0]
     Euc_max = Euc_max[np.logical_and(Euc_max >= Euc_min, Euc_max <= max(E_uncal))][0]
     dEuc = 0.2/roughpars[-2]
@@ -660,7 +663,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
         range_uncal = [(r[0]/d, r[1]/d) if isinstance(r, tuple) else r/d for r, d in zip(range_keV, der)]
         n_bins = [sum(r)/0.5/d if isinstance(r, tuple) else r/0.2/d for r, d in zip(range_keV, der)]
 
-    pk_pars,pk_errors, pk_covs, pk_binws, pk_ranges, pk_pvals = hpge_fit_E_peaks(E_uncal, got_peaks_locs, range_uncal, 
+    pk_pars,pk_errors, pk_covs, pk_binws, pk_ranges, pk_pvals = hpge_fit_E_peaks(E_uncal, got_peaks_locs, range_uncal,
                                     n_bins=n_bins,
                                     funcs=funcs, method=method, gof_funcs = gof_funcs, n_events=n_events,
                                     uncal_is_int=False, simplex=simplex, allowed_p_val=allowed_p_val)
@@ -720,7 +723,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
     derco = np.polyder(np.poly1d(pars)).coefficients
     der = [pgf.poly(Ei, derco) for Ei in fitted_peaks_keV]
 
-    cal_fwhms = uncal_fwhms * der 
+    cal_fwhms = uncal_fwhms * der
     cal_fwhms_errs = uncal_fwhm_errs*der
     results['pk_fwhms'] = np.asarray([(u*d, e*d) for u, e, d in zip(uncal_fwhms, uncal_fwhm_errs, der)])
 
@@ -730,7 +733,7 @@ def hpge_E_calibration(E_uncal, peaks_keV, guess_keV, deg=0, uncal_is_int=False,
         for i, (Ei, fwhm, fwhme) in enumerate(zip(fitted_peaks_keV, cal_fwhms, cal_fwhms_errs)):
             print(f'\t{i}'.ljust(4) + str(Ei).ljust(9) + f'| {fwhm:.2f}+-{fwhme:.2f} keV'.ljust(5))
 
-    return   pars, cov, results  
+    return   pars, cov, results
 
 
 
@@ -988,6 +991,7 @@ def match_peaks(data_pks, cal_pks):
     Match uncalibrated peaks with literature energy values.
     """
     from itertools import combinations
+
     from scipy.stats import linregress
 
     n_pks = len(cal_pks) if len(cal_pks) < len(data_pks) else len(data_pks)
@@ -1013,10 +1017,10 @@ def match_peaks(data_pks, cal_pks):
     print(i, best_err)
     print("cal:",cal)
     print("data:",data)
-    plt.scatter(data, cal, label='min.err:{:.2e}'.format(err))
+    plt.scatter(data, cal, label=f'min.err:{err:.2e}')
     xs = np.linspace(data[0], data[-1], 10)
     plt.plot(xs, best_m * xs + best_b , c="r",
-             label="y = {:.2f} x + {:.2f}".format(best_m,best_b) )
+             label=f"y = {best_m:.2f} x + {best_b:.2f}" )
     plt.xlabel("Energy (uncal)", ha='right', x=1)
     plt.ylabel("Energy (keV)", ha='right', y=1)
     plt.legend()
@@ -1103,8 +1107,8 @@ def calibrate_tl208(energy_series, cal_peaks=None, plotFigure=None):
         # plt.figure()
         # plt.plot(bin_centers,peak_hist,  color="k", ls="steps")
 
-        # inpu = input("q to quit...")
-        # if inpu == "q": exit()
+        # inp = input("q to quit...")
+        # if inp == "q": exit()
 
         try:
             guess_e, guess_sigma, guess_area = get_gaussian_guess(
@@ -1127,8 +1131,8 @@ def calibrate_tl208(energy_series, cal_peaks=None, plotFigure=None):
             gauss(bin_centers, guess_e, guess_sigma, guess_area),
             color="b")
 
-        # inpu = input("q to quit...")
-        # if inpu == "q": exit()
+        # inp = input("q to quit...")
+        # if inp == "q": exit()
 
         bounds = ([0.9 * guess_e, 0.5 * guess_sigma, 0, 0, 0, 0, 0], [
             1.1 * guess_e, 2 * guess_sigma, 0.1, 0.75, window_width_in_adc, 10,
@@ -1143,8 +1147,8 @@ def calibrate_tl208(energy_series, cal_peaks=None, plotFigure=None):
 
         plt.plot(bin_centers, radford_peak(bin_centers, *params), color="r")
 
-        # inpu = input("q to quit...")
-        # if inpu == "q": exit()
+        # inp = input("q to quit...")
+        # if inp == "q": exit()
 
         fit_result_map[energy] = params
         centers[i] = params[0]
diff --git a/pygama/pargen/mse_psd.py b/pygama/pargen/mse_psd.py
index 2a2578776..ecaabb75e 100644
--- a/pygama/pargen/mse_psd.py
+++ b/pygama/pargen/mse_psd.py
@@ -3,9 +3,9 @@
 - get_ae_cut (does A/E)
 """
 
-import numpy as np
-import matplotlib.pyplot as plt
 import matplotlib.gridspec as gs
+import matplotlib.pyplot as plt
+import numpy as np
 from matplotlib.colors import LogNorm
 
 from pygama.analysis.peak_fitting import *
diff --git a/pygama/raw/build_raw.py b/pygama/raw/build_raw.py
index 785b50b04..7b2a2dd14 100644
--- a/pygama/raw/build_raw.py
+++ b/pygama/raw/build_raw.py
@@ -1,10 +1,18 @@
-import os, time, sys, glob, json, tqdm
+import glob
+import json
+import os
+import sys
+import time
+
 import numpy as np
+import tqdm
 
 from pygama import lgdo
 from pygama.math.utils import sizeof_fmt
-from .raw_buffer import RawBuffer, RawBufferList, RawBufferLibrary, write_to_lh5_and_clear
+
 from .fc.fc_streamer import FCStreamer
+from .raw_buffer import (RawBuffer, RawBufferLibrary, RawBufferList,
+                         write_to_lh5_and_clear)
 
 #from orca.stream_orca import *
 #from stream_llama import *
@@ -12,9 +20,9 @@
 #from stream_fc import *
 
 
-def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192, 
+def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
               n_max=np.inf, overwrite=True, verbosity=2, **kwargs):
-    """ Convert data into LEGEND hdf5 `raw` format.  
+    """ Convert data into LEGEND hdf5 `raw` format.
 
     Takes an input stream (in_stream) of a given type (in_stream_type) and
     writes to output file(s) according to the user's a specification (out_spec).
@@ -29,11 +37,11 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
         Type of stream used to write the input file.
         Options are 'ORCA', 'FlashCams', 'LlamaDaq', 'Compass', 'MGDO'
     out_spec : str or json dict or RawBufferLibrary or None
-        Specification for the output stream. 
+        Specification for the output stream.
         - If None, uses '{in_stream}.hdf5' as the output filename.
-        - If a str not ending in '.json', interpretted as the output filename.
-        - If a str ending in '.json', interpretted as a filename containing
-          json-shorthand for the output sepcification (see raw_buffer.py)
+        - If a str not ending in '.json', interpreted as the output filename.
+        - If a str ending in '.json', interpreted as a filename containing
+          json-shorthand for the output specification (see raw_buffer.py)
         - If a json dict, should be a dict loaded from the json shorthand
           notation for RawBufferLibraries (see raw_buffer.py), which is then
           used to build a RawBufferLibrary
@@ -53,12 +61,12 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
 
     # convert any environment variables in in_stream so that we can check for readability
     in_stream = os.path.expandvars(in_stream)
-    # later: fix if in_stream is not a file 
+    # later: fix if in_stream is not a file
     # (e.g. a socket, which exists if open and has size = np.inf)
     if not os.path.exists(in_stream):
         print(f'Error: file {in_stream} not found')
         return
-    in_stream_size = os.stat(in_stream).st_size 
+    in_stream_size = os.stat(in_stream).st_size
 
     # try to guess the input stream type if it's not provided
     if in_stream_type is None:
@@ -76,11 +84,11 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
     rb_lib = None
     if isinstance(out_spec, str) and out_spec.endswith('.json'):
         with open(out_spec) as json_file: out_spec = json.load(json_file)
-    if isinstance(out_spec, dict): 
+    if isinstance(out_spec, dict):
         out_spec = RawBufferLibrary(json_dict=out_spec, kw_dict=kwargs)
     if isinstance(out_spec, RawBufferLibrary): rb_lib = out_spec
     # if no rb_lib, write all data to file
-    if out_spec is None: 
+    if out_spec is None:
         out_spec = in_stream
         i_ext = out_spec.rfind('.')
         if i_ext != -1: out_spec = out_spec[:i_ext]
@@ -96,27 +104,27 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
         return
     if buffer_size > n_max: buffer_size = n_max
 
-    # ouput start of processing info if verbosity > 0
+    # output start of processing info if verbosity > 0
     if verbosity > 0:
         print( 'Starting build_raw processing.')
-        print(f'  Input: {in_stream}')        
+        print(f'  Input: {in_stream}')
         out_files = [out_spec]
-        if isinstance(out_spec, RawBufferLibrary): 
+        if isinstance(out_spec, RawBufferLibrary):
             out_files = out_spec.get_list_of('out_stream')
         if len(out_files) == 1: print(f'  Output: {out_files[0]}')
-        else: 
+        else:
             print(f'  Output:')
             for out_file in out_files: print(f'- {out_file}')
-        print(f'  Buffer size: {buffer_size}')  
-        print(f'  Max num. events: {n_max}')    
-        if verbosity > 1: 
+        print(f'  Buffer size: {buffer_size}')
+        print(f'  Max num. events: {n_max}')
+        if verbosity > 1:
             if n_max < np.inf: progress_bar = tqdm.tqdm(total=n_max, unit='rows')
             else: progress_bar = tqdm.tqdm(total=in_stream_size, unit='B', unit_scale=True)
 
     # start a timer and a byte counter
     t_start = time.time()
 
-    # select the approprate streamer for in_stream
+    # select the appropriate streamer for in_stream
     streamer = None
     if in_stream_type == 'ORCA':
         print(f'Error: ORCA streaming not yet implemented')
@@ -139,8 +147,8 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
     # initialize the stream and read header. Also initializes rb_lib
     if verbosity > 1: progress_bar.update(0)
     out_stream = out_spec if isinstance(out_spec, str) else ''
-    header_data = streamer.open_stream(in_stream, rb_lib=rb_lib, buffer_size=buffer_size, 
-                                       chunk_mode='full_only', out_stream=out_stream, 
+    header_data = streamer.open_stream(in_stream, rb_lib=rb_lib, buffer_size=buffer_size,
+                                       chunk_mode='full_only', out_stream=out_stream,
                                        verbosity=verbosity)
     rb_lib = streamer.rb_lib
     if verbosity > 1 and n_max == np.inf: progress_bar.update(streamer.n_bytes_read)
@@ -169,7 +177,7 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
     n_bytes_last = streamer.n_bytes_read
     while True:
         chunk_list = streamer.read_chunk(verbosity=verbosity-2)
-        if verbosity > 1 and n_max == np.inf: 
+        if verbosity > 1 and n_max == np.inf:
             progress_bar.update(streamer.n_bytes_read-n_bytes_last)
             n_bytes_last = streamer.n_bytes_read
         if len(chunk_list) == 0: break
@@ -186,9 +194,9 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
 
     if verbosity > 0:
         elapsed = time.time() - t_start
-        print("Time elapsed: {:.2f} sec".format(elapsed))
+        print(f"Time elapsed: {elapsed:.2f} sec")
         out_files = rb_lib.get_list_of('out_stream')
-        if len(out_files) == 1: 
+        if len(out_files) == 1:
             out_file = out_files[0]
             colpos = out_file.find(':')
             if colpos != -1: out_file = out_file[:colpos]
@@ -196,7 +204,7 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
                 file_size = os.stat(out_file).st_size
                 print(f"Output file: {out_file} ({sizeof_fmt(file_size)})")
             else: print("Output file: {out_file} (not written)")
-        else: 
+        else:
             print("Output files:")
             for out_file in out_files:
                 colpos = out_file.find(':')
@@ -209,4 +217,3 @@ def build_raw(in_stream, in_stream_type=None, out_spec=None, buffer_size=8192,
         print(f"Conversion speed: {sizeof_fmt(streamer.n_bytes_read/elapsed)}ps")
 
         print('Done.\n')
-
diff --git a/pygama/raw/data_decoder.py b/pygama/raw/data_decoder.py
index ddea13402..49c20d8dd 100644
--- a/pygama/raw/data_decoder.py
+++ b/pygama/raw/data_decoder.py
@@ -1,11 +1,13 @@
 """
 data_decoder.py
 
-base classe for decoding data into raw lgdo tables / files
+base classes for decoding data into raw lgdo tables / files
 """
 
 from abc import ABC
+
 import numpy as np
+
 from pygama import lgdo
 
 
@@ -45,7 +47,7 @@ def __init__(self, garbage_length=256, packet_size_guess=1024):
                                      lgdo.Array(shape=garbage_length, dtype='uint32'))
 
 
-    def get_key_list(self): 
+    def get_key_list(self):
         """ overload with list of keys for this decoder.
         e.g. `return range(n_channels)`
         The default version works for decoders with single / no keys
diff --git a/pygama/raw/data_streamer.py b/pygama/raw/data_streamer.py
index 41ed6926f..935cd12c0 100644
--- a/pygama/raw/data_streamer.py
+++ b/pygama/raw/data_streamer.py
@@ -1,8 +1,10 @@
 from abc import ABC
-from .raw_buffer import RawBuffer, RawBufferList, RawBufferLibrary
+
+from .raw_buffer import RawBuffer, RawBufferLibrary, RawBufferList
+
 
 class DataStreamer(ABC):
-    """ Base clase for data streams 
+    """ Base clase for data streams
 
     Provides a uniform interface for streaming, e.g.
 
@@ -93,7 +95,7 @@ def open_stream(self, stream_name, rb_lib=None, buffer_size=8192,
 
             # dec_name is in rb_lib: store the name, and initialize its buffer lgdos
             dec_names.append(dec_name)
-            for rb in rb_lib[dec_name]: 
+            for rb in rb_lib[dec_name]:
                 # use the first available key
                 key = rb.key_list[0] if len(rb.key_list) > 0 else None
                 rb.lgdo = decoder.make_lgdo(key=key, size=buffer_size)
@@ -173,7 +175,7 @@ def read_chunk(self, chunk_mode_override=None, rp_max=1000000, clear_full_buffer
         if clear_full_buffers: self.rb_lib.clear_full()
         self.any_full = False
 
-        chunk_mode = self.chunk_mode if chunk_mode_override is None else chunk_mode_override 
+        chunk_mode = self.chunk_mode if chunk_mode_override is None else chunk_mode_override
         if verbosity>1: print(f'reading chunk with chunk_mode {chunk_mode}')
 
         read_one_packet = (chunk_mode == 'single_packet')
@@ -199,7 +201,7 @@ def read_chunk(self, chunk_mode_override=None, rp_max=1000000, clear_full_buffer
                 elif rb.is_full(): list_of_rbs.append(rb)
         return list_of_rbs
 
-        
+
 
 
     def get_decoder_list(self):
@@ -230,4 +232,3 @@ def build_default_rb_lib(self, out_stream=''):
             rb_lib[dec_name] = RawBufferList()
             rb_lib[dec_name].append(rb)
         return rb_lib
-
diff --git a/pygama/raw/fc/fc_config_decoder.py b/pygama/raw/fc/fc_config_decoder.py
index ebef8b834..71551f93f 100644
--- a/pygama/raw/fc/fc_config_decoder.py
+++ b/pygama/raw/fc/fc_config_decoder.py
@@ -1,6 +1,8 @@
-from ..data_decoder import *
 from pygama import lgdo
 
+from ..data_decoder import *
+
+
 class FCConfigDecoder(DataDecoder):
     """
     Decode FlashCam config data
diff --git a/pygama/raw/fc/fc_event_decoder.py b/pygama/raw/fc/fc_event_decoder.py
index 9186a9758..08f34ab42 100644
--- a/pygama/raw/fc/fc_event_decoder.py
+++ b/pygama/raw/fc/fc_event_decoder.py
@@ -1,5 +1,5 @@
-from pygama.raw.data_decoder import *
 from pygama import lgdo
+from pygama.raw.data_decoder import *
 
 
 class FCEventDecoder(DataDecoder):
@@ -174,7 +174,7 @@ def decode_packet(self, fcio, evt_rbkd, packet_id, verbosity=0):
         n_bytes : int
             (estimated) number of bytes in the packet that was just decoded.
         """
-        if fcio.numtraces > self.max_numtraces: 
+        if fcio.numtraces > self.max_numtraces:
             self.max_numtraces = fcio.numtraces
             # The buffer might be storing all channels' data, so set the
             # fill_safety to the max number of traces we've seen so far.
diff --git a/pygama/raw/fc/fc_status_decoder.py b/pygama/raw/fc/fc_status_decoder.py
index 8e60b363e..1dc835e29 100644
--- a/pygama/raw/fc/fc_status_decoder.py
+++ b/pygama/raw/fc/fc_status_decoder.py
@@ -1,7 +1,10 @@
 import os
+
 import numpy as np
+
 from ..data_decoder import *
 
+
 class FCStatusDecoder(DataDecoder):
     """
     decode FlashCam digitizer status data.
@@ -10,7 +13,7 @@ def __init__(self, *args, **kwargs):
         """
         """
         self.decoded_values = {
-            'status': { # 0: Errors occured, 1: no errors
+            'status': { # 0: Errors occurred, 1: no errors
               'dtype': 'int32',
             },
             'statustime': { # fc250 seconds, microseconds, dummy, startsec startusec
@@ -73,7 +76,7 @@ def decode_packet(self, fcio, status_rb, packet_id, verbosity=0):
         tbl = status_rb.lgdo
         ii = status_rb.loc
 
-        # status -- 0: Errors occured, 1: no errors
+        # status -- 0: Errors occurred, 1: no errors
         tbl['status'].nda[ii] = fcio.status
 
         # times
@@ -100,4 +103,3 @@ def decode_packet(self, fcio, status_rb, packet_id, verbosity=0):
         status_rb.loc += 1
 
         return len(tbl) - status_rb.loc > 1
-
diff --git a/pygama/raw/fc/fc_streamer.py b/pygama/raw/fc/fc_streamer.py
index ee6ce6f3f..b9609fa92 100644
--- a/pygama/raw/fc/fc_streamer.py
+++ b/pygama/raw/fc/fc_streamer.py
@@ -1,12 +1,16 @@
 import os
+
+import fcutils
 import numpy as np
+
 from pygama import lgdo
-from ..raw_buffer import *
+
 from ..data_streamer import DataStreamer
+from ..raw_buffer import *
 from .fc_config_decoder import FCConfigDecoder
 from .fc_event_decoder import FCEventDecoder
 from .fc_status_decoder import FCStatusDecoder
-import fcutils
+
 
 class FCStreamer(DataStreamer):
     """
@@ -39,14 +43,14 @@ def open_stream(self, fcio_filename, rb_lib=None, buffer_size=8192,
         Parameters
         ----------
         fcio_filename : str
-            the FCIO filename 
+            the FCIO filename
         rb_lib : RawBufferLibrary
             library of buffers for this stream
         buffer_size : int
             length of tables to be read out in read_chunk
         verbosity : int
             verbosity level for the initialize function
- 
+
         Returns
         -------
         header_data : list(RawBuffer)
@@ -78,7 +82,7 @@ def open_stream(self, fcio_filename, rb_lib=None, buffer_size=8192,
         self.packet_id = 0 # for storing packet order in output tables
         #self.max_numtraces = 0 # for checking that the tables are large enough
 
-        if 'FCConfigDecoder' in rb_lib: 
+        if 'FCConfigDecoder' in rb_lib:
             config_rb_list = rb_lib['FCConfigDecoder']
             if len(config_rb_list) != 1:
                 print(f'warning! config_rb_list had length {len(config_rb_list)}, ignoring all but the first')
@@ -119,7 +123,7 @@ def read_packet(self, verbosity=0):
         elif rc == 5: # recevent
             print(f'warning: got a RecEvent packet -- skipping?')
             print(f'         n_bytes_read = {self.n_bytes_read}')
-            # sizeof(fcio_recevent): (6 + 3*10 + 1*2304 + 3*4000)*4 
+            # sizeof(fcio_recevent): (6 + 3*10 + 1*2304 + 3*4000)*4
             self.n_bytes_read += 57360
             return True
 
@@ -144,4 +148,3 @@ def read_packet(self, verbosity=0):
             # sizeof(fcio_event): (5 + 3*10 + 1)*4 + numtraces*(1 + nsamples+2)*2
             self.n_bytes_read += 144 + self.fcio.numtraces*(self.fcio.nsamples + 3)*2
             return True
-
diff --git a/pygama/raw/orca/orca_digitizers.py b/pygama/raw/orca/orca_digitizers.py
index a15888138..eb130c5e6 100644
--- a/pygama/raw/orca/orca_digitizers.py
+++ b/pygama/raw/orca/orca_digitizers.py
@@ -1,9 +1,12 @@
 import sys
+
 import numpy as np
 
-from .orcadaq import OrcaDecoder, get_ccc
 from pygama.lh5 import Table
 
+from .orcadaq import OrcaDecoder, get_ccc
+
+
 class ORCAStruck3302(OrcaDecoder):
     """
     ORCA decoder for Struck 3302 digitizer data
@@ -43,7 +46,7 @@ def __init__(self, *args, **kwargs):
             'waveform': {
               'dtype': 'uint16',
               'datatype': 'waveform',
-              'length': 65532, # max value. override this before initalizing buffers to save RAM
+              'length': 65532, # max value. override this before initializing buffers to save RAM
               'sample_period': 10, # override if a different clock rate is used
               'sample_period_units': 'ns',
               'units': 'adc',
@@ -488,7 +491,7 @@ def __init__(self, *args, **kwargs):
             'waveform': {
               'dtype': 'uint16',
               'datatype': 'waveform',
-              'length': 65532, # max value. override this before initalizing buffers to save RAM
+              'length': 65532, # max value. override this before initializing buffers to save RAM
               'sample_period': 8, # override if a different clock rate is used
               'sample_period_units': 'ns',
               'units': 'adc',
diff --git a/pygama/raw/orca/orca_flashcam.py b/pygama/raw/orca/orca_flashcam.py
index fee9e1750..452dd83b3 100644
--- a/pygama/raw/orca/orca_flashcam.py
+++ b/pygama/raw/orca/orca_flashcam.py
@@ -1,17 +1,18 @@
 import numpy as np
 
-from .orcadaq import OrcaDecoder, get_ccc, get_readout_info, get_auxhw_info
 from .fcdaq import FlashCamEventDecoder
+from .orcadaq import OrcaDecoder, get_auxhw_info, get_ccc, get_readout_info
+
 
 class ORCAFlashCamListenerConfigDecoder(OrcaDecoder):
     '''
     Decoder for FlashCam listener config written by ORCA
     '''
     def __init__(self, *args, **kwargs):
-        
+
         self.decoder_name    = 'ORFlashCamListenerConfigDecoder'
         self.orca_class_name = 'ORFlashCamListenerModel'
-        
+
         # up through ch_inputnum, these are in order of the fcio data format
         # for similicity.  append any additional values after this.
         self.decoded_values = {
@@ -37,7 +38,7 @@ def __init__(self, *args, **kwargs):
                               'array_of_equalsized_arrays<1,1>{real}',
                               'length': 2400, },
             }
-        
+
         super().__init__(args, kwargs)
 
 
@@ -48,17 +49,17 @@ def get_decoded_values(self, channel=None):
     def max_n_rows_per_packet(self):
         return 1
 
-    
+
     def decode_packet(self, packet, lh5_tables,
                       packet_id, header_dict, verbose=False):
-        
+
         data = np.frombuffer(packet, dtype=np.int32)
         tbl  = lh5_tables
         ii   = tbl.loc
 
         tbl['readout_id'].nda[ii]  = (data[0] & 0xffff0000) >> 16
         tbl['listener_id'].nda[ii] =  data[0] & 0x0000ffff
-        
+
         for i,k in enumerate(self.decoded_values):
             if i < 2: continue
             tbl[k].nda[ii] = data[i-1]
@@ -72,7 +73,7 @@ def decode_packet(self, packet, lh5_tables,
 
         tbl.push_row()
 
-        
+
 class ORCAFlashCamListenerStatusDecoder(OrcaDecoder):
     '''
     Decoder for FlashCam status packets written by ORCA
@@ -141,7 +142,7 @@ def __init__(self, *args, **kwargs):
             'card_link_errors': {
                 'dtype':        'uint32',
                 'datatype':     'array<1>{array<1>{real}}',
-                'length_guess':  self.nCards, },        
+                'length_guess':  self.nCards, },
             'card_other_errors': {
                 'dtype':        'uint32',
                 'datatype':     'array<1>{array<1>{real}}',
@@ -184,7 +185,7 @@ def __init__(self, *args, **kwargs):
         # arrays to temporarily store card-level decoded data
         self.cdata = {}
         self.resize_card_data(ncards=self.nCards)
-            
+
         super().__init__(args, kwargs)
 
 
@@ -222,18 +223,18 @@ def resize_card_data(self, ncards):
         # set nCards to allow for not calling this function during decoding
         self.nCards = ncards
 
-        
+
     def get_decoded_values(self, channel=None):
         return self.decoded_values
 
-    
+
     def max_n_rows_per_packet(self):
         return 1
 
-    
+
     def decode_packet(self, packet, lh5_tables,
                       packet_id, header_dict, verbose=False):
-        
+
         data = np.frombuffer(packet, dtype=np.uint32)
         tbl  = lh5_tables
         ii   = tbl.loc
@@ -287,16 +288,16 @@ def decode_packet(self, packet, lh5_tables,
         # populate the card level data with the flattened local data, then push
         for key in self.cdata:
             tbl[key].set_vector(ii, self.cdata[key].flatten())
-            
+
         tbl.push_row()
-        
+
 
 class ORCAFlashCamADCWaveformDecoder(OrcaDecoder):
     """
     Decoder for FlashCam ADC data written by ORCA
     """
     def __init__(self, *args, **kwargs):
-        
+
         self.decoder_name    = 'ORFlashCamADCWaveformDecoder'
         self.orca_class_name = 'ORFlashCamADCModel'
 
@@ -308,7 +309,7 @@ def __init__(self, *args, **kwargs):
             'fcio_id':  { 'dtype': 'uint16', }  }
         fc = FlashCamEventDecoder()
         self.decoded_values_template.update(fc.decoded_values)
-        
+
         super().__init__(*args, **kwargs)
         self.decoded_values = {}
         self.skipped_channels = {}
@@ -332,13 +333,13 @@ def max_n_rows_per_packet(self):
 
     def set_object_info(self, object_info):
         self.object_info = object_info
-        
+
         # get the readout list for looking up the waveform length.
         # catch AttributeError for when the header_dict is not yet set.
         roi = []
         try: roi=get_readout_info(self.header_dict, 'ORFlashCamListenerModel')
         except AttributeError: pass
-        
+
         for card_dict in self.object_info:
             crate   = card_dict['Crate']
             card    = card_dict['Card']
@@ -369,7 +370,7 @@ def set_object_info(self, object_info):
                 if samples > 0:
                     self.decoded_values[ccc]['waveform']['length'] = samples
 
-                    
+
     def decode_packet(self, packet, lh5_tables,
                       packet_id, header_dict, verbose=False):
 
diff --git a/pygama/raw/orca/orca_pollers.py b/pygama/raw/orca/orca_pollers.py
index d2d68f94f..960930a30 100644
--- a/pygama/raw/orca/orca_pollers.py
+++ b/pygama/raw/orca/orca_pollers.py
@@ -1,9 +1,11 @@
+import sys
+
 import numpy as np
 import pandas as pd
-import sys
 
 from .io_base import DataTaker
 
+
 class MJDPreampDecoder(DataTaker):
     """
     members:
@@ -128,8 +130,8 @@ def __init__(self, *args, **kwargs):
 
     def decode_event(self, event_data_bytes, event_number, header_dict,
                      verbose=False):
-        """ 
-        see README for the 32-bit data word diagram 
+        """
+        see README for the 32-bit data word diagram
         """
         event_data_int = np.fromstring(event_data_bytes, dtype=np.uint32)
         event_data_float = np.fromstring(event_data_bytes, dtype=np.float32)
diff --git a/pygama/raw/orca/stream_orca.py b/pygama/raw/orca/stream_orca.py
index ca2ea47b4..d73abbe97 100644
--- a/pygama/raw/orca/stream_orca.py
+++ b/pygama/raw/orca/stream_orca.py
@@ -1,12 +1,15 @@
-import sys, gzip
-import numpy as np
+import gzip
 import plistlib
+import sys
 
+import numpy as np
 from tqdm.std import tqdm
-from ..utils import tqdm_range, update_progress
-from .io_base import DataDecoder
+
 from pygama import lh5
+
+from ..utils import tqdm_range, update_progress
 from .ch_group import *
+from .io_base import DataDecoder
 
 
 class OrcaDecoder(DataDecoder):
@@ -51,32 +54,32 @@ def parse_header(orca_filename):
     """
     with open_orca(orca_filename) as xmlfile_handle:
         #read the first word:
-        ba = bytearray(xmlfile_handle.read(8))
+        barr = bytearray(xmlfile_handle.read(8))
 
         #Replacing this to be python2 friendly
         # #first 4 bytes: header length in long words
-        # i = int.from_bytes(ba[:4], byteorder=sys.byteorder)
+        # i = int.from_bytes(barr[:4], byteorder=sys.byteorder)
         # #second 4 bytes: header length in bytes
-        # j = int.from_bytes(ba[4:], byteorder=sys.byteorder)
+        # j = int.from_bytes(barr[4:], byteorder=sys.byteorder)
 
         big_endian = False if sys.byteorder == "little" else True
-        i = from_bytes(ba[:4], big_endian=big_endian)
-        j = from_bytes(ba[4:], big_endian=big_endian)
+        i = from_bytes(barr[:4], big_endian=big_endian)
+        j = from_bytes(barr[4:], big_endian=big_endian)
         if (np.ceil(j/4) != i-2) and (np.ceil(j/4) != i-3):
             print('Error: header byte length = %d is the wrong size to fit into %d header packet words' % (j, i-2))
             return i, j, {}
 
         #read in the next that-many bytes that occupy the plist header
         as_bytes = xmlfile_handle.read(j)
-        ba = bytearray(as_bytes)
+        barr = bytearray(as_bytes)
 
         #convert to string
         #the readPlistFromBytes method doesn't exist in 2.7
         if sys.version_info[0] < 3:
-            header_string = ba.decode("utf-8")
+            header_string = barr.decode("utf-8")
             header_dict = plistlib.readPlistFromString(header_string)
         elif sys.version_info[1] < 9:
-            header_dict = plistlib.readPlistFromBytes(ba)
+            header_dict = plistlib.readPlistFromBytes(barr)
         else:
             header_dict = plistlib.loads(as_bytes, fmt=plistlib.FMT_XML)
         return i, j, header_dict
@@ -180,7 +183,7 @@ def get_object_info(header_dict, orca_class_name):
 
 def get_readout_info(header_dict, orca_class_name, unique_id=-1):
     """
-    retunrs a list with all the readout list info from the header with name
+    returns a list with all the readout list info from the header with name
     orca_class_name.  optionally, if unique_id >= 0 only return the list for
     that Orca unique id number.
     """
@@ -252,8 +255,8 @@ def get_next_packet(f_in):
     # reserved        = (head[4] + (head[5]<<8))
 
     # Using an array of uint32
-    record_length = int((head[0] & 0x3FFFF))
-    data_id = int((head[0] >> 18))
+    record_length = int(head[0] & 0x3FFFF)
+    data_id = int(head[0] >> 18)
     # reserved =int( (head[1] &0xFFFF))
 
     # /* ========== read in the rest of the event data ========== */
@@ -288,6 +291,7 @@ def get_channel(ccc):
 # Do it here so that orca_digitizers can import the functions above here
 from . import orca_digitizers, orca_flashcam
 
+
 def process_orca(daq_filename, raw_file_pattern, n_max=np.inf, ch_groups_dict=None, verbose=False, buffer_size=1024):
     """
     convert ORCA DAQ data to "raw" lh5
@@ -310,7 +314,7 @@ def process_orca(daq_filename, raw_file_pattern, n_max=np.inf, ch_groups_dict=No
     file_size = float(f_in.tell())
     f_in.seek(0, 0)  # rewind
     file_size_MB = file_size / 1e6
-    print("Total file size: {:.3f} MB".format(file_size_MB))
+    print(f"Total file size: {file_size_MB:.3f} MB")
     print("Run number:", get_run_number(header_dict))
 
 
@@ -469,10 +473,9 @@ def process_orca(daq_filename, raw_file_pattern, n_max=np.inf, ch_groups_dict=No
         print("WARNING, Found the following unknown data IDs:")
         for data_id in unrecognized_data_ids:
             try:
-                print("  {}: {}".format(data_id, id2dn_dict[data_id]))
+                print(f"  {data_id}: {id2dn_dict[data_id]}")
             except KeyError:
-                print("  {}: Unknown".format(data_id))
+                print(f"  {data_id}: Unknown")
         print("hopefully they weren't important!\n")
 
-    print("Wrote RAW File:\n    {}\nFILE INFO:".format(raw_file_pattern))
-
+    print(f"Wrote RAW File:\n    {raw_file_pattern}\nFILE INFO:")
diff --git a/pygama/raw/raw_buffer.py b/pygama/raw/raw_buffer.py
index 5c07d2a72..f312083b3 100644
--- a/pygama/raw/raw_buffer.py
+++ b/pygama/raw/raw_buffer.py
@@ -60,8 +60,10 @@
 """
 
 import os
+
 from pygama import lgdo
 
+
 class RawBuffer:
     '''
     A RawBuffer is in essence a an lgdo object (typically a Table) to which
@@ -71,7 +73,7 @@ class RawBuffer:
 
     Attributes
     ----------
-    lgdo : lgdo 
+    lgdo : lgdo
         the lgdo used as the actual buffer. Typically a table. Set to None upon
         creation so that the user or a decoder can initialize it later.
     key_list : list
@@ -86,7 +88,7 @@ class RawBuffer:
         File example: '/path/filename.lh5:/group'
         Socket example: '198.0.0.100:8000'
     out_name : str (optional)
-        the name / identifier of the object in the ouput stream 
+        the name / identifier of the object in the output stream
     '''
 
 
@@ -120,7 +122,7 @@ def __repr__(self): return str(self)
 class RawBufferList(list):
     '''
     A RawBufferList holds a collection of RawBuffers of identical structure
-    (same format lgdo's with the same fields). 
+    (same format lgdo's with the same fields).
     '''
 
 
@@ -146,11 +148,11 @@ def set_from_json_dict(self, json_dict, kw_dict={}):
         expand_rblist_json_dict(json_dict, kw_dict)
         for name in json_dict:
             rb = RawBuffer()
-            if 'key_list' in json_dict[name]: 
+            if 'key_list' in json_dict[name]:
                 rb.key_list = json_dict[name]['key_list']
-            if 'out_stream' in json_dict[name]: 
+            if 'out_stream' in json_dict[name]:
                 rb.out_stream = json_dict[name]['out_stream']
-            if 'out_name' in json_dict[name]: 
+            if 'out_name' in json_dict[name]:
                 rb.out_name = json_dict[name]['out_name']
             else: rb.out_name = name
             self.append(rb);
@@ -182,7 +184,7 @@ def get_list_of(self, attribute):
         return values
 
     def clear_full(self):
-        for rb in self: 
+        for rb in self:
             if rb.is_full(): rb.loc = 0
 
 
@@ -192,7 +194,7 @@ class RawBufferLibrary(dict):
     names of decoders that can write to them
     '''
     def __init__(self, json_dict=None, kw_dict={}):
-        if json_dict is not None: 
+        if json_dict is not None:
             self.set_from_json_dict(json_dict, kw_dict)
 
 
@@ -210,11 +212,11 @@ def set_from_json_dict(self, json_dict, kw_dict={}):
         }
 
         By default "name" is used for the RawBuffer's "out_name" attribute, but
-        this can be overrided if desired by providing an explicit "out_name"
+        this can be overridden if desired by providing an explicit "out_name"
 
-        Allowed shorthands, in order of exapansion:
+        Allowed shorthands, in order of expansion:
         * key_list may have entries that are 2-integer lists corresponding to
-          the first and last integer keys in a continguous range (e.g. of
+          the first and last integer keys in a contiguous range (e.g. of
           channels) that get stored to the same buffer. These simply get
           replaced with the explicit list of integers in the range. We use lists
           not tuples for json compliance.
@@ -241,7 +243,7 @@ def set_from_json_dict(self, json_dict, kw_dict={}):
         kw_dict : dict
             dict of keyword-value pairs for substitutions into the out_stream
             and out_name fields
-        '''                
+        '''
         for list_name in json_dict:
             if list_name not in self: self[list_name] = RawBufferList()
             self[list_name].set_from_json_dict(json_dict[list_name], kw_dict)
@@ -266,7 +268,7 @@ def get_list_of(self, attribute, unique=True):
         output_file_list = rbl.get_list_of('out_stream')
         """
         values = []
-        for rb_list in self.values(): 
+        for rb_list in self.values():
             values += rb_list.get_list_of(attribute)
         if unique: values = list(set(values))
         return values
@@ -294,7 +296,7 @@ def expand_rblist_json_dict(json_dict, kw_dict):
             return
         info = json_dict[name] # changes to info will change json_dict[name]
         # make sure we have a key list
-        if 'key_list' not in info: 
+        if 'key_list' not in info:
             print(f'expand_json_dict: {name} is missing key_list')
             continue
 
@@ -309,18 +311,18 @@ def expand_rblist_json_dict(json_dict, kw_dict):
                 info['key_list'][i:i+1] = range(key_range[0], key_range[1]+1)
                 i += key_range[1]-key_range[0]
             i += 1
-        
+
         # Expand list_names if name contains a key-based formatter
-        if '{key' in name: 
+        if '{key' in name:
             for key in info['key_list']:
                 expanded_name = name.format(key=key)
                 json_dict[expanded_name] = info.copy()
                 json_dict[expanded_name]['key_list'] = [key];
             json_dict.pop(name)
 
-    # now re-iterate and exand out_paths
+    # now re-iterate and expand out_paths
     for name, info in json_dict.items():
-        if len(info['key_list']) == 1 and info['key_list'][0] != "*": 
+        if len(info['key_list']) == 1 and info['key_list'][0] != "*":
             kw_dict['key'] = info['key_list'][0]
         if 'out_stream' in info:
             if name != '*' and '{name' in info['out_stream']: kw_dict['name'] = name
@@ -344,10 +346,10 @@ def write_to_lh5_and_clear(raw_buffers, lh5_store=None, wo_mode='append', verbos
     for rb in raw_buffers:
         if rb.lgdo is None or rb.loc == 0: continue # no data to write
         ii = rb.out_stream.find(':')
-        if ii == -1: 
+        if ii == -1:
             filename = rb.out_stream
             group = '/'
-        else: 
+        else:
             filename = rb.out_stream[:ii]
             group = rb.out_stream[ii+1:]
             if len(group) == 0: group = '/' # in case out_stream ends with :
@@ -356,5 +358,3 @@ def write_to_lh5_and_clear(raw_buffers, lh5_store=None, wo_mode='append', verbos
                                n_rows=rb.loc, wo_mode=wo_mode, verbosity=verbosity)
         # and clear
         rb.loc = 0
-
-
diff --git a/pygama/vis/mpl/clint.mpl b/pygama/vis/mpl/clint.mpl
index febd0189e..3116dfb10 100644
--- a/pygama/vis/mpl/clint.mpl
+++ b/pygama/vis/mpl/clint.mpl
@@ -45,4 +45,3 @@ ytick.minor.visible  : True
 font.family : STIXGeneral
 font.serif : STIX
 mathtext.fontset : stix
-
diff --git a/pygama/vis/mpl/root.mpl b/pygama/vis/mpl/root.mpl
index d0fdccc24..157c5de2d 100644
--- a/pygama/vis/mpl/root.mpl
+++ b/pygama/vis/mpl/root.mpl
@@ -33,4 +33,3 @@ ytick.minor.visible  : True
 # font.family : Helvetica Neue Light
 # font.style : light
 #font.family : Optima
-
diff --git a/pygama/vis/waveform_browser.py b/pygama/vis/waveform_browser.py
index ad3e23454..a2e324528 100644
--- a/pygama/vis/waveform_browser.py
+++ b/pygama/vis/waveform_browser.py
@@ -1,29 +1,30 @@
-import sys
-import os
 import glob
 import itertools
-import string
 import math
+import os
+import string
+import sys
+
+import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
+from cycler import cycler
+from matplotlib.lines import Line2D
 
-from pygama.dsp.processing_chain import build_processing_chain
 import pygama.lgdo.lh5_store as lh5
+from pygama.dsp.processing_chain import build_processing_chain
 from pygama.math.units import unit_registry as ureg
 
-import matplotlib.pyplot as plt
-from matplotlib.lines import Line2D
-from cycler import cycler
 
 class WaveformBrowser:
     """
     The Waveform Browser is a tool meant for interacting with waveforms from
     LH5 files. This defines an interface for drawing waveforms from a file,
     drawing transformed waveforms defined using build_dsp style json files,
-    drawing horizontal and verticle lines at the values of calculated
+    drawing horizontal and vertical lines at the values of calculated
     parameters, and filling a legend with calculated parameters.
-    
-    
+
+
     """
 
     def __init__(self, files_in, lh5_group, base_path = '',
@@ -102,12 +103,12 @@ def __init__(self, files_in, lh5_group, base_path = '',
         block_width (default 16): block width for processing chain
         """
         self.verbosity = verbosity
-        
+
         self.norm_par = norm
         self.align_par = align
         self.n_drawn = n_drawn
         self.next_entry = 0
-        
+
         # data i/o initialization
         self.lh5_it = lh5.LH5Iterator(files_in,
                                       lh5_group,
@@ -115,7 +116,7 @@ def __init__(self, files_in, lh5_group, base_path = '',
                                       entry_list=entry_list,
                                       entry_mask=entry_mask,
                                       buffer_len=buffer_len )
-        
+
 
         # Get the input buffer and read the first chunk
         self.lh5_in, _ = self.lh5_it.read(0)
@@ -128,12 +129,12 @@ def __init__(self, files_in, lh5_group, base_path = '',
                 entry_offset = self.lh5_it.file_map[i-1] if i>0 else 0
                 entries += [ entry_offset + entry for entry in f_entries ]
             self.aux_vals = self.aux_vals.iloc[entries].reset_index()
-        
+
         # initialize objects to draw: dict from name to list of 2DLines
         if isinstance(lines, str): self.lines = { lines:[] }
         elif lines is None: self.lines = {}
         else: self.lines = { l:[] for l in lines }
-        
+
         # styles
         if isinstance(styles, (list, tuple)):
             self.styles = [ None for _ in self.lines ]
@@ -161,10 +162,10 @@ def __init__(self, files_in, lh5_group, base_path = '',
 
         self.legend_format = [] # list of formatter strings
         self.legend_vals = {}  # Set up dict from names to lists of values
-        
+
         if legend is None: legend = []
         elif isinstance(legend, str): legend = [legend]
-        
+
         for entry in legend:
             legend_format = ""
             for st, name, form, cv in string.Formatter().parse(entry):
@@ -172,7 +173,7 @@ def __init__(self, files_in, lh5_group, base_path = '',
                     # If name is none, this is the last batch of characters
                     legend_format += st
                     break
-                
+
                 if name=='':
                     raise KeyError("Cannot use empty formatter in "+entry)
                 self.legend_vals[name] = []
@@ -180,35 +181,35 @@ def __init__(self, files_in, lh5_group, base_path = '',
                 if form is None or form=='':
                     form = '~0.3P'
                 cv = '' if cv is None or cv=='' else "!"+cv
-                legend_format += "{}{{{}:{}{}}}".format(st, name, form, cv)
+                legend_format += f"{st}{{{name}:{form}{cv}}}"
             self.legend_format.append(legend_format)
 
         self.legend_kwargs = legend_opts if isinstance(legend_opts, dict) else {}
-        
+
         # make processing chain and output buffer
         outputs = list(self.lines) + list(self.legend_vals)
         if isinstance(self.norm_par, str): outputs += [self.norm_par]
         if isinstance(self.align_par, str): outputs += [self.align_par]
 
-            
+
         # Remove any values not found in aux_vals
         if self.aux_vals is not None:
             outputs = [ o for o in outputs if o not in self.aux_vals ]
-        
+
         self.proc_chain, self.lh5_it.field_mask, self.lh5_out = build_processing_chain(self.lh5_in, dsp_config, db_dict=database, outputs=outputs, verbosity=self.verbosity, block_width=block_width)
         self.proc_chain.execute()
-        
+
         # Check if all of our outputs can be found
         for name in outputs:
             if not name in self.lh5_out:
                 raise KeyError("Could not find "+name+" in input lh5 file, DSP config file, or aux values")
-        
+
         self.x_unit = ureg[x_unit] if x_unit else None
         self.x_lim = x_lim
         self.y_lim = y_lim
         self.auto_x_lim = [np.inf, -np.inf]
         self.auto_y_lim = [np.inf, -np.inf]
-        
+
         # Set limit and convert to x-unit if needed; also set x_unit if needed
         if self.x_lim is not None:
             self.x_lim = list(self.x_lim)
@@ -221,16 +222,16 @@ def __init__(self, files_in, lh5_group, base_path = '',
                     else:
                         self.x_unit = self.x_lim[i].u
                         self.x_lim[i] = self.x_lim[i].m
-        
+
         # If we still have no x_unit get it from the first waveform we can find
         if self.x_unit is None:
             for wf in self.lh5_out.values():
                 if not isinstance(wf, lh5.WaveformTable): continue
                 self.x_unit = ureg[wf.dt_units]
-        
+
         self.fig = None
-        self.ax = None        
-    
+        self.ax = None
+
     def new_figure(self):
         """Create a new figure and draw in it"""
         self.fig, self.ax = plt.subplots(1)
@@ -265,12 +266,12 @@ def clear_data(self):
         self.auto_x_lim = [np.inf, -np.inf]
         self.auto_y_lim = [np.inf, -np.inf]
         self.n_stored = 0
-        
+
     def find_entry(self, entry, append=True, safe=False):
         """
         Find the requested data associated with entry in input files and
         place store it internally without drawing it.
-        
+
         Parameters
         ----------
         entry : int or [ints]
@@ -295,7 +296,7 @@ def find_entry(self, entry, append=True, safe=False):
             self.lh5_it.read(entry)
             self.proc_chain.execute()
             i_tb = 0
-        
+
         # get scaling factor/time shift if used
         if self.norm_par is None:
             norm = 1.
@@ -317,7 +318,7 @@ def find_entry(self, entry, append=True, safe=False):
                 ref_time = self.aux_vals[self.align_par][entry]
             else:
                 raise
-            
+
         leg_handle = None
 
         # lines
@@ -332,20 +333,20 @@ def find_entry(self, entry, append=True, safe=False):
                 x = np.linspace(t0, t0+dt*(data.wf_len-1), data.wf_len)
                 lines.append(Line2D(x, y))
                 self._update_auto_limit(x, y)
-                
+
             elif isinstance(data, lh5.Array):
                 val = data.nda[i_tb]
                 unit = data.attrs.get('units', None)
                 if unit and unit in ureg and ureg.is_compatible_with(unit, self.x_unit):
                     # Vertical line
-                    val = val*float(ureg[unit]/self.x_unit) - ref_time 
+                    val = val*float(ureg[unit]/self.x_unit) - ref_time
                     lines.append(Line2D([val]*2, [-lim, lim]))
                     self._update_auto_limit(val, None)
                 else:
                     # Horizontal line
                     lines.append(Line2D([-lim, lim], [val/norm]*2))
                     self._update_auto_limit(None, val)
-                    
+
             elif data is None:
                 # Check for data in auxiliary table. It's unitless so I guess just do an hline...
                 val = self.aux_vals[name][entry]/norm
@@ -369,12 +370,12 @@ def find_entry(self, entry, append=True, safe=False):
                     data = ureg.Quantity(data.nda[i_tb])
             else:
                 raise TypeError("WaveformBrowser does not adding legend entries for data of type " + data.__class__)
-            
+
             vals.append(data)
 
         self.n_stored += 1
         self.next_entry = entry + 1
-    
+
     def draw_current(self, clear=True):
         """
         Draw the waveforms and data currently held internally by this class.
@@ -382,7 +383,7 @@ def draw_current(self, clear=True):
         # Make figure/axis if needed
         if not (self.ax and self.fig and plt.fignum_exists(self.fig.number)):
             self.new_figure()
-        
+
         if clear:
             self.ax.clear()
 
@@ -398,20 +399,20 @@ def draw_current(self, clear=True):
         leg_labels = []
         if not isinstance(self.styles, list):
             styles = self.styles
-            
+
         # draw lines
         for i, lines in enumerate(self.lines.values()):
             if isinstance(self.styles, list):
                 styles = self.styles[i]
             if styles is None:
                 styles = cycler(plt.rcparams)
-            
+
             for line, sty in zip(lines, styles):
                 if sty is not None: line.update(sty)
                 if line.get_figure() is not None: line.remove()
                 line.set_transform(self.ax.transData)
                 self.ax.add_line(line)
-                
+
                 leg_handles.append(line)
 
         # Get legend entries
@@ -448,16 +449,16 @@ def _update_auto_limit(self, x, y):
         if isinstance(x, np.ndarray) and self.x_lim is not None:
             x_where['where'] = ((x>=self.x_lim[0]) & (x<=self.x_lim[1]))
         if x is not None:
-            self.auto_x_lim[0] = np.amin(x, **y_where, initial=self.auto_x_lim[0]) 
+            self.auto_x_lim[0] = np.amin(x, **y_where, initial=self.auto_x_lim[0])
             self.auto_x_lim[1] = np.amax(x, **y_where, initial=self.auto_x_lim[1])
         if y is not None:
             self.auto_y_lim[0] = np.amin(y, **y_where, initial=self.auto_y_lim[0])
             self.auto_y_lim[1] = np.amax(y, **y_where, initial=self.auto_y_lim[1])
-                
+
     def draw_entry(self, entry, append=False, clear=True, safe=False):
         """
         Draw specified entry in the current figure/axes
-        
+
         Parameters
         ----------
         entry : int or [ints]
@@ -478,19 +479,18 @@ def find_next(self, n_wfs = None, append = False):
         entries = (self.next_entry, self.next_entry+n_wfs)
         self.find_entry(range(*entries), append, safe=True)
         return entries
-        
+
     def draw_next(self, n_wfs = None, append = False, clear = True):
         """Draw the next n_wfs waveforms (default self.n_drawn). See draw_next"""
         entries = self.find_next(append)
         self.draw_current(clear)
         return entries
-            
+
     def reset(self):
         """ Reset to the start of the file for draw_next """
         self.clear_data()
         self.next_entry = 0
-        
+
     def __iter__(self):
         while self.next_entry < len(self.lh5_it):
             yield self.draw_next()
-            
diff --git a/setup.py b/setup.py
index e4b25b4ed..fc04d872e 100755
--- a/setup.py
+++ b/setup.py
@@ -1,19 +1,20 @@
- #!/usr/bin/env python3
+#!/usr/bin/env python3
 import os
+import platform
 import re
+import subprocess
 import sys
 import sysconfig
-import platform
-import subprocess
-
 from distutils.version import LooseVersion
-from setuptools import setup, Extension, find_packages
-from setuptools.command.egg_info import egg_info
+from shutil import copyfile, copymode
+
+from setuptools import Extension, find_packages, setup
 from setuptools.command.build_ext import build_ext
 from setuptools.command.build_py import build_py
 from setuptools.command.develop import develop
+from setuptools.command.egg_info import egg_info
 from setuptools.command.test import test as TestCommand
-from shutil import copyfile, copymode
+
 
 class CMakeExtension(Extension):
     def __init__(self, name, sourcedir=''):
@@ -90,7 +91,7 @@ def make_git_file():
         print(ex)
         print('continuing...')
 
-#Add a git hook to clean jupyter notebooks before commiting
+#Add a git hook to clean jupyter notebooks before committing
 def clean_jupyter_notebooks():
     import git
     repo = git.Repo(os.path.dirname(os.path.realpath(__file__)))
diff --git a/tutorials/DataProcessing.ipynb b/tutorials/DataProcessing.ipynb
index 3a8640f0c..9d47a0606 100644
--- a/tutorials/DataProcessing.ipynb
+++ b/tutorials/DataProcessing.ipynb
@@ -143,21 +143,9 @@
     "\n",
     "raw_df.hist('trapEmax', bins=2000)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "94f59e52-73e4-4de4-8999-d2f5cba1da1a",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 5
 }
diff --git a/tutorials/HADESDataTutorial.ipynb b/tutorials/HADESDataTutorial.ipynb
index cd1e5861c..1d00541d1 100644
--- a/tutorials/HADESDataTutorial.ipynb
+++ b/tutorials/HADESDataTutorial.ipynb
@@ -49,7 +49,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This gives access to all the prodenv commands a full desciption of which are found above."
+    "This gives access to all the prodenv commands a full description of which are found above."
    ]
   },
   {
@@ -779,11 +779,7 @@
    ]
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 4
 }
diff --git a/tutorials/IntroToDSP.ipynb b/tutorials/IntroToDSP.ipynb
index a6a3362ff..ae2eaa30c 100644
--- a/tutorials/IntroToDSP.ipynb
+++ b/tutorials/IntroToDSP.ipynb
@@ -8,7 +8,7 @@
     "\n",
     "Ian Guinn, UNC.  Presented at [LEGEND Software Tutorial, Nov. 2021](https://indico.legend-exp.org/event/561/)\n",
     "\n",
-    "This notebook will provide a brief introduction to a few of LEGEND's bread and butter digital signal processing filters. These processors will be applied in sequence to digitized waveforms that are recorded to disk, and a brief explanation of why we use them is provided. The parameters of each processor may be varied as well in order to see the affects that they have. This goes along with C. Wiseman and I. Guinn's slides shown at the Early Career Showcase: https://indico.legend-exp.org/event/699/. Reccomended reading is Knoll chapter 17, which introduces many pulse-shaping concepts (although, from a hardware-perspective more than a digital signal processing perspective; even so, the concepts are similar).\n",
+    "This notebook will provide a brief introduction to a few of LEGEND's bread and butter digital signal processing filters. These processors will be applied in sequence to digitized waveforms that are recorded to disk, and a brief explanation of why we use them is provided. The parameters of each processor may be varied as well in order to see the affects that they have. This goes along with C. Wiseman and I. Guinn's slides shown at the Early Career Showcase: https://indico.legend-exp.org/event/699/. Recommended reading is Knoll chapter 17, which introduces many pulse-shaping concepts (although, from a hardware-perspective more than a digital signal processing perspective; even so, the concepts are similar).\n",
     "\n",
     "### Set up the python environment\n",
     "It is recommended that you use a file from the LEGEND test-data repository: https://github.com/legend-exp/legend-testdata; you may need to change data_file if you want to run this notebook yourself."
@@ -261,7 +261,7 @@
    "metadata": {},
    "source": [
     "### Energy estimation: trap-max\n",
-    "The most common way to measure uncalibrated energy is to use the maximum of the trapezoidal filter. This works fairly well, assuming you have a nice flat top; however, due to our imperfect pole-zero correction, this can result in a degredation in our energy resolution"
+    "The most common way to measure uncalibrated energy is to use the maximum of the trapezoidal filter. This works fairly well, assuming you have a nice flat top; however, due to our imperfect pole-zero correction, this can result in a degradation in our energy resolution"
    ]
   },
   {
@@ -386,20 +386,9 @@
     "This notebook has hopefully provided an informative introduction to some of the signal processing filters used by LEGEND to perform our analysis. The filters shown often have tradeoffs between reducing the impact of noise and preserving the information in the original waveform; for this reason, we typically optimize the input parameters of each processor. In addition, the parameters measured here require calibration to be physically meaningful. The slides in https://indico.legend-exp.org/event/699/ will introduce how each of these processes are done.\n",
     "\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 4
 }
diff --git a/tutorials/LH5ParGridSearch.ipynb b/tutorials/LH5ParGridSearch.ipynb
index 395fdc7b0..d30410a49 100644
--- a/tutorials/LH5ParGridSearch.ipynb
+++ b/tutorials/LH5ParGridSearch.ipynb
@@ -94,7 +94,7 @@
    "outputs": [],
    "source": [
     "# load data from the raw files.\n",
-    "# we need cycle number and packet index information to retreive waveforms, so make a dataframe\n",
+    "# we need cycle number and packet index information to retrieve waveforms, so make a dataframe\n",
     "\n",
     "dg = DataGroup('./metadata/cage.json', load=True)\n",
     "dg.fileDB.query(query, inplace=True)\n",
@@ -1536,21 +1536,9 @@
    "source": [
     "# lol TBD  :-)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "08237f13-6e97-4060-a13e-13c91d5681cf",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 5
 }
diff --git a/tutorials/WaveformBrowserTutorial.ipynb b/tutorials/WaveformBrowserTutorial.ipynb
index 258ff3c12..3c308c150 100644
--- a/tutorials/WaveformBrowserTutorial.ipynb
+++ b/tutorials/WaveformBrowserTutorial.ipynb
@@ -206,7 +206,7 @@
     "aoe_accept = (AoE>0.045) & energy_selection\n",
     "\n",
     "AoE[aoe_accept].hist(bins=200, range=(-0, 0.1))\n",
-    "AoE[aoe_cut].hist(bins=200, range=(-0, 0.1))"
+    "AoE[aoe_cut].hist(bins=200, range=(-0, 0.1))\n",
     "# Use the lpgta dsp json file. TODO: get this from DataGroup\n",
     "# dsp_config_file = os.path.expandvars(\"$HOME/pygama/experiments/lpgta/LPGTA_dsp.json\")\n",
     "dsp_config_file = os.path.expandvars(\"./metadata/LPGTA_dsp.json\")\n"
@@ -260,7 +260,7 @@
    "source": [
     "## Example 5\n",
     "\n",
-    "The waveforms, lines and legend entries are all stored inside of the waveform browser. Sometimes you want to acess these directly; maybe you want to access the raw data, or do control the lines in a way not enabled by the WaveformBrowser interface. It is possible to access them quickly and easily. Waveforms and legend values are stored as a dict from the parameter name to a list of stored values.\n",
+    "The waveforms, lines and legend entries are all stored inside of the waveform browser. Sometimes you want to access these directly; maybe you want to access the raw data, or do control the lines in a way not enabled by the WaveformBrowser interface. It is possible to access them quickly and easily. Waveforms and legend values are stored as a dict from the parameter name to a list of stored values.\n",
     "- The waveforms are as a list of matplotlib Line2D artists\n",
     "- Horizontal and vertical lines are also stored as Line2D artists\n",
     "- Legend entries are stored as pint Quanitities\n",
@@ -300,20 +300,9 @@
     "    print(\"TrapEmax:\", en)\n",
     "    print()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 4
 }
diff --git a/tutorials/WriteProcessors.ipynb b/tutorials/WriteProcessors.ipynb
index 665e556dd..8b9b2dd6b 100644
--- a/tutorials/WriteProcessors.ipynb
+++ b/tutorials/WriteProcessors.ipynb
@@ -133,21 +133,9 @@
     "    for i, x in enumerate(w_in):\n",
     "        w_out[i] = x + a_in\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "b592af1c-8f04-4578-aa5f-35ba3a976188",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
- "metadata": {
-  "language_info": {
-   "name": "python"
-  }
- },
+ "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 5
 }
diff --git a/tutorials/metadata/cage.json b/tutorials/metadata/cage.json
old mode 100755
new mode 100644
diff --git a/tutorials/metadata/dsp_config.json b/tutorials/metadata/dsp_config.json
index fc752af59..509dab4a2 100644
--- a/tutorials/metadata/dsp_config.json
+++ b/tutorials/metadata/dsp_config.json
@@ -189,6 +189,6 @@
       "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 82149dba4..fdcbc39d1 100644
--- a/tutorials/metadata/optimizer_results.json
+++ b/tutorials/metadata/optimizer_results.json
@@ -1,1 +1,1 @@
-{"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"}}
\ No newline at end of file
+{"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"}}


Test patch:
diff --git a/tests/functional/test_l200_dataproc.py b/tests/functional/test_l200_dataproc.py
index 8b1378917..e69de29bb 100644
--- a/tests/functional/test_l200_dataproc.py
+++ b/tests/functional/test_l200_dataproc.py
@@ -1,1 +0,0 @@
-
diff --git a/tests/functional/test_teststand_dataproc.py b/tests/functional/test_teststand_dataproc.py
index 8b1378917..e69de29bb 100644
--- a/tests/functional/test_teststand_dataproc.py
+++ b/tests/functional/test_teststand_dataproc.py
@@ -1,1 +0,0 @@
-
diff --git a/tests/integration/test_integration.py b/tests/integration/test_integration.py
index 8b1378917..e69de29bb 100644
--- a/tests/integration/test_integration.py
+++ b/tests/integration/test_integration.py
@@ -1,1 +0,0 @@
-
diff --git a/tests/lgdo/test_lgdo_utils.py b/tests/lgdo/test_lgdo_utils.py
index ad8146c94..5e21e7d3d 100644
--- a/tests/lgdo/test_lgdo_utils.py
+++ b/tests/lgdo/test_lgdo_utils.py
@@ -1,6 +1,7 @@
-import pygama.lgdo.lgdo_utils as lgdo_utils
 import numpy as np
 
+import pygama.lgdo.lgdo_utils as lgdo_utils
+
 
 def test_get_element_type():
 
diff --git a/tests/lgdo/test_struct.py b/tests/lgdo/test_struct.py
index 03b3fd99f..c8ff3d231 100644
--- a/tests/lgdo/test_struct.py
+++ b/tests/lgdo/test_struct.py
@@ -1,7 +1,8 @@
 import pytest
-import pygama.lgdo as lgdo
 from numpy.testing import assert_
 
+import pygama.lgdo as lgdo
+
 
 @pytest.fixture()
 def struct():
diff --git a/tests/math/test_fwhm.py b/tests/math/test_fwhm.py
index 526886936..0f4340336 100644
--- a/tests/math/test_fwhm.py
+++ b/tests/math/test_fwhm.py
@@ -1,5 +1,6 @@
-from pytest import approx
 import numpy as np
+from pytest import approx
+
 import pygama.math.peak_fitting as pgf
 
 
diff --git a/tests/pargen/test_ecal.py b/tests/pargen/test_ecal.py
index 580bbc390..4a9e4b053 100644
--- a/tests/pargen/test_ecal.py
+++ b/tests/pargen/test_ecal.py
@@ -1,5 +1,6 @@
-import pytest
 import numpy as np
+import pytest
+
 import pygama.math.histogram as pgh
 import pygama.math.peak_fitting as pf
 from pygama.pargen import energy_cal
diff --git a/tests/raw/test_raw_buffer.py b/tests/raw/test_raw_buffer.py
index ea31665ba..96e6911f6 100644
--- a/tests/raw/test_raw_buffer.py
+++ b/tests/raw/test_raw_buffer.py
@@ -1,5 +1,7 @@
-import pytest
 import json
+
+import pytest
+
 import pygama.raw.raw_buffer as prb
 from pygama.raw.fc.fc_event_decoder import FCEventDecoder