diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index aed6bc445..c4c74804e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -86,19 +86,15 @@ jobs:
       - name: Running mypy
         shell: bash -l {0}
         run: |
-          mypy -p firecrown --ignore-missing-imports
-          mypy -p examples --ignore-missing-imports
-          mypy -p tests --ignore-missing-imports
+          mypy -p firecrown
+          mypy -p examples
+          mypy -p tests
       - name: Running pylint
         shell: bash -l {0}
         run: |
           pylint --rcfile tests/pylintrc  tests
           pylint --rcfile firecrown/models/pylintrc firecrown/models
-          pylint firecrown/connector
-          pylint firecrown/*.py
-          pylint firecrown/likelihood/*.py
-          pylint firecrown/likelihood/gauss_family/*.py
-
+          pylint firecrown
       - name: Running pytest
         shell: bash -l {0}
         run: python -m pytest -vv
diff --git a/docs/contrib.rst b/docs/contrib.rst
index 7db2368d9..52c0c7a35 100644
--- a/docs/contrib.rst
+++ b/docs/contrib.rst
@@ -26,7 +26,7 @@ Please run:
 
 .. code:: bash
 
-    mypy firecrown/ --ignore-missing-imports 
+    mypy -p firecrown -p examples -p tests
 
 and fix any errors reported before pushing commits to the GitHub repository.
 
diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py
index 9b7e521e5..7849acccb 100755
--- a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py
+++ b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py
@@ -1,11 +1,12 @@
 #!/usr/bin/env python
 
 """Example factory function for DES Y1 3x2pt likelihood."""
-
+from dataclasses import dataclass
 import os
 
 from typing import Dict, Union, Tuple
 
+import numpy as np
 import sacc
 import pyccl as ccl
 import pyccl.nl_pt
@@ -26,6 +27,48 @@
 )
 
 
+@dataclass
+class CclSetup:
+    """A package of related CCL parameters, to reduce the number of variables
+    used in the :python:`run_likelihood` method."""
+
+    a_1: float = 1.0
+    a_2: float = 0.5
+    a_d: float = 0.5
+    b_1: float = 2.0
+    b_2: float = 1.0
+    b_s: float = 1.0
+    mag_bias: float = 1.0
+
+
+@dataclass
+class CElls:
+    GG: np.ndarray
+    GI: np.ndarray
+    II: np.ndarray
+    cs_total: np.ndarray
+    gG: np.ndarray
+    gI: np.ndarray
+    mI: np.ndarray
+    gg: np.ndarray
+    gm: np.ndarray
+    gg_total: np.ndarray
+
+    def __init__(self, stat0: TwoPoint, stat2: TwoPoint, stat3: TwoPoint):
+        self.GG = stat0.cells[("shear", "shear")]
+        self.GI = stat0.cells[("shear", "intrinsic_pt")]
+        self.II = stat0.cells[("intrinsic_pt", "intrinsic_pt")]
+        self.cs_total = stat0.cells["total"]
+
+        self.gG = stat2.cells[("galaxies", "shear")]
+        self.gI = stat2.cells[("galaxies", "intrinsic_pt")]
+        self.mI = stat2.cells[("magnification+rsd", "intrinsic_pt")]
+
+        self.gg = stat3.cells[("galaxies", "galaxies")]
+        self.gm = stat3.cells[("galaxies", "magnification+rsd")]
+        self.gg_total = stat3.cells["total"]
+
+
 def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
     """Likelihood factory function for DES Y1 3x2pt analysis."""
 
@@ -115,11 +158,6 @@ def run_likelihood() -> None:
     """Produce some plots using the likelihood function built by
     :python:`build_likelihood`.
     """
-    # We do imports here to save a bit of time when importing this module but
-    # not using the run_likelihood function.
-    # pylint: disable=import-outside-toplevel
-    import numpy as np
-    import matplotlib.pyplot as plt
 
     # pylint: enable=import-outside-toplevel
 
@@ -137,19 +175,9 @@ def run_likelihood() -> None:
     ccl_cosmo = ccl.CosmologyVanillaLCDM()
     ccl_cosmo.compute_nonlin_power()
 
-    # Bare CCL setup
-    a_1 = 1.0
-    a_2 = 0.5
-    a_d = 0.5
-
-    b_1 = 2.0
-    b_2 = 1.0
-    b_s = 1.0
-
-    mag_bias = 1.0
-
+    cs = CclSetup()
     c_1, c_d, c_2 = pyccl.nl_pt.translate_IA_norm(
-        ccl_cosmo, z=z, a1=a_1, a1delta=a_d, a2=a_2, Om_m2_for_c2=False
+        ccl_cosmo, z=z, a1=cs.a_1, a1delta=cs.a_d, a2=cs.a_2, Om_m2_for_c2=False
     )
 
     # Code that creates a Pk2D object:
@@ -165,27 +193,27 @@ def run_likelihood() -> None:
         c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d)
     )
     ptt_m = pyccl.nl_pt.PTMatterTracer()
-    ptt_g = pyccl.nl_pt.PTNumberCountsTracer(b1=b_1, b2=b_2, bs=b_s)
+    ptt_g = pyccl.nl_pt.PTNumberCountsTracer(b1=cs.b_1, b2=cs.b_2, bs=cs.b_s)
     # IA
-    pk_im = ptc.get_biased_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m)
-    pk_ii = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_i, tracer2=ptt_i)
-    pk_gi = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_i)
+    pk_im = ptc.get_biased_pk2d(ptt_i, tracer2=ptt_m)
+    pk_ii = ptc.get_biased_pk2d(ptt_i, tracer2=ptt_i)
+    pk_gi = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_i)
     # Galaxies
-    pk_gm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_m)
-    pk_gg = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_g)
+    pk_gm = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_m)
+    pk_gg = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_g)
     # Magnification: just a matter-matter P(k)
-    pk_mm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_m, tracer2=ptt_m)
+    pk_mm = ptc.get_biased_pk2d(ptt_m, tracer2=ptt_m)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(
         {
-            "ia_a_1": a_1,
-            "ia_a_2": a_2,
-            "ia_a_d": a_d,
-            "lens0_bias": b_1,
-            "lens0_b_2": b_2,
-            "lens0_b_s": b_s,
-            "lens0_mag_bias": mag_bias,
+            "ia_a_1": cs.a_1,
+            "ia_a_2": cs.a_2,
+            "ia_a_d": cs.a_d,
+            "lens0_bias": cs.b_1,
+            "lens0_b_2": cs.b_2,
+            "lens0_b_s": cs.b_s,
+            "lens0_mag_bias": cs.mag_bias,
             "src0_delta_z": 0.000,
             "lens0_delta_z": 0.000,
         }
@@ -202,33 +230,70 @@ def run_likelihood() -> None:
 
     print(f"Log-like = {log_like:.1f}")
 
-    # Plot the predicted and measured statistic
-    # x = likelihood.statistics[0].ell_or_theta_
-    # y_data = likelihood.statistics[0].measured_statistic_
-
     assert isinstance(likelihood, ConstGaussian)
     assert likelihood.cov is not None
 
+    stat0 = likelihood.statistics[0].statistic
+
+    # x = likelihood.statistics[0].ell_or_theta_
+    # y_data = likelihood.statistics[0].measured_statistic_
+
     # y_err = np.sqrt(np.diag(likelihood.cov))[: len(x)]
     # y_theory = likelihood.statistics[0].predicted_statistic_
 
-    print(list(likelihood.statistics[0].cells.keys()))
+    print(list(stat0.cells.keys()))
 
-    ells = likelihood.statistics[0].ells
-    cells_GG = likelihood.statistics[0].cells[("shear", "shear")]
-    cells_GI = likelihood.statistics[0].cells[("shear", "intrinsic_pt")]
-    cells_II = likelihood.statistics[0].cells[("intrinsic_pt", "intrinsic_pt")]
-    cells_cs_total = likelihood.statistics[0].cells["total"]
+    stat2 = likelihood.statistics[2].statistic
+    assert isinstance(stat2, TwoPoint)
+    print(list(stat2.cells.keys()))
 
-    print(list(likelihood.statistics[2].cells.keys()))
-    cells_gG = likelihood.statistics[2].cells[("galaxies", "shear")]
-    cells_gI = likelihood.statistics[2].cells[("galaxies", "intrinsic_pt")]
-    cells_mI = likelihood.statistics[2].cells[("magnification+rsd", "intrinsic_pt")]
+    stat3 = likelihood.statistics[3].statistic
+    print(list(stat3.cells.keys()))
 
-    print(list(likelihood.statistics[3].cells.keys()))
-    cells_gg = likelihood.statistics[3].cells[("galaxies", "galaxies")]
-    cells_gm = likelihood.statistics[3].cells[("galaxies", "magnification+rsd")]
-    cells_gg_total = likelihood.statistics[3].cells["total"]
+    plot_predicted_and_measured_statistics(
+        ccl_cosmo,
+        cs,
+        lens_nz,
+        lens_z,
+        nz,
+        pk_gg,
+        pk_gi,
+        pk_gm,
+        pk_ii,
+        pk_im,
+        pk_mm,
+        stat0,
+        stat2,
+        stat3,
+        z,
+    )
+
+
+def plot_predicted_and_measured_statistics(
+    ccl_cosmo,
+    cs,
+    lens_nz,
+    lens_z,
+    nz,
+    pk_gg,
+    pk_gi,
+    pk_gm,
+    pk_ii,
+    pk_im,
+    pk_mm,
+    stat0,
+    stat2,
+    stat3,
+    z,
+):
+    """Plot the predictions and measurements."""
+    # We do imports here to save a bit of time when importing this module but
+    # not using the run_likelihood function.
+    # pylint: disable=import-outside-toplevel
+    import matplotlib.pyplot as plt
+
+    ells = stat0.ells
+    cells = CElls(stat0, stat2, stat3)
 
     # Code that computes effect from IA using that Pk2D object
     t_lens = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))
@@ -250,13 +315,12 @@ def run_likelihood() -> None:
         has_rsd=True,
         dndz=(lens_z, lens_nz),
         bias=None,
-        mag_bias=(lens_z, mag_bias * np.ones_like(lens_z)),
+        mag_bias=(lens_z, cs.mag_bias * np.ones_like(lens_z)),
     )
     cl_GI = ccl.angular_cl(ccl_cosmo, t_lens, t_ia, ells, p_of_k_a=pk_im)
     cl_II = ccl.angular_cl(ccl_cosmo, t_ia, t_ia, ells, p_of_k_a=pk_ii)
     # The weak gravitational lensing power spectrum
     cl_GG = ccl.angular_cl(ccl_cosmo, t_lens, t_lens, ells)
-
     # Galaxies
     cl_gG = ccl.angular_cl(ccl_cosmo, t_g, t_lens, ells, p_of_k_a=pk_gm)
     cl_gI = ccl.angular_cl(ccl_cosmo, t_g, t_ia, ells, p_of_k_a=pk_gi)
@@ -265,36 +329,32 @@ def run_likelihood() -> None:
     cl_mI = ccl.angular_cl(ccl_cosmo, t_m, t_ia, ells, p_of_k_a=pk_im)
     cl_gm = ccl.angular_cl(ccl_cosmo, t_g, t_m, ells, p_of_k_a=pk_gm)
     cl_mm = ccl.angular_cl(ccl_cosmo, t_m, t_m, ells, p_of_k_a=pk_mm)
-
     # The observed angular power spectrum is the sum of the two.
     cl_cs_theory = cl_GG + 2 * cl_GI + cl_II
     cl_gg_theory = cl_gg + 2 * cl_gm + cl_mm
-
     fig, ax = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
     fig.subplots_adjust(hspace=0)
     # ax[0].plot(x, y_theory, label="Total")
-    ax[0].plot(ells, cells_GG, label="GG firecrown")
+    ax[0].plot(ells, cells.GG, label="GG firecrown")
     ax[0].plot(ells, cl_GG, ls="--", label="GG CCL")
-    ax[0].plot(ells, -cells_GI, label="-GI firecrown")
+    ax[0].plot(ells, -cells.GI, label="-GI firecrown")
     ax[0].plot(ells, -cl_GI, ls="--", label="-GI CCL")
-    ax[0].plot(ells, cells_II, label="II firecrown")
+    ax[0].plot(ells, cells.II, label="II firecrown")
     ax[0].plot(ells, cl_II, ls="--", label="II CCL")
-    ax[0].plot(ells, -cells_gI, label="-Ig firecrown")
+    ax[0].plot(ells, -cells.gI, label="-Ig firecrown")
     ax[0].plot(ells, -cl_gI, ls="--", label="-Ig CCL")
-    ax[0].plot(ells, cells_cs_total, label="total CS firecrown")
+    ax[0].plot(ells, cells.cs_total, label="total CS firecrown")
     ax[0].plot(ells, cl_cs_theory, ls="--", label="total CS CCL")
-
-    ax[1].plot(ells, cells_gG, label="Gg firecrown")
+    ax[1].plot(ells, cells.gG, label="Gg firecrown")
     ax[1].plot(ells, cl_gG, ls="--", label="Gg CCL")
-    ax[1].plot(ells, cells_gg, label="gg firecrown")
+    ax[1].plot(ells, cells.gg, label="gg firecrown")
     ax[1].plot(ells, cl_gg, ls="--", label="gg CCL")
-    ax[1].plot(ells, -cells_mI, label="-mI firecrown")
+    ax[1].plot(ells, -cells.mI, label="-mI firecrown")
     ax[1].plot(ells, -cl_mI, ls="--", label="-mI CCL")
-    ax[1].plot(ells, cells_gm, label="gm firecrown")
+    ax[1].plot(ells, cells.gm, label="gm firecrown")
     ax[1].plot(ells, cl_gm, ls="--", label="gm CCL")
-    ax[1].plot(ells, cells_gg_total, label="total gg firecrown")
+    ax[1].plot(ells, cells.gg_total, label="total gg firecrown")
     ax[1].plot(ells, cl_gg_theory, ls="--", label="total gg CCL")
-
     # ax[0].errorbar(x, y_data, y_err, ls="none", marker="o")
     ax[0].set_xscale("log")
     ax[1].set_xlabel(r"$\ell$")
@@ -303,10 +363,8 @@ def run_likelihood() -> None:
         a.set_yscale("log")
         a.set_ylabel(r"$C_\ell$")
         a.legend(fontsize="small")
-
     fig.suptitle("PT Cls, including IA, galaxy bias, magnification")
     fig.savefig("pt_cls.png", facecolor="white", dpi=300)
-
     plt.show()
 
 
diff --git a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py
index 16544c27f..e918a2676 100644
--- a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py
+++ b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py
@@ -14,7 +14,7 @@
 from firecrown.likelihood.likelihood import Likelihood
 
 
-saccfile = os.path.expanduser(
+SACCFILE = os.path.expanduser(
     os.path.expandvars(
         "${FIRECROWN_DIR}/examples/des_y1_3x2pt/des_y1_3x2pt_sacc_data.fits"
     )
@@ -24,41 +24,10 @@
 def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
     """Build the likelihood for the DES Y1 cosmic shear data TATT."""
     # Load sacc file
-    sacc_data = sacc.Sacc.load_fits(saccfile)
+    sacc_data = sacc.Sacc.load_fits(SACCFILE)
 
-    # Define sources
     n_source = 1
-    sources = {}
-
-    # Define the intrinsic alignment systematic. This will be added to the
-    # lensing sources later
-    ia_systematic = wl.TattAlignmentSystematic()
-
-    for i in range(n_source):
-        # Define the photo-z shift systematic.
-        pzshift = wl.PhotoZShift(sacc_tracer=f"src{i}")
-
-        # Create the weak lensing source, specifying the name of the tracer in the
-        # sacc file and a list of systematics
-        sources[f"src{i}"] = wl.WeakLensing(
-            sacc_tracer=f"src{i}", systematics=[pzshift, ia_systematic]
-        )
-
-    # Define the statistics we like to include in the likelihood
-    stats = {}
-    for stat, sacc_stat in [
-        ("xip", "galaxy_shear_xi_plus"),
-        ("xim", "galaxy_shear_xi_minus"),
-    ]:
-        for i in range(n_source):
-            for j in range(i, n_source):
-                # Define two-point statistics, given two sources (from above) and
-                # the type of statistic.
-                stats[f"{stat}_src{i}_src{j}"] = TwoPoint(
-                    source0=sources[f"src{i}"],
-                    source1=sources[f"src{j}"],
-                    sacc_data_type=sacc_stat,
-                )
+    stats = define_stats(n_source)
 
     # Create the likelihood from the statistics
     pt_calculator = pyccl.nl_pt.EulerianPTCalculator(
@@ -85,19 +54,52 @@ def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
     return likelihood, modeling_tools
 
 
+def define_stats(n_source):
+    """Define the TwoPoint objects to be returned by this factory furnciton."""
+    sources = define_sources(n_source)
+    stats = {}
+    for stat, sacc_stat in [
+        ("xip", "galaxy_shear_xi_plus"),
+        ("xim", "galaxy_shear_xi_minus"),
+    ]:
+        for i in range(n_source):
+            for j in range(i, n_source):
+                # Define two-point statistics, given two sources (from above) and
+                # the type of statistic.
+                stats[f"{stat}_src{i}_src{j}"] = TwoPoint(
+                    source0=sources[f"src{i}"],
+                    source1=sources[f"src{j}"],
+                    sacc_data_type=sacc_stat,
+                )
+    return stats
+
+
+def define_sources(n_source):
+    """Return the sources to be used by the factory function."""
+    result = {}
+    # Define the intrinsic alignment systematic. This will be added to the
+    # lensing restult later
+    ia_systematic = wl.TattAlignmentSystematic()
+    for i in range(n_source):
+        # Define the photo-z shift systematic.
+        pzshift = wl.PhotoZShift(sacc_tracer=f"src{i}")
+
+        # Create the weak lensing source, specifying the name of the tracer in the
+        # sacc file and a list of systematics
+        result[f"src{i}"] = wl.WeakLensing(
+            sacc_tracer=f"src{i}", systematics=[pzshift, ia_systematic]
+        )
+    return result
+
+
 # We can also run the likelihood directly
 def run_likelihood() -> None:
     """Run the likelihood."""
 
-    import numpy as np  # pylint: disable-msg=import-outside-toplevel
-    import matplotlib.pyplot as plt  # pylint: disable-msg=import-outside-toplevel
-
-    like_and_tools = build_likelihood(None)
-    likelihood: Likelihood = like_and_tools[0]
-    tools: ModelingTools = like_and_tools[1]
+    likelihood, tools = build_likelihood(None)
 
     # Load sacc file
-    sacc_data = sacc.Sacc.load_fits(saccfile)
+    sacc_data = sacc.Sacc.load_fits(SACCFILE)
 
     src0_tracer = sacc_data.get_tracer("src0")
     z, nz = src0_tracer.z, src0_tracer.nz
@@ -156,8 +158,7 @@ def run_likelihood() -> None:
     print(f"Log-like = {log_like:.1f}")
 
     # Plot the predicted and measured statistic
-
-    two_point_0 = likelihood.statistics[0]
+    two_point_0 = likelihood.statistics[0].statistic
     assert isinstance(two_point_0, TwoPoint)
 
     # x = two_point_0.ell_or_theta_
@@ -172,13 +173,20 @@ def run_likelihood() -> None:
     # pylint: disable=no-member
     print(list(two_point_0.cells.keys()))
 
+    make_plot(ccl_cosmo, nz, pk_ii, pk_im, two_point_0, z)
+
+
+def make_plot(ccl_cosmo, nz, pk_ii, pk_im, two_point_0, z):
+    """Create and show a diagnostic plot."""
+    import numpy as np  # pylint: disable-msg=import-outside-toplevel
+    import matplotlib.pyplot as plt  # pylint: disable-msg=import-outside-toplevel
+
     ells = two_point_0.ells
     cells_gg = two_point_0.cells[("shear", "shear")]
     cells_gi = two_point_0.cells[("shear", "intrinsic_pt")]
     cells_ii = two_point_0.cells[("intrinsic_pt", "intrinsic_pt")]
     cells_total = two_point_0.cells["total"]
     # pylint: enable=no-member
-
     # Code that computes effect from IA using that Pk2D object
     t_lens = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))
     t_ia = ccl.WeakLensingTracer(
@@ -196,7 +204,6 @@ def run_likelihood() -> None:
     cl_theory = (
         cl_GG + 2 * cl_GI + cl_II
     )  # normally we would also have a third term, +cl_II).
-
     # plt.plot(x, y_theory, label="Total")
     plt.plot(ells, cells_gg, label="GG firecrown")
     plt.plot(ells, cl_GG, ls="--", label="GG CCL")
@@ -206,7 +213,6 @@ def run_likelihood() -> None:
     plt.plot(ells, cl_II, ls="--", label="II CCL")
     plt.plot(ells, cells_total, label="total firecrown")
     plt.plot(ells, cl_theory, ls="--", label="total CCL")
-
     # plt.errorbar(x, y_data, y_err, ls="none", marker="o")
     plt.xscale("log")
     plt.yscale("log")
@@ -217,7 +223,6 @@ def run_likelihood() -> None:
     # plt.ylim(bottom=1e-12)
     plt.title("TATT IA")
     plt.savefig("tatt.png", facecolor="white", dpi=300)
-
     plt.show()
 
 
diff --git a/examples/des_y1_3x2pt/numcosmo_run.py b/examples/des_y1_3x2pt/numcosmo_run.py
index 96755ad43..b4c26b02b 100755
--- a/examples/des_y1_3x2pt/numcosmo_run.py
+++ b/examples/des_y1_3x2pt/numcosmo_run.py
@@ -14,7 +14,7 @@
 
 Ncm.cfg_init()
 
-with open(r"numcosmo_firecrown_model.yml", "r", encoding="utf8") as modelfile:
+with open(r"numcosmo_firecrown_model.yml", "r", encoding="utf-8") as modelfile:
     ncmodel = yaml.load(modelfile, Loader=yaml.Loader)
 
 NcFirecrown = define_numcosmo_model(ncmodel)
diff --git a/examples/des_y1_3x2pt/numcosmo_run_PT.py b/examples/des_y1_3x2pt/numcosmo_run_PT.py
index 3b505e40f..66cea4802 100755
--- a/examples/des_y1_3x2pt/numcosmo_run_PT.py
+++ b/examples/des_y1_3x2pt/numcosmo_run_PT.py
@@ -15,7 +15,7 @@
 
 Ncm.cfg_init()
 
-with open(r"numcosmo_firecrown_model_PT.yml", "r", encoding="utf8") as modelfile:
+with open(r"numcosmo_firecrown_model_PT.yml", "r", encoding="utf-8") as modelfile:
     ncmodel = yaml.load(modelfile, Loader=yaml.Loader)
 
 NcFirecrownPT = define_numcosmo_model(ncmodel)
diff --git a/examples/srd_sn/generate_sn_data.py b/examples/srd_sn/generate_sn_data.py
index f74e3ad72..3bfde3b08 100644
--- a/examples/srd_sn/generate_sn_data.py
+++ b/examples/srd_sn/generate_sn_data.py
@@ -1,6 +1,5 @@
 """Generate SACC data into file srd-y1-converted.sacc.
 """
-from typing import Any, Optional
 import os
 import tarfile
 import urllib.request
@@ -142,17 +141,15 @@ def ensure_local_data_files(dirname_year10, url):
 def main(args):
     """Generate the file srd-y1-converted.sacc."""
     y1cov, y1dat = read_hubble_data(args)
-    #  set up the sacc data name for the astrophysical sources involved.
-    sources = ["supernova"]
-    properties = ["distance"]
-    # The statistc
-    statistic = "mu"
+
     # There is no futher specified needed here - everything is scalar.
-    subtype: Optional[Any] = None
-    sndata_type = sacc.build_data_type_name(sources, properties, statistic, subtype)
+    sndata_type = sacc.build_data_type_name(
+        sources=["supernova"], properties=["distance"], statistic="mu", subtype=None
+    )
     type_details = sacc.parse_data_type_name(sndata_type)
     print(
-        "type_details.sources, type_details.properties, type_details.statistic, type_details.subtype"
+        "type_details.sources, type_details.properties, "
+        "type_details.statistic, type_details.subtype"
     )
     print(
         type_details.sources,
diff --git a/examples/srd_sn/numcosmo_run.py b/examples/srd_sn/numcosmo_run.py
index ca551a2dc..7644d0e44 100755
--- a/examples/srd_sn/numcosmo_run.py
+++ b/examples/srd_sn/numcosmo_run.py
@@ -21,7 +21,7 @@
 
 Ncm.cfg_init()
 
-with open(r"numcosmo_firecrown_model_snia.yml", "r", encoding="utf8") as modelfile:
+with open(r"numcosmo_firecrown_model_snia.yml", "r", encoding="utf-8") as modelfile:
     ncmodel = yaml.load(modelfile, Loader=yaml.Loader)
 
 NcFirecrownSNIa = define_numcosmo_model(ncmodel)
@@ -60,9 +60,7 @@
     require_nonlinear_pk=False, dist=dist, model_list=["NcFirecrownSNIa"]
 )
 
-sacc_file = os.path.expandvars(
-    "${FIRECROWN_DIR}" "/examples/srd_sn/srd-y1-converted.sacc"
-)
+sacc_file = os.path.expandvars("${FIRECROWN_DIR}/examples/srd_sn/srd-y1-converted.sacc")
 
 nc_factory = NumCosmoFactory(
     "sn_srd.py", NamedParameters({"sacc_file": sacc_file}), map_cosmo
diff --git a/firecrown/connector/cosmosis/likelihood.py b/firecrown/connector/cosmosis/likelihood.py
index 6ddb1c8cd..d9ffa0f98 100644
--- a/firecrown/connector/cosmosis/likelihood.py
+++ b/firecrown/connector/cosmosis/likelihood.py
@@ -13,10 +13,13 @@
 from cosmosis.datablock import option_section
 from cosmosis.datablock import names as section_names
 import pyccl as ccl
+
 from firecrown.connector.mapping import mapping_builder, MappingCosmoSIS
+from firecrown.likelihood.gauss_family.gauss_family import GaussFamily
 from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint
 from firecrown.likelihood.likelihood import load_likelihood, Likelihood, NamedParameters
 from firecrown.parameters import ParamsMap
+from firecrown.updatable import MissingSamplerParameterError
 
 
 def extract_section(sample: cosmosis.datablock, section: str) -> NamedParameters:
@@ -68,11 +71,28 @@ def __init__(self, config: cosmosis.datablock):
             print(f"The Firecrown likelihood needs a required parameter: {err}")
             print("*" * 30)
             raise
-
         self.map = mapping_builder(
             input_style="CosmoSIS", require_nonlinear_pk=require_nonlinear_pk
         )
 
+        # If sampling_sections  is empty, but we have required parameters, then
+        # we have a configuration problem, and ParamsMap can never be built
+        # correctly.
+        if len(self.sampling_sections) == 0:
+            required_parameters = self.likelihood.required_parameters()
+            if len(required_parameters) != 0:
+                msg = (
+                    f"The configured likelihood has required "
+                    f"parameters, but CosmoSIS is not providing them.\n"
+                    f"The required parameters are:\n"
+                    f"{list(required_parameters.get_params_names())}\n"
+                    f"You need to provide the names of the DataBlock "
+                    f"sections where these parameters are to be found\n"
+                    f"in the `sampling_parameters_sections` parameter in the "
+                    f"likelihood configuration."
+                )
+                raise RuntimeError(msg)
+
     def execute(self, sample: cosmosis.datablock) -> int:
         """This is the method called for each sample generated by the sampler."""
         cosmological_params: NamedParameters = extract_section(
@@ -90,8 +110,12 @@ def execute(self, sample: cosmosis.datablock) -> int:
         # calculations. e.g., data_vector/firecrown_theory  data_vector/firecrown_data
 
         firecrown_params = self.calculate_firecrown_params(sample)
+        try:
+            self.likelihood.update(firecrown_params)
+        except MissingSamplerParameterError as exc:
+            msg = self.form_error_message(exc)
+            raise RuntimeError(msg) from exc
 
-        self.likelihood.update(firecrown_params)
         self.tools.prepare(ccl_cosmo)
         loglike = self.likelihood.compute_loglike(self.tools)
 
@@ -105,16 +129,19 @@ def execute(self, sample: cosmosis.datablock) -> int:
         self.tools.reset()
 
         # Save concatenated data vector and inverse covariance to enable support
-        # for the CosmoSIS fisher sampler.
-        sample.put(
-            "data_vector", "firecrown_theory", self.likelihood.predicted_data_vector
-        )
-        sample.put(
-            "data_vector", "firecrown_data", self.likelihood.measured_data_vector
-        )
-        sample.put(
-            "data_vector", "firecrown_inverse_covariance", self.likelihood.inv_cov
-        )
+        # for the CosmoSIS Fisher sampler. This can only work for likelihoods
+        # that have these quantities. Currently, this is only GaussFamily.
+
+        if isinstance(self.likelihood, GaussFamily):
+            sample.put(
+                "data_vector", "firecrown_theory", self.likelihood.predicted_data_vector
+            )
+            sample.put(
+                "data_vector", "firecrown_data", self.likelihood.measured_data_vector
+            )
+            sample.put(
+                "data_vector", "firecrown_inverse_covariance", self.likelihood.inv_cov
+            )
 
         # Write out theory and data vectors to the data block the ease
         # debugging.
@@ -145,6 +172,25 @@ def execute(self, sample: cosmosis.datablock) -> int:
 
         return 0
 
+    def form_error_message(self, exc):
+        """Form the error message that will be used to report a missing
+        parameter, when that parameter should have been supplied by the
+        sampler."""
+        msg = (
+            "A required parameter was not found in any of the "
+            "sections searched on DataBlock.\n"
+            "These are specified by the space-separated string "
+            "`sampling_parameter_sections`.\n"
+            "The supplied value was"
+        )
+        sampling_parameters_sections = " ".join(self.sampling_sections)
+        if sampling_parameters_sections:
+            msg += f": `{sampling_parameters_sections}`"
+        else:
+            msg += " an empty string."
+        msg += f"\nThe missing parameter is named: `{exc.parameter}`\n"
+        return msg
+
     def calculate_firecrown_params(self, sample: cosmosis.datablock) -> ParamsMap:
         """Calculate the ParamsMap for this sample."""
 
diff --git a/firecrown/likelihood/gauss_family/gauss_family.py b/firecrown/likelihood/gauss_family/gauss_family.py
index 884739b79..c17862fd1 100644
--- a/firecrown/likelihood/gauss_family/gauss_family.py
+++ b/firecrown/likelihood/gauss_family/gauss_family.py
@@ -21,7 +21,7 @@
 from ..likelihood import Likelihood
 from ...modeling_tools import ModelingTools
 from ...updatable import UpdatableCollection
-from .statistic.statistic import Statistic
+from .statistic.statistic import Statistic, GuardedStatistic
 
 
 class GaussFamily(Likelihood):
@@ -38,7 +38,9 @@ def __init__(
         super().__init__()
         if len(statistics) == 0:
             raise ValueError("GaussFamily requires at least one statistic")
-        self.statistics = UpdatableCollection(statistics)
+        self.statistics: UpdatableCollection = UpdatableCollection(
+            GuardedStatistic(s) for s in statistics
+        )
         self.cov: Optional[npt.NDArray[np.float64]] = None
         self.cholesky: Optional[npt.NDArray[np.float64]] = None
         self.inv_cov: Optional[npt.NDArray[np.float64]] = None
@@ -46,11 +48,18 @@ def __init__(
     def read(self, sacc_data: sacc.Sacc) -> None:
         """Read the covariance matrix for this likelihood from the SACC file."""
 
+        if sacc_data.covariance is None:
+            msg = (
+                f"The {type(self).__name__} likelihood requires a covariance, "
+                f"but the SACC data object being read does not have one."
+            )
+            raise RuntimeError(msg)
+
         covariance = sacc_data.covariance.dense
         for stat in self.statistics:
             stat.read(sacc_data)
 
-        indices_list = [stat.sacc_indices.copy() for stat in self.statistics]
+        indices_list = [s.statistic.sacc_indices.copy() for s in self.statistics]
         indices = np.concatenate(indices_list)
         cov = np.zeros((len(indices), len(indices)))
 
@@ -117,6 +126,8 @@ def compute_chisq(self, tools: ModelingTools) -> float:
             data_vector = self.get_data_vector()
         except NotImplementedError:
             data_vector, theory_vector = self.compute(tools)
+
+        assert len(data_vector) == len(theory_vector)
         residuals = data_vector - theory_vector
 
         self.predicted_data_vector: npt.NDArray[np.float64] = theory_vector
diff --git a/firecrown/likelihood/gauss_family/statistic/cluster_number_counts.py b/firecrown/likelihood/gauss_family/statistic/cluster_number_counts.py
index 36f01fc48..0f3246cbc 100644
--- a/firecrown/likelihood/gauss_family/statistic/cluster_number_counts.py
+++ b/firecrown/likelihood/gauss_family/statistic/cluster_number_counts.py
@@ -172,6 +172,7 @@ def read(self, sacc_data: sacc.Sacc):
 
         self.data_vector = DataVector.from_list(data_vector_list)
         self.sacc_indices = np.array(sacc_indices_list)
+        super().read(sacc_data)
 
     def get_data_vector(self) -> DataVector:
         """Return the data vector; raise exception if there is none."""
diff --git a/firecrown/likelihood/gauss_family/statistic/statistic.py b/firecrown/likelihood/gauss_family/statistic/statistic.py
index aa842c906..111be335b 100644
--- a/firecrown/likelihood/gauss_family/statistic/statistic.py
+++ b/firecrown/likelihood/gauss_family/statistic/statistic.py
@@ -9,16 +9,18 @@
 """
 
 from __future__ import annotations
-from typing import List
+from typing import List, Optional, final
 from dataclasses import dataclass
+from abc import abstractmethod
 import warnings
 import numpy as np
 import numpy.typing as npt
 import sacc
 
-from ....modeling_tools import ModelingTools
-from ....updatable import Updatable
-from .source.source import SourceSystematic
+import firecrown.parameters
+from firecrown.parameters import DerivedParameterCollection, RequiredParameters
+from firecrown.modeling_tools import ModelingTools
+from firecrown.updatable import Updatable
 
 
 class DataVector(npt.NDArray[np.float64]):
@@ -90,31 +92,158 @@ def __iter__(self):
         yield self.theory
 
 
+class StatisticUnreadError(RuntimeError):
+    """Run-time error indicating an attempt has been made to use a statistic
+    that has not had `read` called in it."""
+
+    def __init__(self, stat: Statistic):
+        msg = (
+            f"The statistic {stat} was used for calculation before `read` "
+            f"was called.\nIt may be that a likelihood factory function did not"
+            f"call `read` before returning the likelihood."
+        )
+        super().__init__(msg)
+        self.statstic = stat
+
+
 class Statistic(Updatable):
-    """An abstract statistic class.
+    """The abstract base class for all physics-related statistics.
 
     Statistics read data from a SACC object as part of a multi-phase
     initialization. The manage a :python:`DataVector` and, given a
     :python:`ModelingTools` object, can compute a :python:`TheoryVector`.
 
-    Statistics represent things like two-point functions and mass functinos.
+    Statistics represent things like two-point functions and mass functions.
     """
 
-    systematics: List[SourceSystematic]
-    sacc_indices: npt.NDArray[np.int64]
+    def __init__(self):
+        super().__init__()
+        self.sacc_indices: Optional[npt.NDArray[np.int64]]
+        self.ready = False
+
+    def read(self, _: sacc.Sacc) -> None:
+        """Read the data for this statistic from the SACC data, and mark it
+        as ready for use. Derived classes that override this function
+        should make sure to call the base class method using:
+            super().read(sacc_data)
+        as the last thing they do in `__init__`.
+
+        Note that currently the argument is not used; it is present so that this
+        method will have the correct argument type for the override.
+        """
+        self.ready = True
+        if len(self.get_data_vector()) == 0:
+            raise RuntimeError(
+                f"the statistic {self} has read a data vector "
+                f"of length 0; the length must be positive"
+            )
+
+    @abstractmethod
+    def get_data_vector(self) -> DataVector:
+        """Gets the statistic data vector."""
+
+    @abstractmethod
+    def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector:
+        """Compute a statistic from sources, applying any systematics."""
+
+
+class GuardedStatistic(Updatable):
+    """:python:`GuardedStatistic` is used by the framework to maintain and
+    validate the state of instances of classes derived from
+    :python:`Statistic`."""
+
+    def __init__(self, stat: Statistic):
+        """Initialize the GuardedStatistic to contain the given
+        :python:`Statistic`."""
+        super().__init__()
+        self.statistic = stat
 
     def read(self, sacc_data: sacc.Sacc) -> None:
-        """Read the data for this statistic from the SACC file."""
+        """Read whatever data is needed from the given :python:`sacc.Sacc
+        object.
+
+        After this function is called, the object should be prepared for the
+        calling of the methods :python:`get_data_vector` and
+        :python:`compute_theory_vector`.
+        """
+        if self.statistic.ready:
+            raise RuntimeError("Firecrown has called read twice on a GuardedStatistic")
+        try:
+            self.statistic.read(sacc_data)
+        except TypeError as exc:
+            msg = (
+                f"A statistic of type {type(self.statistic).__name__} has raised "
+                f"an exception during `read`.\nThe problem may be a malformed "
+                f"SACC data object."
+            )
+            raise RuntimeError(msg) from exc
 
     def get_data_vector(self) -> DataVector:
-        """Gets the statistic data vector."""
-        raise NotImplementedError("Method `get_data_vector` is not implemented!")
+        """Return the contained :python:`Statistic`'s data vector.
+
+        :python:`GuardedStatistic` ensures that :python:`read` has been called.
+        first."""
+        if not self.statistic.ready:
+            raise StatisticUnreadError(self.statistic)
+        return self.statistic.get_data_vector()
 
     def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector:
-        """Compute a statistic from sources, applying any systematics."""
-        raise NotImplementedError("Method `compute_theory_vector` is not implemented!")
+        """Return the contained :python:`Statistic`'s computed theory vector.
+
+        :python:`GuardedStatistic` ensures that :python:`read` has been called.
+        first."""
+        if not self.statistic.ready:
+            raise StatisticUnreadError(self.statistic)
+        return self.statistic.compute_theory_vector(tools)
 
-    def compute(self, tools: ModelingTools) -> StatisticsResult:
-        """Compute a statistic from sources, applying any systematics."""
 
-        raise NotImplementedError("Method `compute` is not implemented!")
+class TrivialStatistic(Statistic):
+    """A minimal statistic only to be used for testing Gaussian likelihoods.
+
+    It returns a :python:`DataVector` and :python:`TheoryVector` that is three
+    elements long. The SACC data provided to :python:`TrivialStatistic.read`
+    must supply the necessary values.
+    """
+
+    def __init__(self) -> None:
+        """Initialize this statistic."""
+        super().__init__()
+        # Data and theory will both be of length self.count
+        self.count = 3
+        self.data_vector: Optional[DataVector] = None
+        self.mean = firecrown.parameters.create()
+        self.computed_theory_vector = False
+
+    def read(self, sacc_data: sacc.Sacc):
+        """Read the necessary items from the sacc data."""
+
+        our_data = sacc_data.get_mean(data_type="count")
+        assert len(our_data) == self.count
+        self.data_vector = DataVector.from_list(our_data)
+        self.sacc_indices = np.arange(len(self.data_vector))
+        super().read(sacc_data)
+
+    @final
+    def _reset(self):
+        """Reset this statistic."""
+        self.computed_theory_vector = False
+
+    @final
+    def _required_parameters(self) -> RequiredParameters:
+        """Return an empty RequiredParameters."""
+        return RequiredParameters([])
+
+    @final
+    def _get_derived_parameters(self) -> DerivedParameterCollection:
+        """Return an empty DerivedParameterCollection."""
+        return DerivedParameterCollection([])
+
+    def get_data_vector(self) -> DataVector:
+        """Return the data vector; raise exception if there is none."""
+        assert self.data_vector is not None
+        return self.data_vector
+
+    def compute_theory_vector(self, _: ModelingTools) -> TheoryVector:
+        """Return a fixed theory vector."""
+        self.computed_theory_vector = True
+        return TheoryVector.from_list([self.mean] * self.count)
diff --git a/firecrown/likelihood/gauss_family/statistic/supernova.py b/firecrown/likelihood/gauss_family/statistic/supernova.py
index ebd03e6dc..9fcf06744 100644
--- a/firecrown/likelihood/gauss_family/statistic/supernova.py
+++ b/firecrown/likelihood/gauss_family/statistic/supernova.py
@@ -9,6 +9,7 @@
 
 import pyccl
 import sacc
+from sacc.tracers import MiscTracer
 
 from ....modeling_tools import ModelingTools
 from .statistic import Statistic, DataVector, TheoryVector
@@ -31,6 +32,22 @@ def __init__(self, sacc_tracer) -> None:
     def read(self, sacc_data: sacc.Sacc):
         """Read the data for this statistic from the SACC file."""
 
+        # We do not actually need the tracer, but we want to make sure the SACC
+        # data contains this tracer.
+        # TODO: remove the work-around when the new version of SACC supporting
+        # sacc.Sacc.has_tracer is available.
+        try:
+            tracer = sacc_data.get_tracer(self.sacc_tracer)
+        except KeyError as exc:
+            # Translate to the error type we want
+            raise ValueError(
+                f"The SACC file does not contain the MiscTracer {self.sacc_tracer}"
+            ) from exc
+        if not isinstance(tracer, MiscTracer):
+            raise ValueError(
+                f"The SACC tracer {self.sacc_tracer} is not a " f"MiscTracer"
+            )
+
         data_points = sacc_data.get_data_points(
             data_type="supernova_distance_mu", tracers=(self.sacc_tracer,)
         )
@@ -38,6 +55,7 @@ def read(self, sacc_data: sacc.Sacc):
         self.a = 1.0 / (1.0 + z)
         self.data_vector = DataVector.from_list([dp.value for dp in data_points])
         self.sacc_indices = np.arange(len(self.data_vector))
+        super().read(sacc_data)
 
     def get_data_vector(self) -> DataVector:
         """Return the data vector; raise exception if there is none."""
diff --git a/firecrown/likelihood/gauss_family/statistic/two_point.py b/firecrown/likelihood/gauss_family/statistic/two_point.py
index 1aee616b9..8cc701ea3 100644
--- a/firecrown/likelihood/gauss_family/statistic/two_point.py
+++ b/firecrown/likelihood/gauss_family/statistic/two_point.py
@@ -298,6 +298,8 @@ def read(self, sacc_data: sacc.Sacc) -> None:
         self.measured_statistic_ = self.data_vector
         self.sacc_tracers = tracers
 
+        super().read(sacc_data)
+
     def calculate_ell_or_theta(self):
         """See _ell_for_xi.
 
diff --git a/firecrown/parameters.py b/firecrown/parameters.py
index dbc25201b..a0f30a242 100644
--- a/firecrown/parameters.py
+++ b/firecrown/parameters.py
@@ -86,6 +86,10 @@ def __init__(self, params_names: Iterable[str]):
         """Construct an instance from an Iterable yielding strings."""
         self.params_names: Set[str] = set(params_names)
 
+    def __len__(self):
+        """Return the number of parameters contained."""
+        return len(self.params_names)
+
     def __add__(self, other: RequiredParameters):
         """Return a new RequiredParameters with the concatenated names.
 
diff --git a/firecrown/updatable.py b/firecrown/updatable.py
index 0b50af02a..e117cb1be 100644
--- a/firecrown/updatable.py
+++ b/firecrown/updatable.py
@@ -29,6 +29,23 @@
 GeneralUpdatable = Union["Updatable", "UpdatableCollection"]
 
 
+class MissingSamplerParameterError(RuntimeError):
+    """Error class raised when an Updatable failes to be updated because the
+    ParamsMap supplied for the update is missing a parameter that should have
+    been provided by the sampler."""
+
+    def __init__(self, parameter: str):
+        """Create the error, with a meaning error message."""
+        self.parameter = parameter
+        msg = (
+            f"The parameter `{parameter}` is required to update "
+            "something in this likelihood.\nIt should have been supplied"
+            "by the sampling framework.\nThe object being updated was:\n"
+            "{context}\n"
+        )
+        super().__init__(msg)
+
+
 class Updatable(ABC):
     """Abstract class Updatable is the base class for Updatable objects in Firecrown.
 
@@ -130,11 +147,8 @@ def update(self, params: ParamsMap) -> None:
             try:
                 value = params.get_from_prefix_param(self.sacc_tracer, parameter)
             except KeyError as exc:
-                raise RuntimeError(
-                    f"Missing required parameter "
-                    f"`{parameter_get_full_name(self.sacc_tracer, parameter)}`,"
-                    f" the sampling framework should provide this parameter."
-                    f" The object requiring this parameter is {self}."
+                raise MissingSamplerParameterError(
+                    parameter_get_full_name(self.sacc_tracer, parameter)
                 ) from exc
             setattr(self, parameter, value)
 
diff --git a/pylintrc b/pylintrc
index 4d471c1b9..c7fefa57c 100644
--- a/pylintrc
+++ b/pylintrc
@@ -61,3 +61,7 @@ ignore-none=false
 
 # Expected format of line ending, e.g. empty (any line ending), LF or CRLF.
 expected-line-ending-format=LF
+
+[DESIGN]
+max-args=20
+max-attributes=20
diff --git a/tests/conftest.py b/tests/conftest.py
index df6a415ac..219fc87b8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,26 +1,15 @@
 """
-pytest configuration additions.
+Pytest configuration additions.
+
+Fixtures defined here are available to any test in Firecrown.
 """
-from typing import final, Optional
 
-import numpy as np
 import pytest
 
 import sacc
 
-from firecrown.likelihood.gauss_family.statistic.statistic import (
-    Statistic,
-    DataVector,
-    TheoryVector,
-)
-from firecrown import parameters
-from firecrown.parameters import (
-    RequiredParameters,
-    DerivedParameterCollection,
-    ParamsMap,
-)
-
-from firecrown.modeling_tools import ModelingTools
+from firecrown.likelihood.gauss_family.statistic.statistic import TrivialStatistic
+from firecrown.parameters import ParamsMap
 
 
 def pytest_addoption(parser):
@@ -57,49 +46,6 @@ def pytest_collection_modifyitems(config, items):
 # Fixtures
 
 
-class TrivialStatistic(Statistic):
-    """A minimal statistic for testing Gaussian likelihoods."""
-
-    def __init__(self) -> None:
-        """Initialize this statistic."""
-        super().__init__()
-        self.data_vector: Optional[DataVector] = None
-        self.mean = parameters.create()
-        self.computed_theory_vector = False
-
-    def read(self, sacc_data: sacc.Sacc):
-        """This trivial class does not actually need to read anything."""
-
-        our_data = sacc_data.get_mean(data_type="count")
-        self.data_vector = DataVector.from_list(our_data)
-        self.sacc_indices = np.arange(len(self.data_vector))
-
-    @final
-    def _reset(self):
-        """Reset this statistic. This implementation has nothing to do."""
-        self.computed_theory_vector = False
-
-    @final
-    def _required_parameters(self) -> RequiredParameters:
-        """Return an empty RequiredParameters."""
-        return RequiredParameters([])
-
-    @final
-    def _get_derived_parameters(self) -> DerivedParameterCollection:
-        """Return an empty DerivedParameterCollection."""
-        return DerivedParameterCollection([])
-
-    def get_data_vector(self) -> DataVector:
-        """Return the data vector; raise exception if there is none."""
-        assert self.data_vector is not None
-        return self.data_vector
-
-    def compute_theory_vector(self, _: ModelingTools) -> TheoryVector:
-        """Return a fixed theory vector."""
-        self.computed_theory_vector = True
-        return TheoryVector.from_list([self.mean, self.mean, self.mean])
-
-
 @pytest.fixture(name="trivial_stats")
 def make_stats():
     """Return a non-empty list of TrivialStatistics."""
@@ -112,12 +58,13 @@ def make_trivial_params() -> ParamsMap:
     return ParamsMap({"mean": 1.0})
 
 
-@pytest.fixture(name="sacc_data")
+@pytest.fixture(name="sacc_data_for_trivial_stat")
 def make_sacc_data():
-    """Create a trivial sacc.Sacc object."""
+    """Create a sacc.Sacc object suitable for configuring a
+    :python:`TrivialStatistic`."""
     result = sacc.Sacc()
     result.add_data_point("count", (), 1.0)
     result.add_data_point("count", (), 4.0)
     result.add_data_point("count", (), -3.0)
-    result.add_covariance(np.diag([4.0, 9.0, 16.0]))
+    result.add_covariance([4.0, 9.0, 16.0])
     return result
diff --git a/tests/connector/cosmosis/test_cosmosis_module.py b/tests/connector/cosmosis/test_cosmosis_module.py
index 3aa75022a..cb276f2d2 100644
--- a/tests/connector/cosmosis/test_cosmosis_module.py
+++ b/tests/connector/cosmosis/test_cosmosis_module.py
@@ -3,11 +3,18 @@
 As a unit test, what this can test is very limited.
 This test do not invoke the `cosmosis` executable.
 """
-
-
-from cosmosis.datablock import DataBlock
+from os.path import expandvars
+import yaml
 import pytest
-from firecrown.connector.cosmosis.likelihood import FirecrownLikelihood
+import numpy as np
+from cosmosis.datablock import DataBlock, option_section, names as section_names
+
+from firecrown.likelihood.likelihood import NamedParameters
+from firecrown.connector.cosmosis.likelihood import (
+    FirecrownLikelihood,
+    extract_section,
+    MissingSamplerParameterError,
+)
 
 
 @pytest.fixture(name="minimal_module_config")
@@ -17,7 +24,7 @@ def fixture_minimal_module_config() -> DataBlock:
     This is the minimal possible configuration."""
     block = DataBlock()
     block.put_string(
-        "module_options", "likelihood_source", "tests/likelihood/lkdir/lkscript.py"
+        option_section, "likelihood_source", "tests/likelihood/lkdir/lkscript.py"
     )
     return block
 
@@ -28,13 +35,125 @@ def fixture_defective_module_config() -> DataBlock:
     parameter to configure a ParameterizedLikelihood."""
     block = DataBlock()
     block.put_string(
-        "module_options",
+        option_section,
         "likelihood_source",
-        "tests/likelihood/lkdir/lk_needing_param.py",
+        expandvars("${FIRECROWN_DIR}/tests/likelihood/lkdir/lk_needing_param.py"),
     )
     return block
 
 
+@pytest.fixture(name="minimal_config")
+def fixture_minimal_config() -> DataBlock:
+    result = DataBlock()
+    result.put_string(
+        option_section,
+        "likelihood_source",
+        expandvars("${FIRECROWN_DIR}/tests/likelihood/lkdir/lkscript.py"),
+    )
+    return result
+
+
+@pytest.fixture(name="config_with_derived_parameters")
+def fixture_config_with_derived_parameters() -> DataBlock:
+    result = DataBlock()
+    result.put_string(
+        option_section,
+        "Likelihood_source",
+        expandvars("${FIRECROWN_DIR}/tests/likelihood/lkdir/lk_derived_parameter.py"),
+    )
+    return result
+
+
+@pytest.fixture(name="config_with_const_gaussian_missing_sampling_parameters_sections")
+def fixture_config_with_const_gaussian_missing_sampling_parameters_sections() -> (
+    DataBlock
+):
+    result = DataBlock()
+    result.put_string(
+        option_section,
+        "Likelihood_source",
+        expandvars(
+            "${FIRECROWN_DIR}/tests/likelihood/gauss_family/lkscript_const_gaussian.py"
+        ),
+    )
+    # result.put_string(option_section, )
+    return result
+
+
+@pytest.fixture(name="minimal_firecrown_mod")
+def fixture_minimal_firecrown_mod(minimal_config: DataBlock) -> FirecrownLikelihood:
+    return FirecrownLikelihood(minimal_config)
+
+
+@pytest.fixture(name="firecrown_mod_with_derived_parameters")
+def fixture_firecrown_mod_with_derived_parameters(
+    config_with_derived_parameters: DataBlock,
+) -> FirecrownLikelihood:
+    return FirecrownLikelihood(config_with_derived_parameters)
+
+
+@pytest.fixture(name="firecrown_mod_with_const_gaussian")
+def fixture_firecrown_mod_with_const_gaussian(
+    working_config_for_const_gaussian: DataBlock,
+) -> FirecrownLikelihood:
+    result = FirecrownLikelihood(working_config_for_const_gaussian)
+    return result
+
+
+@pytest.fixture(name="sample_with_cosmo")
+def fixture_sample_with_cosmo() -> DataBlock:
+    """Return a DataBlock that contains some cosmological parameters."""
+    result = DataBlock()
+    params = {
+        "h0": 0.83,
+        "omega_b": 0.22,
+        "omega_c": 0.120,
+        "omega_k": 0.0,
+        "omega_nu": 0.0,
+        "w": -1.0,
+        "wa": 0.0,
+    }
+    for name, val in params.items():
+        result.put_double("cosmological_parameters", name, val)
+    return result
+
+
+@pytest.fixture(name="minimal_sample")
+def fixture_minimal_sample(sample_with_cosmo: DataBlock) -> DataBlock:
+    with open("tests/distances.yml", encoding="utf-8") as stream:
+        rawdata = yaml.load(stream, yaml.CLoader)
+    sample = sample_with_cosmo
+    for section_name, section_data in rawdata.items():
+        for parameter_name, value in section_data.items():
+            sample.put(section_name, parameter_name, np.array(value))
+    return sample
+
+
+@pytest.fixture(name="sample_with_M")
+def fixture_sample_with_M(minimal_sample: DataBlock) -> DataBlock:
+    minimal_sample.put("supernova_parameters", "pantheon_M", 4.5)
+    return minimal_sample
+
+
+@pytest.fixture(name="sample_without_M")
+def fixture_sample_without_M(minimal_sample: DataBlock) -> DataBlock:
+    minimal_sample.put("supernova_parameters", "nobody_loves_me", 4.5)
+    return minimal_sample
+
+
+def test_extract_section_gets_named_parameters(defective_module_config: DataBlock):
+    params = extract_section(defective_module_config, "module_options")
+    assert isinstance(params, NamedParameters)
+    assert params.get_string("likelihood_source") == expandvars(
+        "${FIRECROWN_DIR}/tests/likelihood/lkdir/lk_needing_param.py"
+    )
+
+
+def test_extract_section_raises_on_missing_section(defective_module_config: DataBlock):
+    with pytest.raises(RuntimeError, match="Datablock section `xxx' does not exist"):
+        _ = extract_section(defective_module_config, "xxx")
+
+
 def test_parameterless_module_construction(minimal_module_config):
     """Make sure we can create a CosmoSIS likelihood modules that does not
     introduce any new parameters."""
@@ -47,3 +166,98 @@ def test_missing_required_parameter(defective_module_config):
     failure."""
     with pytest.raises(KeyError):
         _ = FirecrownLikelihood(defective_module_config)
+
+
+def test_initialize_minimal_module(minimal_firecrown_mod: FirecrownLikelihood):
+    assert isinstance(minimal_firecrown_mod, FirecrownLikelihood)
+
+
+def test_execute_missing_cosmological_parameters(
+    minimal_firecrown_mod: FirecrownLikelihood,
+):
+    no_cosmo_params = DataBlock()
+    with pytest.raises(
+        RuntimeError,
+        match="Datablock section `cosmological_parameters' does not exist.",
+    ):
+        _ = minimal_firecrown_mod.execute(no_cosmo_params)
+
+
+def test_execute_with_cosmo(
+    minimal_firecrown_mod: FirecrownLikelihood, minimal_sample: DataBlock
+):
+    assert minimal_firecrown_mod.execute(minimal_sample) == 0
+    assert minimal_sample[section_names.likelihoods, "firecrown_like"] == -3.0
+
+
+def test_execute_with_derived_parameters(
+    firecrown_mod_with_derived_parameters: FirecrownLikelihood,
+    minimal_sample: DataBlock,
+):
+    assert firecrown_mod_with_derived_parameters.execute(minimal_sample) == 0
+    assert minimal_sample.get_double("derived_section", "derived_param0") == 1.0
+
+
+def test_module_init_with_missing_sampling_sections(
+    config_with_const_gaussian_missing_sampling_parameters_sections: DataBlock,
+):
+    with pytest.raises(RuntimeError, match=r"\['pantheon_M'\]"):
+        s = config_with_const_gaussian_missing_sampling_parameters_sections.to_string()
+        assert s is not None
+        _ = FirecrownLikelihood(
+            config_with_const_gaussian_missing_sampling_parameters_sections
+        )
+
+
+@pytest.fixture(name="config_with_wrong_sampling_sections")
+def fixture_config_with_wrong_sampling_sections(
+    config_with_const_gaussian_missing_sampling_parameters_sections: DataBlock,
+) -> DataBlock:
+    # The giant name is good documentation, but I can't type that correctly
+    # twice.
+    config = config_with_const_gaussian_missing_sampling_parameters_sections
+    config[option_section, "sampling_parameters_sections"] = "this_is_wrong"
+    return config
+
+
+@pytest.fixture(name="working_config_for_const_gaussian")
+def fixture_working_config_for_const_gaussian(
+    config_with_const_gaussian_missing_sampling_parameters_sections: DataBlock,
+) -> DataBlock:
+    config = config_with_const_gaussian_missing_sampling_parameters_sections
+    config[option_section, "sampling_parameters_sections"] = "supernova_parameters"
+    return config
+
+
+def test_module_init_with_wrong_sampling_sections(
+    config_with_wrong_sampling_sections: DataBlock,
+):
+    mod = FirecrownLikelihood(config_with_wrong_sampling_sections)
+    assert isinstance(mod, FirecrownLikelihood)
+
+
+def test_module_exec_with_wrong_sampling_sections(
+    config_with_wrong_sampling_sections: DataBlock, sample_with_M: DataBlock
+):
+    mod = FirecrownLikelihood(config_with_wrong_sampling_sections)
+    with pytest.raises(
+        RuntimeError, match="Datablock section `this_is_wrong' does not exist"
+    ):
+        _ = mod.execute(sample_with_M)
+
+
+def test_module_exec_missing_parameter_in_sampling_sections(
+    firecrown_mod_with_const_gaussian: FirecrownLikelihood, sample_without_M: DataBlock
+):
+    with pytest.raises(RuntimeError, match="`supernova_parameters`") as exc:
+        _ = firecrown_mod_with_const_gaussian.execute(sample_without_M)
+    outer_execption = exc.value
+    inner_exception = outer_execption.__cause__
+    assert isinstance(inner_exception, MissingSamplerParameterError)
+    assert inner_exception.parameter == "pantheon_M"
+
+
+def test_module_exec_working(
+    firecrown_mod_with_const_gaussian: FirecrownLikelihood, sample_with_M: DataBlock
+):
+    assert firecrown_mod_with_const_gaussian.execute(sample_with_M) == 0
diff --git a/tests/connector/numcosmo/test_numcosmo_mapping.py b/tests/connector/numcosmo/test_numcosmo_mapping.py
index 97901447f..30a35fde9 100644
--- a/tests/connector/numcosmo/test_numcosmo_mapping.py
+++ b/tests/connector/numcosmo/test_numcosmo_mapping.py
@@ -266,7 +266,11 @@ def test_numcosmo_mapping(numcosmo_cosmo_fixture, request):
     ["numcosmo_cosmo_xcdm", "numcosmo_cosmo_xcdm_no_nu", "numcosmo_cosmo_cpl"],
 )
 def test_numcosmo_data(
-    numcosmo_cosmo_fixture, trivial_stats, sacc_data, nc_model_trivial, request
+    numcosmo_cosmo_fixture,
+    trivial_stats,
+    sacc_data_for_trivial_stat,
+    nc_model_trivial,
+    request,
 ):
     """Test the NumCosmo data connector for NcmData."""
 
@@ -282,7 +286,7 @@ def test_numcosmo_data(
     )
 
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
 
     fc_data = NumCosmoData(
         likelihood=likelihood,
@@ -314,7 +318,11 @@ def test_numcosmo_data(
     ["numcosmo_cosmo_xcdm", "numcosmo_cosmo_xcdm_no_nu", "numcosmo_cosmo_cpl"],
 )
 def test_numcosmo_gauss_cov(
-    numcosmo_cosmo_fixture, trivial_stats, sacc_data, nc_model_trivial, request
+    numcosmo_cosmo_fixture,
+    trivial_stats,
+    sacc_data_for_trivial_stat,
+    nc_model_trivial,
+    request,
 ):
     """Test the NumCosmo data connector for NcmDataGaussCov."""
 
@@ -330,7 +338,7 @@ def test_numcosmo_gauss_cov(
     )
 
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
 
     fc_data = NumCosmoGaussCov(
         likelihood=likelihood,
diff --git a/tests/connector/numcosmo/test_numcosmo_model.py b/tests/connector/numcosmo/test_numcosmo_model.py
index 9fda446cf..fdd8270f0 100644
--- a/tests/connector/numcosmo/test_numcosmo_model.py
+++ b/tests/connector/numcosmo/test_numcosmo_model.py
@@ -42,12 +42,12 @@ def test_model_save_load(tmp_path, nc_model):
     """Test saving and loading a NumCosmo model."""
 
     with open(
-        tmp_path / r"numcosmo_firecrown_model.yml", "w", encoding="utf8"
+        tmp_path / r"numcosmo_firecrown_model.yml", "w", encoding="utf-8"
     ) as modelfile:
         yaml.dump(nc_model, modelfile, Dumper=yaml.Dumper)
 
     with open(
-        tmp_path / r"numcosmo_firecrown_model.yml", "r", encoding="utf8"
+        tmp_path / r"numcosmo_firecrown_model.yml", "r", encoding="utf-8"
     ) as modelfile:
         model_copy = yaml.load(modelfile, Loader=yaml.Loader)
 
diff --git a/tests/distances.yml b/tests/distances.yml
new file mode 100644
index 000000000..7ba923b2c
--- /dev/null
+++ b/tests/distances.yml
@@ -0,0 +1,4 @@
+distances:
+  d_m: [0.0, 30.599016773238837, 61.164969838209146, 91.69730457953247, 122.19546576621251, 152.65889770161618, 183.08704437377514, 213.47934960593844, 243.83525720727596, 274.1542111236613, 304.43565558844614, 334.67903527312603, 364.8837954378314, 395.04938208154766, 425.17524209196944, 455.26082339491893, 485.30557510322814, 515.3089476650088, 545.2703930112107, 575.1893647023988, 605.0653180746485, 634.8977103844911, 664.6860009528034, 694.4296513075827, 724.1281253255095, 753.780889372219, 783.3874124412047, 812.9471662912787, 842.4596255824987, 871.9242680104995, 901.3405744391514, 930.7080290314592, 960.0261193786538, 989.2943366273881, 1018.512175604976, 1047.6791349426137, 1076.7947171965068, 1105.8584289668547, 1134.8697810146282, 1163.8282883760776, 1192.7334704749217, 1221.5848512321663, 1250.381959173494, 1279.124327534182, 1307.811494361507, 1336.4430026145776, 1365.0184002615756, 1393.5372403743443, 1421.9990812203082, 1450.4034863516767, 1478.7500246919053, 1507.0382706193889, 1535.26780404836, 1563.4382105069606, 1591.549081212481, 1619.600013143741, 1647.5906091105878, 1675.52047782052, 1703.389233942409, 1731.1964981673086, 1758.9418972663686, 1786.6250641458269, 1814.2456378990812, 1841.8032638558682, 1869.297593628516, 1896.7282851553073, 1924.095002740947, 1951.3974170941488, 1978.6352053623566, 2005.808051163614, 2032.9156446155973, 2059.957682361837, 2086.9338675951485, 2113.843910078286, 2140.6875261618634, 2167.4644387995504, 2194.1743775605864, 2220.81707863964, 2247.3922848640364, 2273.8997456984016, 2300.339217246755, 2326.7104622520724, 2353.013250093383, 2379.247356780422, 2405.4125649458747, 2431.5086638352805, 2457.535449294613, 2483.492723755584, 2509.3802962187356, 2535.1979822343474, 2560.945603881203, 2586.622989743288, 2612.2299748844357, 2637.766400820986, 2663.2321154925203, 2688.6269732306846, 2713.9508347261917, 2739.203566994024, 2764.385043336902, 2789.495143307066]
+  h: [0.0003299327135551882, 0.00033028642429886265, 0.00033064689311028655, 0.0003310141692726998, 0.0003313883014348929, 0.00033176933760470134, 0.0003321573251426722, 0.00033255231075590697, 0.0003329543404920849, 0.0003333634597336719, 0.000333779713192319, 0.00033420314490345386, 0.00033463379822107017, 0.0003350717158127181, 0.0003355169396546992, 0.00033596951102746967, 0.0003364294705112553, 0.00033689685798188004, 0.00033737171260681264, 0.0003378540728414326, 0.000338343976425519, 0.0003388414603799634, 0.00033934656100371036, 0.00033985931387092586, 0.0003403797538283965, 0.0003409079149931599, 0.00034144383075036957, 0.00034198753375139245, 0.0003425390559121421, 0.00034309842841164797, 0.0003436656816908607, 0.0003442408454516933, 0.0003448239486562996, 0.0003454150195265885, 0.00034601408554397507, 0.00034662117344936585, 0.00034723630924338015, 0.0003478595181868038, 0.0003484908248012773, 0.0003491302528702132, 0.00034977782543994503, 0.0003504335648211034, 0.0003510974925902181, 0.0003517696295915456, 0.00035244999593911667, 0.00035313861101900494, 0.00035383549349181175, 0.0003545406612953651, 0.00035525413164763013, 0.0003559759210498276, 0.0003567060452897582, 0.0003574445194453269, 0.0003581913578882677, 0.00035894657428806097, 0.00035971018161604334, 0.000360482192149704, 0.0003612626174771633, 0.0003620514685018326, 0.00036284875544724657, 0.00036365448786206757, 0.00036446867462525617, 0.0003652913239514028, 0.00036612244339621805, 0.0003669620398621739, 0.00036781011960429515, 0.0003686666882360922, 0.00036953175073563445, 0.00037040531145175553, 0.00037128737411038935, 0.00037217794182102905, 0.0003730770170833059, 0.0003739846017936815, 0.000374900697252251, 0.0003758253041696485, 0.00037675842267405384, 0.0003777000523182926, 0.0003786501920870262, 0.0003796088404040263, 0.00038057599513952853, 0.0003815516536176613, 0.00038253581262394495, 0.0003835284684128543, 0.00038452961671544254, 0.0003855392527470207, 0.00038655737121488646, 0.0003875839663260999, 0.0003886190317953005, 0.0003896625608525598, 0.000390714546251268, 0.0003917749802760464, 0.00039284385475068433, 0.0003939211610460946, 0.00039500689008828416, 0.00039610103236633395, 0.0003972035779403871, 0.00039831451644963803, 0.00039943383712032053, 0.00040056152877369016, 0.0004016975798339982, 0.0004028419783364516]
+  z: [0.0, 0.010101010101010102, 0.020202020202020204, 0.030303030303030304, 0.04040404040404041, 0.05050505050505051, 0.06060606060606061, 0.07070707070707072, 0.08080808080808081, 0.09090909090909091, 0.10101010101010102, 0.11111111111111112, 0.12121212121212122, 0.13131313131313133, 0.14141414141414144, 0.15151515151515152, 0.16161616161616163, 0.17171717171717174, 0.18181818181818182, 0.19191919191919193, 0.20202020202020204, 0.21212121212121213, 0.22222222222222224, 0.23232323232323235, 0.24242424242424243, 0.25252525252525254, 0.26262626262626265, 0.27272727272727276, 0.2828282828282829, 0.29292929292929293, 0.30303030303030304, 0.31313131313131315, 0.32323232323232326, 0.33333333333333337, 0.3434343434343435, 0.3535353535353536, 0.36363636363636365, 0.37373737373737376, 0.38383838383838387, 0.393939393939394, 0.4040404040404041, 0.4141414141414142, 0.42424242424242425, 0.43434343434343436, 0.4444444444444445, 0.4545454545454546, 0.4646464646464647, 0.4747474747474748, 0.48484848484848486, 0.494949494949495, 0.5050505050505051, 0.5151515151515152, 0.5252525252525253, 0.5353535353535354, 0.5454545454545455, 0.5555555555555556, 0.5656565656565657, 0.5757575757575758, 0.5858585858585859, 0.595959595959596, 0.6060606060606061, 0.6161616161616162, 0.6262626262626263, 0.6363636363636365, 0.6464646464646465, 0.6565656565656566, 0.6666666666666667, 0.6767676767676768, 0.686868686868687, 0.696969696969697, 0.7070707070707072, 0.7171717171717172, 0.7272727272727273, 0.7373737373737375, 0.7474747474747475, 0.7575757575757577, 0.7676767676767677, 0.7777777777777778, 0.787878787878788, 0.797979797979798, 0.8080808080808082, 0.8181818181818182, 0.8282828282828284, 0.8383838383838385, 0.8484848484848485, 0.8585858585858587, 0.8686868686868687, 0.8787878787878789, 0.888888888888889, 0.8989898989898991, 0.9090909090909092, 0.9191919191919192, 0.9292929292929294, 0.9393939393939394, 0.9494949494949496, 0.9595959595959597, 0.9696969696969697, 0.9797979797979799, 0.98989898989899, 1.0]
diff --git a/tests/likelihood/gauss_family/lkscript_const_gaussian.py b/tests/likelihood/gauss_family/lkscript_const_gaussian.py
index f02312fb8..99d592f67 100644
--- a/tests/likelihood/gauss_family/lkscript_const_gaussian.py
+++ b/tests/likelihood/gauss_family/lkscript_const_gaussian.py
@@ -3,19 +3,25 @@
 for testing purposes.
 """
 
-import numpy as np
+import sacc
 
 from firecrown.likelihood.gauss_family.gaussian import ConstGaussian
 from firecrown.likelihood.gauss_family.statistic.supernova import Supernova
-from firecrown.likelihood.gauss_family.statistic.statistic import DataVector
 
 
 def build_likelihood(_):
-    """Return an EmptyLikelihood object."""
-    statistic = Supernova(sacc_tracer="no-tracer")
+    """Return a ConstGaussian (likelihood) object."""
+    statistic = Supernova(sacc_tracer="pantheon")
     likelihood = ConstGaussian(statistics=[statistic])
-    likelihood.cov = np.array([1.0])
-    likelihood.cov.shape = (1, 1)
-    statistic.data_vector = DataVector.create(np.array([0.0]))
 
+    # We need the Sacc object to contain everything that the Supernova stat wil
+    # read.
+    sacc_data = sacc.Sacc()
+    tracer_tuple = ("pantheon",)
+    sacc_data.add_tracer("misc", tracer_tuple[0])
+    sacc_data.add_data_point("supernova_distance_mu", tracer_tuple, 1.0, z=0.1)
+    sacc_data.add_data_point("supernova_distance_mu", tracer_tuple, 4.0, z=0.2)
+    sacc_data.add_data_point("supernova_distance_mu", tracer_tuple, -3.0, z=0.3)
+    sacc_data.add_covariance([4.0, 9.0, 16.0])
+    likelihood.read(sacc_data)
     return likelihood
diff --git a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py
index ecea740a5..4b11b3302 100644
--- a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py
+++ b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py
@@ -1,4 +1,4 @@
-"""Tests for ClusterNumberCounts.
+"""Tests for the ClusterNumberCounts statistic.
 """
 import pytest
 
diff --git a/tests/likelihood/gauss_family/statistic/test_statistic.py b/tests/likelihood/gauss_family/statistic/test_statistic.py
index 2030be1e2..e4153c0a7 100644
--- a/tests/likelihood/gauss_family/statistic/test_statistic.py
+++ b/tests/likelihood/gauss_family/statistic/test_statistic.py
@@ -1,9 +1,14 @@
 """
 Tests for the module firecrown.likelihood.gauss_family.statistic.statistic.
 """
+from typing import List
 import numpy as np
+import pytest
+import sacc
+
 import firecrown.likelihood.gauss_family.statistic.statistic as stat
 
+
 VECTOR_CLASSES = (stat.TheoryVector, stat.DataVector)
 
 
@@ -75,3 +80,27 @@ def test_vector_residuals():
     assert isinstance(difference, np.ndarray)
     for cls in VECTOR_CLASSES:
         assert not isinstance(difference, cls)
+
+
+def test_guarded_statistic_read_only_once(
+    sacc_data_for_trivial_stat: sacc.Sacc, trivial_stats: List[stat.TrivialStatistic]
+):
+    gs = stat.GuardedStatistic(trivial_stats.pop())
+    assert not gs.statistic.ready
+    gs.read(sacc_data_for_trivial_stat)
+    assert gs.statistic.ready
+    with pytest.raises(
+        RuntimeError, match="Firecrown has called read twice on a GuardedStatistic"
+    ):
+        gs.read(sacc_data_for_trivial_stat)
+
+
+def test_guarded_statistic_get_data_before_read(trivial_stats):
+    s = trivial_stats.pop()
+    with pytest.raises(
+        stat.StatisticUnreadError,
+        match=f"The statistic {s} was used for "
+        f"calculation before `read` was called.",
+    ):
+        g = stat.GuardedStatistic(s)
+        _ = g.get_data_vector()
diff --git a/tests/likelihood/gauss_family/statistic/test_supernova.py b/tests/likelihood/gauss_family/statistic/test_supernova.py
new file mode 100644
index 000000000..8c14a51a5
--- /dev/null
+++ b/tests/likelihood/gauss_family/statistic/test_supernova.py
@@ -0,0 +1,86 @@
+"""Tests for the Supernova statistic.
+"""
+import pytest
+
+import sacc
+import pyccl
+
+from firecrown.likelihood.gauss_family.statistic.supernova import Supernova
+from firecrown.modeling_tools import ModelingTools
+from firecrown.parameters import ParamsMap
+
+
+@pytest.fixture(name="minimal_stat")
+def fixture_minimal_stat() -> Supernova:
+    """Return a correctly initialized :python:`Supernova` object."""
+    stat = Supernova(sacc_tracer="sn_fake_sample")
+    return stat
+
+
+@pytest.fixture(name="missing_sacc_tracer")
+def fixture_missing_tracer() -> sacc.Sacc:
+    """Return a sacc.Sacc object that lacks a sacc_tracer."""
+    return sacc.Sacc()
+
+
+@pytest.fixture(name="wrong_tracer_type")
+def fixture_wrong_tracer_type() -> sacc.Sacc:
+    data = sacc.Sacc()
+    data.add_tracer("NZ", "sn_fake_sample", 1.0, 5.0)
+    return data
+
+
+@pytest.fixture(name="good_sacc_data")
+def fixture_sacc_data() -> sacc.Sacc:
+    """Return a sacc.Sacc object sufficient to correctly read a
+    :python:`Supernova` object.
+    """
+    data = sacc.Sacc()
+    data.add_tracer("Misc", "sn_fake_sample")
+    # The value of 16.95 supplied here is what will be recovered as entry 0
+    # of the statistic's data vector.
+    data.add_data_point("supernova_distance_mu", ("sn_fake_sample",), 16.95, z=0.0413)
+    # TODO: fill in the right stuff.
+
+    return data
+
+
+def test_missing_sacc_tracer_fails_read(
+    minimal_stat: Supernova, missing_sacc_tracer: sacc.Sacc
+):
+    with pytest.raises(
+        ValueError,
+        match="The SACC file does not contain the MiscTracer sn_fake_sample",
+    ):
+        minimal_stat.read(missing_sacc_tracer)
+
+
+def test_wrong_tracer_type_fails_read(
+    minimal_stat: Supernova, wrong_tracer_type: sacc.Sacc
+):
+    with pytest.raises(
+        ValueError,
+        match=f"The SACC tracer {minimal_stat.sacc_tracer} is not a MiscTracer",
+    ):
+        minimal_stat.read(wrong_tracer_type)
+
+
+def test_read_works(minimal_stat: Supernova, good_sacc_data: sacc.Sacc):
+    """After read() is called, we should be able to get the statistic's
+
+    :python:`DataVector` and also should be able to call
+    :python:`compute_theory_vector`.
+    """
+    minimal_stat.read(good_sacc_data)
+    data_vector = minimal_stat.get_data_vector()
+    assert len(data_vector) == 1
+    assert data_vector[0] == 16.95
+
+    tools = ModelingTools()
+    tools.prepare(pyccl.CosmologyVanillaLCDM())
+    params = ParamsMap({"sn_fake_sample_M": 1.1})
+    minimal_stat.update(params)
+    theory_vector = minimal_stat.compute_theory_vector(tools)
+    assert len(theory_vector) == 1
+    # How do we test the implementation of compute_theory_vector() without
+    # repeating the implementation here?
diff --git a/tests/likelihood/gauss_family/test_const_gaussian.py b/tests/likelihood/gauss_family/test_const_gaussian.py
index 40d860321..e868a572b 100644
--- a/tests/likelihood/gauss_family/test_const_gaussian.py
+++ b/tests/likelihood/gauss_family/test_const_gaussian.py
@@ -1,10 +1,14 @@
 """Unit testsing for ConstGaussian
 """
-
+import pyccl
 import pytest
 import numpy as np
 
+import sacc
+
+import firecrown.parameters
 from firecrown.likelihood.gauss_family.gaussian import ConstGaussian
+from firecrown.likelihood.gauss_family.gauss_family import Statistic
 from firecrown.modeling_tools import ModelingTools
 from firecrown.parameters import (
     RequiredParameters,
@@ -12,6 +16,33 @@
 )
 
 
+@pytest.fixture(name="sacc_with_data_points")
+def fixture_sass_missing_covariance() -> sacc.Sacc:
+    """Return a Sacc object for configuring a ConstGaussian likelihood,
+    but which is missing a covariance matrix."""
+    result = sacc.Sacc()
+    result.add_tracer("misc", "sn_fake_sample")
+    for cnt in [7.0, 4.0]:
+        result.add_data_point("misc", ("sn_fake_sample",), cnt)
+    return result
+
+
+@pytest.fixture(name="sacc_with_covariance")
+def fixture_sacc_with_covariance(sacc_with_data_points: sacc.Sacc) -> sacc.Sacc:
+    result = sacc_with_data_points
+    cov = np.array([[1.0, -0.5], [-0.5, 1.0]])
+    result.add_covariance(cov)
+    return result
+
+
+@pytest.fixture(name="tools_with_vanilla_cosmology")
+def fixture_tools_with_vanilla_cosmology():
+    """Return a ModelingTools object containing the LCDM cosmology from
+    pyccl."""
+    result = ModelingTools()
+    result.prepare(pyccl.CosmologyVanillaLCDM())
+
+
 def test_require_nonempty_statistics():
     with pytest.raises(ValueError):
         _ = ConstGaussian(statistics=[])
@@ -23,41 +54,77 @@ def test_get_cov_fails_before_read(trivial_stats):
         _ = likelihood.get_cov()
 
 
-def test_get_cov_works_after_read(trivial_stats, sacc_data):
+def test_get_cov_works_after_read(trivial_stats, sacc_data_for_trivial_stat):
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
     assert np.all(likelihood.get_cov() == np.diag([4.0, 9.0, 16.0]))
 
 
-def test_chisquared(trivial_stats, sacc_data, trivial_params):
+def test_chisquared(trivial_stats, sacc_data_for_trivial_stat, trivial_params):
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
     likelihood.update(trivial_params)
     assert likelihood.compute_chisq(ModelingTools()) == 2.0
 
 
-def test_required_parameters(trivial_stats, sacc_data, trivial_params):
+def test_required_parameters(trivial_stats, sacc_data_for_trivial_stat, trivial_params):
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
     likelihood.update(trivial_params)
     expected_params = RequiredParameters(params_names=["mean"])
     assert likelihood.required_parameters() == expected_params
 
 
-def test_derived_parameters(trivial_stats, sacc_data, trivial_params):
+def test_derived_parameters(trivial_stats, sacc_data_for_trivial_stat, trivial_params):
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
     likelihood.update(trivial_params)
     expected_params = DerivedParameterCollection([])
     assert likelihood.get_derived_parameters() == expected_params
 
 
-def test_reset(trivial_stats, sacc_data, trivial_params):
+def test_reset(trivial_stats, sacc_data_for_trivial_stat, trivial_params):
     likelihood = ConstGaussian(statistics=trivial_stats)
-    likelihood.read(sacc_data)
+    likelihood.read(sacc_data_for_trivial_stat)
     likelihood.update(trivial_params)
     assert not trivial_stats[0].computed_theory_vector
     assert likelihood.compute_chisq(ModelingTools()) == 2.0
     assert trivial_stats[0].computed_theory_vector
     likelihood.reset()
     assert not trivial_stats[0].computed_theory_vector
+
+
+def test_missing_covariance(trivial_stats, sacc_with_data_points: sacc.Sacc):
+    likelihood = ConstGaussian(statistics=trivial_stats)
+    with pytest.raises(
+        RuntimeError,
+        match="The ConstGaussian likelihood "
+        "requires a covariance, but the "
+        "SACC data object being read does "
+        "not have one.",
+    ):
+        likelihood.read(sacc_with_data_points)
+
+
+def test_using_good_sacc(
+    trivial_stats,
+    sacc_data_for_trivial_stat: sacc.Sacc,
+    tools_with_vanilla_cosmology: ModelingTools,
+):
+    likelihood = ConstGaussian(statistics=trivial_stats)
+    likelihood.read(sacc_data_for_trivial_stat)
+    params = firecrown.parameters.ParamsMap(mean=10.5)
+    likelihood.update(params)
+    chisq = likelihood.compute_chisq(tools_with_vanilla_cosmology)
+    assert isinstance(chisq, float)
+    assert chisq > 0.0
+
+
+def test_after_read_all_statistics_are_ready(
+    trivial_stats, sacc_data_for_trivial_stat: sacc.Sacc
+):
+    likelihood = ConstGaussian(statistics=trivial_stats)
+    likelihood.read(sacc_data_for_trivial_stat)
+    for gs in likelihood.statistics:
+        stat: Statistic = gs.statistic
+        assert stat.ready
diff --git a/tests/test_parameters.py b/tests/test_parameters.py
index 394eaf510..ec1339da0 100644
--- a/tests/test_parameters.py
+++ b/tests/test_parameters.py
@@ -35,6 +35,15 @@ def test_create_with_wrong_arg():
         _ = create("cow")  # type: ignore
 
 
+def test_required_parameters_length():
+    empty = RequiredParameters([])
+    assert len(empty) == 0
+    a = RequiredParameters(["a"])
+    assert len(a) == 1
+    b = RequiredParameters(["a", "b"])
+    assert len(b) == 2
+
+
 def test_get_params_names_does_not_allow_mutation():
     """The caller of RequiredParameters.get_params_names should not be able to modify
     the state of the object on which the call was made."""
diff --git a/tests/test_pt_systematics.py b/tests/test_pt_systematics.py
index a13bfb4f7..d6816a454 100644
--- a/tests/test_pt_systematics.py
+++ b/tests/test_pt_systematics.py
@@ -137,24 +137,28 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data):
 
     # print(list(likelihood.statistics[0].cells.keys()))
     # pylint: disable=no-member
-    ells = likelihood.statistics[0].ells
-    cells_GG = likelihood.statistics[0].cells[("shear", "shear")]
-    cells_GI = likelihood.statistics[0].cells[("shear", "intrinsic_pt")]
-    cells_II = likelihood.statistics[0].cells[("intrinsic_pt", "intrinsic_pt")]
-    cells_cs_total = likelihood.statistics[0].cells["total"]
+
+    # TODO:  We need to have a way to test systematics without requiring
+    #  digging  into the innards of a likelihood object.
+    s0 = likelihood.statistics[0].statistic
+    ells = s0.ells
+    cells_GG = s0.cells[("shear", "shear")]
+    cells_GI = s0.cells[("shear", "intrinsic_pt")]
+    cells_II = s0.cells[("intrinsic_pt", "intrinsic_pt")]
+    cells_cs_total = s0.cells["total"]
 
     # print(list(likelihood.statistics[2].cells.keys()))
-    cells_gG = likelihood.statistics[2].cells[("galaxies", "shear")]
-    cells_gI = likelihood.statistics[2].cells[("galaxies", "intrinsic_pt")]
-    cells_mI = likelihood.statistics[2].cells[("magnification+rsd", "intrinsic_pt")]
+    s2 = likelihood.statistics[2].statistic
+    cells_gG = s2.cells[("galaxies", "shear")]
+    cells_gI = s2.cells[("galaxies", "intrinsic_pt")]
+    cells_mI = s2.cells[("magnification+rsd", "intrinsic_pt")]
 
     # print(list(likelihood.statistics[3].cells.keys()))
-    cells_gg = likelihood.statistics[3].cells[("galaxies", "galaxies")]
-    cells_gm = likelihood.statistics[3].cells[("galaxies", "magnification+rsd")]
-    cells_mm = likelihood.statistics[3].cells[
-        ("magnification+rsd", "magnification+rsd")
-    ]
-    cells_gg_total = likelihood.statistics[3].cells["total"]
+    s3 = likelihood.statistics[3].statistic
+    cells_gg = s3.cells[("galaxies", "galaxies")]
+    cells_gm = s3.cells[("galaxies", "magnification+rsd")]
+    cells_mm = s3.cells[("magnification+rsd", "magnification+rsd")]
+    cells_gg_total = s3.cells["total"]
     # pylint: enable=no-member
     # Code that computes effect from IA using that Pk2D object
     t_lens = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))
@@ -295,13 +299,12 @@ def test_pt_mixed_systematics(sacc_data):
     # print(list(likelihood.statistics[0].cells.keys()))
     # pylint: disable=no-member
 
-    ells = likelihood.statistics[0].ells
+    s0 = likelihood.statistics[0].statistic
+    ells = s0.ells
 
     # print(list(likelihood.statistics[2].cells.keys()))
-    cells_gG = likelihood.statistics[0].cells[("galaxies+magnification+rsd", "shear")]
-    cells_gI = likelihood.statistics[0].cells[
-        ("galaxies+magnification+rsd", "intrinsic_pt")
-    ]
+    cells_gG = s0.cells[("galaxies+magnification+rsd", "shear")]
+    cells_gI = s0.cells[("galaxies+magnification+rsd", "intrinsic_pt")]
     # pylint: enable=no-member
 
     # Code that computes effect from IA using that Pk2D object
diff --git a/tutorial/introduction_to_firecrown.qmd b/tutorial/introduction_to_firecrown.qmd
index ad90fefc1..75b43d6af 100644
--- a/tutorial/introduction_to_firecrown.qmd
+++ b/tutorial/introduction_to_firecrown.qmd
@@ -206,7 +206,7 @@ python setup.py build
 ```{.bash code-copy=true}
 black --check firecrown/ examples/ tests/
 flake8 firecrown examples tests
-mypy --ignore-missing-imports -p firecrown -p examples -p tests
+mypy -p firecrown -p examples -p tests
 pylint --rcfile pylintrc_for_tests tests
 pylint --rcfile firecrown/models/pylintrc firecrown/models
 pylint firecrown/connector