diff --git a/crow/__init__.py b/crow/__init__.py
index 72b4e73..3867c39 100644
--- a/crow/__init__.py
+++ b/crow/__init__.py
@@ -8,4 +8,4 @@
 from .recipes.binned_grid import GridBinnedClusterRecipe
 from .recipes.binned_parent import BinnedClusterRecipe
 
-__version__ = "1.0.0"
+__version__ = "1.0.1"
diff --git a/crow/recipes/binned_exact.py b/crow/recipes/binned_exact.py
index 435e3ec..a4a73b6 100644
--- a/crow/recipes/binned_exact.py
+++ b/crow/recipes/binned_exact.py
@@ -60,64 +60,43 @@ def _setup_with_completeness(self):
             self._completeness_distribution = self.completeness.distribution
 
     def _setup_with_purity(self):
-        """Makes mass distribution use additional integral with completeness"""
+        """Makes mass distribution use additional integral with purity"""
         if self.purity is None:
             self._mass_distribution_distribution = self.mass_distribution.distribution
         else:
-            self._mass_distribution_distribution = self._impure_mass_distribution
+            self._mass_distribution_distribution = self._mass_distribution_purity
 
-    def _impure_mass_distribution(self, log_mass, z, log_mass_proxy_limits):
+    def _mass_distribution_purity(self, log_mass, z, log_mass_proxy):
 
-        ##############################
-        # Fix this function, Henrique
-        # Good luck!
-        ##############################
-
-        integrator = NumCosmoIntegrator(
-            relative_tolerance=1e-6,
-            absolute_tolerance=1e-12,
+        return (
+            self.mass_distribution.gaussian_kernel(log_mass, z, log_mass_proxy)
+            / self.purity.distribution(log_mass_proxy, z)
+            * np.log(10)
         )
 
-        def integration_func(int_args, extra_args):
-            ln_mass_proxy = int_args[:, 0]
-            log_mass_proxy = ln_mass_proxy / np.log(10.0)
-            return np.array(
-                [
-                    self.mass_distribution.gaussian_kernel(
-                        log_mass, z, np.array([_log_mass_proxy])
-                    )
-                    / self.purity.distribution(np.array([_log_mass_proxy]), z)
-                    for _log_mass_proxy in log_mass_proxy
-                ]
-            )
-
-        integrator.integral_bounds = [
-            (
-                np.log(10.0) * log_mass_proxy_limits[0],
-                np.log(10.0) * log_mass_proxy_limits[1],
-            )
-        ]
-
-        return integrator.integrate(integration_func)
-
     def _get_theory_prediction_counts(
         self,
         average_on: None | ClusterProperty = None,
     ) -> Callable[
-        [npt.NDArray[np.float64], npt.NDArray[np.float64], tuple[float, float], float],
+        [
+            npt.NDArray[np.float64],
+            npt.NDArray[np.float64],
+            npt.NDArray[np.float64],
+            float,
+        ],
         npt.NDArray[np.float64],
     ]:
         """Get a callable that evaluates a cluster theory prediction.
 
-        Returns a callable function that accepts mass, redshift, mass proxy limits,
+        Returns a callable function that accepts mass, redshift, mass_proxy,
         and the sky area of your survey and returns the theoretical prediction for the
-        expected number of clusters.
+        expected number of clusters including false detections.
         """
 
         def theory_prediction(
             mass: npt.NDArray[np.float64],
             z: npt.NDArray[np.float64],
-            mass_proxy_limits: tuple[float, float],
+            mass_proxy: npt.NDArray[np.float64],
             sky_area: float,
         ):
             prediction = (
@@ -125,9 +104,18 @@ def theory_prediction(
                 * self.cluster_theory.mass_function(mass, z)
                 * self._completeness_distribution(mass, z)
                 * self.redshift_distribution.distribution()
-                * self._mass_distribution_distribution(mass, z, mass_proxy_limits)
             )
 
+            if self.purity == None:
+                assert (
+                    len(mass_proxy) == 2
+                ), "mass_proxy with no purity should be size 2"
+                prediction *= self._mass_distribution_distribution(
+                    mass, z, (mass_proxy[0], mass_proxy[1])
+                )
+            else:
+                prediction *= self._mass_distribution_distribution(mass, z, mass_proxy)
+
             if average_on is None:
                 return prediction
 
@@ -149,7 +137,7 @@ def _get_function_to_integrate_counts(
             [
                 npt.NDArray[np.float64],
                 npt.NDArray[np.float64],
-                tuple[float, float],
+                npt.NDArray[np.float64],
                 float,
             ],
             npt.NDArray[np.float64],
@@ -167,11 +155,14 @@ def function_mapper(
             mass = int_args[:, 0]
             z = int_args[:, 1]
 
-            mass_proxy_low = extra_args[0]
-            mass_proxy_high = extra_args[1]
-            sky_area = extra_args[2]
+            if self.purity == None:
+                mass_proxy = np.array([extra_args[0], extra_args[1]])
+                sky_area = extra_args[2]
+            else:
+                mass_proxy = int_args[:, 2]
+                sky_area = extra_args[0]
 
-            return prediction(mass, z, (mass_proxy_low, mass_proxy_high), sky_area)
+            return prediction(mass, z, mass_proxy, sky_area)
 
         return function_mapper
 
@@ -188,11 +179,22 @@ def evaluate_theory_prediction_counts(
         using the Murata 2019 binned mass-richness relation and assuming perfectly
         measured redshifts.
         """
-        self.integrator.integral_bounds = [
-            self.mass_interval,
-            z_edges,
-        ]
-        self.integrator.extra_args = np.array([*log_proxy_edges, sky_area])
+        assert len(log_proxy_edges) == 2, "log_proxy_edges should be size 2"
+        assert len(z_edges) == 2, "z_edges should be size 2"
+
+        if self.purity == None:
+            self.integrator.integral_bounds = [
+                self.mass_interval,
+                z_edges,
+            ]
+            self.integrator.extra_args = np.array([*log_proxy_edges, sky_area])
+        else:
+            self.integrator.integral_bounds = [
+                self.mass_interval,
+                z_edges,
+                log_proxy_edges,
+            ]
+            self.integrator.extra_args = np.array([sky_area])
 
         theory_prediction = self._get_theory_prediction_counts(average_on)
         prediction_wrapper = self._get_function_to_integrate_counts(theory_prediction)
@@ -208,7 +210,7 @@ def _get_theory_prediction_shear_profile(
         [
             npt.NDArray[np.float64],
             npt.NDArray[np.float64],
-            tuple[float, float],
+            npt.NDArray[np.float64],
             float,
             float,
         ],
@@ -224,7 +226,7 @@ def _get_theory_prediction_shear_profile(
         def theory_prediction(
             mass: npt.NDArray[np.float64],
             z: npt.NDArray[np.float64],
-            mass_proxy_limits: tuple[float, float],
+            mass_proxy: npt.NDArray[np.float64],
             sky_area: float,
             radius_center: float,
         ):
@@ -233,7 +235,6 @@ def theory_prediction(
                 * self.cluster_theory.mass_function(mass, z)
                 * self.redshift_distribution.distribution()
                 * self._completeness_distribution(mass, z)
-                * self.mass_distribution.distribution(mass, z, mass_proxy_limits)
             )
             if average_on is None:
                 # pylint: disable=no-member
@@ -248,6 +249,18 @@ def theory_prediction(
                     z=z,
                     radius_center=radius_center,
                 )
+                if self.purity == None:
+                    assert (
+                        len(mass_proxy) == 2
+                    ), "mass_proxy with no purity should be size 2"
+                    prediction *= self._mass_distribution_distribution(
+                        mass, z, (mass_proxy[0], mass_proxy[1])
+                    )
+                else:
+                    prediction *= self._mass_distribution_distribution(
+                        mass, z, mass_proxy
+                    )
+
             return prediction
 
         return theory_prediction
@@ -258,7 +271,7 @@ def _get_function_to_integrate_shear_profile(
             [
                 npt.NDArray[np.float64],
                 npt.NDArray[np.float64],
-                tuple[float, float],
+                npt.NDArray[np.float64],
                 float,
                 float,
             ],
@@ -274,16 +287,20 @@ def _get_function_to_integrate_shear_profile(
         def function_mapper(
             int_args: npt.NDArray, extra_args: npt.NDArray
         ) -> npt.NDArray[np.float64]:
+
             mass = int_args[:, 0]
             z = int_args[:, 1]
 
-            mass_proxy_low = extra_args[0]
-            mass_proxy_high = extra_args[1]
-            sky_area = extra_args[2]
-            radius_center = extra_args[3]
-            return prediction(
-                mass, z, (mass_proxy_low, mass_proxy_high), sky_area, radius_center
-            )
+            if self.purity == None:
+                mass_proxy = np.array([extra_args[0], extra_args[1]])
+                sky_area = extra_args[2]
+                radius_center = extra_args[3]
+            else:
+                mass_proxy = int_args[:, 2]
+                sky_area = extra_args[0]
+                radius_center = extra_args[1]
+
+            return prediction(mass, z, mass_proxy, sky_area, radius_center)
 
         return function_mapper
 
@@ -301,15 +318,28 @@ def evaluate_theory_prediction_lensing_profile(
         using the Murata 2019 binned mass-richness relation and assuming perfectly
         measured redshifts.
         """
-        self.integrator.integral_bounds = [
-            self.mass_interval,
-            z_edges,
-        ]
+        assert len(log_proxy_edges) == 2, "log_proxy_edges should be size 2"
+        assert len(z_edges) == 2, "z_edges should be size 2"
+
+        if self.purity == None:
+            self.integrator.integral_bounds = [
+                self.mass_interval,
+                z_edges,
+            ]
+            extra_args = [*log_proxy_edges]
+        else:
+            self.integrator.integral_bounds = [
+                self.mass_interval,
+                z_edges,
+                log_proxy_edges,
+            ]
+            extra_args = []
+
         deltasigma_list = []
+
         for radius_center in radius_centers:
-            self.integrator.extra_args = np.array(
-                [*log_proxy_edges, sky_area, radius_center]
-            )
+            extra_args.extend([sky_area, radius_center])
+            self.integrator.extra_args = np.array(extra_args)
             if self.cluster_theory._beta_parameters is not None:
                 self.cluster_theory.set_beta_s_interp(*z_edges)
             theory_prediction = self._get_theory_prediction_shear_profile(average_on)
diff --git a/crow/recipes/binned_grid.py b/crow/recipes/binned_grid.py
index ee28b6a..20e1ebd 100644
--- a/crow/recipes/binned_grid.py
+++ b/crow/recipes/binned_grid.py
@@ -111,7 +111,6 @@ def _get_hmf_grid(
             vol = self.cluster_theory.comoving_volume(z, sky_area)
             # assign
             self._hmf_grid[key] = vol[:, np.newaxis] * mass_function_2d
-
         return self._hmf_grid[key]
 
     def _get_mass_richness_grid(
@@ -147,7 +146,6 @@ def _get_completeness_grid(self, z: npt.NDArray[np.float64], key):
             self._completeness_grid[key] = self._completeness_distribution(
                 self.log_mass_grid[np.newaxis, :], z[:, np.newaxis]
             )
-
         return self._completeness_grid[key]
 
     def _get_purity_grid(
diff --git a/notebooks/grid_new.ipynb b/notebooks/grid_new.ipynb
index f34fbe9..286ed0b 100644
--- a/notebooks/grid_new.ipynb
+++ b/notebooks/grid_new.ipynb
@@ -65,7 +65,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 7,
    "id": "0b708fef-7cdb-4239-84af-3f1520364076",
    "metadata": {},
    "outputs": [],
@@ -140,6 +140,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "ola\n",
       "Simpson Integral 87221.31147840535, dblquad integral 87220.6489541402\n",
       "Abs error 0.6625242651498411, rel error 7.595956612371779e-06\n"
      ]
@@ -353,8 +354,8 @@
       "\n",
       "\n",
       " (12, 10)\n",
-      "Simpson Integral 0.13497569756262148, dblquad integral 0.1349756810357927\n",
-      "Abs error 1.6526828794383164e-08, rel error 1.224430110635666e-07\n"
+      "Simpson Integral 0.13497569756262148, dblquad integral 0.06724380521649297\n",
+      "Abs error 0.06773189234612852, rel error 1.0072584698034883\n"
      ]
     }
    ],
@@ -505,8 +506,8 @@
       "(0.2, 0.4) (1.0, 1.3)\n",
       "Simpson Integral 465.86054457432704, dblquad integral 465.88615450982513\n",
       "Abs error 0.025609935498096092, rel error 5.4970372590390504e-05\n",
-      "First eval took: 0.0015347003936767578, second eval took: 0.0007061958312988281, integral took: 0.0427401065826416\n",
-      "After reset: 0.0011837482452392578\n",
+      "First eval took: 0.0014820098876953125, second eval took: 0.0008060932159423828, integral took: 0.33889269828796387\n",
+      "After reset: 0.0011088848114013672\n",
       " Counts 1, 2 ,3: (np.float64(465.86054457432704), np.float64(465.86054457432704), np.float64(465.86054457432704))\n"
      ]
     }
@@ -555,8 +556,8 @@
       "[3.67451479e+15]\n",
       "Simpson Integral [3.67407585e+15], dblquad integral [3.67451479e+15]\n",
       "Abs error [4.38941366e+11], rel error [0.00011946]\n",
-      "First eval took: 0.0026946067810058594, second eval took: 0.0008661746978759766, integral took: 0.11332583427429199\n",
-      "After reset: 0.0024514198303222656\n",
+      "First eval took: 0.002621173858642578, second eval took: 0.0007369518280029297, integral took: 0.09407925605773926\n",
+      "After reset: 0.0021927356719970703\n",
       " Counts 1, 2 ,3: (array([3.67407585e+15]), array([3.67407585e+15]), array([3.67407585e+15]))\n"
      ]
     }
@@ -636,8 +637,8 @@
      "text": [
       "Simpson Integral [1.22802065], dblquad integral [1.36087356]\n",
       "Abs error [0.13285291], rel error [0.09762325]\n",
-      "First eval took: 0.04190349578857422, second eval took: 0.0007469654083251953, integral took: 0.28732752799987793\n",
-      "After reset: 0.0007152557373046875\n",
+      "First eval took: 0.02831244468688965, second eval took: 0.0010771751403808594, integral took: 0.22215914726257324\n",
+      "After reset: 0.0007371902465820312\n",
       " Values 1, 2 ,3: (array([1.22802065]), array([1.22802065]), array([1.22802065]))\n"
      ]
     }
@@ -686,8 +687,8 @@
       "(0.2, 0.4) (1.0, 1.3)\n",
       "Simpson Integral 6643.852449550316, dblquad integral 6644.228454470479\n",
       "Abs error 0.3760049201628135, rel error 5.659120885737057e-05\n",
-      "First eval took: 0.0013675689697265625, second eval took: 0.0006885528564453125, integral took: 0.0020668506622314453\n",
-      "After reset: 0.0012717247009277344\n",
+      "First eval took: 0.0015234947204589844, second eval took: 0.0005614757537841797, integral took: 0.002446413040161133\n",
+      "After reset: 0.0008502006530761719\n",
       " Counts 1, 2 ,3: (np.float64(6643.852449550316), np.float64(6643.852449550316), np.float64(6643.852449550316))\n"
      ]
     }
@@ -738,8 +739,8 @@
       "(0.2, 0.4) (1.0, 1.3)\n",
       "Simpson Integral 145.5788864174151, dblquad integral 145.5957901462481\n",
       "Abs error 0.016903728833000287, rel error 0.00011610039559539764\n",
-      "First eval took: 0.0013070106506347656, second eval took: 0.0006926059722900391, integral took: 0.0016515254974365234\n",
-      "After reset: 0.0012030601501464844\n",
+      "First eval took: 0.0014867782592773438, second eval took: 0.0005896091461181641, integral took: 0.0023458003997802734\n",
+      "After reset: 0.0009424686431884766\n",
       " Counts 1, 2 ,3: (np.float64(145.5788864174151), np.float64(145.5788864174151), np.float64(145.5788864174151))\n"
      ]
     }
@@ -841,12 +842,12 @@
      "text": [
       "\n",
       "--- One-shot counts ---\n",
-      "Grid time:    0.0523 s\n",
-      "Integral time: 0.0023 s\n",
+      "Grid time:    0.0346 s\n",
+      "Integral time: 0.0019 s\n",
       "\n",
       "--- One-shot shear ---\n",
-      "Grid time:    0.0021 s\n",
-      "Integral time: 0.4202 s\n"
+      "Grid time:    0.0020 s\n",
+      "Integral time: 0.3056 s\n"
      ]
     }
    ],
@@ -936,8 +937,8 @@
      "output_type": "stream",
      "text": [
       "\n",
-      "TOTAL GRID TIME = 0.23 s\n",
-      "TOTAL INTEGRAL TIME = 24.13 s\n",
+      "TOTAL GRID TIME = 0.16 s\n",
+      "TOTAL INTEGRAL TIME = 17.28 s\n",
       "\n",
       "--- RESULTS COMPARISON ---\n",
       "Total number of bins computed: 20\n",
@@ -1014,7 +1015,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 20,
    "id": "80383885-0945-4e1c-bf2e-ed64072d9130",
    "metadata": {
     "scrolled": true
@@ -1076,7 +1077,16 @@
       "  shear diff  = 0.000007\n",
       "  counts diff = 0.000010\n",
       "\n",
-      "Trying mass_grid_size = 80\n"
+      "Trying mass_grid_size = 80\n",
+      "  shear diff  = 0.000007\n",
+      "  counts diff = 0.000010\n",
+      "\n",
+      "Trying mass_grid_size = 120\n",
+      "  shear diff  = 0.000007\n",
+      "  counts diff = 0.000010\n",
+      "\n",
+      "===== END ACCURACY SWEEP WITH COUNTS =====\n",
+      "\n"
      ]
     }
    ],
@@ -1253,9 +1263,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python (firecrown20.0)",
+   "display_name": "crow",
    "language": "python",
-   "name": "firecrown_20"
+   "name": "crow"
   },
   "language_info": {
    "codemirror_mode": {
@@ -1267,7 +1277,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.0"
+   "version": "3.14.1"
   }
  },
  "nbformat": 4,
diff --git a/tests/test_recipe_binned_counts.py b/tests/test_recipe_binned_counts.py
index 42138b1..641e7f0 100644
--- a/tests/test_recipe_binned_counts.py
+++ b/tests/test_recipe_binned_counts.py
@@ -8,7 +8,7 @@
 import pytest
 from hypothesis import HealthCheck, given, settings
 from hypothesis.strategies import floats
-from scipy.integrate import dblquad, simpson
+from scipy.integrate import dblquad, quad, simpson
 
 from crow import (
     ClusterAbundance,
@@ -127,6 +127,7 @@ def test_binned_grid_init(
 def test_get_theory_prediction_returns_value(
     binned_exact: ExactBinnedClusterRecipe,
 ):
+
     prediction = binned_exact._get_theory_prediction_counts(ClusterProperty.COUNTS)
 
     assert prediction is not None
@@ -350,6 +351,47 @@ def test_evaluates_theory_prediction_with_completeness(
     assert np.abs(prediction_grid_w_comp / prediction_w_comp - 1.0) <= 1.0e-4
 
 
+def test_evaluates_theory_prediction_assertions(
+    binned_exact: ExactBinnedClusterRecipe,
+):
+    mass_proxy_edges = (2, 5)
+    z_edges = (0.5, 1)
+    sky_area = 360**2
+
+    mass_proxy_edges_err = (2, 5, 6)
+    z_edges_err = (0.5, 1, 1.2)
+
+    with pytest.raises(AssertionError):
+        binned_exact.evaluate_theory_prediction_counts(
+            z_edges=z_edges,
+            log_proxy_edges=mass_proxy_edges_err,
+            sky_area=sky_area,
+        )
+
+    with pytest.raises(AssertionError):
+        binned_exact.evaluate_theory_prediction_counts(
+            z_edges=z_edges_err,
+            log_proxy_edges=mass_proxy_edges,
+            sky_area=sky_area,
+        )
+
+    binned_exact.purity = None
+    prediction_exact = binned_exact.evaluate_theory_prediction_counts(
+        z_edges, mass_proxy_edges, sky_area
+    )
+
+    assert len(binned_exact.integrator.extra_args) == 3
+    assert len(binned_exact.integrator.integral_bounds) == 2
+
+    binned_exact.purity = purity_models.PurityAguena16()
+    prediction_exact_w_pur = binned_exact.evaluate_theory_prediction_counts(
+        z_edges, mass_proxy_edges, sky_area
+    )
+
+    assert len(binned_exact.integrator.extra_args) == 1
+    assert len(binned_exact.integrator.integral_bounds) == 3
+
+
 def test_evaluates_theory_prediction_with_purity(
     binned_exact: ExactBinnedClusterRecipe,
     binned_grid: GridBinnedClusterRecipe,
@@ -357,9 +399,17 @@ def test_evaluates_theory_prediction_with_purity(
     mass_proxy_edges = (2, 5)
     z_edges = (0.5, 1)
     sky_area = 360**2
-    ######################################
-    # Henrique, also set up the tests here
-    ######################################
+
+    prediction_exact = binned_exact.evaluate_theory_prediction_counts(
+        z_edges, mass_proxy_edges, sky_area
+    )
+
+    binned_exact.purity = purity_models.PurityAguena16()
+
+    prediction_exact_w_pur = binned_exact.evaluate_theory_prediction_counts(
+        z_edges, mass_proxy_edges, sky_area
+    )
+
     binned_grid_w_pur = get_base_binned_grid(  # Create grid recipe with purity
         None, purity_models.PurityAguena16()
     )
@@ -375,6 +425,10 @@ def test_evaluates_theory_prediction_with_purity(
     )
 
     assert prediction_grid <= prediction_grid_w_pur
+    assert prediction_exact <= prediction_exact_w_pur
+
+    assert np.abs(prediction_grid / prediction_exact - 1.0) <= 1.0e-4
+    assert np.abs(prediction_grid_w_pur / prediction_exact_w_pur - 1.0) <= 1.0e-4
 
 
 @given(
@@ -425,9 +479,21 @@ def test_evaluates_theory_mass_distribution_with_purity(
     binned_exact_w_pur = ExactBinnedClusterRecipe(
         **_kwargs, purity=purity_models.PurityAguena16()
     )
-    probability_w_pur = binned_exact_w_pur._mass_distribution_distribution(
-        mass_array, z_array, mass_proxy_limits
-    )
+
+    def mass_distribtuion_purity_integrand(mass_proxy, mass, z):
+        return binned_exact_w_pur._mass_distribution_distribution(
+            np.array([mass]), np.array([z]), mass_proxy
+        ).item()
+
+    probability_w_pur = [
+        quad(
+            mass_distribtuion_purity_integrand,
+            mass_proxy_limits[0],
+            mass_proxy_limits[1],
+            args=(mass, z),
+        )[0]
+        for mass, z in zip(mass_array, z_array)
+    ]
 
     assert (probability < probability_w_pur).all()
 
@@ -562,7 +628,7 @@ def integrand(ln_proxy_scalar, z_scalar):
         log_proxy_scalar = ln_proxy_scalar / np.log(10.0)
         z_array = np.array([z_scalar])
         log_proxy_array = np.array([log_proxy_scalar])
-        return binned_grid_w_pur._purity_distribution(log_proxy_array, z_array)
+        return binned_grid_w_pur._purity_distribution(log_proxy_array, z_array).item()
 
     z_bin = (z_points[0], z_points[-1])
     proxy_bin = (proxy_grid_size[0], proxy_grid_size[-1])
diff --git a/tests/test_recipe_binned_shear_profile.py b/tests/test_recipe_binned_shear_profile.py
index 7884905..8437f5d 100644
--- a/tests/test_recipe_binned_shear_profile.py
+++ b/tests/test_recipe_binned_shear_profile.py
@@ -227,16 +227,30 @@ def test_evaluates_theory_prediction_returns_value(
 ):
 
     mass_proxy_edges = (2, 5)
+    mass_proxy_edges_err = (2, 5, 6)
     z_edges = (0.5, 1)
+    z_edges_err = (0.5, 1.0, 1.2)
     radius_center = np.atleast_1d(1.5)
+    sky_area = 360**2
     average_on = ClusterProperty.DELTASIGMA
 
+    with pytest.raises(AssertionError):
+        binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges_err, radius_center, sky_area, average_on
+        )
+
+    with pytest.raises(AssertionError):
+        binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
+            z_edges_err, mass_proxy_edges, radius_center, sky_area, average_on
+        )
+
     prediction = binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
-        z_edges, mass_proxy_edges, radius_center, 360**2, average_on
+        z_edges, mass_proxy_edges, radius_center, sky_area, average_on
     )
     prediction_c = binned_exact_deltasigma.evaluate_theory_prediction_counts(
-        z_edges, mass_proxy_edges, 360**2
+        z_edges, mass_proxy_edges, sky_area
     )
+
     assert prediction > 0
     assert prediction_c > 0
 
@@ -252,6 +266,9 @@ def test_grid_shear_matches_exact_within_tolerance(
     sky_area = 360**2
     average_on = ClusterProperty.DELTASIGMA
 
+    binned_exact_deltasigma.completeness = None
+    binned_exact_deltasigma.purity = None
+
     pred_exact = binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
         z_edges, mass_proxy_edges, radii, sky_area, average_on
     )
@@ -259,6 +276,57 @@ def test_grid_shear_matches_exact_within_tolerance(
         z_edges, mass_proxy_edges, radii, sky_area, average_on
     )
 
+    binned_exact_deltasigma.completeness = completeness_models.CompletenessAguena16()
+    binned_exact_deltasigma.purity = None
+    binned_grid_deltasigma_w_c = get_base_binned_grid(
+        completeness_models.CompletenessAguena16(), None, True
+    )
+
+    pred_exact_w_comp = (
+        binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+    pred_grid_w_comp = (
+        binned_grid_deltasigma_w_c.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+
+    binned_exact_deltasigma.completeness = None
+    binned_exact_deltasigma.purity = purity_models.PurityAguena16()
+    binned_grid_deltasigma_w_p = get_base_binned_grid(
+        None, purity_models.PurityAguena16(), True
+    )
+
+    pred_exact_w_pur = (
+        binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+    pred_grid_w_pur = (
+        binned_grid_deltasigma_w_p.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+
+    binned_exact_deltasigma.completeness = completeness_models.CompletenessAguena16()
+    binned_exact_deltasigma.purity = purity_models.PurityAguena16()
+    binned_grid_deltasigma_w_cp = get_base_binned_grid(
+        completeness_models.CompletenessAguena16(), purity_models.PurityAguena16(), True
+    )
+
+    pred_exact_w_cp = (
+        binned_exact_deltasigma.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+    pred_grid_w_cp = (
+        binned_grid_deltasigma_w_cp.evaluate_theory_prediction_lensing_profile(
+            z_edges, mass_proxy_edges, radii, sky_area, average_on
+        )
+    )
+
     # Allow a modest relative tolerance for grid approximation
     rel_tol = 1.0e-4
     assert pred_exact.shape == pred_grid.shape
@@ -266,6 +334,24 @@ def test_grid_shear_matches_exact_within_tolerance(
     denom = np.where(pred_exact == 0, 1.0, pred_exact)
     assert np.all(np.abs((pred_grid - pred_exact) / denom) < rel_tol)
 
+    rel_tol = 1.0e-4
+    assert pred_exact_w_comp.shape == pred_grid_w_comp.shape
+    # avoid division by zero
+    denom = np.where(pred_exact_w_comp == 0, 1.0, pred_grid_w_comp)
+    assert np.all(np.abs((pred_grid_w_comp - pred_exact_w_comp) / denom) < rel_tol)
+
+    rel_tol = 1.0e-4
+    assert pred_exact_w_pur.shape == pred_grid_w_pur.shape
+    # avoid division by zero
+    denom = np.where(pred_exact_w_pur == 0, 1.0, pred_grid_w_pur)
+    assert np.all(np.abs((pred_grid_w_pur - pred_exact_w_pur) / denom) < rel_tol)
+
+    rel_tol = 1.0e-4
+    assert pred_exact_w_cp.shape == pred_grid_w_cp.shape
+    # avoid division by zero
+    denom = np.where(pred_exact_w_cp == 0, 1.0, pred_grid_w_cp)
+    assert np.all(np.abs((pred_grid_w_cp - pred_exact_w_cp) / denom) < rel_tol)
+
 
 def test_grid_reduced_shear_matches_exact_within_tolerance(
     binned_exact_gt: ExactBinnedClusterRecipe,
@@ -401,5 +487,6 @@ def test_shear_respects_completeness_and_purity_effects(
 
     # purity model usually acts to reduce contamination; behavior depends on model details,
     # but it should produce finite, non-negative numbers (sanity)
+
     assert np.isfinite(pur_val)
     assert pur_val >= 0.0