From 8ddb0e27db2b45a760aaa94f00bc05a0d52abe29 Mon Sep 17 00:00:00 2001
From: Kevin Meagher <11620178+kjmeagher@users.noreply.github.com>
Date: Fri, 9 Feb 2024 18:37:28 -0600
Subject: [PATCH] Switch to using pytest-regressions for test_fluxes (#36)

* Remove restriction on tables version by switching runner to macos-13
* Add gh action to add a test report to PR
* fix boolean trap in corsika_weighter
* remove asfarray, switch to asarray(dtype=np.float64)
* switch to pytest.parallelization from using a bunch of unittest functions in dataset test
* sum all of the species for global spline fit test
* move tox ini to pyproject
---
 .github/workflows/tests.yml                   |  34 ++-
 .pre-commit-config.yaml                       |   2 +-
 .reuse/dep5                                   |   6 +-
 pyproject.toml                                |  23 +-
 src/simweights/__init__.py                    |   2 +
 src/simweights/_corsika_weighter.py           |  10 +-
 src/simweights/_fluxes.py                     |   6 +-
 src/simweights/_generation_surface.py         |   4 +-
 src/simweights/_powerlaw.py                   |  18 +-
 src/simweights/_spatial.py                    |  15 +-
 src/simweights/_utils.py                      |  10 +-
 src/simweights/_weighter.py                   |   2 +-
 tests/test_corsika_datasets.py                | 222 ++++++++----------
 tests/test_fluxes.py                          | 157 +++++--------
 .../test_flux_model_FixedFractionFlux0_.npz   | Bin 0 -> 1085 bytes
 .../test_flux_model_FixedFractionFlux1_.npz   | Bin 0 -> 1085 bytes
 .../test_flux_model_GaisserH3a_.npz           | Bin 0 -> 1355 bytes
 .../test_flux_model_GaisserH4a_.npz           | Bin 0 -> 1355 bytes
 .../test_flux_model_GaisserH4a_IT_.npz        | Bin 0 -> 1086 bytes
 .../test_flux_model_GaisserHillas_.npz        | Bin 0 -> 1342 bytes
 .../test_flux_model_GlobalFitGST_.npz         | Bin 0 -> 1348 bytes
 .../test_flux_model_GlobalFitGST_IT_.npz      | Bin 0 -> 1079 bytes
 .../test_flux_model_GlobalSplineFit5Comp_.npz | Bin 0 -> 1354 bytes
 .../test_flux_model_GlobalSplineFit_.npz      | Bin 0 -> 7536 bytes
 .../test_flux_model_GlobalSplineFit_IT_.npz   | Bin 0 -> 1084 bytes
 .../test_flux_model_Hoerandel5_.npz           | Bin 0 -> 1356 bytes
 .../test_flux_model_Hoerandel_.npz            | Bin 0 -> 6995 bytes
 .../test_flux_model_Hoerandel_IT_.npz         | Bin 0 -> 1084 bytes
 .../test_flux_model_Honda2004_.npz            | Bin 0 -> 1354 bytes
 .../test_fluxes/test_flux_model_TIG1996_.npz  | Bin 0 -> 280 bytes
 tests/test_genie_datasets.py                  |  97 ++++----
 tests/test_icetop_datasets.py                 |  89 ++++---
 tests/test_nugen_datasets.py                  | 185 ++++++---------
 tox.ini                                       |  12 -
 34 files changed, 412 insertions(+), 482 deletions(-)
 create mode 100644 tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GaisserH3a_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GaisserH4a_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GaisserHillas_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GlobalFitGST_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_Hoerandel5_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_Hoerandel_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_Honda2004_.npz
 create mode 100644 tests/test_fluxes/test_flux_model_TIG1996_.npz
 delete mode 100644 tox.ini

diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 0b23826..28da8bf 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -20,10 +20,8 @@ jobs:
         python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
         os: [ubuntu-22.04]
         include:
-          - python-version: "3.11"
-            os: macos-12
           - python-version: "3.12"
-            os: macos-12
+            os: macos-13
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python ${{ matrix.python-version }}
@@ -41,7 +39,14 @@ jobs:
       - name: Run Tests
         env:
           SIMWEIGHTS_TESTDATA: .
-        run: python3 -m pytest  --cov-report=xml
+        run: python3 -m pytest  --cov-report=xml --junit-xml=test-results-${{matrix.os}}-${{matrix.python-version}}.junit.xml
+      - name: Upload Test Results
+        uses: actions/upload-artifact@v4
+        if: always()
+        with:
+          if-no-files-found: error
+          name: test-results-${{matrix.os}}-${{matrix.python-version}}
+          path: test-results-${{matrix.os}}-${{matrix.python-version}}.junit.xml
       - name: Upload Coverage to Codecov
         if: ${{ !github.event.act }}
         uses: codecov/codecov-action@v4
@@ -50,3 +55,24 @@ jobs:
         with:
           fail_ci_if_error: false
           verbose: true
+  publish-test-results:
+    name: "Publish Tests Results"
+    needs: Tests
+    runs-on: ubuntu-latest
+    permissions:
+      checks: write
+      pull-requests: write
+      contents: read
+    if: always()
+    steps:
+      - name: Download Artifacts
+        uses: actions/download-artifact@v4
+        with:
+          path: .
+          pattern: test-results-*
+          merge-multiple: true
+      - name: Publish Test Results
+        uses: EnricoMi/publish-unit-test-result-action@v2
+        with:
+          files: "*.xml"
+          deduplicate_classes_by_file_name: true
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index de7e6ee..6e70819 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -44,7 +44,7 @@ repos:
         exclude: ^contrib/
         additional_dependencies: [numpy, pandas]
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.1.14
+    rev: v0.2.1
     hooks:
       - id: ruff
         args: [--fix, --show-fixes]
diff --git a/.reuse/dep5 b/.reuse/dep5
index c51c9ff..75aa168 100644
--- a/.reuse/dep5
+++ b/.reuse/dep5
@@ -3,10 +3,10 @@ Upstream-Name: SimWeights
 Upstream-Contact: the SimWeights contributors
 Source: https://github.com/icecube/simweights
 
-Files: docs/*.svg
+Files: docs/*.svg tests/flux_values.json
 Copyright: 2022 the SimWeights contributors
 License: BSD-2-Clause
 
-Files: tests/flux_values.json
-Copyright: 2022 the SimWeights contributors
+Files: tests/test_fluxes/*.npz
+Copyright: 2024 the SimWeights contributors
 License: BSD-2-Clause
diff --git a/pyproject.toml b/pyproject.toml
index e9d46af..ab70d5c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,7 +34,15 @@ requires-python = "~=3.7"
 [project.optional-dependencies]
 dev = ["pytest", "pre-commit", "reuse", "black", "ruff", "pylint", "mypy"]
 docs = ["sphinx", "sphinx-rtd-theme", "pandas"]
-test = ["h5py", "tables < 3.8; python_version < '3.9'", "tables <= 3.9.2", "pandas", "uproot", "pytest-cov"]
+test = [
+  "h5py",
+  "tables < 3.8; python_version < '3.9'",
+  "tables",
+  "pandas",
+  "uproot",
+  "pytest-cov",
+  'pytest-regressions'
+]
 
 [project.scripts]
 simweights = "simweights.cmdline:main"
@@ -89,7 +97,6 @@ target-version = "py38"
 fixable = ["I"]
 ignore = [
   "ANN401",  # any-type
-  "FBT",  # flake8-boolean-trap
   "S101",  # assert-used
   "COM812",  # confilcts with ruff formatter
   "ISC001",  # confilcts with ruff formatter
@@ -118,3 +125,15 @@ select = ["ALL"]
 
 [tool.ruff.lint.pydocstyle]
 convention = "google"
+
+[tool.tox]
+legacy_tox_ini = """
+[tox]
+envlist = py3{8,9,10,11,12}
+isolated_build = True
+
+[testenv]
+passenv = SIMWEIGHTS_TESTDATA, HDF5_DIR
+deps = .[test]
+commands = pytest
+"""
diff --git a/src/simweights/__init__.py b/src/simweights/__init__.py
index 66c631d..ed0cdcf 100644
--- a/src/simweights/__init__.py
+++ b/src/simweights/__init__.py
@@ -28,6 +28,7 @@
     "GaisserHillas",
     "GlobalFitGST",
     "GlobalFitGST_IT",
+    "GlobalSplineFit",
     "GlobalSplineFit5Comp",
     "GlobalSplineFit_IT",
     "Hoerandel",
@@ -54,6 +55,7 @@
     GaisserHillas,
     GlobalFitGST,
     GlobalFitGST_IT,
+    GlobalSplineFit,
     GlobalSplineFit5Comp,
     GlobalSplineFit_IT,
     Hoerandel,
diff --git a/src/simweights/_corsika_weighter.py b/src/simweights/_corsika_weighter.py
index e90c969..430619b 100644
--- a/src/simweights/_corsika_weighter.py
+++ b/src/simweights/_corsika_weighter.py
@@ -13,11 +13,11 @@
 from ._generation_surface import GenerationSurface, generation_surface
 from ._powerlaw import PowerLaw
 from ._spatial import NaturalRateCylinder
-from ._utils import Column, constcol, get_column, get_table, has_table
+from ._utils import Column, constcol, get_column, get_table, has_column, has_table
 from ._weighter import Weighter
 
 
-def sframe_corsika_surface(table: Any, oversampling: bool) -> GenerationSurface:
+def sframe_corsika_surface(table: Any) -> GenerationSurface:
     """Inspect the rows of a CORSIKA S-Frame table object to generate a surface object.
 
     This function works on files generated with either triggered CORSIKA or corsika-reader because
@@ -40,7 +40,7 @@ def sframe_corsika_surface(table: Any, oversampling: bool) -> GenerationSurface:
             get_column(table, "max_energy")[i],
             "energy",
         )
-        oversampling_val = get_column(table, "oversampling")[i] if oversampling else 1
+        oversampling_val = get_column(table, "oversampling")[i] if has_column(table, "oversampling") else 1
         pdgid = int(get_column(table, "primary_type")[i])
         surfaces.append(
             get_column(table, "n_events")[i]
@@ -108,7 +108,7 @@ def CorsikaWeighter(file_obj: Any, nfiles: float | None = None) -> Weighter:  #
             )
             raise RuntimeError(mesg)
 
-        surface = sframe_corsika_surface(get_table(file_obj, info_obj), oversampling=False)
+        surface = sframe_corsika_surface(get_table(file_obj, info_obj))
         triggered = True
 
     elif nfiles is None:
@@ -119,7 +119,7 @@ def CorsikaWeighter(file_obj: Any, nfiles: float | None = None) -> Weighter:  #
                 "nfiles and no I3CorsikaInfo table was found."
             )
             raise RuntimeError(msg)
-        surface = sframe_corsika_surface(get_table(file_obj, info_obj), oversampling=True)
+        surface = sframe_corsika_surface(get_table(file_obj, info_obj))
         triggered = False
 
     else:
diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py
index 5207960..dc44b15 100644
--- a/src/simweights/_fluxes.py
+++ b/src/simweights/_fluxes.py
@@ -21,7 +21,7 @@
 from pathlib import Path
 from typing import TYPE_CHECKING, Callable, Mapping, Sequence
 
-from numpy import asfarray, bool_, broadcast_arrays, exp, float64, genfromtxt, int32, piecewise, sqrt
+from numpy import asarray, bool_, broadcast_arrays, exp, float64, genfromtxt, int32, piecewise, sqrt
 from numpy import sum as nsum
 from scipy.interpolate import CubicSpline  # pylint: disable=import-error
 
@@ -398,7 +398,7 @@ def __init__(
         self: FixedFractionFlux,
         fractions: Mapping[PDGCode, float],
         basis: CosmicRayFlux | None = None,
-        normalized: bool = True,
+        normalized: bool = True,  # noqa: FBT001, FBT002
     ) -> None:
         """Flux that is a fixed fraction of another flux.
 
@@ -422,7 +422,7 @@ def __call__(self: FixedFractionFlux, energy: ArrayLike, pdgid: ArrayLike) -> ND
         energy_arr, pdgid_arr = broadcast_arrays(energy, pdgid)
         fluxsum = sum(self.flux(energy_arr, p) for p in self.pdgids)
         cond = self._condition(energy_arr, pdgid_arr)
-        return asfarray(fluxsum * piecewise(energy, cond, self.fracs))
+        return asarray(fluxsum * piecewise(energy, cond, self.fracs), dtype=float64)
 
 
 class _references:
diff --git a/src/simweights/_generation_surface.py b/src/simweights/_generation_surface.py
index 58dd694..24bbffb 100644
--- a/src/simweights/_generation_surface.py
+++ b/src/simweights/_generation_surface.py
@@ -101,13 +101,13 @@ def get_epdf(
         cols = {}
         shape = None
         for key, value in kwargs.items():
-            cols[key] = np.asfarray(value)
+            cols[key] = np.asarray(value, dtype=np.float64)
             if shape is None:
                 shape = cols[key].shape
             else:
                 assert shape == cols[key].shape  # type: ignore[unreachable]
         assert shape is not None
-        count = np.zeros(shape, dtype=np.float_)
+        count = np.zeros(shape, dtype=np.float64)
         # loop over particle type
         for ptype in np.unique(pdgid):
             mask = ptype == pdgid
diff --git a/src/simweights/_powerlaw.py b/src/simweights/_powerlaw.py
index 16e78f4..7786ef7 100644
--- a/src/simweights/_powerlaw.py
+++ b/src/simweights/_powerlaw.py
@@ -52,17 +52,17 @@ def __init__(self: PowerLaw, g: float, a: float, b: float, colname: str | None =
         self.columns = (colname,)
 
     def _pdf(self: PowerLaw, x: NDArray[np.float64]) -> NDArray[np.float64]:
-        return np.asfarray(x**self.g / self.integral)
+        return np.asarray(x**self.g / self.integral, dtype=np.float64)
 
     def _cdf(self: PowerLaw, x: NDArray[np.float64]) -> NDArray[np.float64]:
         if self.G == 0:
-            return np.asfarray(np.log(x / self.a) / self.integral)
-        return np.asfarray((x**self.G - self.a**self.G) / self.G / self.integral)
+            return np.asarray(np.log(x / self.a) / self.integral, dtype=np.float64)
+        return np.asarray((x**self.G - self.a**self.G) / self.G / self.integral, dtype=np.float64)
 
     def _ppf(self: PowerLaw, q: NDArray[np.float64]) -> NDArray[np.float64]:
         if self.G == 0:
-            return np.asfarray(self.a * np.exp(q * self.integral))
-        return np.asfarray((q * self.G * self.integral + self.a**self.G) ** (1 / self.G))
+            return np.asarray(self.a * np.exp(q * self.integral), dtype=np.float64)
+        return np.asarray((q * self.G * self.integral + self.a**self.G) ** (1 / self.G), np.float64)
 
     def pdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]:
         r"""Probability density function.
@@ -74,7 +74,7 @@ def pdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]:
         Returns:
             array_like: Probability density function evaluated at `x`
         """
-        xa = np.asfarray(x)
+        xa = np.asarray(x, dtype=np.float64)
         return np.piecewise(xa, [(xa >= self.a) & (xa <= self.b)], [self._pdf])
 
     def cdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]:
@@ -86,7 +86,7 @@ def cdf(self: PowerLaw, x: ArrayLike) -> NDArray[np.float64]:
         Returns:
             array_like: Cumulative distribution function evaluated at `x`
         """
-        qa = np.asfarray(x)
+        qa = np.asarray(x, dtype=np.float64)
         return np.piecewise(qa, [qa < self.a, qa > self.b], [0, 1, self._cdf])
 
     def ppf(self: PowerLaw, q: ArrayLike) -> NDArray[np.float64]:
@@ -98,7 +98,7 @@ def ppf(self: PowerLaw, q: ArrayLike) -> NDArray[np.float64]:
         Returns:
             array_like: quantile corresponding to the lower tail probability `q`.
         """
-        qa = np.asfarray(q)
+        qa = np.asarray(q, dtype=np.float64)
         return np.piecewise(qa, [(qa >= 0) & (qa <= 1)], [self._ppf, np.nan])
 
     def rvs(self: PowerLaw, size: Any = None, random_state: SeedType = None) -> NDArray[np.float64]:
@@ -116,7 +116,7 @@ def rvs(self: PowerLaw, size: Any = None, random_state: SeedType = None) -> NDAr
                 Default is None.
         """
         rand_state = check_random_state(random_state)
-        return self._ppf(np.asfarray(rand_state.uniform(0, 1, size)))
+        return self._ppf(np.asarray(rand_state.uniform(0, 1, size), dtype=np.float64))
 
     def __repr__(self: PowerLaw) -> str:
         return f"{self.__class__.__name__}({self.g} ,{self.a}, {self.b})"
diff --git a/src/simweights/_spatial.py b/src/simweights/_spatial.py
index 3ecc6a1..e58df28 100644
--- a/src/simweights/_spatial.py
+++ b/src/simweights/_spatial.py
@@ -42,17 +42,18 @@ def projected_area(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64
 
         As seen from the angle described by cos_zen.
         """
-        cosz = np.asfarray(cos_zen)
+        cosz = np.asarray(cos_zen, dtype=np.float64)
         assert np.all(cosz >= -1)
         assert np.all(cosz <= +1)
-        return np.asfarray(self._cap * np.abs(cosz) + self._side * np.sqrt(1 - cosz**2))
+        return np.asarray(self._cap * np.abs(cosz) + self._side * np.sqrt(1 - cosz**2), dtype=np.float64)
 
     def _diff_etendue(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64]:
-        cosz = np.asfarray(cos_zen)
+        cosz = np.asarray(cos_zen, dtype=np.float64)
         assert np.all(cosz >= -1)
         assert np.all(cosz <= +1)
-        return np.asfarray(
+        return np.asarray(
             np.pi * (self._cap * cosz * np.abs(cosz) + self._side * (cosz * np.sqrt(1 - cosz**2) - np.arccos(cosz))),
+            dtype=np.float64,
         )
 
     def pdf(self: CylinderBase, cos_zen: ArrayLike) -> NDArray[np.float64]:
@@ -91,7 +92,7 @@ def _pdf(self: UniformSolidAngleCylinder, cos_zen: NDArray[np.float64]) -> NDArr
         return 1 / (2 * np.pi * (self.cos_zen_max - self.cos_zen_min) * self.projected_area(cos_zen))
 
     def pdf(self: UniformSolidAngleCylinder, cos_zen: ArrayLike) -> NDArray[np.float64]:
-        cosz = np.asfarray(cos_zen)
+        cosz = np.asarray(cos_zen, dtype=np.float64)
         return np.piecewise(cosz, [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)], [self._pdf])
 
 
@@ -123,7 +124,7 @@ def __init__(
         self._normalization = 1 / self.etendue
 
     def pdf(self: NaturalRateCylinder, cos_zen: ArrayLike) -> NDArray[np.float64]:
-        cosz = np.asfarray(cos_zen)
+        cosz = np.asarray(cos_zen, dtype=np.float64)
         return np.piecewise(
             cosz,
             [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)],
@@ -157,7 +158,7 @@ def projected_area(self: CircleInjector, cos_zen: float) -> float:  # noqa: ARG0
 
     def pdf(self: CircleInjector, cos_zen: ArrayLike) -> NDArray[np.float64]:
         """The probability density function for the given zenith angle."""
-        cosz = np.asfarray(cos_zen)
+        cosz = np.asarray(cos_zen, dtype=np.float64)
         return np.piecewise(
             cosz,
             [(cosz >= self.cos_zen_min) & (cosz <= self.cos_zen_max)],
diff --git a/src/simweights/_utils.py b/src/simweights/_utils.py
index 6759951..c63cbe4 100644
--- a/src/simweights/_utils.py
+++ b/src/simweights/_utils.py
@@ -29,7 +29,7 @@ def __init__(self: Column, colname: str | None = None) -> None:
 
     def pdf(self: Column, value: ArrayLike) -> NDArray[np.float64]:
         r"""Probability density function."""
-        return 1 / np.asfarray(value)
+        return 1 / np.asarray(value, dtype=np.float64)
 
     def __eq__(self: Column, other: object) -> bool:
         return isinstance(other, Column) and self.columns == other.columns
@@ -47,7 +47,7 @@ def __init__(self: Const, v: float) -> None:
 
     def pdf(self: Const) -> NDArray[np.float64]:
         r"""Probability density function."""
-        return np.asfarray(self.v)
+        return np.asarray(self.v, dtype=np.float64)
 
     def __eq__(self: Const, other: object) -> bool:
         return isinstance(other, Const) and self.v == other.v
@@ -84,11 +84,11 @@ def has_column(table: Any, name: str) -> bool:
 def get_column(table: Any, name: str) -> NDArray[np.float64]:
     """Helper function getting a column from a table, works with h5py, pytables, and pandas."""
     if hasattr(table, "cols"):
-        return np.asfarray(getattr(table.cols, name)[:])
+        return np.asarray(getattr(table.cols, name)[:], dtype=np.float64)
     column = table[name]
     if hasattr(column, "array") and callable(column.array):
-        return np.asfarray(column.array(library="np"))
-    return np.asfarray(column)
+        return np.asarray(column.array(library="np"), dtype=np.float64)
+    return np.asarray(column, dtype=np.float64)
 
 
 def constcol(table: Any, colname: str, mask: NDArray[np.bool_] | None = None) -> float:
diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py
index 36c76db..2896109 100644
--- a/src/simweights/_weighter.py
+++ b/src/simweights/_weighter.py
@@ -199,7 +199,7 @@ def effective_area(
         assert np.array_equal(enbin, energy_bins)
         assert np.array_equal(czbin, cos_zenith_bins)
         e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin))
-        return np.asfarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies))
+        return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64)
 
     def __add__(self: Weighter, other: Weighter | int) -> Weighter:
         if other == 0:
diff --git a/tests/test_corsika_datasets.py b/tests/test_corsika_datasets.py
index f949de2..77996ac 100755
--- a/tests/test_corsika_datasets.py
+++ b/tests/test_corsika_datasets.py
@@ -5,135 +5,119 @@
 # SPDX-License-Identifier: BSD-2-Clause
 
 import os
-import unittest
+import sys
+from pathlib import Path
 
 import h5py
 import numpy as np
 import pandas as pd
+import pytest
 import tables
 import uproot
 
 from simweights import CorsikaWeighter, GaisserH4a
 from simweights._utils import constcol
 
-
-class TestCorsikaDatasets(unittest.TestCase):
-    @classmethod
-    def setUpClass(cls):
-        cls.flux = GaisserH4a()
-        datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
-        if not datadir:
-            cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set")
-        cls.datadir = datadir + "/"
-
-    def untriggered_weights(self, f):
-        cwm = f["CorsikaWeightMap"]
-        pflux = self.flux(cwm["PrimaryEnergy"], cwm["PrimaryType"])
-        if (cwm["PrimarySpectralIndex"] == -1).any():
-            assert (cwm["PrimarySpectralIndex"] == -1).all()
-            energy_integral = np.log(cwm["EnergyPrimaryMax"] / cwm["EnergyPrimaryMin"])
-        else:
-            energy_integral = (
-                cwm["EnergyPrimaryMax"] ** (cwm["PrimarySpectralIndex"] + 1)
-                - cwm["EnergyPrimaryMin"] ** (cwm["PrimarySpectralIndex"] + 1)
-            ) / (cwm["PrimarySpectralIndex"] + 1)
-
-        energy_weight = cwm["PrimaryEnergy"] ** cwm["PrimarySpectralIndex"]
-        return 1e4 * pflux * energy_integral / energy_weight * cwm["AreaSum"] / (cwm["NEvents"] * cwm["OverSampling"])
-
-    def triggered_weights(self, f):
-        i3cw = f["I3CorsikaWeight"]
-        flux_val = self.flux(i3cw["energy"], i3cw["type"])
-        info = f["I3PrimaryInjectorInfo"]
-        energy = i3cw["energy"]
-        epdf = np.zeros_like(energy, dtype=float)
-
-        for ptype in np.unique(info["primary_type"]):
-            info_mask = info["primary_type"] == ptype
-            n_events = info["n_events"][info_mask].sum()
-            min_energy = constcol(info, "min_energy", info_mask)
-            max_energy = constcol(info, "max_energy", info_mask)
-            min_zenith = constcol(info, "min_zenith", info_mask)
-            max_zenith = constcol(info, "max_zenith", info_mask)
-            cylinder_height = constcol(info, "cylinder_height", info_mask)
-            cylinder_radius = constcol(info, "cylinder_radius", info_mask)
-            power_law_index = constcol(info, "power_law_index", info_mask)
-
-            G = power_law_index + 1
-            side = 2e4 * cylinder_radius * cylinder_height
-            cap = 1e4 * np.pi * cylinder_radius**2
-            cos_minz = np.cos(min_zenith)
-            cos_maxz = np.cos(max_zenith)
-            ET1 = cap * cos_minz * np.abs(cos_minz) + side * (cos_minz * np.sqrt(1 - cos_minz**2) - min_zenith)
-            ET2 = cap * cos_maxz * np.abs(cos_maxz) + side * (cos_maxz * np.sqrt(1 - cos_maxz**2) - max_zenith)
-            etendue = np.pi * (ET1 - ET2)
-
-            mask = ptype == i3cw["type"]
-            energy_term = energy[mask] ** power_law_index * G / (max_energy**G - min_energy**G)
-            epdf[mask] += n_events * energy_term / etendue
-
-        return i3cw["weight"] * flux_val / epdf
-
-    def cmp_dataset(self, triggered, fname, rate):
-        fname = self.datadir + fname
-
-        if triggered:
-            nfiles = None
-            refweight = self.triggered_weights
-        else:
-            nfiles = 1
-            refweight = self.untriggered_weights
-
-        reffile = h5py.File(fname + ".hdf5", "r")
-        w0 = refweight(reffile)
-
-        inputfiles = [
-            ("h5py", reffile),
-            ("uproot", uproot.open(fname + ".root")),
-            ("tables", tables.open_file(fname + ".hdf5", "r")),
-            ("pandas", pd.HDFStore(fname + ".hdf5", "r")),
-        ]
-
-        for name, infile in inputfiles:
-            with self.subTest(name):
-                wobj = CorsikaWeighter(infile, nfiles)
-                w = wobj.get_weights(self.flux)
-                self.assertAlmostEqual(w.sum(), rate)
-                np.testing.assert_allclose(w0, w, 1e-6)
-
-        for _, infile in inputfiles:
-            infile.close()
-
-    def test_012602(self):
-        self.cmp_dataset(False, "Level2_IC86.2015_corsika.012602.000000", 102.01712611701736)
-
-    def test_020014(self):
-        self.cmp_dataset(False, "Level2_IC86.2015_corsika.020014.000000", 23.015500214424705)
-
-    def test_020021(self):
-        self.cmp_dataset(False, "Level2_IC86.2015_corsika.020021.000000", 69.75465614509928)
-
-    def test_020208(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020208.000001", 22.622983704306385)
-
-    def test_020243(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020243.000001", 4.590586137762489)
-
-    def test_020263(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020263.000000", 10.183937153798436)
-
-    def test_020777(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020777.000000", 362.94284441826704)
-
-    def test_020778(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020778.000000", 6.2654796956603)
-
-    def test_020780(self):
-        self.cmp_dataset(False, "Level2_IC86.2016_corsika.020780.000000", 14.215947086098588)
-
-    def test_21889(self):
-        self.cmp_dataset(True, "Level2_IC86.2016_corsika.021889.000000", 122.83809329321922)
+flux = GaisserH4a()
+datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
+if datadir:
+    datadir = Path(datadir)
+
+datasets = [
+    (False, "Level2_IC86.2015_corsika.012602.000000", 102.01712611701736),
+    (False, "Level2_IC86.2015_corsika.020014.000000", 23.015500214424705),
+    (False, "Level2_IC86.2015_corsika.020021.000000", 69.75465614509928),
+    (False, "Level2_IC86.2016_corsika.020208.000001", 22.622983704306385),
+    (False, "Level2_IC86.2016_corsika.020243.000001", 4.590586137762489),
+    (False, "Level2_IC86.2016_corsika.020263.000000", 10.183937153798436),
+    (False, "Level2_IC86.2016_corsika.020777.000000", 362.94284441826704),
+    (False, "Level2_IC86.2016_corsika.020778.000000", 6.2654796956603),
+    (False, "Level2_IC86.2016_corsika.020780.000000", 14.215947086098588),
+    (True, "Level2_IC86.2016_corsika.021889.000000", 122.83809329321922),
+]
+
+
+def untriggered_weights(f):
+    cwm = f["CorsikaWeightMap"]
+    pflux = flux(cwm["PrimaryEnergy"], cwm["PrimaryType"])
+    if (cwm["PrimarySpectralIndex"] == -1).any():
+        assert (cwm["PrimarySpectralIndex"] == -1).all()
+        energy_integral = np.log(cwm["EnergyPrimaryMax"] / cwm["EnergyPrimaryMin"])
+    else:
+        energy_integral = (
+            cwm["EnergyPrimaryMax"] ** (cwm["PrimarySpectralIndex"] + 1)
+            - cwm["EnergyPrimaryMin"] ** (cwm["PrimarySpectralIndex"] + 1)
+        ) / (cwm["PrimarySpectralIndex"] + 1)
+
+    energy_weight = cwm["PrimaryEnergy"] ** cwm["PrimarySpectralIndex"]
+    return 1e4 * pflux * energy_integral / energy_weight * cwm["AreaSum"] / (cwm["NEvents"] * cwm["OverSampling"])
+
+
+def triggered_weights(f):
+    i3cw = f["I3CorsikaWeight"]
+    flux_val = flux(i3cw["energy"], i3cw["type"])
+    info = f["I3PrimaryInjectorInfo"]
+    energy = i3cw["energy"]
+    epdf = np.zeros_like(energy, dtype=float)
+
+    for ptype in np.unique(info["primary_type"]):
+        info_mask = info["primary_type"] == ptype
+        n_events = info["n_events"][info_mask].sum()
+        min_energy = constcol(info, "min_energy", info_mask)
+        max_energy = constcol(info, "max_energy", info_mask)
+        min_zenith = constcol(info, "min_zenith", info_mask)
+        max_zenith = constcol(info, "max_zenith", info_mask)
+        cylinder_height = constcol(info, "cylinder_height", info_mask)
+        cylinder_radius = constcol(info, "cylinder_radius", info_mask)
+        power_law_index = constcol(info, "power_law_index", info_mask)
+
+        G = power_law_index + 1
+        side = 2e4 * cylinder_radius * cylinder_height
+        cap = 1e4 * np.pi * cylinder_radius**2
+        cos_minz = np.cos(min_zenith)
+        cos_maxz = np.cos(max_zenith)
+        ET1 = cap * cos_minz * np.abs(cos_minz) + side * (cos_minz * np.sqrt(1 - cos_minz**2) - min_zenith)
+        ET2 = cap * cos_maxz * np.abs(cos_maxz) + side * (cos_maxz * np.sqrt(1 - cos_maxz**2) - max_zenith)
+        etendue = np.pi * (ET1 - ET2)
+
+        mask = ptype == i3cw["type"]
+        energy_term = energy[mask] ** power_law_index * G / (max_energy**G - min_energy**G)
+        epdf[mask] += n_events * energy_term / etendue
+
+    return i3cw["weight"] * flux_val / epdf
+
+
+@pytest.mark.parametrize(("triggered", "fname", "rate"), datasets)
+@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set")
+def test_dataset(triggered, fname, rate):
+    fname = datadir / fname
+
+    if triggered:
+        nfiles = None
+        refweight = triggered_weights
+    else:
+        nfiles = 1
+        refweight = untriggered_weights
+
+    reffile = h5py.File(str(fname) + ".hdf5", "r")
+    w0 = refweight(reffile)
+
+    inputfiles = [
+        ("h5py", reffile),
+        ("uproot", uproot.open(str(fname) + ".root")),
+        ("tables", tables.open_file(str(fname) + ".hdf5", "r")),
+        ("pandas", pd.HDFStore(str(fname) + ".hdf5", "r")),
+    ]
+
+    for _, infile in inputfiles:
+        wobj = CorsikaWeighter(infile, nfiles)
+        w = wobj.get_weights(flux)
+        assert w.sum() == pytest.approx(rate)
+        assert w0 == pytest.approx(w, 1e-6)
+
+    for _, infile in inputfiles:
+        infile.close()
 
 
 if __name__ == "__main__":
-    unittest.main()
+    sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]]))
diff --git a/tests/test_fluxes.py b/tests/test_fluxes.py
index 2b18fdc..2c161e3 100755
--- a/tests/test_fluxes.py
+++ b/tests/test_fluxes.py
@@ -5,123 +5,74 @@
 # SPDX-License-Identifier: BSD-2-Clause
 
 import json
-import unittest
+import sys
 from pathlib import Path
 
 import numpy as np
+import pytest
 
-from simweights import _fluxes
+import simweights
 
+with (Path(__file__).parent / "flux_values.json").open() as f:
+    flux_values = json.load(f)
 E = np.logspace(2, 10, 9)
 
-
-class TestCosmicRayModels(unittest.TestCase):
-    @classmethod
-    def setUpClass(cls):
-        with (Path(__file__).parent / "flux_values.json").open() as f:
-            cls.flux_values = json.load(f)
-
-    def flux_cmp(self, name, *args):
-        flux = getattr(_fluxes, name)(*args)
-        v1 = flux(*np.meshgrid(E, [int(i) for i in self.flux_values[name]]))
-        v2 = np.array(list(self.flux_values[name].values())) / 1e4
-        np.testing.assert_allclose(v1, v2, 1e-13)
-
-        # make sure you get zero for non CR primaries
-        np.testing.assert_allclose(flux(E, 22), 0)
-
-    def test_Hoerandel(self):
-        self.flux_cmp("Hoerandel")
-
-    def test_Hoerandel5(self):
-        self.flux_cmp("Hoerandel5")
-
-    def test_Hoerandel_IT(self):
-        self.flux_cmp("Hoerandel_IT")
-
-    def test_GaisserHillas(self):
-        self.flux_cmp("GaisserHillas")
-
-    def test_GaisserH3a(self):
-        self.flux_cmp("GaisserH3a")
-
-    def test_GaisserH4a(self):
-        self.flux_cmp("GaisserH4a")
-
-    def test_GaisserH4a_IT(self):
-        self.flux_cmp("GaisserH4a_IT")
-
-    def test_Honda2004(self):
-        self.flux_cmp("Honda2004")
-
-    def test_TIG1996(self):
-        self.flux_cmp("TIG1996")
-
-    def test_GlobalFitGST(self):
-        self.flux_cmp("GlobalFitGST")
-
-    def test_GlobalFitGST_IT(self):
-        self.flux_cmp("GlobalFitGST_IT")
-
-    def test_GlobalSplineFit(self):
-        self.flux_cmp("GlobalSplineFit")
-
-    def test_GlobalSplineFit5Comp(self):
-        self.flux_cmp("GlobalSplineFit5Comp")
-
-    def test_GlobalSplineFit_IT(self):
-        self.flux_cmp("GlobalSplineFit_IT")
-
-    def test_FixedFractionFlux(self):
-        self.flux_cmp("FixedFractionFlux", {2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4})
-        self.flux_cmp(
-            "FixedFractionFlux",
-            {2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4},
-            _fluxes.GaisserH4a_IT(),
-        )
-
-
-def test_GlobalSplineFit_similar():
-    """
-    Test if GlobalSplineFit is similar to to other models within a factor of 500,
-    mainly to check if the units provided by the GST files match.
-    This can be not transparent in the file and web-interface.
-    If the units mismatch, we should expect at least a deviation of 1000 (one prefix)
-    or most likely a mismatch of 10 000 (cm^2 vs m^2).
-    """
-    gsf = _fluxes.GlobalSplineFit()
-
-    for name in ("GlobalFitGST", "GaisserH3a", "Hoerandel5"):
-        model = getattr(_fluxes, name)()
-
-        same_pdgids = [pdgid for pdgid in gsf.pdgids if pdgid in model.pdgids]
-
-        f_gsf = gsf(*np.meshgrid(E, [int(i) for i in same_pdgids]))
-        f = model(*np.meshgrid(E, [int(i) for i in same_pdgids]))
-
-        assert 0.2 < np.mean(f / f_gsf) < 500
-
-
-def test_GlobalSplineFit5Comp_similar():
+flux_models = [
+    simweights.Hoerandel(),
+    simweights.Hoerandel5(),
+    simweights.Hoerandel_IT(),
+    simweights.GaisserHillas(),
+    simweights.GaisserH3a(),
+    simweights.GaisserH4a(),
+    simweights.GaisserH4a_IT(),
+    simweights.Honda2004(),
+    simweights.TIG1996(),
+    simweights.GlobalFitGST(),
+    simweights.GlobalFitGST_IT(),
+    simweights.GlobalSplineFit(),
+    simweights.GlobalSplineFit5Comp(),
+    simweights.GlobalSplineFit_IT(),
+    simweights.FixedFractionFlux({2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}),
+    simweights.FixedFractionFlux({2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4}, simweights.GaisserH4a_IT()),
+]
+
+
+@pytest.mark.parametrize("flux", flux_models, ids=[x.__class__.__name__ for x in flux_models])
+def test_flux_model(flux, ndarrays_regression):
+    # this is the old regression test it can stick around for a bit but will be deleted at a certain point
+    for pdgid in flux.pdgids:
+        v1 = flux(E, pdgid)
+        v2 = np.array(flux_values[flux.__class__.__name__][str(int(pdgid))]) / 1e4
+    assert v1 == pytest.approx(v2, rel=1e-13)
+
+    ndarrays_regression.check({pdgid.name: flux(E, pdgid) for pdgid in flux.pdgids}, default_tolerance={"rtol": 1e-13})
+    # make sure you get zero for non CR primaries
+    assert flux(E, 22) == pytest.approx(0)
+
+
+gsfmodels = [
+    simweights.GlobalSplineFit(),
+    simweights.GlobalSplineFit5Comp(),
+    simweights.GlobalSplineFit_IT(),
+]
+
+
+@pytest.mark.parametrize("gsf", gsfmodels)
+def test_GlobalSplineFit_similar(gsf):
     """
-    Test if GlobalSplineFit is similar to to other models within a factor of 500,
+    Test if GlobalSplineFit is similar to to other models within a factor of 0.4,
     mainly to check if the units provided by the GST files match.
+    Sum all species before check because different models use different particles.
     This can be not transparent in the file and web-interface.
     If the units mismatch, we should expect at least a deviation of 1000 (one prefix)
     or most likely a mismatch of 10 000 (cm^2 vs m^2).
     """
-    gsf = _fluxes.GlobalSplineFit5Comp()
-
     for name in ("GlobalFitGST", "GaisserH3a", "Hoerandel5"):
-        model = getattr(_fluxes, name)()
-
-        same_pdgids = [pdgid for pdgid in gsf.pdgids if pdgid in model.pdgids]
-
-        f_gsf = gsf(*np.meshgrid(E, [int(i) for i in same_pdgids]))
-        f = model(*np.meshgrid(E, [int(i) for i in same_pdgids]))
-
-        assert 0.2 < np.mean(f / f_gsf) < 500
+        model = getattr(simweights, name)()
+        f_gsf = gsf(*np.meshgrid(E, gsf.pdgids)).sum(axis=0)
+        f = model(*np.meshgrid(E, model.pdgids)).sum(axis=0)
+        assert f == pytest.approx(f_gsf, rel=0.4)
 
 
 if __name__ == "__main__":
-    unittest.main()
+    sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]]))
diff --git a/tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz b/tests/test_fluxes/test_flux_model_FixedFractionFlux0_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..e8be2b4dd051e714dc14ce5ce5fc3b0a4426a6c4
GIT binary patch
literal 1085
zcmWIWW@gc4U|`??Vnv24xe-79Ljfm)2tz<XPHC}TUO^=zg8;(`pb`+60HGLWzY%{W
zEpYNgz_Ea}^Cm^jiCL1jbV1(aw1w*yFPJwiJ|t+$4DtE#Qzmf<m3yDj^keN3PgfCJ
z!t_&Wn$;{TuB}`m4D0^HM&0jvzkd7fKZ_PV*|&bv+*Faz*_Tr<=c@b8eQV~}dwrT`
z!O`lM%DiR^FU2x;=i6+4mzmV7HTQ8+Z9@Cr8;#z^Yu-q&eWG%DUlbEV0E$-<uYGtA
z^9mo(D;}vPex=DdsdVtqdzgPVLj7|Z>L1^0Q2z+$I4#}HBzFHt+TF`Xq(ZIKkNi5o
zSG|#WdRW9w^Qe=a(>6viqWQ;{HM0%bKmLYhwD-^HJKS|X)pd8yS$$h@@29?5+3Hzy
zxx(Ju%vinqVfV?|)~mMpPBJ}jdNu9V6@Sn5D}$1@?0j(M;J$~UbH26&zh-55mb<Cy
zf5O+To<dAvX#V;CpgR@jAAVr`xTTun@=!mu6Wn}Q1Z{>!&}~=*!4e!dG=dT~zXT_^
z$*=@x4NY(^Yk~e@WYT3u%{-vI0m>(wP$z?^21XDG&r0apKnWG5jSr*}3K|$s16ioa
z7F{nWQDE0Q8)zwN5<%AsiU3f^f}G0_wG%`&FlI9`phri5H!B-Rk_8Cg0BLJx5Dx&u
C4UO~w

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz b/tests/test_fluxes/test_flux_model_FixedFractionFlux1_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..e8be2b4dd051e714dc14ce5ce5fc3b0a4426a6c4
GIT binary patch
literal 1085
zcmWIWW@gc4U|`??Vnv24xe-79Ljfm)2tz<XPHC}TUO^=zg8;(`pb`+60HGLWzY%{W
zEpYNgz_Ea}^Cm^jiCL1jbV1(aw1w*yFPJwiJ|t+$4DtE#Qzmf<m3yDj^keN3PgfCJ
z!t_&Wn$;{TuB}`m4D0^HM&0jvzkd7fKZ_PV*|&bv+*Faz*_Tr<=c@b8eQV~}dwrT`
z!O`lM%DiR^FU2x;=i6+4mzmV7HTQ8+Z9@Cr8;#z^Yu-q&eWG%DUlbEV0E$-<uYGtA
z^9mo(D;}vPex=DdsdVtqdzgPVLj7|Z>L1^0Q2z+$I4#}HBzFHt+TF`Xq(ZIKkNi5o
zSG|#WdRW9w^Qe=a(>6viqWQ;{HM0%bKmLYhwD-^HJKS|X)pd8yS$$h@@29?5+3Hzy
zxx(Ju%vinqVfV?|)~mMpPBJ}jdNu9V6@Sn5D}$1@?0j(M;J$~UbH26&zh-55mb<Cy
zf5O+To<dAvX#V;CpgR@jAAVr`xTTun@=!mu6Wn}Q1Z{>!&}~=*!4e!dG=dT~zXT_^
z$*=@x4NY(^Yk~e@WYT3u%{-vI0m>(wP$z?^21XDG&r0apKnWG5jSr*}3K|$s16ioa
z7F{nWQDE0Q8)zwN5<%AsiU3f^f}G0_wG%`&FlI9`phri5H!B-Rk_8Cg0BLJx5Dx&u
C4UO~w

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GaisserH3a_.npz b/tests/test_fluxes/test_flux_model_GaisserH3a_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..495eed6f25c979e2f03af1773c1ca8b32fc544f5
GIT binary patch
literal 1355
zcmWIWW@gc4U|`??Vnqgr&5_;zp@5S?gdrdxr?glvub`5VL4e@|PzeZ3fKUvx--y4G
z7C3n#;8?)gd6S~%#4O2Mx*%_I+QM~<7tEU$9}+ZWhWPyWDU-N_%DvBM`muJ2r>lr9
zVfraG&1#kv*H$hOhIP*J_vXCXU3gFL!P%d0cHftId+4TV&P}E5H}axx3)jABo29#Y
z)s<~$+n3#X^L*8$s{s{>dxLJN{km-E{Cdm%L#lIRxYwku+Pr0JRQ`4rh5!_=Z12Bi
z1M><W&?_FPCVr*KIjJ=DkN@+Y=|;AeN#Co&PH(>>AZ@MN%_-eF$LzM%+1`JLbXKHR
zf4jbM@~U0jITrKWkKI3&x@G5~6F;^TZtvc@Z2r`<$5wi;xs;a6zw*|Hi$ZArG4=J{
zi0mIfLtOspr(XR0n||rr+tu4|yiDEv#V+8=meVf-R#a3Z&RhH4Puh=9GBNrf^Y$;Z
z#J1g*VCUR%<JJESdhzPHGCS07-1_Z*bLUgJ?ln61Q_e>1jJv>v<{i-%qjs2g_<_;m
zm}6v4`}kRP__mFBUS!(xvh0Izdsa>=dw5o{+yDBPTXwO}gLSt}xL7v%<Gt7Wf>O=I
zlXpyyFyFV$?CytPos{5~uWhT}nXj!fi`bf7x<YQ2ULq@+hmJQg?M3#GTdJuU?LD;q
zZuzt~729u~+w*%}{j$Ancm3{8`MsV0?AFjvGXn~?xSj}m{NvSi?ajAyqax1zz8&*7
zdGEK<GpBa#DtUdocT$-4dAUtdVqc#N+^&q&26~8*NtYS5Gyr9GP{!wkCMposzz8DY
zWdpi4P|kyC;{&OLf(FLZKo)8~MAr*S<)HKnvJzGAY@ns6X&zlKC`rNW<%ik{q8b>p
mm>AF#7`kRq^kLU*&CI}n7MTIwtZX2YS%FX!=q?@>5Dx(VhV&-@

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GaisserH4a_.npz b/tests/test_fluxes/test_flux_model_GaisserH4a_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..3ae51478c601b3f32b9b9d001d4810d71c70752a
GIT binary patch
literal 1355
zcmWIWW@gc4U|`??Vnv1pg4Le?p@5S?gdrdxr?glvub`5VL4e@|PzeZ3fKUvx--y4G
z7C3n#;8?)gd6S~%#4O2Mx*%_I+QM~<7tEU$9}+ZWhWPyWDU-N_%DvBM`muJ2r>lr9
zVfraG&1#kv*H$hOhIO2*<`#0IF&lR8Y^m?NTBE+ZVtZ+B(XQy+e#v!_H-0$W{q~D#
zbF}E0%eh5yd3$%fdLRATxcKeu?)+U*r&ZQ=|4H9__|(dTqYF=+jtXN4K=I1?pL$#{
zukZoA;*o0NSDKuYis_$z>iOs0w8)*)zW3b--4|}X^Tgp}=d2VfmDk@)yRG|r{*3K^
zE^LVn|78&spVH`iFa308qfCz(>z>&8i}Y?qOkW}8V_&Iz(&+K~=KVH08*ia`=lhK*
zX2{<0GsNW`Dn`$Y^WU=jbnYd6Z}>c&fB%GKbMM6UM|Yc_+5B|o<q0c3N1yoiOm463
z+`_kJ`lrqn#z%d-`)hOc(Z$S@nb&@INtwQ8<@t?Wm-pO!E1{RbhUTA#uXSW#{^18k
zk7JIJIql<T!u|TGsnhf}f1VjATDvCXQs3djin<4uUH^9XT<fNfAA7c)&NGbbnHQAx
zS@YSZZSe=SUmsgECsh8{;WgoT8y^`xh)O7ZdHVig+1ERD(ERh}vcL>v|G1@^n&I%z
zMCzrt;Pd~v>#DcieAo6r@So%_sWR!ET{6#fXJ!Ac_`7AV{^Zsi>G}3^_c$)Tn-={0
z)$QLk|00YnLtJG~dTjppW7DU5FQ-?pS;m!<`f6`R3MbG<j7+-BsHFiYHGwidC)Cj(
zs(}$i!pjD9ZJ?Y7%2Ob1d?1xj(7<>a$U@DB=z2lP9A+h|-dR9PQ4>A7UQm*P*~<^J
w6i7EPW-&3KCopu)p!fqt49L-_nyr`_7|>!fz?+o~WC1G>Y69KG!vf*~06*CE2><{9

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz b/tests/test_fluxes/test_flux_model_GaisserH4a_IT_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..3639bc7c9e1156b3164b975408c736900400d850
GIT binary patch
literal 1086
zcmWIWW@gc4U|`??Vnv1pg4Le?p@5S?gdrdxr?glvub`5VL4e@|PzeZ3fKUvx--y4G
z7C3n#;8?)gd6S~%#4O2Mx*%_I+QM~<7tEU$9}+ZWhWPyWDU-N_%DvBM`muJ2r>lr9
zVfraG&1#kv*H$hOhIO2*<`#0IF&lR8Y^m?NTBE+ZVtZ+B(XQy+e#v!_H-0$W{q~D#
zbF}E0%eh5yd3$%fdLRATxcKeu?)+U*r&ZQ=|4H9__|(dTqYF=+jtXN4K=I1?pL$#{
zukZoA;*o0NSDKuYis_$z>iOs0w8)*)zW3b--4|}X^Tgp}=d2VfmDk@)yRG|r{*3K^
zE^LVn|78&spVH`iFa308qfCz(>z>&8i}Y?qOkW}8V_&Iz(&+K~=KVH08*ia`r!=xF
z7uh@hhGsatL%rzH+ty$9yJ_Fs?a7g{xzRt*oXAT$D7J0O)6+RKwVCbj#Hr?Lmzn+v
z+H%5WcCS;p+Lb59HzSUU$?iTlJ1KY4)r!E^mz-8#I9lXxb*(s(0nI;eE(^?n`G+4E
zJ#ML{xI8qGdhru{{y%qJ^|qVu+WrUrll&!BCcU#u=9%uS?B5lCx9ruQ+?pdj-+t~M
z$HjNkf`7ld{k!I0gt29atL#aS&HsLE`gHH*^y)RsxN=fo?afHx1p0`PNtYQl_kcnl
zlu<aLjs{T;j35%8m(aC=QYt77g0%60R6;=m<7prZHQl1?1tklZm8g1W0WC#MB<OlU
lQ2+{CkiGm+J3&+fV-^zwdV~acv$BCCS%B~jkhWn4@c><?h^qhq

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GaisserHillas_.npz b/tests/test_fluxes/test_flux_model_GaisserHillas_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..958a1825cbd5cac51acfae37d29d80d3882de270
GIT binary patch
literal 1342
zcmWIWW@gc4U|`??VnqhVohFw5p@5S?gdrdxr?glvub`5VL4e@|PzeZhK`4gVZ^U0o
z3!FR=a4cZ$yh%}WVwU7BU6409ZQ;7b3+7FW4+)wwLwtVxlu2Ad<=$sB{aCxi(^bTl
zF#VL8W;M%-Yb%!s!@7+5o93OZeg5`a{;x~c$y*oSP1EJyy!5yEZM)2_u9TP+S=Eb-
zqBg20ZR(xX`@XAn-SjovZ7OF~SI%GibZ+*dXPf1ht}<|7WC%d<gqL2j9?TPbKu>t2
zn)sC_=cLlqC(;w|eSEdK?}n?cNW`Ct^;_rGemuguIy`5S?W(2v<xkq$-cA<2Rj+10
zJ3Bb*jr6;<zfSLZ`{3nF@fn}j@U`pD%9vK+FrSeN%_pKeoR%W{#Lp0yPx`4BC-UFc
zmw&Cw*%*8!Odzc9!iLq?c}?ynSY6*1oqMRb)ywhrrjzWoH~vKDEW098SD16w{q^?S
z=blX4<sB{0Kl4(~%e|Lu^6T{`M{Sw}Og|_QV|e8|AIv-az=(0oF*3&yGZU!ipQ-w0
zx7Sv0zq0Plj-Si6_Xw3S9n{WSvZM57UWnGY?RUN%fA+J_zqsG#+-A$%nrn*|Jm@QR
zI&!>Bu+B7FJnd!W=CFc@xR@vFDhgjR^P+iZVx;MQWDmKenwsJ85EUb+!SwQzZ{_*h
zxOd;$R8_iiS8VHoxnGOET5dCL3Vr{!P9x#8+}(c*1J)cpruQ&Dc=Z|c*Sec7HU4>+
zcunv0w7Nel59X9w-oL+iZ}w6}pnn*dbeT~L08su0<#kSI9Dt|>Mi2=v4$!rMvK&ks
zA4nw>G%y|ovQRT1x?WHs2ZavEN>si5Kub}RJi1;`iUK7SkY0WuhJgmgI3@=4^o6b&
i6nij>Q8nu`GccgVWq>y;8_0B4Ak+l9>pxHl0|NjFtK>rf

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GlobalFitGST_.npz b/tests/test_fluxes/test_flux_model_GlobalFitGST_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..7e94810604c3f150f9300567bcc88211ca585dab
GIT binary patch
literal 1348
zcmWIWW@gc4U|`??Vnv2Dt=e<{p@5S?gdrdxr?glvub`5VL4e@|Pzeb1LnwyXZ^U0o
z3!FR=a4cZ$yh%}WVwU7BU6409ZQ;7b3+7FW4+)wwLwtVxlu2Ad<=$sB{aCxi(^bTl
zF#VL8W;M%-Yb%!s!@7TSzaHW`_x#PJP^s>HF@;;K*G44lE`Kj}yKLj?-JjC5Z@trv
z{`gv_dj46l-z#UW)74+|_PfW|D{nhCtUJRuEyTn=Gqq;H*F9fXPo2gPfZ~&x-<-=}
zKH&rU#3R+juQWL)71KN2)bmd7Yu~)PwdFU?l^+eg-<Bh7y|(GMESL7S8~a|}e_x~S
zc={Xf;s(caw<67+-?yy~E^od1Ec?JQYd`*5YxUCJmhZ3m*=olu!hq%#R_TTkWUu%c
z;_?a=V`g&D^P_oHTW_y@yyJA)fq93|rCoeFqw1{RZQa~S4@1AnY&-9N)$CMB*bSFi
zXUo6D{l9VRM_9>?U}@IZE2X-{Z+puayyyC69C_r8CM%kE_7}fC5AzN`FlHQcjLdOF
z4Hf+($@TU7&Y8d8Z1dj}6O(o#Ki%x{)k5{noZGmUn3ZVF=Du-`Q-Af;|F^Ry*WZZC
z-26NH^>K0gyxRun4@d2s`;u$@qUNR9Wy`MZHS_gF^Us+{hux6<<Cbb_hQmJ-sF&1a
z%gT-|6kS`^Yjo!Hx7ughtOKHc7B2d!a$Ek?-Jj|;M|0)hMVS^ZKfyEm=69~3@{<3X
z=S??C4&IWs@auK&h1ciT#^zmJe4_B3`_6TuKo2o8=`y1h1E8=6<#<l0n?Y0qBZ!0-
z3h3HEc@C7DK-%~~Dxsi(@f47Sngh}Gf)Y8*N>sg5ftI2sd33#?6a~}E53>|VH!!9!
qF`%a}bj_gX1H}u-(WshDm>C$*A~V37l?`M8D-db|-NniR;sF5bqwVGZ

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz b/tests/test_fluxes/test_flux_model_GlobalFitGST_IT_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..461896200e7ac32fcf4457e0c066094f43b7ecc1
GIT binary patch
literal 1079
zcmWIWW@gc4U|`??Vnv2Dt=e<{p@5S?gdrdxr?glvub`5VL4e@|Pzeb1LnwyXZ^U0o
z3!FR=a4cZ$yh%}WVwU7BU6409ZQ;7b3+7FW4+)wwLwtVxlu2Ad<=$sB{aCxi(^bTl
zF#VL8W;M%-Yb%!s!@7TSzaHW`_x#PJP^s>HF@;;K*G44lE`Kj}yKLj?-JjC5Z@trv
z{`gv_dj46l-z#UW)74+|_PfW|D{nhCtUJRuEyTn=Gqq;H*F9fXPo2gPfZ~&x-<-=}
zKH&rU#3R+juQWL)71KN2)bmd7Yu~)PwdFU?l^+eg-<Bh7y|(GMESL7S8~a|}e_x~S
zc={Xf;s(caw<67+-?yy~E^od1Ec?JQYd`*5YxUCJmhZ3m*=olu!hq%#lRauS$X@X`
zG{fN)D#lFjrdL<zPdmHG_*3tOxqF>&>8Thm>3R7@=5~B}?lwD)$s6C!P%r&I>yXQA
zb<_8}uUGGn7N7ItLFLwK$AA4WT&cdNi*MSxn=_-n-!d^}M2(p<lMcJVyu%NS8MjnZ
zT>hCry{M5bD?7GObZuF$(V5fVYM*Vh4v6|$xag<KZTVAof2!9U&6R%_Wm>rW1kdc7
z-?@UyOa5=3H{B>XcuU&Cuh+d7UY}nZn|F2biNbg8JJ*Q<<A;$+ml-wNfWjV>LpY&s
z22l-+AQGO5(6xcmC@A%TwDEydLO}!LDIg0q#iHv4B?_38sCuUYEk#Wt=z2jh018@=
gz5Gx+K~w`{3KIi*d<1y2vVkO7fbb2FHf09!0KgTCs{jB1

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit5Comp_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..ab596b4ce7935770b43649b8d433cd43ed38d960
GIT binary patch
literal 1354
zcmWIWW@gc4U|`??Vnv3h)hS>9Ljfm)2tz<XPHC}TUO^=zg8;(`pb`+60HGLWzY%{W
zEpYNgz_Ea}^Cm^jiCL1jbV1(aw1w*yFPJwiJ|t+$4DtE#Qzmf<m3yDj^keN3PgfCJ
z!t_&Wn$;{TuB}`m4D0sX-T%#R^S+$rAM)QnSryAyCii;UakkpKN!D{eCaNp9&wAUq
zXH~_#KapmJ@AdQ^pZS%y*XvMlOhVO!ucvs9%U_S=i!|-JJM;Ub`G>g}0#Lm2ao2BU
zm{<6KUhzmZ@heTvNyYR}KlS`0$`u^)dB(XLE7Qb1rm5e!p?i0GePhY>wA<yIYO7<;
zue^P4=H2r<wxsUzoqBQU?0H)CYnhFje#_72+VC$rJ1}m4S@YG_pP7HMGFOG6d8dDx
zNdvNX{0!;f9l!7G{qv>XY>k#oR(kON!i9HV4j0VMKh&M;D$Mg*t?=}%Te+42yBsSP
zosE2W-8+3o?dHCy%WF1xDITlZml@jT^Q1_-y7W4Guj_hUM%383*mw39%sc$R*m2A;
zGRF};)Jtx9ac@6->7Ld7c82IX+nAqm=ltA%zdzIKdUsn{?bFxg6(4WDtKZ)J{r;?=
z*~Q|`U&8k7&Oe(~b2(_mt<%3&9BmUgd3<V^tNH3>ai8udaiYaflv2GSvWMJKP0eWT
zp>@4i-%qd6i^w^?;L4_Dt9>u~UCN3+#?&9Vjmyi}_f|N^t#t)k_f1m>OEtcGp(uIR
zZ|CgGn<qclv3Znl``Y2y;`YzWBk$IHwtjbTc9J*HLySzi%&3I{DDXk~o)el!KvV-G
zh=dmn=-NOT50s@q+W0^!p`d~BG?0aw4bk<2QaN_Lvw)VOrg?O|pd<xzE<emtAl<;2
m$;5!3z|b{=q7S=fD`o};w8#wbW@Q6ezzT$#KzDJofOr5DN%On_

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..fbe3310e17e930cd251d8ffd10318191eba777cf
GIT binary patch
literal 7536
zcmchbc~n!^7RG}^6;!BzP=qH{3n-(^G6)r4(F!EU93nEvAdrM2GRW*mkx3Du5+;>F
zfkeU(h5)fzPzaz;1PKvPMkRz4Frb9t<@Pmi<zC5usbN9FTKU#Fzq`-b``g>rYW)Tc
z7!0<V`@0)97Uj)e{re>flZM&aqC-P8vA>)bfqe^m#{CNLy8gd67;Z{-QvTbvXSS`j
z_`G!2oRfvvq5|wq-#dhR1$mk74t5#8$>zIfq>D+0TXpP*Z@eb!_k&E~y3ac^jkAr#
zYQ&^rgvJ5NJiMwIUGirTYqZo&JZu*u^XWG+170S%GoC7Y6SA2X8cQScmXHS8EKSXk
zjzw(K@yGTjJtD8O+qXVde2o-G>R%hUy^x;&L=0vNLOEV?iO53{=b{|*(SwJcM*Dop
zlO{AC&3u{gk*81Lv_wXTtfrAKW+J^@n{)WRFP&8^Zb&Wd!60&@FivGg1LKr}rN>m0
zRhKxD#^Q>Z=A-@$dP`>6`Ad#<^CsHv=8kF*9O>4f9)28ae?tK{H-v^mxSTEL-dy!Z
zAKqyG2u5uLj>@`)!&)y9yl{IQb<cAeuNYwI4>e-(`_EJ{W3#>5NM4L5A)n0Kr;YtW
zoui&-sQF#Sr>iN~&mnTWMIk`$r=FeW2eR-n5&$H;^QOsa+=4KiO04egy4|Tz5t6{{
zx2ItBXGQyV>SZmRZ?5i8Rtk$;+;~T}S3<=)4UYLpGF$Y%@_TWK8eRrxL!L7lsZwGg
z6vtlBh?x+cn^u>0f_SvG1@H(lKnWt|SDFtv(3+M=vv0QaqDb<u_Z`nAGrZ|!J*_bE
zbS$o=R-(`Jjo+Em^#Q}rOR`E_%#)FgcYm0|-Ae5mlp;IimVb2Gf0a@9oU{P~bR9XU
zzz=jpTSowp@B&0Iy-HVjKA~P2GTeO5c>dx2+-bwIw;wyi>AFFIaDDW0P<PgiSjO8X
z3!IXCS!##r*`CEAsqrR5F?IYhYkvSo0aN4U7L?t5eXddo!lSp>YBxU*Tw72=5F(F;
z32VzPdVhf7a$NHDrfGYW)ZTpAT`514!Xs2>J}t$ZxmRUeKFq#)0GFqzr2FsrY2sjU
z!O5A)g1c-I>w$#{XAXaIf{5)juWFNo@Z3lYxy8?e&^GwT<D*I5)P2b@w!BhzE9Gu<
zzXc6`@$iooZ%4i9ihgl9as5yHa1HQB)wdqDjAI#!hAEewd&6i;b6RHN6)Se1oqzp2
zXR1bF9<Db$HAsY>w<pIeX8Cz6wFM;v;U#Rs$TpvaqQOnH{pi@n{(G@0%&O7r^E;)G
z^p4QgJagY%-<#A39Shcuhn>n-hc-30X1G)#%|<0%!l<fB%L@i-+ilcBYxVRA>C#q^
z@XQ&1Lh^7Vxb+9_qjNw&e1ynzOWr{|cwaSbAg!*(d$pTFvT;NwKQ#2EISq)mMPom;
zx8~n#@<N2h+r)SvgB=)f26b{q)`_T!3S1eXB5nlN=;Sqrwsp15X^4R9sT7?x3F7h8
z(f#LnvdtUZF{{dS(FvQ0TH~-cPJX0@W1h7|_!joXYpGd*7z=i;eb@4iO4=_)qdC3!
zos9mQ(TDt8^N_*rs2h=1=|jh(oiZ%ros*&IAo*za7(b7VpUzjv2j91yvwG?y#Y|5t
zJd^)=tP)jX7}VR*vx5|9y&e(XwHOqFh<R}L4aWh;a4(r;!W%u7BQ9i)Ey$=bROEIg
zttJ&B<&^yMQMme-(0rg#FTH^u=rCGGkVQxpo<(qL_+T!lSIHFr^d0f!z6cK!&YzZ?
zJ!L=AElASDgz}@9KUyZ_pOXFQeOm_YGAD3e%ySm|Y|m`D#CFQ+OU5UD<J8fWR-VzI
zs}4DDa9sFHejIy$orC}Ud6Pou8b{ET`~H>f**`$O^1Q1{9UpcTU4$PcI}!a1WBXs_
zM<8qW?WS3>C>UbtlpN!>xIFbc^WRWGv(<ehgfU&!m-EfU$G9pekDs~S0FMXoqqe$&
z{DBb5&-$H_^nfPQDYIBNTE;%&Qv>Tq1;NA0HYft+d|sKpOyB!vVpOh6dq4>G!BoPn
zsM+#-)6w=6R46r|^+WM}jme8nQ%EhDyzv9I&~*JWXW5dU$6i;EJrJVb(W2&aA}MV(
zQB7BkGAVOcn1fL!Wse6QV;(Nk8G)3y<(-Me2nCfqIeJIkYORucR4SWIua(`LJ62(~
z;xSu1_xxY@y2oj`x!Z4hKn-b^ghv(+2uRpR(7O5p4xA7pv{yvTLUvvBRH00cS%&%N
zc>Fs;(NQfGMqg#&%1BZ;(tB=lrKM+&vk}?_bz~@_P})^DdrB6y%E{{yug4-j1WG!!
zN~lxBmR?msBeaNNYw-geKC36F0tpdl$d#_yb?jPAlIv7g&A@y%=0gknT1x0NDcf$3
zPqP~b9Vvl0U-a3%#W&UebE;{}DMrsx_4zo1C7+s6iHY^Qe-0_5Okk<^3Bg}D(uByK
zXdV!7g{*YHf<vgaXU~NNi4(1EtF(qUt!7}zCzzr4+$>2;suMCpXr!mT&Y0u1jrj$g
z^AXG2s7dX3L+t8zs+hS~$8VA=jLNNBRPXB0CYK^Zp$=g&dA%JE=dT0j`BmZ*WY2<Y
zc2#pKhP%mQt<TZ22{A4Wugw0KkN=31OT@ux=YBnqYaZ9}OnJ&UF#P4z?OwNaH=CJT
zeQSF;H@0`)e(%XnQ8~p>AXh<4Dz@x;2S1SgX+1$VNO+b<(|qRgAxX*UFjdN!5vyIH
zdTc1x?)1aT3mEze!ugb+W=$gFrflkiw^c^}A?!LAd@7N>`xf4~d-_miK!EWH4)v#=
z_FVjZHyeu%XhpRZwK~Yp<KVCN6-tl^IocsDQ;BY~;{TbZzS4}3d6j4%-MWjk(xi;;
zRNaI&o6f!;(P}p^J{-;(@M&G{eTf}EyCrjp7<pF<{YrlPj8CRpTo5X$0U955f>aBS
z2e>~@`hrT35W_Qb;NDL~<V(!HFy*cVLW4IoC!OTGf>A>yh%1F3YPB&4G4*0HwiTkI
zyS!CN&IYbpr44oO6W?zttDTGeFq^AOj4X|#9-P7zKyi-SC#&#q{<=SB^#zq6A@bZy
z$~<dP{2(k?sxrK2!mP$-CTex%_00El8Or!oH-ZZ0WmQ&Z-yzyFOU-c4_|;*#-qOqt
ztCrd-B*p8`G<@(}ZRy1Fo?#>B-=W7XUHkO`ejXdFzM$(P#C*LFE=w6QRqLO=JJeQY
zh2BM9>95E8Q1BFHi=(B<R;ASb==K5HpmcFicicRI)8lAge*RqMl~`$CQYz(#PTOcI
zz1`aM#nOqKF9H<zV_6tKkfo3QS7=dtd~9B?iXuqe(%u7gD;XCDVu&SqCa<k4F%B;K
zpG~*<(-ls<;C%ACVdLXAGnbGS<dQGF7_7XIh6{Z{=$vFNz^m;>jP3+-6CKndz^`<;
z^8+0T(*Ft-NWcDM{?7EmsaHd&ab+wL9o6NWa>q-N!W^~>WtTMJ4dw~GkDELd%Y)iG
zw!2Cy1YkzASTAbdRrEB9o+bG8Z74TDFgBa)ukBn1=~2n=Qt13VaDRb3LOgM5a!Ch|
z^nRk{27GX7Dx1B)j43k-G$}iT)~}0_&Q~|`ctU(Pb0W#r28Uy9iw$5SDoC^uRAcpO
zzT?aBpQ(-djpdP@YU6FzbwMIrAdz(%>%kLjV6VWiTJ-NW2Y8DY0sj7fo(+`&(_$b=
z05Wmln}5IJMX20Q;Mp-$3yfv~kNc%Io4Y7@Pz%)pQ&8S^!RrRDgTcTvQK%M}IIg{}
z$$G99G=GF@f#KfT>u%iugMmkVP%SV>TU+~vs|AnMpju#PwYIi=BMb%}VL`RPfN1R)
zMSlZ>fyYHqEijGYX(f2U<-W&@SZ(6&89bYTih<F`+G6?5TrqeM0u=*u1|BhZWd>Ia
znL0qlKu-@e>%h?hiyOIz8`9xJ#XtwW_O_G6U@%BO4HX0JDsS81gNfS0y+@$t6{-bV
z!?h2l+E(tuphgj@1-dZax?tiZt`^dhLB&AFv$nW%8~1U6`yZ$n=q`AV3rw8H{TJM8
zK(&CCUwdH;7vc!m*h9sDZRRZuMmGDF`?$a+8Y%`n>e^!Q?JyX`u|mawA>@g{$b(!l
z*d9W)fcxTU!Mnywa*@Ga3@Qe!&f3o_Mv8kWz{UqE25iOJVz@Ln4}eVuR1B2xwYROf
ngS%~TL5GTgdb#$>EUp+*J#DQva?=tFrojEl=6-r41N`@2{4(J}

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz b/tests/test_fluxes/test_flux_model_GlobalSplineFit_IT_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..ffe04c14980b5b685b864ae1afcaf758e0e98d4d
GIT binary patch
literal 1084
zcmWIWW@gc4U|`??Vnv3h)hS>9Ljfm)2tz<XPHC}TUO^=zg8;(`pb`+60HGLWzY%{W
zEpYNgz_Ea}^Cm^jiCL1jbV1(aw1w*yFPJwiJ|t+$4DtE#Qzmf<m3yDj^keN3PgfCJ
z!t_&Wn$;{TuB}`m4D0sX-T%#R^S+$rAM)QnSryAyCii;UakkpKN!D{eCaNp9&wAUq
zXH~_#KapmJ@AdQ^pZS%y*XvMlOhVO!ucvs9%U_S=i!|-JJM;Ub`G>g}0#Ll-wY+&F
z%qx6AuXv=I_?0H-q+<G~pKAWmn=e-W`1o<Zo1L$V#q})Ha-UCHWBA+gMcVE6%lOPs
zuezG+KeswJY_;XxXZsiYIdt0ke?<K@-CzIg^bXvba{K7C)s@YmJHA;P#@mN#qj~30
zsGbk9cl-^_XzrbL7LDEJESBGH-#xjeXy&Wj6!CYG)1_YRykU0beP3bD$~9}7WSuXr
zebH&A`ulvoS?15Qwc9Sc*Ll_SeBI8g9xOXecY9Xhoy7mo9!h(od8h7CmLbeL{J_|8
zOEtyiAL=DHy~IBDP2O9#-7x$(!>dGc=B>X0_qWa{yd`v7b<wQvpZm0Lu)V(j&%%ED
z;b&g2C&#P)w~Nr<rKG#G#OB(sNk``&`!v6JtyR{m6dQ?ckxW1jF*4~gqvjouTS3``
z6PlzzR0AW3gy$r5ZJ=}tN`WA4d?1xj(7<>a$U;rE=z2kk0=wQ>Kub}R2)bTS48WYr
f53>|VH!x;0F`&msfHx}}NDT`Rz5&u!%pe{BkO_(e

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_Hoerandel5_.npz b/tests/test_fluxes/test_flux_model_Hoerandel5_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..029d652ee3cdda7a0d862edbbb5862f258980d37
GIT binary patch
literal 1356
zcmWIWW@gc4U|`??Vnqhq9qgt5p@5S?gdrdxr?glvub`5VL4e@|PzeZ3fKUvx--y4G
z7C3n#;8?)gd6S~%#4O2Mx*%_I+QM~<7tEU$9}+ZWhWPyWDU-N_%DvBM`muJ2r>lr9
zVfraG&1#kv*H$hO2EF>+sS6H0p8iH=^Z%2_-j;9L{A}i%_giWwzZIDo(0%8s*EaRK
z=_{X?_-uK(?nlViE01SwJ^MDADZS!-#L2eeRGnA*lPyzUeYc$c|K3Ihh5!_=2<`q+
z3G)gc&?_FPCVr*KIjNZb>8F~1)+OAGURWG2v+;4yj}!UV(h}37<u-jg_Rw_O+{jqX
z+3&S($R+>#W4I<|lka+$=MVFx=C3`^7MmKf_f(eb^uroT0*6vdKe(0H$6M`JMe~l#
ztJ+)0-tjZU<(-Mt^G;dwYCBP#-EaK%)fWAIC7Nq}uu@$@U4Cx%>76fsOqp!{Bloho
zyj5}E#*52d&U;<wx0nCg!{Esh_n%}%*xz{ir{t$(?CR(H>#}UN8Es%e^U%7t%dKD@
z;s?f$V~&wIjtHV+f;)TkknXDVvKx2v>k2+E?b<f?^u5>VR@=^QU0d?Y*7J7cfp4#0
zJvH7pKWod%^pq{7_Fqm``8}D~C0smr?fO`=9~rk^q+gwM?Budj@)?|H9(uA~st?&i
zZmFhbH22WDwfaTpMd!xnmTNvgd%4<r>YIC0HJ92M&d;2^d9!QJ|LqNH-}YU&9$9`n
z=Ue^5lHI{y&;LqWbGhd48I#wWlD1jO+s<BlE@oBTpU0WgK}n90NtYS5H~<AcDCcuR
zlNE?+U<8ry!U0_yDC>bT6-XN&NF@|BFrEgoP%|RBUQjv*rC?0Gvw)VOrh0U}pd<xz
uE<emtAl<;2&BTD7z|b{=q7S=f8)gOuw8#wbW@Q6ezzT$#KzH%7fOr6niuWD>

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_Hoerandel_.npz b/tests/test_fluxes/test_flux_model_Hoerandel_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..ed13c42ca8b516460a410b472faebd0f21dfe5bd
GIT binary patch
literal 6995
zcmchbd00~U8pqLcO3furEw#GNj4fuCYg(FZW)q?2j<_q3iddSOCE1;3s)=q^DlV9m
znNA@%N?<O?jAd#>w5ceX<5+lI&}@+mjp584=JA}1|7rNhf#>1-{LbgR@9+11-`5Mi
zVBtm>47OPQZw7O@yC_-q_e&kN4Cdv9OHQ)FA5K$&>B3s%KLK8s{&#_qM~sF{bX!}z
z9(d7;t^^eZm*8(+!)IZx({EfWB8B;SXXhDF!m_Vu=%vCRt@ltnZ-o8Uutaszn6$NE
zs|HhJ8O$mpKId9vf9{AWL)Q90oLa>Y<&TWrA!Lb{<$Dpk+xD~8p3J>9aP+{PLHbAE
z4=?%)m^?8`6?<$fu=Pw5+H!~!<%AxRIOWP-R>ELjAQZj2k7)`NZ8^$5w4FzC6b}9W
zJehOGLl5snT~7)#t?a+>u_d7$6M;pos+xGvXV1Ed3|UVevk<I``aES_9$kf~*fY?V
zU`(O$zY0P7vLD8-%)PhnYu!e4@|(R!UBk9=%po{uFX~Ru#5oYNWfq({3=grS36&Fh
zUOeJf%=o&DUCHd5dK6)HHJ0imbcm`D5lkdIe$2kIA;N*-DkF3-KAlT0+V@79Z^fq>
zaxz}jO2*ycj71VPeD5?`1;W#GetO|do?U3CdFP>T<Vdqd9}_F<>?+2DMh!kJ{judO
zt6#?8$<(lil~E#ZlYXkI^jO(q;qu_6mrs@+M2Fni^lq}G4@d20{|BMw&#<i69LoGS
zYJX1^!V`e#sR!{`ThGEXhv`=4V>FnIC3X2(3|EBbCErU)xNqvR)5t`qPO>O&#u*B5
zPj36O9kN>Y29>=(o#|ejEOO7e=sBh3pJqh-<7D~<ma}o)1#?U19w<(ScLH`M&Tea)
zdB?G0Z4Aj@ZE7U2j}^HP;w++bpO+qan;J``wVpVg6|{xZS<ilU3|X<^IcER$|I!xl
zwv@k4*tW($y^pD#&*V6zlP(#!$CeJDd=A1P(Q)*y?3#(=VLjVnFn7C$ZZx@-vsU6D
zNb?kY*ySTWEmj$?+1Y~_XOZcIM<aq0k1<{BOpJJSe0+WG;Pxq5+)d^ucd36_vkt95
zaHpc#^#l&NHIpxUce&XR!jpf_^VgX?p4JZY9v-^THZh8lPrE}b3`>94qqDxw)x6cJ
zulhL7$T|9uO}&)A20`u>ArSxOsjb=>gRflFUb)w{ZB0`vQpVTwnqGGk=Vy`UDwDby
zK*M8P8oFU7&R**;a2VtzU8htC7j+?No5yT^cZ}Y^xzRs3Mc&J@A*?MoVOg?bP({DA
zB%eQrupU0Rg-_wNhZK+}j7wHAhHeK+ooj|Jye`$fc-p;@kA&79mmRDY1&@xr_IRLe
zHqR2DIn1{ejLFg;GDfTPe>uV~+&;jeiQFng%~#!6hE2oCWQhy2l03m33P~>Eoz^Vw
zyqJhB50P0kgid<i4@%$?26ly}QTdsLo{J86S1p9(o9&L25fD&>jqUsc8E9^@SslBt
zDy~{R^mxjEfzVTOy1JmCPH1yjeatVi8=Vz>GQBe(CdB0IGvwJNCwYv|RV~jr@$cYN
zln`T-&i>Hqy|&{!-|G^p8U*NM=+bR7fe=`mFK{9_y?j)ez%3QhMwbS$YIB76#`>2%
zmKIv9%f-SD2evetk%A4`aVDZbt9Lt(G+2=khz(8T6dvAlufQ$kS;HXq0u7UHPh<9C
z@ahntD>3*(Gl6#EZ06q}&`DO~rN@b;+a_7a+~)(`)vw3zzJHJezea1MmY1*tw1t65
z``BrbIE-x^m*n_#6PK31;0L=2%c&T@yKhi}$z_Jr1-7D>i8|P7Xv$qHFkCki$S20e
zah8;u!-^Vm#VjDYM${>`TNlJKyCTF>Qp%gNum&u0qabsrJPmV0xA6CsoOSmfBv49v
zR<L)sZz6QmO8AC%Fu`}%h{X{{s4JW%mey-Tk{|&Rtz>}=f|u=V24N15gjRT-v8?2i
zAmmj6KP_!K6jhCKL25d_kLP+U7gYyZ{_>RSlBG9={1sy>9h21YOlYaZ$J{jY#I~WS
zai&?2s?J8Vm345&(3=*hLvSQ9LKN`<{I!qmYz|=#d5A)42Hz<xHo#%>jVke>gtfvy
z^P64OrJW6l_!uqt@$jUs?T6jD#fMEC##;3=d7SX+&C?&Xv(U%eZ&qKRbyJ*z+)Qv?
zgT8;Dp&RJifabRf9^lO0jkC3%CEw-_$SPtmeK$tab=w$vn-TrFb9fzoTtYSuJ;Ne%
zE2>B7sOTGhlrFuDU;(D=xkl6O;h=H@A8mL>UOL}ma=(5+Uug$vN!(|d!-Whi)L~P#
zf_-NK?M$@$0!2tjEkzu<T7;`^65nHSIrPq#cRvP)iHN!45Bf>%clkXD{n+w{wC6F=
z?}m9jJqH`$w2qboFVsgPs~W!XAqxjAa7`)Kf2;2F>ehqopo%!VB?=%Q<KVWlxrDhp
z5Ugx9GwtrtE~il4Tv>N=DNOseoZ#0V1i$%*jvq=r@*YE$6{*Lbc>e<D;pEaj@!8?t
zCob<vniIDH*|sKdiyHArAi<G}(>LHnK@%>JhFYw^`77ZfzCe7e24f^)y{`Va)Sug4
zz#qQ)N!;F_yNVuW8aeFczlpK`{m;C!>LzKf$$4F-<dNkEYYyFFa^F5FW)Kdv-|C9g
zcg3h9#rcB~Hw?^hP@vICAL2|PpD4T8cF-IK=sNw@Cc-dbB;|l3h0GS85*N7yPd#Q_
zh|d^w=nj`wKPshgo-(J&**tYySxlCnl-5M5Z`<D{3uj*Kh%DLzWB9pQ44g%(j6hFp
zn<b_HoNz#Zd}HjsKm}S<#$VgrNs5p-jI5W&ewMOsv3Ri()c5g1?PhlFsAGhWVPZFX
zvg#_2bLrVO4{tJ!a;Rak<xjedqEgfBH*a$21_=-SG-=18pe`G=?3{rjKEVF*x1awC
zWJM$QtQw;<2|k#j+QfQCo=&{*&1zm@6oS)9>vrD$fv;_%{#{v<VaCyT%V9+H*vimD
z?1s*+zE5U`airhBx#{Zp9SYI&xG43}o~O|B;4E_AB?Zr48z|9!wh}aV`9>EBn%(t(
z)D-X6u>1DME!(TKE<Pbl-nP5X`r$|~2YIU*Cn~Q^>SlH}o{`pc`kp5FwFjq<VO-bC
zo76xFy#H5obLznH8bbQ2+t7Ub>7vs=BN2cHfw!Ok{X@^q$vOFq?=zCEW>@4gN{Ftj
zz1ElO=Jvh#{6s>rB`Y=l6_MY2S!Nm&RruD=eDJ<~TpfRdPSwJipe;8;d%M1-;F-_c
zqInCr-q3tINV^ZRq<hi!^WUXf(JBriJ2eA3#lD4;g_Ve>c+}V@Z1cx5VXgDNR}}qL
ztjUPcxTr~jF=_8cBsu40AcsHsRUEH1YYl3AeY7?<g55?maLpDb25nUFf~%-*TmT-A
z0>vMggsT5NH3TlFRDkz?hNVy$Fp2~+5RhpDzx?}kN`))`1RhL6wZPm@u_RbqAYT+b
z1%zsW8Jlu#tSSr!p147^z&uO&J#Sqg*Mg>7P%SV#Qm$2B2!nw~Nl+~?sZoB<0=X7E
zw}EPbsfhA*e^!IRz%vr478pS&UpMtD7z{j+fNFtuU!m1e1Xq4NrQ*Iwe$U|MA1Ve~
zYvtnQ@|nB`xY34+fi_mTcu+2eG__DM&@w6)*UFO{(m+DRKzpWm)?j9$1`Gyi)}Ugb
z7gDSYPHL2<d>?^2Bd8YWJCtiTERin^>O!Dez~d{l;G34pwGhV-6$5@+x%iQm{BeQZ
zHB=0ESH<H36JM2o1v^=&7BGs+7skmUc7tsrR1BCh#lm1@7hU<|0$Vhw7;r$!#ac^Y
zP(K6}19n4k+hF7#xfpCdpjx1`D_@wVCr1VsdZ-wvqsqkuefg#US5>GOC|b(Jh-LC~
Z09?Yn;A-;20R}Uc|NTq;X{7=1?N7(F{Y3x(

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz b/tests/test_fluxes/test_flux_model_Hoerandel_IT_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..ebb94c5bf3fc5486023bf3ac57ba53b3ecef0051
GIT binary patch
literal 1084
zcmWIWW@gc4U|`??Vnqhq9qgt5p@5S?gdrdxr?glvub`5VL4e@|PzeZ3fKUvx--y4G
z7C3n#;8?)gd6S~%#4O2Mx*%_I+QM~<7tEU$9}+ZWhWPyWDU-N_%DvBM`muJ2r>lr9
zVfraG&1#kv*H$hO2EF>+sS6H0p8iH=^Z%2_-j;9L{A}i%_giWwzZIDo(0%8s*EaRK
z=_{X?_-uK(?nlViE01SwJ^MDADZS!-#L2eeRGnA*lPyzUeYc$c|K3Ihh5!_=2<`q+
z3G)gc&?_FPCVr*KIjNZb>8F~1)+OAGURWG2v+;4yj}!UV(h}37<u-jg_Rw_O+{jqX
z+3&S($R+>#W4I<|lka+$=MVFx=C3`^7MmKf_f(eb^uroT0*6vdKe(0H$6M`JMf1*r
zi$&{^z2k3aMtko({oQ)+R-RSn@l)IX9{Vemy*vAIZprG)*Kck68ouwc-o)*1e@N8i
zeY(kMc0Xy==JMqRH!`XsUux8C+Z=X!Lb>1Sf^DbPnR4$;c$dYA=A9?&rTSpr;RnW!
zTdFB8|4=Wvty`;KbY66Bd~Uhs^Rt($t*5@ZH&t_~o#Fh<*_$`J_Wa-8u=Z`=h3k>!
zw{yPLKP=fD{Pp~=v^AG&{+=;;y(wv%rM&IzwdZ11)%|&#IUN)~j7+-BsCftER!}zK
zgeEBv)xZcM;W-Ii8z`NEQXohhA4nw>G%%h9vQSekx?WJCz^-=|&{EVSg02@712E_E
f!z=~T4UCyg4CwI@;LXYgQo{m-Z-BHFGl&NO?_85&

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_Honda2004_.npz b/tests/test_fluxes/test_flux_model_Honda2004_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..f4973a90c6db54dd3aab481e51ef170625fbd2ad
GIT binary patch
literal 1354
zcmWIWW@gc4U|`??Vnqf|-Xk^tp@5S?gdrdxr?glvub`5VL4e@|Pzeb1LnwyXZ^U0o
z3!FR=a4cZ$yh%}WVwU7BU6409ZQ;7b3+7FW4+)wwLwtVxlu2Ad<=$sB{aCxi(^bTl
zF#VL8W;M%-Yb%!s!@8M$RX$F(;cvX2dgp#$U$Akt_uV`HYW|7cuHKz{$YS~G8(C|s
zzPPT>?>hUo?AVfyFtOR$iTBgCUz-x~dhXQd9czDdUOhQs_Rsa3oz^e}p!j6MqHjxJ
zKH&rU#3R+juQWL)71KKtsOOz6rR8SRZqC`zJ*DPk(v#%$<-4D6jy&?RcN?qEoYOAv
zR^K@H;DT4e`gw_)lh5Dz|9kt#8x<KE?`rDT9<*Aj{w42L<N2_spDtPVWiE|oMDx$-
zZt-4Z|M(f=@(&fGXO&R!;mc2@-<J8V%`6hSx_P$JCZD}8-zG(Wzh*2p-OTG&onF<)
z-Ai(te0irQJhx*HoBcCx&FMMLtJiL|=u$89{g?6eqx|#c%bPEuN6)6bOFl5~@B^dA
zF~`W9_VLruTwZi!*}CoAd;hHo__lCkY`4cewd0pvw>`gB`61WlRQ~FO(m%)lebJiz
z)aiM#&po@X?}X#)dXJi%h^qLvcz1N}x5&E|Wj`-7W%_EP#m|K~x{}Bqa!WNeqrHd9
zUxuo0dp|GZ`iw`PqMCh=uglC_`;qxK>o)(dyJz>Sao*yZS6LHZkv5t4`Lwy-k@M=d
z{rmj!XSrv8$m{H)B!7FEeQ8VIJb3a@=4m1;&_j$&y3DAB0Vpwn@;xUslYpoOMi2=v
z8ql?Yavms8fwb{~R6;=m<0&8uH6Nnu1*LMBm8g1W0WC#M^XPg(NeX5!Kg?1f-N2a1
l#DJc_&^3di54&b7W(Eed$PDmiWdm8j3WS<KcX6|TcmQR)_YwdA

literal 0
HcmV?d00001

diff --git a/tests/test_fluxes/test_flux_model_TIG1996_.npz b/tests/test_fluxes/test_flux_model_TIG1996_.npz
new file mode 100644
index 0000000000000000000000000000000000000000..c989a4b2d7499aef8b5316bea25b06b44a2ae92c
GIT binary patch
literal 280
zcmWIWW@gc4U|`??Vnv480gopAhXPIp5r%+(oYG>wyn;$b1_6c>KqVkB0YWj%ek1-$
zTHxe~fMWq`=S_;56SE|5>4Ln;X$#jaUNCQ3d`Qrg8RGNfr%d7!D)&C4>Brh7o~|Oc
zgz2Z$G^<%wTwA$B81$siw%?wwa(Clo(cFsdzxMavwVZibCj9yPt$+9WmSj)dJnyaK
zZgv^txU-ul?%!LLe>}=>_Vq4)>zX~A*G#X>iz)ehR&%wdschZ5$y?Sk1b8zt=`y2w
j2IN+dyE&n522l-+AQJAi0B=?{kN_hPngi+6APxfnh*4nk

literal 0
HcmV?d00001

diff --git a/tests/test_genie_datasets.py b/tests/test_genie_datasets.py
index 56f51b5..ba74c64 100755
--- a/tests/test_genie_datasets.py
+++ b/tests/test_genie_datasets.py
@@ -5,81 +5,74 @@
 # SPDX-License-Identifier: BSD-2-Clause
 
 import os
-import unittest
+import sys
 from pathlib import Path
 
 import h5py
 import numpy as np
 import pandas as pd
+import pytest
 import tables
 import uproot
 
 from simweights import GenieWeighter
 from simweights._utils import get_column, get_table
 
+datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
+datasets = [
+    "upgrade_genie_step3_140021_000000",
+    "upgrade_genie_step3_141828_000000",
+    "GENIE_NuMu_IceCubeUpgrade_v58.22590.000000",
+]
+approx = pytest.approx
 
-class TestNugenDatasets(unittest.TestCase):
-    @classmethod
-    def setUpClass(cls):
-        datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
-        if not datadir:
-            cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set")
-        cls.datadir = datadir + "/"
 
-    def cmp_dataset(self, fname):
-        filename = Path(self.datadir) / fname
-        reffile = h5py.File(str(filename) + ".hdf5", "r")
-        wd = reffile["I3MCWeightDict"]
+@pytest.mark.parametrize("fname", datasets)
+@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set")
+def test_dataset(fname):
+    filename = Path(datadir) / fname
+    reffile = h5py.File(str(filename) + ".hdf5", "r")
+    wd = reffile["I3MCWeightDict"]
 
-        solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"]))
-        injection_area = np.pi * (wd["InjectionSurfaceR"] * 1e2) ** 2
-        global_probability_scale = wd["GlobalProbabilityScale"]
-        genie_weight = wd["GENIEWeight"]
+    solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"]))
+    injection_area = np.pi * (wd["InjectionSurfaceR"] * 1e2) ** 2
+    global_probability_scale = wd["GlobalProbabilityScale"]
+    genie_weight = wd["GENIEWeight"]
 
-        pli = -wd["PowerLawIndex"][0]
-        energy_integral = ((10 ** wd["MaxEnergyLog"][0]) ** (pli + 1) - (10 ** wd["MinEnergyLog"][0]) ** (pli + 1)) / (pli + 1)
-        energy_factor = 1 / (wd["PrimaryNeutrinoEnergy"] ** pli / energy_integral)
-        one_weight = global_probability_scale * genie_weight * energy_factor * solid_angle * injection_area
-        np.testing.assert_allclose(one_weight, wd["OneWeight"])
-        final_weight = wd["OneWeight"] / (get_column(get_table(reffile, "I3GenieInfo"), "n_flux_events")[0])
+    pli = -wd["PowerLawIndex"][0]
+    energy_integral = ((10 ** wd["MaxEnergyLog"][0]) ** (pli + 1) - (10 ** wd["MinEnergyLog"][0]) ** (pli + 1)) / (pli + 1)
+    energy_factor = 1 / (wd["PrimaryNeutrinoEnergy"] ** pli / energy_integral)
+    one_weight = global_probability_scale * genie_weight * energy_factor * solid_angle * injection_area
+    np.testing.assert_allclose(one_weight, wd["OneWeight"])
+    final_weight = wd["OneWeight"] / (get_column(get_table(reffile, "I3GenieInfo"), "n_flux_events")[0])
 
-        fobjs = [
-            reffile,
-            uproot.open(str(filename) + ".root"),
-            tables.open_file(str(filename) + ".hdf5", "r"),
-            pd.HDFStore(str(filename) + ".hdf5", "r"),
-        ]
+    fobjs = [
+        reffile,
+        uproot.open(str(filename) + ".root"),
+        tables.open_file(str(filename) + ".hdf5", "r"),
+        pd.HDFStore(str(filename) + ".hdf5", "r"),
+    ]
 
-        for fobj in fobjs:
-            with self.subTest(lib=str(fobj)):
-                w = GenieWeighter(fobj)
+    for fobj in fobjs:
+        w = GenieWeighter(fobj)
 
-                pdf0 = next(iter(w.surface.spectra.values()))[0].dists[0]
-                np.testing.assert_allclose(1 / pdf0.v, global_probability_scale * solid_angle * injection_area, 1e-5)
+        pdf0 = next(iter(w.surface.spectra.values()))[0].dists[0]
+        assert 1 / pdf0.v == approx(global_probability_scale * solid_angle * injection_area, rel=1e-5)
 
-                np.testing.assert_allclose(w.get_weight_column("wght"), genie_weight)
+        assert w.get_weight_column("wght") == approx(genie_weight)
 
-                power_law = next(iter(w.surface.spectra.values()))[0].dists[2]
-                energy_term = 1 / power_law.pdf(w.get_weight_column("energy"))
-                np.testing.assert_allclose(energy_term, energy_factor)
+        power_law = next(iter(w.surface.spectra.values()))[0].dists[2]
+        energy_term = 1 / power_law.pdf(w.get_weight_column("energy"))
+        assert energy_term == approx(energy_factor)
 
-                one_weight = w.get_weight_column("wght") * energy_term * solid_angle * injection_area * global_probability_scale
-                np.testing.assert_allclose(one_weight, wd["OneWeight"], 1e-5)
+        one_weight = w.get_weight_column("wght") * energy_term * solid_angle * injection_area * global_probability_scale
+        assert one_weight == approx(wd["OneWeight"], rel=1e-5)
 
-                np.testing.assert_allclose(w.get_weights(1), final_weight, 1e-5)
+        assert w.get_weights(1) == approx(final_weight, rel=1e-5)
 
-        for fobj in fobjs:
-            fobj.close()
-
-    def test_140021(self):
-        self.cmp_dataset("upgrade_genie_step3_140021_000000")
-
-    def test_141828(self):
-        self.cmp_dataset("upgrade_genie_step3_141828_000000")
-
-    def test_22590(self):
-        self.cmp_dataset("GENIE_NuMu_IceCubeUpgrade_v58.22590.000000")
+    for fobj in fobjs:
+        fobj.close()
 
 
 if __name__ == "__main__":
-    unittest.main()
+    sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]]))
diff --git a/tests/test_icetop_datasets.py b/tests/test_icetop_datasets.py
index 3248e5a..aabfaaf 100755
--- a/tests/test_icetop_datasets.py
+++ b/tests/test_icetop_datasets.py
@@ -5,63 +5,58 @@
 # SPDX-License-Identifier: BSD-2-Clause
 
 import os
-import unittest
+import sys
 from pathlib import Path
 
 import h5py
 import numpy as np
 import pandas as pd
+import pytest
 import tables
 import uproot
 
 from simweights import IceTopWeighter
 
-
-class TestIceTopDatasets(unittest.TestCase):
-    @classmethod
-    def setUpClass(cls):
-        datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
-        if not datadir:
-            cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set")
-        cls.datadir = datadir + "/"
-
-    def cmp_dataset(self, fname):
-        filename = Path(self.datadir) / fname
-        reffile = h5py.File(str(filename) + ".hdf5", "r")
-
-        assert len(reffile["I3TopInjectorInfo"]) == 1
-        si = reffile["I3TopInjectorInfo"][0]
-        pri = reffile["MCPrimary"]
-        solid_angle = np.pi * (np.cos(si["min_zenith"]) ** 2 - np.cos(si["max_zenith"]) ** 2)
-        injection_area = np.pi * (si["sampling_radius"] * 1e2) ** 2
-        energy_integral = np.log(si["max_energy"] / si["min_energy"])  # assuming E^-1
-        energy_factor = energy_integral * pri["energy"]
-        final_weight = energy_factor * solid_angle * injection_area / si["n_events"]
-
-        fobjs = [
-            reffile,
-            uproot.open(str(filename) + ".root"),
-            tables.open_file(str(filename) + ".hdf5", "r"),
-            pd.HDFStore(str(filename) + ".hdf5", "r"),
-        ]
-
-        for fobj in fobjs:
-            with self.subTest(lib=str(fobj)):
-                w = IceTopWeighter(fobj)
-                spatial = w.surface.spectra[2212][0].dists[1]
-                proj_area = spatial.projected_area(1)
-                np.testing.assert_allclose(proj_area, injection_area)
-                sw_etendue = 1 / spatial.pdf(1)
-                np.testing.assert_allclose(sw_etendue, solid_angle * injection_area, 1e-5)
-                np.testing.assert_allclose(energy_integral * w.get_weight_column("energy"), energy_factor)
-                np.testing.assert_allclose(w.get_weights(1), final_weight, 1e-5)
-
-        for fobj in fobjs:
-            fobj.close()
-
-    def test_12360(self):
-        self.cmp_dataset("Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0_0")
+datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
+datasets = ["Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0_0"]
+approx = pytest.approx
+
+
+@pytest.mark.parametrize("fname", datasets)
+@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set")
+def test_dataset(fname):
+    filename = Path(datadir) / fname
+    reffile = h5py.File(str(filename) + ".hdf5", "r")
+
+    assert len(reffile["I3TopInjectorInfo"]) == 1
+    si = reffile["I3TopInjectorInfo"][0]
+    pri = reffile["MCPrimary"]
+    solid_angle = np.pi * (np.cos(si["min_zenith"]) ** 2 - np.cos(si["max_zenith"]) ** 2)
+    injection_area = np.pi * (si["sampling_radius"] * 1e2) ** 2
+    energy_integral = np.log(si["max_energy"] / si["min_energy"])  # assuming E^-1
+    energy_factor = energy_integral * pri["energy"]
+    final_weight = energy_factor * solid_angle * injection_area / si["n_events"]
+
+    fobjs = [
+        reffile,
+        uproot.open(str(filename) + ".root"),
+        tables.open_file(str(filename) + ".hdf5", "r"),
+        pd.HDFStore(str(filename) + ".hdf5", "r"),
+    ]
+
+    for fobj in fobjs:
+        w = IceTopWeighter(fobj)
+        spatial = w.surface.spectra[2212][0].dists[1]
+        proj_area = spatial.projected_area(1)
+        assert proj_area == approx(injection_area)
+        sw_etendue = 1 / spatial.pdf(1)
+        assert sw_etendue == approx(solid_angle * injection_area, 1e-5)
+        assert energy_factor == approx(energy_integral * w.get_weight_column("energy"))
+        assert final_weight == approx(w.get_weights(1), 1e-5)
+
+    for fobj in fobjs:
+        fobj.close()
 
 
 if __name__ == "__main__":
-    unittest.main()
+    sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]]))
diff --git a/tests/test_nugen_datasets.py b/tests/test_nugen_datasets.py
index f65737b..b1e9a50 100755
--- a/tests/test_nugen_datasets.py
+++ b/tests/test_nugen_datasets.py
@@ -5,123 +5,94 @@
 # SPDX-License-Identifier: BSD-2-Clause
 
 import os
-import unittest
+import sys
 from pathlib import Path
 
 import h5py
 import numpy as np
 import pandas as pd
+import pytest
 import tables
 import uproot
 
 from simweights import NuGenWeighter
 
-
-class TestNugenDatasets(unittest.TestCase):
-    @classmethod
-    def setUpClass(cls):
-        datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
-        if not datadir:
-            cls.skipTest(None, "environment variable SIMWEIGHTS_TESTDATA not set")
-        cls.datadir = datadir + "/"
-
-    def cmp_dataset(self, fname):
-        filename = Path(self.datadir) / fname
-        reffile = h5py.File(str(filename) + ".hdf5", "r")
-
-        wd = reffile["I3MCWeightDict"]
-        pdgid = wd["PrimaryNeutrinoType"][0]
-
-        solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"]))
-        if "SolidAngle" in wd.dtype.names:
-            np.testing.assert_allclose(solid_angle, wd["SolidAngle"])
-
-        if "InjectionAreaCGS" in wd.dtype.names:
-            injection_area = wd["InjectionAreaCGS"]
-        if "InjectionAreaNormCGS" in wd.dtype.names:
-            injection_area = wd["InjectionAreaNormCGS"]
-
-        if "TotalWeight" in wd.dtype.names:
-            total_weight = wd["TotalWeight"]
-        elif "TotalInteractionProbabilityWeight" in wd.dtype.names:
-            total_weight = wd["TotalInteractionProbabilityWeight"]
-
-        type_weight = wd["TypeWeight"] if "TypeWeight" in wd.dtype.names else 0.5
-        w0 = wd["OneWeight"] / (wd["NEvents"] * type_weight)
-
-        fobjs = [
-            reffile,
-            uproot.open(str(filename) + ".root"),
-            tables.open_file(str(filename) + ".hdf5", "r"),
-            pd.HDFStore(str(filename) + ".hdf5", "r"),
-        ]
-
-        for fobj in fobjs:
-            with self.subTest(lib=str(fobj)):
-                w = NuGenWeighter(fobj, nfiles=1)
-
-                event_weight = w.get_weight_column("event_weight")
-                np.testing.assert_allclose(event_weight, total_weight)
-
-                cylinder = w.surface.spectra[pdgid][0].dists[2]
-                proj_area = cylinder.projected_area(w.get_weight_column("cos_zen"))
-                np.testing.assert_allclose(proj_area, injection_area)
-
-                sw_etendue = 1 / cylinder.pdf(w.get_weight_column("cos_zen"))
-                np.testing.assert_allclose(sw_etendue, solid_angle * injection_area, 1e-5)
-
-                power_law = w.surface.spectra[pdgid][0].dists[1]
-                energy_factor = 1 / power_law.pdf(w.get_weight_column("energy"))
-                one_weight = w.get_weight_column("event_weight") * energy_factor * solid_angle * injection_area
-                np.testing.assert_allclose(one_weight, wd["OneWeight"])
-
-                np.testing.assert_allclose(w.get_weights(1), w0, 1e-5)
-
-        for fobj in fobjs:
-            fobj.close()
-
-    def test_20885(self):
-        self.cmp_dataset("Level2_IC86.2016_NuE.020885.000000")
-
-    def test_20878(self):
-        self.cmp_dataset("Level2_IC86.2016_NuMu.020878.000000")
-
-    def test_20895(self):
-        self.cmp_dataset("Level2_IC86.2016_NuTau.020895.000000")
-
-    def test_12646(self):
-        self.cmp_dataset("Level2_IC86.2012_nugen_nue.012646.000000.clsim-base-4.0.5.0.99_eff")
-
-    def test_11836(self):
-        self.cmp_dataset("Level2_IC86.2012_nugen_nutau.011836.000000.clsim-base-4.0.3.0.99_eff")
-
-    def test_11477(self):
-        self.cmp_dataset("Level2_IC86.2012_nugen_nutau.011477.000000.clsim-base-4.0.3.0.99_eff")
-
-    def test_11374(self):
-        self.cmp_dataset("Level2_IC86.2012_nugen_numu.011374.000050.clsim-base-4.0.3.0.99_eff")
-
-    def test_11297(self):
-        self.cmp_dataset("Level2_nugen_nutau_IC86.2012.011297.000000")
-
-    def test_11070(self):
-        self.cmp_dataset("Level2_nugen_numu_IC86.2012.011070.000000")
-
-    def test_11069(self):
-        self.cmp_dataset("Level2_nugen_numu_IC86.2012.011069.000000")
-
-    def test_11065(self):
-        self.cmp_dataset("Level2_IC86.2012_nugen_NuTau.011065.000001")
-
-    def test_11029(self):
-        self.cmp_dataset("Level2_nugen_numu_IC86.2012.011029.000000")
-
-    def test_20692(self):
-        self.cmp_dataset("Level2_IC86.2011_nugen_NuE.010692.000000")
-
-    def test_10634(self):
-        self.cmp_dataset("Level2_IC86.2011_nugen_NuMu.010634.000000")
+datasets = [
+    "Level2_IC86.2016_NuE.020885.000000",
+    "Level2_IC86.2016_NuMu.020878.000000",
+    "Level2_IC86.2016_NuTau.020895.000000",
+    "Level2_IC86.2012_nugen_nue.012646.000000.clsim-base-4.0.5.0.99_eff",
+    "Level2_IC86.2012_nugen_nutau.011836.000000.clsim-base-4.0.3.0.99_eff",
+    "Level2_IC86.2012_nugen_nutau.011477.000000.clsim-base-4.0.3.0.99_eff",
+    "Level2_IC86.2012_nugen_numu.011374.000050.clsim-base-4.0.3.0.99_eff",
+    "Level2_nugen_nutau_IC86.2012.011297.000000",
+    "Level2_nugen_numu_IC86.2012.011070.000000",
+    "Level2_nugen_numu_IC86.2012.011069.000000",
+    "Level2_IC86.2012_nugen_NuTau.011065.000001",
+    "Level2_nugen_numu_IC86.2012.011029.000000",
+    "Level2_IC86.2011_nugen_NuE.010692.000000",
+    "Level2_IC86.2011_nugen_NuMu.010634.000000",
+]
+approx = pytest.approx
+datadir = os.environ.get("SIMWEIGHTS_TESTDATA", None)
+
+
+@pytest.mark.parametrize("fname", datasets)
+@pytest.mark.skipif(not datadir, reason="environment variable SIMWEIGHTS_TESTDATA not set")
+def test_dataset(fname):
+    filename = Path(datadir) / fname
+    reffile = h5py.File(str(filename) + ".hdf5", "r")
+
+    wd = reffile["I3MCWeightDict"]
+    pdgid = wd["PrimaryNeutrinoType"][0]
+
+    solid_angle = 2 * np.pi * (np.cos(wd["MinZenith"]) - np.cos(wd["MaxZenith"]))
+    if "SolidAngle" in wd.dtype.names:
+        np.testing.assert_allclose(solid_angle, wd["SolidAngle"])
+
+    if "InjectionAreaCGS" in wd.dtype.names:
+        injection_area = wd["InjectionAreaCGS"]
+    if "InjectionAreaNormCGS" in wd.dtype.names:
+        injection_area = wd["InjectionAreaNormCGS"]
+
+    if "TotalWeight" in wd.dtype.names:
+        total_weight = wd["TotalWeight"]
+    elif "TotalInteractionProbabilityWeight" in wd.dtype.names:
+        total_weight = wd["TotalInteractionProbabilityWeight"]
+
+    type_weight = wd["TypeWeight"] if "TypeWeight" in wd.dtype.names else 0.5
+    w0 = wd["OneWeight"] / (wd["NEvents"] * type_weight)
+
+    fobjs = [
+        reffile,
+        uproot.open(str(filename) + ".root"),
+        tables.open_file(str(filename) + ".hdf5", "r"),
+        pd.HDFStore(str(filename) + ".hdf5", "r"),
+    ]
+
+    for fobj in fobjs:
+        w = NuGenWeighter(fobj, nfiles=1)
+
+        event_weight = w.get_weight_column("event_weight")
+        assert event_weight == approx(total_weight)
+
+        cylinder = w.surface.spectra[pdgid][0].dists[2]
+        proj_area = cylinder.projected_area(w.get_weight_column("cos_zen"))
+        assert proj_area == approx(injection_area)
+
+        sw_etendue = 1 / cylinder.pdf(w.get_weight_column("cos_zen"))
+        assert sw_etendue == approx(solid_angle * injection_area, rel=1e-5)
+
+        power_law = w.surface.spectra[pdgid][0].dists[1]
+        energy_factor = 1 / power_law.pdf(w.get_weight_column("energy"))
+        one_weight = w.get_weight_column("event_weight") * energy_factor * solid_angle * injection_area
+        assert one_weight == approx(wd["OneWeight"])
+
+        assert w0 == approx(w.get_weights(1), rel=1e-5)
+
+    for fobj in fobjs:
+        fobj.close()
 
 
 if __name__ == "__main__":
-    unittest.main()
+    sys.exit(pytest.main(["-v", __file__, *sys.argv[1:]]))
diff --git a/tox.ini b/tox.ini
deleted file mode 100644
index 6942f1a..0000000
--- a/tox.ini
+++ /dev/null
@@ -1,12 +0,0 @@
-; SPDX-FileCopyrightText: © 2022 the SimWeights contributors
-;
-; SPDX-License-Identifier: BSD-2-Clause
-
-[tox]
-envlist = py3{8,9,10,11,12}
-isolated_build = True
-
-[testenv]
-passenv = SIMWEIGHTS_TESTDATA, HDF5_DIR
-deps = .[test]
-commands = pytest