diff --git a/clmm/__init__.py b/clmm/__init__.py
index 4071ee026..9689bbd54 100644
--- a/clmm/__init__.py
+++ b/clmm/__init__.py
@@ -26,4 +26,4 @@
 )
 from . import support
 
-__version__ = "1.9.0"
+__version__ = "1.10.0"
diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py
index 4a2e3de44..33adb2dd9 100644
--- a/clmm/dataops/__init__.py
+++ b/clmm/dataops/__init__.py
@@ -13,6 +13,7 @@
     validate_argument,
     _validate_ra,
     _validate_dec,
+    _validate_is_deltasigma_sigma_c,
 )
 from ..redshift import (
     _integ_pzfuncs,
@@ -30,13 +31,7 @@ def compute_tangential_and_cross_components(
     shear2,
     geometry="curve",
     is_deltasigma=False,
-    cosmo=None,
-    z_lens=None,
-    z_source=None,
     sigma_c=None,
-    use_pdz=False,
-    pzbins=None,
-    pzpdf=None,
     validate_input=True,
 ):
     r"""Computes tangential- and cross- components for shear or ellipticity
@@ -76,17 +71,15 @@ def compute_tangential_and_cross_components(
         g_t =& -\left( g_1\cos\left(2\phi\right)+g_2\sin\left(2\phi\right)\right)\\
         g_x =& g_1 \sin\left(2\phi\right)-g_2\cos\left(2\phi\right)
 
-    Finally, and if requested by the user throught the `is_deltasigma` flag, an estimate of the
-    excess surface density :math:`\widehat{\Delta\Sigma}` is obtained from
+    Finally, if  the critical surface density (:math:`\Sigma_\text{crit}`) is provided, an estimate
+    of the excess surface density :math:`\widehat{\Delta\Sigma}` is obtained from
 
     .. math::
-        \widehat{\Delta\Sigma_{t,x}} = g_{t,x} \times \Sigma_c(cosmo, z_l, z_{\text{src}})
+        \widehat{\Delta\Sigma_{t,x}} = g_{t,x} \times \Sigma_\text{crit}(z_l, z_{\text{src}})
 
-    where :math:`\Sigma_c` is the critical surface density that depends on the cosmology and on
-    the lens and source redshifts. If :math:`g_{t,x}` correspond to the shear, the above
-    expression is an accurate. However, if :math:`g_{t,x}` correspond to ellipticities or reduced
-    shear, this expression only gives an estimate :math:`\widehat{\Delta\Sigma_{t,x}}`, valid only
-    in the weak lensing regime.
+    If :math:`g_{t,x}` correspond to the shear, the above expression is an accurate. However, if
+    :math:`g_{t,x}` correspond to ellipticities or reduced shear, this expression only gives an
+    estimate :math:`\widehat{\Delta\Sigma_{t,x}}`, valid only in the weak lensing regime.
 
     Parameters
     ----------
@@ -107,28 +100,10 @@ def compute_tangential_and_cross_components(
         Options are curve (uses astropy) or flat.
     is_deltasigma: bool
         If `True`, the tangential and cross components returned are multiplied by Sigma_crit.
-        Results in units of :math:`M_\odot\ Mpc^{-2}`
-    cosmo: clmm.Cosmology, optional
-        Required if `is_deltasigma` is True and `sigma_c` not provided.
-        Not used if `sigma_c` is provided.
-    z_lens: float, optional
-        Redshift of the lens, required if `is_deltasigma` is True and `sigma_c` not provided.
-        Not used if `sigma_c` is provided.
-    z_source: array, optional
-        Redshift of the source, required if `is_deltasigma` is True and `sigma_c` not provided.
-        Not used if `sigma_c` is provided or `use_pdz=True`.
-    sigma_c : float, optional
-        Critical surface density in units of :math:`M_\odot\ Mpc^{-2}`,
-        if provided, `cosmo`, `z_lens` and `z_source` are not used.
-    use_pdz: bool
-        Flag to use or not the source redshift p(z), required if `is_deltasigma` is True
-        and `sigma_c` not provided. If `False`, the point estimate provided by `z_source` is used.
-    pzpdf : array, optional
-        Photometric probablility density functions of the source galaxies, required if
-        `is_deltasigma=True` and `use_pdz=True` and `sigma_c` not provided.
-    pzbins : array, optional
-        Redshift axis on which the individual photoz pdf is tabulated, required if
-        `is_deltasigma=True` and `use_pdz=True` and `sigma_c` not provided.
+        It requires `sigma_c` argument. Results in units of :math:`M_\odot\ Mpc^{-2}`
+    sigma_c : None, array_like
+        Critical (effective) surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+        Used only when is_deltasigma=True.
     validate_input: bool
         Validade each input argument
 
@@ -152,15 +127,13 @@ def compute_tangential_and_cross_components(
         validate_argument(locals(), "shear1", "float_array")
         validate_argument(locals(), "shear2", "float_array")
         validate_argument(locals(), "geometry", str)
-        validate_argument(locals(), "is_deltasigma", bool)
-        validate_argument(locals(), "z_lens", float, argmin=0, eqmin=True, none_ok=True)
-        validate_argument(locals(), "z_source", "float_array", argmin=0, eqmin=True, none_ok=True)
         validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
         ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency(
             [ra_source, dec_source, shear1, shear2],
             names=("Ra", "Dec", "Shear1", "Shear2"),
             prefix="Tangential- and Cross- shape components sources",
         )
+        _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c)
     elif np.iterable(ra_source):
         ra_source_, dec_source_, shear1_, shear2_ = (
             np.array(col) for col in [ra_source, dec_source, shear1, shear2]
@@ -182,58 +155,33 @@ def compute_tangential_and_cross_components(
     # Compute the tangential and cross shears
     tangential_comp = _compute_tangential_shear(shear1_, shear2_, phi)
     cross_comp = _compute_cross_shear(shear1_, shear2_, phi)
-    # If the is_deltasigma flag is True, multiply the results by Sigma_crit.
-
-    if is_deltasigma:
-        if sigma_c is None and use_pdz is False:
-            # Need to verify that cosmology and redshifts are provided
-            if any(t_ is None for t_ in (z_lens, z_source, cosmo)):
-                raise TypeError(
-                    "To compute DeltaSigma, please provide a "
-                    "i) cosmology, ii) redshift of lens and sources"
-                )
-
-            sigma_c = cosmo.eval_sigma_crit(z_lens, z_source)
-
-        elif sigma_c is None:
-            # Need to verify that cosmology, lens redshift, source redshift bins and
-            # source redshift pdf are provided
-            if any(t_ is None for t_ in (z_lens, cosmo, pzbins, pzpdf)):
-                raise TypeError(
-                    "To compute DeltaSigma using the redshift pdz of the sources, "
-                    "please provide a "
-                    "i) cosmology, ii) lens redshift, iii) source redshift bins and"
-                    "iv) source redshift pdf"
-                )
-
-            sigma_c = compute_critical_surface_density_eff(
-                cosmo, z_lens, pzbins=pzbins, pzpdf=pzpdf
-            )
 
-        tangential_comp *= sigma_c
-        cross_comp *= sigma_c
+    if sigma_c is not None:
+        _sigma_c_arr = np.array(sigma_c)
+        tangential_comp *= _sigma_c_arr
+        cross_comp *= _sigma_c_arr
+
     return angsep, tangential_comp, cross_comp
 
 
 def compute_background_probability(
-    z_lens, z_source=None, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True
+    z_lens, z_src=None, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True
 ):
-    r"""Probability for being a background galaxy
+    r"""Probability for being a background galaxy, defined by:
+
+        .. math::
+            P(z_s > z_l) = \int_{z_l}^{+\infty} dz_s \; p_{\text{photoz}}(z_s),
+
+    when the photometric probability density functions (:math:`p_{\text{photoz}}(z_s)`) are
+    provided. In the case of true redshifts, it returns 1 if :math:`z_s > z_l` else returns 0.
+
 
     Parameters
     ----------
     z_lens: float
         Redshift of the lens.
-    z_source: array, optional
+    z_src: array, optional
         Redshift of the source. Used only if pzpdf=pzbins=None.
-    use_pdz: bool
-        Flag to use or not the source redshif. If `False`,
-        the point estimate provided by `z_source` is used.
-    pzpdf : array, optional
-        Photometric probablility density functions of the source galaxies.
-        Used instead of z_source if provided.
-    pzbins : array, optional
-        Redshift axis on which the individual photoz pdf is tabulated.
 
     Returns
     -------
@@ -242,12 +190,12 @@ def compute_background_probability(
     """
     if validate_input:
         validate_argument(locals(), "z_lens", float, argmin=0, eqmin=True)
-        validate_argument(locals(), "z_source", "float_array", argmin=0, eqmin=True, none_ok=True)
+        validate_argument(locals(), "z_src", "float_array", argmin=0, eqmin=True, none_ok=True)
 
     if use_pdz is False:
-        if z_source is None:
-            raise ValueError("z_source must be provided.")
-        p_background = np.array(z_source > z_lens, dtype=float)
+        if z_src is None:
+            raise ValueError("z_src must be provided.")
+        p_background = np.array(z_src > z_lens, dtype=float)
     else:
         if pzpdf is None or pzbins is None:
             raise ValueError("pzbins must be provided with pzpdf.")
@@ -257,12 +205,6 @@ def compute_background_probability(
 
 
 def compute_galaxy_weights(
-    z_lens,
-    cosmo,
-    z_source=None,
-    use_pdz=False,
-    pzpdf=None,
-    pzbins=None,
     use_shape_noise=False,
     shape_component1=None,
     shape_component2=None,
@@ -278,55 +220,39 @@ def compute_galaxy_weights(
     The weights :math:`w_{ls}` express as : :math:`w_{ls} = w_{ls, \text{geo}} \times w_{ls,
     \text{shape}}`, following E. S. Sheldon et al. (2003), arXiv:astro-ph/0312036:
 
-    1. The geometrical weight :math:`w_{ls, \text{geo}}` depends on lens and source redshift
-    information. When considering only redshift point estimates, the weights read
+    1. If computed for shear, the geometrical weights :math:`w_{ls, \text{geo}}` are equal to 1. If
+    computed for :math:`\Delta \Sigma`, it depends on lens and source redshift information via the
+    critical surface density. This component can be expressed as:
 
         .. math::
-            w_{ls, \text{geo}} = \Sigma_c(\text{cosmo}, z_l, z_{\text{src}})^{-2}\;.
+            w_{ls, \text{geo}} = \Sigma_\text{crit}(z_l, z_{\text{src}})^{-2}\;.
+
+        when only redshift point estimates are provided, or as:
 
-        If the redshift pdf of each source, :math:`p_{\text{photoz}}(z_s)`, is known,
-        the weights are computed instead as
 
         .. math::
-            w_{ls, \text{geo}} = \left[\int_{\delta + z_l} dz_s p_{\text{photoz}}(z_s)
-            \Sigma_c(\text{cosmo}, z_l, z_s)^{-1}\right]^2
+            w_{ls, \text{geo}} = \Sigma_\text{crit}^\text{eff}(z_l, z_{\text{src}})^{-2}
+            = \left[\int_{\delta + z_l} dz_s \; p_{\text{photoz}}(z_s)
+            \Sigma_\text{crit}(z_l, z_s)^{-1}\right]^2
 
-        for the tangential shear, the weights :math:`w_{ls, \text{geo}}` are 1.
+        when the redshift pdf of each source, :math:`p_{\text{photoz}}(z_s)`, is known.
 
     2. The shape weight :math:`w_{ls,{\text{shape}}}` depends on shapenoise and/or shape
     measurement errors
 
         .. math::
-            w_{ls, \text{shape}} = 1/(\sigma_{\text{shapenoise}}^2 +
-            \sigma_{\text{measurement}}^2)
-
-
-    3. The probability for a galaxy to be in the background of the cluster is defined by:
-
-        .. math::
-            P(z_s > z_l) = \int_{z_l}^{+\infty} dz_s p_{\text{photoz}}(z_s)
-
-        The function return the probability for a galaxy to be in the background of the cluster;
-        if photometric probability density functions are provoded, the function computes the above
-        integral. In the case of true redshifts, it returns 1 if :math:`z_s > z_l` else returns 0.
+            w_{ls, \text{shape}}^{-1} = \sigma_{\text{shapenoise}}^2 +
+            \sigma_{\text{measurement}}^2
 
 
     Parameters
     ----------
-    z_lens: float
-        Redshift of the lens.
-    z_source: array_like, optional
-        Redshift of the source (point estimate). Used only if `use_pdz=False`.
-    cosmo: clmm.Comology object, None
-        CLMM Cosmology object.
-    use_pdz: bool
-        Flag to use or not the source redshift p(z). If `False` (default) the point estimate
-        provided by `z_source` is used.
-    pzpdf : array_like, optional
-        Photometric probablility density functions of the source galaxies.
-        Used instead of z_source if `use_pdz=True`
-    pzbins : array_like, optional
-        Redshift axis on which the individual photoz pdf is tabulated. Required if `use_pdz=True`
+    is_deltasigma: bool
+        If `False`, weights are computed for shear, else weights are computed for
+        :math:`\Delta \Sigma`.
+    sigma_c : None, array_like
+        Critical (effective) surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+        Used only when is_deltasigma=True.
     use_shape_noise: bool
         If `True` shape noise is included in the weight computation. It then requires
         `shape_componenet{1,2}` to be provided. Default: False.
@@ -345,8 +271,6 @@ def compute_galaxy_weights(
     shape_component2_err: array_like
         The measurement error on the 2nd-component of ellipticity of the source galaxies,
         used if `use_shape_error=True`
-    is_deltasigma: bool
-        Indicates whether it is the excess surface density or the tangential shear
     validate_input: bool
         Validade each input argument
 
@@ -356,54 +280,27 @@ def compute_galaxy_weights(
         Individual lens source pair weights
     """
     if validate_input:
-        validate_argument(locals(), "z_lens", float, argmin=0, eqmin=True)
-        validate_argument(locals(), "z_source", "float_array", argmin=0, eqmin=True, none_ok=True)
-        validate_argument(locals(), "use_pdz", bool)
-        # validate_argument(locals(), 'pzpdf', 'float_array', none_ok=True)
-        # validate_argument(locals(), 'pzbins', 'float_array', none_ok=True)
+        validate_argument(locals(), "sigma_c", "float_array", none_ok=True)
         validate_argument(locals(), "shape_component1", "float_array", none_ok=True)
         validate_argument(locals(), "shape_component2", "float_array", none_ok=True)
         validate_argument(locals(), "shape_component1_err", "float_array", none_ok=True)
         validate_argument(locals(), "shape_component2_err", "float_array", none_ok=True)
         validate_argument(locals(), "use_shape_noise", bool)
-        validate_argument(locals(), "is_deltasigma", bool)
+        validate_argument(locals(), "use_shape_error", bool)
         arguments_consistency(
             [shape_component1, shape_component2],
             names=("shape_component1", "shape_component2"),
             prefix="Shape components sources",
         )
+        _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c)
 
     # computing w_ls_geo
-
-    if is_deltasigma is False:
-        w_ls_geo = 1.0
-    else:
-        if sigma_c is None and use_pdz is False:
-            # Need to verify that cosmology and redshifts are provided
-            if any(t_ is None for t_ in (z_lens, z_source, cosmo)):
-                raise TypeError(
-                    "To compute DeltaSigma, please provide a "
-                    "i) cosmology, ii) redshift of lens and sources"
-                )
-            sigma_c = cosmo.eval_sigma_crit(z_lens, z_source)
-        elif sigma_c is None:
-            # Need to verify that cosmology, lens redshift, source redshift bins and
-            # source redshift pdf are provided
-            if any(t_ is None for t_ in (z_lens, cosmo, pzbins, pzpdf)):
-                raise TypeError(
-                    "To compute DeltaSigma using the redshift pdz of the sources, "
-                    "please provide a "
-                    "i) cosmology, ii) lens redshift, iii) source redshift bins and"
-                    "iv) source redshift pdf"
-                )
-            sigma_c = compute_critical_surface_density_eff(
-                cosmo, z_lens, pzbins=pzbins, pzpdf=pzpdf
-            )
-        w_ls_geo = 1.0 / sigma_c**2
+    w_ls_geo = 1.0
+    if sigma_c is not None:
+        w_ls_geo /= np.array(sigma_c) ** 2
 
     # computing w_ls_shape
-    ngals = len(pzpdf) if use_pdz else len(z_source)
-    err_e2 = np.zeros(ngals)
+    err_e2 = 0
 
     if use_shape_noise:
         if shape_component1 is None or shape_component2 is None:
@@ -420,8 +317,12 @@ def compute_galaxy_weights(
             )
         err_e2 += shape_component1_err**2
         err_e2 += shape_component2_err**2
-    w_ls_shape = np.ones(ngals)
-    w_ls_shape[err_e2 > 0] = 1.0 / err_e2[err_e2 > 0]
+
+    if hasattr(err_e2, "__len__"):
+        w_ls_shape = np.ones(len(err_e2))
+        w_ls_shape[err_e2 > 0] = 1.0 / err_e2[err_e2 > 0]
+    else:
+        w_ls_shape = 1.0 / err_e2 if err_e2 > 0 else 1.0
 
     w_ls = w_ls_shape * w_ls_geo
 
diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py
index cbd6e51e4..b4db5446c 100644
--- a/clmm/galaxycluster.py
+++ b/clmm/galaxycluster.py
@@ -114,7 +114,7 @@ def _repr_html_(self):
             f"<br>{self.galcat._html_table()}"
         )
 
-    def add_critical_surface_density(self, cosmo, use_pdz=False):
+    def add_critical_surface_density(self, cosmo, use_pdz=False, force=False):
         r"""Computes the critical surface density for each galaxy in `galcat`.
         It only runs if input cosmo != galcat cosmo or if `sigma_c` not in `galcat`.
 
@@ -128,6 +128,8 @@ def add_critical_surface_density(self, cosmo, use_pdz=False):
             `sigma_c` is computed as 1/<1/Sigma_crit>, where the average is performed using
             the individual galaxy redshift pdf. In that case, the `galcat` table should have
             pzbins` and `pzpdf` columns.
+        force : bool
+            Force recomputation of sigma_c.
 
         Returns
         -------
@@ -135,7 +137,12 @@ def add_critical_surface_density(self, cosmo, use_pdz=False):
         """
         if cosmo is None:
             raise TypeError("To compute Sigma_crit, please provide a cosmology")
-        if cosmo.get_desc() != self.galcat.meta["cosmo"] or "sigma_c" not in self.galcat.columns:
+        sigmac_colname = "sigma_c_eff" if use_pdz else "sigma_c"
+        if (
+            cosmo.get_desc() != self.galcat.meta["cosmo"]
+            or sigmac_colname not in self.galcat.columns
+            or force
+        ):
             if self.z is None:
                 raise TypeError("Cluster's redshift is None. Cannot compute Sigma_crit")
             if not use_pdz and "z" not in self.galcat.columns:
@@ -149,19 +156,18 @@ def add_critical_surface_density(self, cosmo, use_pdz=False):
                 )
 
             self.galcat.update_cosmo(cosmo, overwrite=True)
-            if use_pdz is False:
-                self.galcat["sigma_c"] = cosmo.eval_sigma_crit(self.z, self.galcat["z"])
-                self.galcat.meta["sigmac_type"] = "standard"
+            if not use_pdz:
+                self.galcat[sigmac_colname] = cosmo.eval_sigma_crit(self.z, self.galcat["z"])
             else:
                 zdata = self._get_input_galdata({"pzpdf": "pzpdf", "pzbins": "pzbins"})
-                self.galcat["sigma_c"] = compute_critical_surface_density_eff(
+                self.galcat[sigmac_colname] = compute_critical_surface_density_eff(
                     cosmo=cosmo,
                     z_cluster=self.z,
                     pzbins=zdata["pzbins"],
                     pzpdf=zdata["pzpdf"],
                     validate_input=self.validate_input,
                 )
-                self.galcat.meta["sigmac_type"] = "effective"
+        return sigmac_colname
 
     def _get_input_galdata(self, col_dict):
         """
@@ -211,15 +217,14 @@ def compute_tangential_and_cross_components(
         r"""Adds a tangential- and cross- components for shear or ellipticity to self
 
         Calls `clmm.dataops.compute_tangential_and_cross_components` with the following arguments:
-        ra_lens: cluster Ra
-        dec_lens: cluster Dec
+        ra_lens: `cluster` Ra
+        dec_lens: `cluster` Dec
         ra_source: `galcat` Ra
         dec_source: `galcat` Dec
         shear1: `galcat` shape_component1
         shear2: `galcat` shape_component2
         geometry: `input` geometry
         is_deltasigma: `input` is_deltasigma
-        sigma_c: `galcat` sigma_c | None
 
         Parameters
         ----------
@@ -258,26 +263,23 @@ def compute_tangential_and_cross_components(
             Cross shear (or assimilated quantity) for each source galaxy
         """
         # Check is all the required data is available
-        cols = self._get_input_galdata(
-            {
-                "ra_source": "ra",
-                "dec_source": "dec",
-                "shear1": shape_component1,
-                "shear2": shape_component2,
-            }
-        )
-
+        col_dict = {
+            "ra_source": "ra",
+            "dec_source": "dec",
+            "shear1": shape_component1,
+            "shear2": shape_component2,
+        }
         if is_deltasigma:
-            self.add_critical_surface_density(cosmo, use_pdz=use_pdz)
-            cols["sigma_c"] = self.galcat["sigma_c"]
+            sigmac_colname = self.add_critical_surface_density(cosmo, use_pdz=use_pdz)
+            col_dict.update({"sigma_c": sigmac_colname})
+        cols = self._get_input_galdata(col_dict)
 
         # compute shears
         angsep, tangential_comp, cross_comp = compute_tangential_and_cross_components(
+            is_deltasigma=is_deltasigma,
             ra_lens=self.ra,
             dec_lens=self.dec,
             geometry=geometry,
-            is_deltasigma=is_deltasigma,
-            use_pdz=use_pdz,
             validate_input=self.validate_input,
             **cols,
         )
@@ -285,6 +287,10 @@ def compute_tangential_and_cross_components(
             self.galcat["theta"] = angsep
             self.galcat[tan_component] = tangential_comp
             self.galcat[cross_component] = cross_comp
+            if is_deltasigma:
+                sigmac_type = "effective" if use_pdz else "standard"
+                self.galcat.meta[f"{tan_component}_sigmac_type"] = sigmac_type
+                self.galcat.meta[f"{cross_component}_sigmac_type"] = sigmac_type
         return angsep, tangential_comp, cross_comp
 
     def compute_background_probability(
@@ -308,7 +314,7 @@ def compute_background_probability(
             Probability for being a background galaxy
         """
         cols = self._get_input_galdata(
-            {"pzpdf": "pzpdf", "pzbins": "pzbins"} if use_pdz else {"z_source": "z"}
+            {"pzpdf": "pzpdf", "pzbins": "pzbins"} if use_pdz else {"z_src": "z"}
         )
         p_background = compute_background_probability(
             self.z, use_pdz=use_pdz, validate_input=self.validate_input, **cols
@@ -369,12 +375,9 @@ def compute_galaxy_weights(
         """
         # input cols
         col_dict = {}
-        if use_pdz:
-            col_dict.update({"pzpdf": "pzpdf", "pzbins": "pzbins"})
         if is_deltasigma:
-            if "sigma_c" not in self.galcat.columns:
-                self.add_critical_surface_density(cosmo, use_pdz=use_pdz)
-            col_dict.update({"z_source": "z", "sigma_c": "sigma_c"})
+            sigmac_colname = self.add_critical_surface_density(cosmo, use_pdz=use_pdz)
+            col_dict.update({"sigma_c": sigmac_colname})
         if use_shape_noise:
             col_dict.update(
                 {
@@ -393,17 +396,18 @@ def compute_galaxy_weights(
 
         # computes weights
         w_ls = compute_galaxy_weights(
-            self.z,
-            cosmo,
-            use_pdz=use_pdz,
+            is_deltasigma=is_deltasigma,
             use_shape_noise=use_shape_noise,
             use_shape_error=use_shape_error,
-            is_deltasigma=is_deltasigma,
             validate_input=self.validate_input,
             **cols,
         )
         if add:
             self.galcat[weight_name] = w_ls
+            if is_deltasigma:
+                self.galcat.meta[f"{weight_name}_sigmac_type"] = (
+                    "effective" if use_pdz else "standard"
+                )
         return w_ls
 
     def draw_gal_z_from_pdz(self, zcol_out="z", overwrite=False, nobj=1, xmin=None, xmax=None):
diff --git a/clmm/support/mock_data.py b/clmm/support/mock_data.py
index 8ef515f90..170b7c94b 100644
--- a/clmm/support/mock_data.py
+++ b/clmm/support/mock_data.py
@@ -392,7 +392,7 @@ def _generate_galaxy_catalog(
         mdelta=cluster_m,
         cdelta=cluster_c,
         z_cluster=cluster_z,
-        z_source=galaxy_catalog["ztrue"],
+        z_src=galaxy_catalog["ztrue"],
         cosmo=cosmo,
         delta_mdef=delta_so,
         halo_profile_model=halo_profile_model,
@@ -407,7 +407,7 @@ def _generate_galaxy_catalog(
         mdelta=cluster_m,
         cdelta=cluster_c,
         z_cluster=cluster_z,
-        z_source=galaxy_catalog["ztrue"],
+        z_src=galaxy_catalog["ztrue"],
         cosmo=cosmo,
         delta_mdef=delta_so,
         halo_profile_model=halo_profile_model,
diff --git a/clmm/theory/func_layer.py b/clmm/theory/func_layer.py
index 0eb2d50c5..ac43d69bd 100644
--- a/clmm/theory/func_layer.py
+++ b/clmm/theory/func_layer.py
@@ -532,7 +532,7 @@ def compute_tangential_shear(
     mdelta,
     cdelta,
     z_cluster,
-    z_source,
+    z_src,
     cosmo,
     delta_mdef=200,
     halo_profile_model="nfw",
@@ -563,7 +563,7 @@ def compute_tangential_shear(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float, function
+    z_src : array_like, float, function
         Information on the background source galaxy redshift(s). Value required depends on
         `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
@@ -590,18 +590,18 @@ def compute_tangential_shear(
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
     z_src_info : str, optional
-        Type of redshift information provided by the `z_source` argument.
+        Type of redshift information provided by the `z_src` argument.
         The following supported options are:
 
-            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
-              It can be individual redshifts for each source galaxy when `z_source` is an array
-              or all sources are at the same redshift when `z_source` is a float.
+            * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+              It can be individual redshifts for each source galaxy when `z_src` is an array
+              or all sources are at the same redshift when `z_src` is a float.
 
-            * 'distribution' : A redshift distribution function is provided by `z_source`.
-              `z_source` must be a one dimensional function.
+            * 'distribution' : A redshift distribution function is provided by `z_src`.
+              `z_src` must be a one dimensional function.
 
-            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
-              `z_source` must be a tuple containing
+            * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+              `z_src` must be a tuple containing
               ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
               the lensing efficiency and square of the lensing efficiency averaged over
               the galaxy redshift distribution repectively.
@@ -655,7 +655,7 @@ def compute_tangential_shear(
     tangential_shear = _modeling_object.eval_tangential_shear(
         r_proj,
         z_cluster,
-        z_source,
+        z_src,
         z_src_info=z_src_info,
         beta_kwargs=beta_kwargs,
         verbose=verbose,
@@ -670,7 +670,7 @@ def compute_convergence(
     mdelta,
     cdelta,
     z_cluster,
-    z_source,
+    z_src,
     cosmo,
     delta_mdef=200,
     halo_profile_model="nfw",
@@ -701,7 +701,7 @@ def compute_convergence(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float, function
+    z_src : array_like, float, function
         Information on the background source galaxy redshift(s). Value required depends on
         `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
@@ -728,18 +728,18 @@ def compute_convergence(
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
     z_src_info : str, optional
-        Type of redshift information provided by the `z_source` argument.
+        Type of redshift information provided by the `z_src` argument.
         The following supported options are:
 
-            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
-              It can be individual redshifts for each source galaxy when `z_source` is an array
-              or all sources are at the same redshift when `z_source` is a float.
+            * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+              It can be individual redshifts for each source galaxy when `z_src` is an array
+              or all sources are at the same redshift when `z_src` is a float.
 
-            * 'distribution' : A redshift distribution function is provided by `z_source`.
-              `z_source` must be a one dimensional function.
+            * 'distribution' : A redshift distribution function is provided by `z_src`.
+              `z_src` must be a one dimensional function.
 
-            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
-              `z_source` must be a tuple containing
+            * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+              `z_src` must be a tuple containing
               ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
               the lensing efficiency and square of the lensing efficiency averaged over
               the galaxy redshift distribution repectively.
@@ -789,7 +789,7 @@ def compute_convergence(
     convergence = _modeling_object.eval_convergence(
         r_proj,
         z_cluster,
-        z_source,
+        z_src,
         z_src_info=z_src_info,
         beta_kwargs=beta_kwargs,
         verbose=verbose,
@@ -804,7 +804,7 @@ def compute_reduced_tangential_shear(
     mdelta,
     cdelta,
     z_cluster,
-    z_source,
+    z_src,
     cosmo,
     delta_mdef=200,
     halo_profile_model="nfw",
@@ -831,7 +831,7 @@ def compute_reduced_tangential_shear(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float, function
+    z_src : array_like, float, function
         Information on the background source galaxy redshift(s). Value required depends on
         `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
@@ -858,18 +858,18 @@ def compute_reduced_tangential_shear(
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
     z_src_info : str, optional
-        Type of redshift information provided by the `z_source` argument.
+        Type of redshift information provided by the `z_src` argument.
         The following supported options are:
 
-            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
-              It can be individual redshifts for each source galaxy when `z_source` is an array
-              or all sources are at the same redshift when `z_source` is a float.
+            * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+              It can be individual redshifts for each source galaxy when `z_src` is an array
+              or all sources are at the same redshift when `z_src` is a float.
 
-            * 'distribution' : A redshift distribution function is provided by `z_source`.
-              `z_source` must be a one dimensional function.
+            * 'distribution' : A redshift distribution function is provided by `z_src`.
+              `z_src` must be a one dimensional function.
 
-            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
-              `z_source` must be a tuple containing
+            * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+              `z_src` must be a tuple containing
               ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
               the lensing efficiency and square of the lensing efficiency averaged over
               the galaxy redshift distribution repectively.
@@ -887,7 +887,7 @@ def compute_reduced_tangential_shear(
 
             * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
               If `z_src_info='discrete'`, full computation is made for each
-              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, reduced
+              `r_proj, z_src` pair individually. If `z_src_info='distribution'`, reduced
               tangential shear at each value of `r_proj` is calculated as
 
               .. math::
@@ -959,7 +959,7 @@ def compute_reduced_tangential_shear(
     red_tangential_shear = _modeling_object.eval_reduced_tangential_shear(
         r_proj,
         z_cluster,
-        z_source,
+        z_src,
         z_src_info=z_src_info,
         approx=approx,
         beta_kwargs=beta_kwargs,
@@ -975,7 +975,7 @@ def compute_magnification(
     mdelta,
     cdelta,
     z_cluster,
-    z_source,
+    z_src,
     cosmo,
     delta_mdef=200,
     halo_profile_model="nfw",
@@ -1002,7 +1002,7 @@ def compute_magnification(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float, function
+    z_src : array_like, float, function
         Information on the background source galaxy redshift(s). Value required depends on
         `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
@@ -1029,18 +1029,18 @@ def compute_magnification(
         If None, use the default value of the backend. (0.25 for the NumCosmo backend and a
         cosmology-dependent value for the CCL backend.)
     z_src_info : str, optional
-        Type of redshift information provided by the `z_source` argument.
+        Type of redshift information provided by the `z_src` argument.
         The following supported options are:
 
-            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
-              It can be individual redshifts for each source galaxy when `z_source` is an array
-              or all sources are at the same redshift when `z_source` is a float.
+            * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+              It can be individual redshifts for each source galaxy when `z_src` is an array
+              or all sources are at the same redshift when `z_src` is a float.
 
-            * 'distribution' : A redshift distribution function is provided by `z_source`.
-              `z_source` must be a one dimensional function.
+            * 'distribution' : A redshift distribution function is provided by `z_src`.
+              `z_src` must be a one dimensional function.
 
-            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
-              `z_source` must be a tuple containing
+            * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+              `z_src` must be a tuple containing
               ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
               the lensing efficiency and square of the lensing efficiency averaged over
               the galaxy redshift distribution repectively.
@@ -1058,7 +1058,7 @@ def compute_magnification(
 
             * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
               If `z_src_info='discrete'`, full computation is made for each
-              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, magnification
+              `r_proj, z_src` pair individually. If `z_src_info='distribution'`, magnification
               at each value of `r_proj` is calculated as
 
               .. math::
@@ -1125,7 +1125,7 @@ def compute_magnification(
     magnification = _modeling_object.eval_magnification(
         r_proj,
         z_cluster,
-        z_source,
+        z_src,
         z_src_info=z_src_info,
         approx=approx,
         beta_kwargs=beta_kwargs,
@@ -1142,7 +1142,7 @@ def compute_magnification_bias(
     mdelta,
     cdelta,
     z_cluster,
-    z_source,
+    z_src,
     cosmo,
     delta_mdef=200,
     halo_profile_model="nfw",
@@ -1185,7 +1185,7 @@ def compute_magnification_bias(
         Galaxy cluster NFW concentration.
     z_cluster : float
         Galaxy cluster redshift
-    z_source : array_like, float, function
+    z_src : array_like, float, function
         Information on the background source galaxy redshift(s). Value required depends on
         `z_src_info` (see below).
     cosmo : clmm.cosmology.Cosmology object
@@ -1212,18 +1212,18 @@ def compute_magnification_bias(
             * 'virial'
 
     z_src_info : str, optional
-        Type of redshift information provided by the `z_source` argument.
+        Type of redshift information provided by the `z_src` argument.
         The following supported options are:
 
-            * 'discrete' (default) : The redshift of sources is provided by `z_source`.
-              It can be individual redshifts for each source galaxy when `z_source` is an array
-              or all sources are at the same redshift when `z_source` is a float.
+            * 'discrete' (default) : The redshift of sources is provided by `z_src`.
+              It can be individual redshifts for each source galaxy when `z_src` is an array
+              or all sources are at the same redshift when `z_src` is a float.
 
-            * 'distribution' : A redshift distribution function is provided by `z_source`.
-              `z_source` must be a one dimensional function.
+            * 'distribution' : A redshift distribution function is provided by `z_src`.
+              `z_src` must be a one dimensional function.
 
-            * 'beta' : The averaged lensing efficiency is provided by `z_source`.
-              `z_source` must be a tuple containing
+            * 'beta' : The averaged lensing efficiency is provided by `z_src`.
+              `z_src` must be a tuple containing
               ( :math:`\langle \beta_s \rangle, \langle \beta_s^2 \rangle`),
               the lensing efficiency and square of the lensing efficiency averaged over
               the galaxy redshift distribution repectively.
@@ -1241,7 +1241,7 @@ def compute_magnification_bias(
 
             * None (default): Requires `z_src_info` to be 'discrete' or 'distribution'.
               If `z_src_info='discrete'`, full computation is made for each
-              `r_proj, z_source` pair individually. If `z_src_info='distribution'`, magnification
+              `r_proj, z_src` pair individually. If `z_src_info='distribution'`, magnification
               bias at each value of `r_proj` is calculated as
 
               .. math::
@@ -1310,7 +1310,7 @@ def compute_magnification_bias(
     magnification_bias = _modeling_object.eval_magnification_bias(
         r_proj,
         z_cluster,
-        z_source,
+        z_src,
         alpha,
         z_src_info=z_src_info,
         approx=approx,
diff --git a/clmm/utils/__init__.py b/clmm/utils/__init__.py
index 9f8d1197c..cf15420ee 100644
--- a/clmm/utils/__init__.py
+++ b/clmm/utils/__init__.py
@@ -38,6 +38,7 @@
     validate_argument,
     _validate_ra,
     _validate_dec,
+    _validate_is_deltasigma_sigma_c,
 )
 
 from .units import (
diff --git a/clmm/utils/validation.py b/clmm/utils/validation.py
index 29d4dffac..ed020d183 100644
--- a/clmm/utils/validation.py
+++ b/clmm/utils/validation.py
@@ -207,3 +207,20 @@ def _validate_dec(loc, dec_name, is_array):
     """
     v_type = "float_array" if is_array else (float, str)
     validate_argument(loc, dec_name, v_type, argmin=-90, eqmin=True, argmax=90, eqmax=True)
+
+
+def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c):
+    r""" "Validate the compatibility between is_deltasigma and sigma_c arguments.
+
+
+    Parameters
+    ----------
+    is_deltasigma: bool
+        If `False`, values are computed for shear, else they are computed for :math:`\Delta \Sigma`.
+    sigma_c : None, array_like
+        Critical (effective) surface density in units of :math:`M_\odot\ Mpc^{-2}`.
+    """
+    if is_deltasigma and sigma_c is None:
+        raise TypeError("sigma_c (=None) must be provided when is_deltasigma=True")
+    if not is_deltasigma and sigma_c is not None:
+        raise TypeError(f"sigma_c (={sigma_c}) must be None when is_deltasigma=False")
diff --git a/examples/NumCosmo/Example1_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb b/examples/NumCosmo/Example1_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
index 2d963f0cd..1c6716138 100644
--- a/examples/NumCosmo/Example1_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
+++ b/examples/NumCosmo/Example1_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
@@ -568,21 +568,21 @@
    "source": [
     "class GaussGammaT(Ncm.DataGaussDiag):\n",
     "    z_cluster = GObject.Property(type=float, flags=GObject.PARAM_READWRITE)\n",
-    "    z_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
+    "    z_src = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "    r_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "\n",
-    "    def __init__(self, z_cluster, r_source, z_source, gt_profile, moo=None):\n",
+    "    def __init__(self, z_cluster, r_source, z_src, gt_profile, moo=None):\n",
     "        Ncm.DataGaussDiag.__init__(self, n_points=len(gt_profile))\n",
     "\n",
     "        self.moo = moo if moo else clmm.Modeling()\n",
     "\n",
-    "        assert len(gt_profile) == len(z_source)\n",
+    "        assert len(gt_profile) == len(z_src)\n",
     "        assert len(gt_profile) == len(r_source)\n",
     "\n",
     "        self.set_size(len(gt_profile))\n",
     "\n",
     "        self.props.z_cluster = z_cluster\n",
-    "        self.props.z_source = Ncm.Vector.new_array(z_source)\n",
+    "        self.props.z_src = Ncm.Vector.new_array(z_src)\n",
     "        self.props.r_source = Ncm.Vector.new_array(r_source)\n",
     "\n",
     "        self.y.set_array(gt_profile)\n",
@@ -611,7 +611,7 @@
     "            self.moo.eval_reduced_tangential_shear(\n",
     "                self.props.r_source.dup_array(),\n",
     "                self.props.z_cluster,\n",
-    "                self.props.z_source.dup_array(),\n",
+    "                self.props.z_src.dup_array(),\n",
     "            )\n",
     "        )\n",
     "        return\n",
@@ -643,8 +643,8 @@
     "moo2.set_cosmo(cosmo)\n",
     "moo2.set_concentration(4.0)\n",
     "\n",
-    "ggt1 = GaussGammaT(z_cluster=cluster_z, r_source=r1, z_source=z1, gt_profile=gt_profile1, moo=moo1)\n",
-    "ggt2 = GaussGammaT(z_cluster=cluster_z, r_source=r2, z_source=z2, gt_profile=gt_profile2, moo=moo2)\n",
+    "ggt1 = GaussGammaT(z_cluster=cluster_z, r_source=r1, z_src=z1, gt_profile=gt_profile1, moo=moo1)\n",
+    "ggt2 = GaussGammaT(z_cluster=cluster_z, r_source=r2, z_src=z2, gt_profile=gt_profile2, moo=moo2)\n",
     "\n",
     "mset1 = ggt1.moo.get_mset()\n",
     "mset2 = ggt2.moo.get_mset()\n",
diff --git a/examples/NumCosmo/Example2_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb b/examples/NumCosmo/Example2_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
index 2e1933930..aab95262e 100644
--- a/examples/NumCosmo/Example2_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
+++ b/examples/NumCosmo/Example2_Fit_Halo_Mass_to_Shear_Catalog_NC.ipynb
@@ -444,7 +444,7 @@
    "source": [
     "class GaussGammaTErr(Ncm.DataGaussDiag):\n",
     "    z_cluster = GObject.Property(type=float, flags=GObject.PARAM_READWRITE)\n",
-    "    z_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
+    "    z_src = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "    r_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "    z_err = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "\n",
@@ -452,20 +452,18 @@
     "        Ncm.DataGaussDiag.__init__(self, n_points=0)\n",
     "        self.moo = clmm.Modeling()\n",
     "\n",
-    "    def init_from_data(\n",
-    "        self, z_cluster, r_source, z_source, gt_profile, gt_err, z_err=None, moo=None\n",
-    "    ):\n",
+    "    def init_from_data(self, z_cluster, r_source, z_src, gt_profile, gt_err, z_err=None, moo=None):\n",
     "        if moo:\n",
     "            self.moo = moo\n",
     "\n",
-    "        assert len(gt_profile) == len(z_source)\n",
+    "        assert len(gt_profile) == len(z_src)\n",
     "        assert len(gt_profile) == len(r_source)\n",
     "        assert len(gt_profile) == len(gt_err)\n",
     "\n",
     "        self.set_size(len(gt_profile))\n",
     "\n",
     "        self.props.z_cluster = z_cluster\n",
-    "        self.props.z_source = Ncm.Vector.new_array(z_source)\n",
+    "        self.props.z_src = Ncm.Vector.new_array(z_src)\n",
     "        self.props.r_source = Ncm.Vector.new_array(r_source)\n",
     "        if z_err:\n",
     "            self.props.r_source = Ncm.Vector.new_array(z_err)\n",
@@ -496,7 +494,7 @@
     "            self.moo.eval_reduced_tangential_shear(\n",
     "                self.props.r_source.dup_array(),\n",
     "                self.props.z_cluster,\n",
-    "                self.props.z_source.dup_array(),\n",
+    "                self.props.z_src.dup_array(),\n",
     "            )\n",
     "        )\n",
     "        return\n",
@@ -536,7 +534,7 @@
     "ggt1.init_from_data(\n",
     "    z_cluster=cluster_z,\n",
     "    r_source=profile1[\"radius\"],\n",
-    "    z_source=profile1[\"z\"],\n",
+    "    z_src=profile1[\"z\"],\n",
     "    gt_profile=profile1[\"gt\"],\n",
     "    gt_err=profile1[\"gt_err\"],\n",
     "    moo=moo1,\n",
@@ -544,7 +542,7 @@
     "ggt2.init_from_data(\n",
     "    z_cluster=cluster_z,\n",
     "    r_source=profile2[\"radius\"],\n",
-    "    z_source=profile2[\"z\"],\n",
+    "    z_src=profile2[\"z\"],\n",
     "    gt_profile=profile2[\"gt\"],\n",
     "    gt_err=profile2[\"gt_err\"],\n",
     "    moo=moo2,\n",
@@ -552,7 +550,7 @@
     "ggt3.init_from_data(\n",
     "    z_cluster=cluster_z,\n",
     "    r_source=profile3[\"radius\"],\n",
-    "    z_source=profile3[\"z\"],\n",
+    "    z_src=profile3[\"z\"],\n",
     "    gt_profile=profile3[\"gt\"],\n",
     "    gt_err=profile3[\"gt_err\"],\n",
     "    moo=moo3,\n",
@@ -806,7 +804,7 @@
    "outputs": [],
    "source": [
     "def create_nc_data_cluster_wl(\n",
-    "    theta, g_t, z_source, z_cluster, cosmo, dist, sigma_z=None, sigma_g=None\n",
+    "    theta, g_t, z_src, z_cluster, cosmo, dist, sigma_z=None, sigma_g=None\n",
     "):\n",
     "    r = clmm.convert_units(theta, \"radians\", \"Mpc\", redshift=z_cluster, cosmo=cosmo)\n",
     "    ga = Ncm.ObjArray.new()\n",
@@ -819,11 +817,11 @@
     "\n",
     "    if sigma_z:\n",
     "        gzgs = Nc.GalaxyRedshiftGauss()\n",
-    "        z_obs = np.column_stack((z_source, (1.0 + z_source) * sigma_z))\n",
+    "        z_obs = np.column_stack((z_src, (1.0 + z_src) * sigma_z))\n",
     "        gzgs.set_obs(Ncm.Matrix.new_array(z_obs.flatten(), 2))\n",
     "    else:\n",
     "        gzgs = Nc.GalaxyRedshiftSpec()\n",
-    "        gzgs.set_z(Ncm.Vector.new_array(z_source))\n",
+    "        gzgs.set_z(Ncm.Vector.new_array(z_src))\n",
     "\n",
     "    gwl = Nc.GalaxyWL(wl_dist=grsg, gz_dist=gzgs)\n",
     "    ga.add(gwl)\n",
diff --git a/examples/NumCosmo/modeling_cmp_numcosmo.ipynb b/examples/NumCosmo/modeling_cmp_numcosmo.ipynb
index d7a2b5e3e..1e2accb3c 100644
--- a/examples/NumCosmo/modeling_cmp_numcosmo.ipynb
+++ b/examples/NumCosmo/modeling_cmp_numcosmo.ipynb
@@ -121,7 +121,7 @@
     "cluster_concentration = 4.0\n",
     "z_max = 3.0\n",
     "z_cluster = 1.0\n",
-    "z_source = 2.0"
+    "z_src = 2.0"
    ]
   },
   {
@@ -446,8 +446,8 @@
     }
    ],
    "source": [
-    "Sigmac = pp.compute_critical_surface_density(cosmo_ccl, z_cluster=z_cluster, z_source=z_source)\n",
-    "nc_Sigmac = smd.sigma_critical(cosmo, z_source, z_cluster, z_cluster)\n",
+    "Sigmac = pp.compute_critical_surface_density(cosmo_ccl, z_cluster=z_cluster, z_src=z_src)\n",
+    "nc_Sigmac = smd.sigma_critical(cosmo, z_src, z_cluster, z_cluster)\n",
     "\n",
     "print(\"% 22.15g % 22.15g %e\" % (Sigmac, nc_Sigmac, Sigmac / nc_Sigmac - 1.0))"
    ]
@@ -484,14 +484,14 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cluster,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo_ccl,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=\"nfw\",\n",
     "    z_src_model=\"single_plane\",\n",
     ")\n",
     "\n",
-    "nc_gammat = np.array([smd.shear(dp, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d])\n",
+    "nc_gammat = np.array([smd.shear(dp, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d])\n",
     "\n",
     "fig, axs = compare(r3d, gammat, nc_gammat, x_name=\"radius [mpc]\", y_name=r\"{\\gamma_t}\")\n",
     "axs[0].set_xscale(\"log\")"
@@ -529,7 +529,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cluster,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo_ccl,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=\"nfw\",\n",
@@ -537,7 +537,7 @@
     ")\n",
     "\n",
     "nc_kappa = np.array(\n",
-    "    [smd.convergence(dp, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.convergence(dp, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "\n",
     "fig, axs = compare(r3d, kappa, nc_kappa, x_name=\"radius [mpc]\", y_name=r\"\\kappa\")\n",
@@ -576,14 +576,14 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cluster,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo_ccl,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=\"nfw\",\n",
     "    z_src_model=\"single_plane\",\n",
     ")\n",
     "nc_gt = np.array(\n",
-    "    [smd.reduced_shear(dp, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.reduced_shear(dp, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "\n",
     "fig, axs = compare(r3d, gt, nc_gt, x_name=\"radius [mpc]\", y_name=r\"{g_t}\")\n",
@@ -598,7 +598,7 @@
    "outputs": [],
    "source": [
     "nc_mu = np.array(\n",
-    "    [smd.magnification(dp, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.magnification(dp, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")"
    ]
   },
@@ -639,15 +639,15 @@
    "source": [
     "config_dict = {\n",
     "    \"dl\": dist.angular_diameter(cosmo, z_cluster) * cosmo.RH_Mpc(),\n",
-    "    \"z_source\": z_source,\n",
+    "    \"z_src\": z_src,\n",
     "    \"cluster_mass\": cluster_mass,\n",
     "    \"G[m3/km.s2]\": Ncm.C.G(),\n",
     "    \"aexp_cluster\": 1.0 / (1.0 + z_cluster),\n",
     "    \"density_profile_parametrization\": density_profile_parametrization,\n",
-    "    \"ds\": dist.angular_diameter(cosmo, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"ds\": dist.angular_diameter(cosmo, z_src) * cosmo.RH_Mpc(),\n",
     "    \"cosmo_Ob0\": cosmo_ccl[\"Omega_b0\"],\n",
-    "    \"aexp_source\": 1.0 / (1.0 + z_source),\n",
-    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"aexp_source\": 1.0 / (1.0 + z_src),\n",
+    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_src) * cosmo.RH_Mpc(),\n",
     "    \"z_cluster\": z_cluster,\n",
     "    \"mass_Delta\": mass_Delta,\n",
     "    \"lightspeed[km/s]\": Ncm.C.c() / 1000.0,\n",
diff --git a/examples/Paper_v1.0/gt_and_use_case.ipynb b/examples/Paper_v1.0/gt_and_use_case.ipynb
index 31c94f3ea..ddd60862b 100644
--- a/examples/Paper_v1.0/gt_and_use_case.ipynb
+++ b/examples/Paper_v1.0/gt_and_use_case.ipynb
@@ -325,7 +325,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=np.mean(cluster.galcat[\"z\"]),  # Mean value of source galaxies redshift\n",
+    "        z_src=np.mean(cluster.galcat[\"z\"]),  # Mean value of source galaxies redshift\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -369,7 +369,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=z_inf,  # Redshift value at infinity\n",
+    "        z_src=z_inf,  # Redshift value at infinity\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -379,7 +379,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=z_inf,  # Redshift value at infinity\n",
+    "        z_src=z_inf,  # Redshift value at infinity\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -421,7 +421,7 @@
     "                    cdelta=4,  # Concentration of the cluster\n",
     "                    z_cluster=cluster_z,  # Redshift of the cluster\n",
     "                    # Redshift value of each source galaxy inside the radial bin\n",
-    "                    z_source=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
+    "                    z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
     "                    cosmo=cosmo,\n",
     "                    delta_mdef=200,\n",
     "                    halo_profile_model=\"nfw\",\n",
diff --git a/examples/Paper_v1.0/mcmc.ipynb b/examples/Paper_v1.0/mcmc.ipynb
index d2c65d157..caa0393ca 100644
--- a/examples/Paper_v1.0/mcmc.ipynb
+++ b/examples/Paper_v1.0/mcmc.ipynb
@@ -442,7 +442,7 @@
    "source": [
     "class GaussGammaTErr(Ncm.DataGaussDiag):\n",
     "    z_cluster = GObject.Property(type=float, flags=GObject.PARAM_READWRITE)\n",
-    "    z_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
+    "    z_src = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "    r_source = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "    z_err = GObject.Property(type=Ncm.Vector, flags=GObject.PARAM_READWRITE)\n",
     "\n",
@@ -450,20 +450,18 @@
     "        Ncm.DataGaussDiag.__init__(self, n_points=0)\n",
     "        self.moo = clmm.Modeling()\n",
     "\n",
-    "    def init_from_data(\n",
-    "        self, z_cluster, r_source, z_source, gt_profile, gt_err, z_err=None, moo=None\n",
-    "    ):\n",
+    "    def init_from_data(self, z_cluster, r_source, z_src, gt_profile, gt_err, z_err=None, moo=None):\n",
     "        if moo:\n",
     "            self.moo = moo\n",
     "\n",
-    "        assert len(gt_profile) == len(z_source)\n",
+    "        assert len(gt_profile) == len(z_src)\n",
     "        assert len(gt_profile) == len(r_source)\n",
     "        assert len(gt_profile) == len(gt_err)\n",
     "\n",
     "        self.set_size(len(gt_profile))\n",
     "\n",
     "        self.props.z_cluster = z_cluster\n",
-    "        self.props.z_source = Ncm.Vector.new_array(z_source)\n",
+    "        self.props.z_src = Ncm.Vector.new_array(z_src)\n",
     "        self.props.r_source = Ncm.Vector.new_array(r_source)\n",
     "        if z_err:\n",
     "            self.props.r_source = Ncm.Vector.new_array(z_err)\n",
@@ -494,7 +492,7 @@
     "            self.moo.eval_reduced_tangential_shear(\n",
     "                self.props.r_source.dup_array(),\n",
     "                self.props.z_cluster,\n",
-    "                self.props.z_source.dup_array(),\n",
+    "                self.props.z_src.dup_array(),\n",
     "            )\n",
     "        )\n",
     "        return\n",
@@ -526,7 +524,7 @@
     "ggt.init_from_data(\n",
     "    z_cluster=cluster_z,\n",
     "    r_source=cl.profile[\"radius\"],\n",
-    "    z_source=cl.profile[\"z\"],\n",
+    "    z_src=cl.profile[\"z\"],\n",
     "    gt_profile=cl.profile[\"gt\"],\n",
     "    gt_err=cl.profile[\"gt_err\"],\n",
     "    moo=moo,\n",
@@ -882,7 +880,7 @@
    "outputs": [],
    "source": [
     "def create_nc_data_cluster_wl(\n",
-    "    theta, g_t, z_source, z_cluster, cosmo, dist, sigma_z=None, sigma_g=None\n",
+    "    theta, g_t, z_src, z_cluster, cosmo, dist, sigma_z=None, sigma_g=None\n",
     "):\n",
     "    r = clmm.convert_units(theta, \"radians\", \"Mpc\", redshift=z_cluster, cosmo=cosmo)\n",
     "    ga = Ncm.ObjArray.new()\n",
@@ -895,11 +893,11 @@
     "\n",
     "    if sigma_z:\n",
     "        gzgs = Nc.GalaxyRedshiftGauss()\n",
-    "        z_obs = np.column_stack((z_source, (1.0 + z_source) * sigma_z))\n",
+    "        z_obs = np.column_stack((z_src, (1.0 + z_src) * sigma_z))\n",
     "        gzgs.set_obs(Ncm.Matrix.new_array(z_obs.flatten(), 2))\n",
     "    else:\n",
     "        gzgs = Nc.GalaxyRedshiftSpec()\n",
-    "        gzgs.set_z(Ncm.Vector.new_array(z_source))\n",
+    "        gzgs.set_z(Ncm.Vector.new_array(z_src))\n",
     "\n",
     "    gwl = Nc.GalaxyWL(wl_dist=grsg, gz_dist=gzgs)\n",
     "    ga.add(gwl)\n",
diff --git a/examples/demo_compute_deltasigma_weights.ipynb b/examples/demo_compute_deltasigma_weights.ipynb
index 2ebdaeb7c..d261446b4 100644
--- a/examples/demo_compute_deltasigma_weights.ipynb
+++ b/examples/demo_compute_deltasigma_weights.ipynb
@@ -55,6 +55,7 @@
     "from clmm import GalaxyCluster\n",
     "from clmm.dataops import compute_galaxy_weights, compute_background_probability\n",
     "from clmm.support import mock_data as mock\n",
+    "from clmm.theory import compute_critical_surface_density_eff\n",
     "\n",
     "clmm.__version__"
    ]
@@ -122,20 +123,7 @@
     "z_gal = noisy_data_z[\"ztrue\"]\n",
     "# add redshift dependency on shape measurement error\n",
     "noisy_data_z[\"e_err\"] = noisy_data_z[\"e_err\"] * (1 + 0.4 * noisy_data_z[\"ztrue\"])\n",
-    "cl0 = GalaxyCluster(\"mock_cluster\", cluster_ra, cluster_dec, cluster_z, noisy_data_z)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "71a6f518-f910-4d5e-b5a5-7c9d82366b86",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "args = (\"mock_cluster\", cluster_ra, cluster_dec, cluster_z, noisy_data_z)\n",
-    "cl1 = GalaxyCluster(*args)\n",
-    "cl2 = GalaxyCluster(*args)\n",
-    "cl3 = GalaxyCluster(*args)"
+    "cluster = GalaxyCluster(\"mock_cluster\", cluster_ra, cluster_dec, cluster_z, noisy_data_z)"
    ]
   },
   {
@@ -163,9 +151,19 @@
    "id": "747e99bf",
    "metadata": {},
    "source": [
-    "#### using the functional interface\n",
+    "#### Using the functional interface\n",
     "\n",
-    "Here, we use the true source redshifts stored in `z_gal`"
+    "First, we need to compute $\\Sigma_{\\rm crit}$. Here, we use the true source redshifts stored in `z_gal`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3432e008",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_c = cosmo.eval_sigma_crit(cluster_z, z_gal)"
    ]
   },
   {
@@ -176,18 +174,9 @@
    "outputs": [],
    "source": [
     "w_ls_true = compute_galaxy_weights(\n",
-    "    cluster_z,\n",
-    "    cosmo,\n",
-    "    z_source=z_gal,\n",
-    "    shape_component1=noisy_data_z[\"e1\"],\n",
-    "    shape_component2=noisy_data_z[\"e2\"],\n",
-    "    shape_component1_err=noisy_data_z[\"e_err\"],\n",
-    "    shape_component2_err=noisy_data_z[\"e_err\"],\n",
-    "    use_pdz=False,\n",
-    "    pzpdf=None,\n",
-    "    pzbins=None,\n",
-    "    use_shape_noise=False,\n",
+    "    sigma_c=sigma_c,\n",
     "    is_deltasigma=True,\n",
+    "    use_shape_noise=False,\n",
     ")"
    ]
   },
@@ -196,7 +185,7 @@
    "id": "206a9562",
    "metadata": {},
    "source": [
-    "#### or as a method of the `GalaxyCluster` object\n",
+    "#### As a method of the `GalaxyCluster` object\n",
     "As a method of the `GalaxyCluster` object, `compute_galaxy_weights` uses the content of the `z` column of the object's `galcat` table as the source redshift point estimates. Given that mock data was generated with a photoz error, this point estimate is different from the true redshift used in the cell above."
    ]
   },
@@ -207,7 +196,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "w_ls_point = cl0.compute_galaxy_weights(\n",
+    "w_ls_point = cluster.compute_galaxy_weights(\n",
     "    use_pdz=False, weight_name=\"w_ls_point\", cosmo=cosmo, is_deltasigma=True, add=True\n",
     ")"
    ]
@@ -225,16 +214,65 @@
    "id": "5f288e2c-32cc-4f35-8744-c66c073b3024",
    "metadata": {},
    "source": [
+    "When considering the photo-z distribution, we can compute the weight based on an effective critical surface density:\n",
+    "\n",
+    "$$\n",
+    "w_{\\rm ls} = \\Sigma_{\\rm crit}^{\\rm eff}(z_l)^{-2},\n",
     "$$\n",
-    "w_{\\rm ls} = \\left(\\int_{z_l}^{+\\infty}\\Sigma_{\\rm crit}(z_l, z)^{-1}p(z) dz\\right)^2\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$$\n",
+    "\\Sigma_{\\rm crit}^{\\rm eff}(z_l) = \\left(\\int_{z_l}^{+\\infty}\\Sigma_{\\rm crit}(z_l, z)^{-1}p(z) dz\\right)^{-1}.\n",
     "$$"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "6355d76b",
+   "metadata": {},
+   "source": [
+    "#### Using the functional interface\n",
+    "\n",
+    "First, we need to compute $\\Sigma_{\\rm crit}^{\\rm eff}$. Here, we use the true source redshifts stored in `z_gal`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0dbbfa25",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_c_eff = compute_critical_surface_density_eff(\n",
+    "    cosmo=cosmo,\n",
+    "    z_cluster=cluster_z,\n",
+    "    pzbins=noisy_data_z.pzpdf_info[\"zbins\"],\n",
+    "    pzpdf=noisy_data_z[\"pzpdf\"],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "aee36658",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "w_ls_photoz = compute_galaxy_weights(\n",
+    "    sigma_c=sigma_c_eff,\n",
+    "    is_deltasigma=True,\n",
+    "    use_shape_noise=False,\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "5285a1a7-e990-4d83-a418-15e04d80a948",
    "metadata": {},
    "source": [
+    "#### As a method of the `GalaxyCluster` object\n",
+    "\n",
     "When used as a method of the `GalaxyCluster` object, with `use_pdz=True`, `compute_galaxy_weights` uses the content of the `pzpdf` and `pzbins` columns of the `galcat` table. It will return an error if those do not exist."
    ]
   },
@@ -245,7 +283,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "w_ls_photoz = cl1.compute_galaxy_weights(\n",
+    "w_ls_photoz = cluster.compute_galaxy_weights(\n",
     "    use_pdz=True, weight_name=\"w_ls_photoz\", cosmo=cosmo, is_deltasigma=True, add=True\n",
     ")"
    ]
@@ -275,7 +313,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "w_ls_shape = cl2.compute_galaxy_weights(\n",
+    "w_ls_shape = cluster.compute_galaxy_weights(\n",
     "    use_pdz=False,\n",
     "    use_shape_noise=True,\n",
     "    shape_component1=\"e1\",\n",
@@ -304,7 +342,7 @@
    "metadata": {},
    "source": [
     "$$\n",
-    "w_{\\rm ls} = \\frac{\\left(\\int_{z_l}^{+\\infty}\\Sigma_{\\rm crit}(z_l, z)^{-1}p(z) dz\\right)^2}{\\sigma_{\\rm shapenoise}^2 + \\sigma^2_{\\rm shape\\ error}}\n",
+    "w_{\\rm ls} = \\frac{\\Sigma_{\\rm crit}^{\\rm eff}(z_l, z)^{-2}}{\\sigma_{\\rm shapenoise}^2 + \\sigma^2_{\\rm shape\\ error}}\n",
     "$$"
    ]
   },
@@ -315,7 +353,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "w_ls_photoz_shape = cl3.compute_galaxy_weights(\n",
+    "w_ls_photoz_shape = cluster.compute_galaxy_weights(\n",
     "    use_pdz=True,\n",
     "    use_shape_noise=True,\n",
     "    shape_component1=\"e1\",\n",
@@ -345,7 +383,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cl0.galcat.colnames"
+    "cluster.galcat.colnames"
    ]
   },
   {
@@ -355,15 +393,15 @@
     "tags": []
    },
    "source": [
-    "### $\\Sigma_c$ metadata\n",
+    "### $\\Sigma_{\\rm crit}$ metadata\n",
     "\n",
-    "As seen above, the weights are computed either using the standard definition of $\\Sigma_c^{\\rm std}$, using point estimate redshift, or the effective definition based on the inverse $\\Sigma_c^{\\rm eff}$ averaged over the redshift PDF $p(z)$:\n",
+    "As seen above, the weights are computed either using the standard definition of $\\Sigma_{\\rm crit}^{\\rm std}$, using point estimate redshift, or the effective definition based on the inverse $\\Sigma_c^{\\rm eff}$ averaged over the redshift PDF $p(z)$:\n",
     "\n",
     "$$ \n",
-    "\\Sigma_c^{\\rm std}\\equiv\\Sigma_{\\rm crit}(z_l, z_s) \\;\\;{\\rm and}\\;\\; \\Sigma_c^{\\rm eff}\\equiv\\left( \\int \\Sigma_{\\rm crit}(z_l, z)^{-1} p(z) dz\\right)^{-1}\n",
+    "\\Sigma_{\\rm crit}^{\\rm std}\\equiv\\Sigma_{\\rm crit}(z_l, z_s) \\;\\;{\\rm and}\\;\\; \\Sigma_{\\rm crit}^{\\rm eff}\\equiv\\left( \\int \\Sigma_{\\rm crit}(z_l, z)^{-1} p(z) dz\\right)^{-1}\n",
     "$$\n",
     "\n",
-    "To keep track of what definition of $\\Sigma_c$ was used, `standard` or `effective` the user may look into the `galcat` metadata "
+    "To keep track of what definition of $\\Sigma_{\\rm crit}$ was used, `standard` or `effective` the user may look into the `galcat` metadata "
    ]
   },
   {
@@ -373,8 +411,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(f\"For cluster 0: {cl0.galcat.meta['sigmac_type']}\")\n",
-    "print(f\"For cluster 1: {cl1.galcat.meta['sigmac_type']}\")"
+    "print(\n",
+    "    \"For cluster w_ls_point column        :\",\n",
+    "    cluster.galcat.meta[\"w_ls_point_sigmac_type\"],\n",
+    ")\n",
+    "print(\n",
+    "    \"For cluster w_ls_photoz column       :\",\n",
+    "    cluster.galcat.meta[\"w_ls_photoz_sigmac_type\"],\n",
+    ")\n",
+    "print(\n",
+    "    \"For cluster w_ls_shape column        :\",\n",
+    "    cluster.galcat.meta[\"w_ls_shape_sigmac_type\"],\n",
+    ")\n",
+    "print(\n",
+    "    \"For cluster w_ls_photoz_shape column :\",\n",
+    "    cluster.galcat.meta[\"w_ls_photoz_shape_sigmac_type\"],\n",
+    ")"
    ]
   },
   {
@@ -405,30 +457,30 @@
     "\n",
     "ax[0].scatter(\n",
     "    z_gal,\n",
-    "    1e3 * cl0.galcat[\"w_ls_point\"] / np.sum(cl0.galcat[\"w_ls_point\"]),\n",
+    "    1e3 * cluster.galcat[\"w_ls_point\"] / cluster.galcat[\"w_ls_point\"].sum(),\n",
     "    c=\"r\",\n",
     "    label=\"point estimate z\",\n",
     "    s=5,\n",
     ")\n",
     "ax[0].scatter(\n",
     "    z_gal,\n",
-    "    1e3 * cl1.galcat[\"w_ls_photoz\"] / np.sum(cl1.galcat[\"w_ls_photoz\"]),\n",
+    "    1e3 * cluster.galcat[\"w_ls_photoz\"] / cluster.galcat[\"w_ls_photoz\"].sum(),\n",
     "    c=\"orange\",\n",
     "    label=\"photoz pdf\",\n",
     "    marker=\"+\",\n",
     "    s=5,\n",
     ")\n",
-    "ax[0].scatter(z_gal, 1e3 * w_ls_true / np.sum(w_ls_true), c=\"g\", label=\"true z\", marker=\"x\", s=5)\n",
+    "ax[0].scatter(z_gal, 1e3 * w_ls_true / w_ls_true.sum(), c=\"g\", label=\"true z\", marker=\"x\", s=5)\n",
     "ax[1].scatter(\n",
     "    z_gal,\n",
-    "    1e3 * cl2.galcat[\"w_ls_shape\"] / np.sum(cl2.galcat[\"w_ls_shape\"]),\n",
+    "    1e3 * cluster.galcat[\"w_ls_shape\"] / cluster.galcat[\"w_ls_shape\"].sum(),\n",
     "    c=\"blue\",\n",
     "    label=r\"shape_error, point z\",\n",
     "    s=5,\n",
     ")\n",
     "ax[1].scatter(\n",
     "    z_gal,\n",
-    "    1e3 * cl3.galcat[\"w_ls_photoz_shape\"] / np.sum(cl3.galcat[\"w_ls_photoz_shape\"]),\n",
+    "    1e3 * cluster.galcat[\"w_ls_photoz_shape\"] / cluster.galcat[\"w_ls_photoz_shape\"].sum(),\n",
     "    c=\"cyan\",\n",
     "    label=\"shape_error + photoz pdf\",\n",
     "    s=5,\n",
@@ -500,7 +552,7 @@
     "    field_size=10.0,\n",
     "    ngals=3000,\n",
     ")\n",
-    "cl0_new = clmm.GalaxyCluster(\"mock_cluster\", 0, 0, cluster_z, noisy_data_z_2)"
+    "cluster_new = clmm.GalaxyCluster(\"mock_cluster\", 0, 0, cluster_z, noisy_data_z_2)"
    ]
   },
   {
@@ -520,7 +572,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "p_background_photoz_1 = cl0.compute_background_probability(\n",
+    "p_background_photoz_1 = cluster.compute_background_probability(\n",
     "    use_pdz=True, p_background_name=\"p_background_photoz_1\", add=True\n",
     ")"
    ]
@@ -532,7 +584,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "p_background_photoz_2 = cl0_new.compute_background_probability(\n",
+    "p_background_photoz_2 = cluster_new.compute_background_probability(\n",
     "    use_pdz=True, p_background_name=\"p_background_photoz_2\", add=True\n",
     ")"
    ]
@@ -549,7 +601,7 @@
     "\n",
     "plt.tick_params(axis=\"both\", which=\"major\", labelsize=10)\n",
     "plt.scatter(\n",
-    "    cl0_new.galcat[\"ztrue\"],\n",
+    "    cluster_new.galcat[\"ztrue\"],\n",
     "    p_background_photoz_2,\n",
     "    c=\"darkorange\",\n",
     "    label=r\"$\\sigma_z = 0.02$\",\n",
@@ -558,7 +610,7 @@
     "    alpha=0.5,\n",
     ")\n",
     "plt.scatter(\n",
-    "    cl0.galcat[\"ztrue\"],\n",
+    "    cluster.galcat[\"ztrue\"],\n",
     "    p_background_photoz_1,\n",
     "    c=\"blue\",\n",
     "    label=r\"$\\sigma_z = 0.05$\",\n",
@@ -574,12 +626,12 @@
     "\n",
     "plt.xlabel(\"true source redshift\", fontsize=15)\n",
     "plt.ylabel(r\"$P(z_s > z_l)$\", fontsize=15)\n",
-    "plt.legend(frameon=True, fontsize=15)\n",
+    "plt.legend(frameon=True, fontsize=15, loc=4)\n",
     "\n",
     "\n",
     "ax2 = fig.add_axes([0.155, 0.65, 0.17, 0.18])\n",
     "ax2.scatter(\n",
-    "    cl0_new.galcat[\"z\"],\n",
+    "    cluster_new.galcat[\"z\"],\n",
     "    p_background_photoz_2,\n",
     "    c=\"darkorange\",\n",
     "    label=r\"$\\sigma_z = 0.02$\",\n",
@@ -589,7 +641,7 @@
     "    lw=0.1,\n",
     ")\n",
     "ax2.scatter(\n",
-    "    cl0.galcat[\"z\"],\n",
+    "    cluster.galcat[\"z\"],\n",
     "    p_background_photoz_1,\n",
     "    c=\"blue\",\n",
     "    label=r\"$\\sigma_z = 0.05$\",\n",
diff --git a/examples/demo_theory_functionality.ipynb b/examples/demo_theory_functionality.ipynb
index 5bb640cce..31d720432 100644
--- a/examples/demo_theory_functionality.ipynb
+++ b/examples/demo_theory_functionality.ipynb
@@ -110,7 +110,7 @@
     "z_cl = 1.0\n",
     "\n",
     "# source properties\n",
-    "z_source = 2.0  # all sources in the same plan\n",
+    "z_src = 2.0  # all sources in the same plan\n",
     "z_distrib_func = clmm.redshift.distributions.chang2013  # sources redshift following a distribution\n",
     "alpha = [2, -0.5]"
    ]
@@ -187,7 +187,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
@@ -206,7 +206,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
@@ -225,7 +225,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
@@ -244,7 +244,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
@@ -264,7 +264,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_source,\n",
+    "    z_src=z_src,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
@@ -290,7 +290,7 @@
     "    mdelta=cluster_mass,\n",
     "    cdelta=cluster_concentration,\n",
     "    z_cluster=z_cl,\n",
-    "    z_source=z_distrib_func,\n",
+    "    z_src=z_distrib_func,\n",
     "    cosmo=cosmo,\n",
     "    delta_mdef=mass_Delta,\n",
     "    halo_profile_model=density_profile_parametrization,\n",
diff --git a/examples/demo_theory_functionality_oo.ipynb b/examples/demo_theory_functionality_oo.ipynb
index 3e8083c3a..2274307df 100644
--- a/examples/demo_theory_functionality_oo.ipynb
+++ b/examples/demo_theory_functionality_oo.ipynb
@@ -105,7 +105,7 @@
     "z_cl = 1.0\n",
     "\n",
     "# source properties\n",
-    "z_source = 2.0  # all sources in the same plane\n",
+    "z_src = 2.0  # all sources in the same plane\n",
     "z_distrib_func = clmm.redshift.distributions.chang2013  # sources redshift following a distribution\n",
     "alpha = [2, -0.5]"
    ]
@@ -127,17 +127,17 @@
     "rho = moo.eval_3d_density(r3d, z_cl)\n",
     "Sigma = moo.eval_surface_density(r3d, z_cl)\n",
     "DeltaSigma = moo.eval_excess_surface_density(r3d, z_cl)\n",
-    "gammat = moo.eval_tangential_shear(r3d, z_cl, z_source)\n",
-    "kappa = moo.eval_convergence(r3d, z_cl, z_source)\n",
+    "gammat = moo.eval_tangential_shear(r3d, z_cl, z_src)\n",
+    "kappa = moo.eval_convergence(r3d, z_cl, z_src)\n",
     "\n",
-    "gt = moo.eval_reduced_tangential_shear(r3d, z_cl, z_source)\n",
+    "gt = moo.eval_reduced_tangential_shear(r3d, z_cl, z_src)\n",
     "# Lensing quantities assuming sources follow a given redshift distribution.\n",
     "gt_z = moo.eval_reduced_tangential_shear(\n",
     "    r3d, z_cl, z_distrib_func, z_src_info=\"distribution\", approx=\"order2\"\n",
     ")\n",
     "\n",
-    "mu = moo.eval_magnification(r3d, z_cl, z_source)\n",
-    "mu_bias = moo.eval_magnification_bias(r3d, z_cl, z_source, alpha)"
+    "mu = moo.eval_magnification(r3d, z_cl, z_src)\n",
+    "mu_bias = moo.eval_magnification_bias(r3d, z_cl, z_src, alpha)"
    ]
   },
   {
diff --git a/examples/mass_fitting/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb b/examples/mass_fitting/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
index e3efcfc9b..0cbce9bb9 100644
--- a/examples/mass_fitting/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
+++ b/examples/mass_fitting/Example3_Fit_Halo_Mass_to_Shear_Catalog.ipynb
@@ -452,7 +452,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=concentration,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=clmm.redshift.distributions.chang2013,  # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
+    "        z_src=clmm.redshift.distributions.chang2013,  # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        # massdef='critical',\n",
@@ -473,7 +473,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=concentration,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=clmm.redshift.distributions.chang2013,  # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
+    "        z_src=clmm.redshift.distributions.chang2013,  # We provide the redshift distribution (default: Chang et al. 2013) for calculating the beta_s statistics\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        # massdef='critical',\n",
diff --git a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb
index 45a76d5ec..0c423d25b 100644
--- a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb
+++ b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb
@@ -476,7 +476,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=concentration,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=(bs_mean, bs2_mean),  # tuple of (bs_mean, bs2_mean)\n",
+    "        z_src=(bs_mean, bs2_mean),  # tuple of (bs_mean, bs2_mean)\n",
     "        z_src_info=\"beta\",\n",
     "        approx=\"order1\",\n",
     "        cosmo=cosmo,\n",
@@ -504,7 +504,7 @@
     "                    cdelta=concentration,  # Concentration of the cluster\n",
     "                    z_cluster=cluster_z,  # Redshift of the cluster\n",
     "                    # Redshift value of each source galaxy inside the radial bin\n",
-    "                    z_source=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
+    "                    z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
     "                    cosmo=cosmo,\n",
     "                    delta_mdef=200,\n",
     "                    massdef=\"critical\",\n",
diff --git a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb
index a1439d1f7..3e24d6f69 100644
--- a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb
+++ b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb
@@ -404,7 +404,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=concentration,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=(bs_mean, bs2_mean),  # tuple of (bs_mean, bs2_mean)\n",
+    "        z_src=(bs_mean, bs2_mean),  # tuple of (bs_mean, bs2_mean)\n",
     "        z_src_info=\"beta\",\n",
     "        approx=\"order1\",\n",
     "        cosmo=cosmo,\n",
@@ -432,7 +432,7 @@
     "                    cdelta=concentration,  # Concentration of the cluster\n",
     "                    z_cluster=cluster_z,  # Redshift of the cluster\n",
     "                    # Redshift value of each source galaxy inside the radial bin\n",
-    "                    z_source=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
+    "                    z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
     "                    cosmo=cosmo,\n",
     "                    delta_mdef=200,\n",
     "                    massdef=\"critical\",\n",
diff --git a/examples/other_flat_sky_limitations.ipynb b/examples/other_flat_sky_limitations.ipynb
index 75a9b2ea6..4447c592b 100644
--- a/examples/other_flat_sky_limitations.ipynb
+++ b/examples/other_flat_sky_limitations.ipynb
@@ -371,7 +371,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=np.mean(cluster.galcat[\"z\"]),  # Mean value of source galaxies redshift\n",
+    "        z_src=np.mean(cluster.galcat[\"z\"]),  # Mean value of source galaxies redshift\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -415,7 +415,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=z_inf,  # Redshift value at infinity\n",
+    "        z_src=z_inf,  # Redshift value at infinity\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -425,7 +425,7 @@
     "        mdelta=10**logm,  # Mass of the cluster [M_sun]\n",
     "        cdelta=4,  # Concentration of the cluster\n",
     "        z_cluster=cluster_z,  # Redshift of the cluster\n",
-    "        z_source=z_inf,  # Redshift value at infinity\n",
+    "        z_src=z_inf,  # Redshift value at infinity\n",
     "        cosmo=cosmo,\n",
     "        delta_mdef=200,\n",
     "        halo_profile_model=\"nfw\",\n",
@@ -467,7 +467,7 @@
     "                    cdelta=4,  # Concentration of the cluster\n",
     "                    z_cluster=cluster_z,  # Redshift of the cluster\n",
     "                    # Redshift value of each source galaxy inside the radial bin\n",
-    "                    z_source=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
+    "                    z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n",
     "                    cosmo=cosmo,\n",
     "                    delta_mdef=200,\n",
     "                    halo_profile_model=\"nfw\",\n",
diff --git a/tests/data/numcosmo/config.json b/tests/data/numcosmo/config.json
index c81716df1..c38ce5e87 100644
--- a/tests/data/numcosmo/config.json
+++ b/tests/data/numcosmo/config.json
@@ -1 +1 @@
-{"dl": 1696.50282336234, "z_source": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "nfw", "ds": 1767.4306833378118, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962517, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435031.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
\ No newline at end of file
+{"dl": 1696.50282336234, "z_src": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "nfw", "ds": 1767.4306833378118, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962517, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435031.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
diff --git a/tests/data/numcosmo/config_einasto_benchmarks.json b/tests/data/numcosmo/config_einasto_benchmarks.json
index 0d6f949fb..dd201b9aa 100644
--- a/tests/data/numcosmo/config_einasto_benchmarks.json
+++ b/tests/data/numcosmo/config_einasto_benchmarks.json
@@ -1 +1 @@
-{"dl": 1696.502823362356, "z_source": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "einasto", "ds": 1767.4306833377823, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962118, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435132.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16, "alpha_einasto": 0.3683247081598754}
\ No newline at end of file
+{"dl": 1696.502823362356, "z_src": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "einasto", "ds": 1767.4306833377823, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962118, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435132.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16, "alpha_einasto": 0.3683247081598754}
diff --git a/tests/data/numcosmo/config_hernquist_benchmarks.json b/tests/data/numcosmo/config_hernquist_benchmarks.json
index bc425840c..59eb97c68 100644
--- a/tests/data/numcosmo/config_hernquist_benchmarks.json
+++ b/tests/data/numcosmo/config_hernquist_benchmarks.json
@@ -1 +1 @@
-{"dl": 1696.502823362356, "z_source": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "hernquist", "ds": 1767.4306833377823, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962118, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435132.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
\ No newline at end of file
+{"dl": 1696.502823362356, "z_src": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "aexp_cluster": 0.5, "density_profile_parametrization": "hernquist", "ds": 1767.4306833377823, "cosmo_Ob0": 0.049, "aexp_source": 0.3333333333333333, "dsl": 636.4288010962118, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 67.66, "cluster_concentration": 4.0, "cosmo_Om0": 0.31240715689016263, "cosmo_Odm0": 0.262, "nc_Sigmac": 2722127302435132.5, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
diff --git a/tests/data/numcosmo/config_old_units.json b/tests/data/numcosmo/config_old_units.json
index 42ae6f974..b483a268d 100644
--- a/tests/data/numcosmo/config_old_units.json
+++ b/tests/data/numcosmo/config_old_units.json
@@ -1 +1 @@
-{"dl": 1177679548.8532429, "z_source": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "nc_cor_factor": 1.0, "aexp_cluster": 0.5, "density_profile_parametrization": "nfw", "ds": 1241373844.9905102, "cosmo_Ob0": 0.045, "aexp_source": 0.3333333333333333, "dsl": 456254145.7550621, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 70.0, "cluster_concentration": 4.0, "cosmo_Om0": 0.27, "nc_Sigmac": 3841.8383445955924, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
\ No newline at end of file
+{"dl": 1177679548.8532429, "z_src": 2.0, "cluster_mass": 1000000000000000.0, "G[m3/km.s2]": 6.6743e-11, "nc_cor_factor": 1.0, "aexp_cluster": 0.5, "density_profile_parametrization": "nfw", "ds": 1241373844.9905102, "cosmo_Ob0": 0.045, "aexp_source": 0.3333333333333333, "dsl": 456254145.7550621, "z_cluster": 1.0, "mass_Delta": 200, "lightspeed[km/s]": 299792.458, "cosmo_H0": 70.0, "cluster_concentration": 4.0, "cosmo_Om0": 0.27, "nc_Sigmac": 3841.8383445955924, "Msun[kg]": 1.988409870698051e+30, "pc_to_m": 3.085677581491367e+16}
diff --git a/tests/data/numcosmo/generate_einasto_hernquist_benchmarks_fromNC.ipynb b/tests/data/numcosmo/generate_einasto_hernquist_benchmarks_fromNC.ipynb
index aeaf5e279..fbd7f9f4d 100644
--- a/tests/data/numcosmo/generate_einasto_hernquist_benchmarks_fromNC.ipynb
+++ b/tests/data/numcosmo/generate_einasto_hernquist_benchmarks_fromNC.ipynb
@@ -107,7 +107,7 @@
     "cluster_concentration = 4.0\n",
     "z_max = 3.0\n",
     "z_cluster = 1.0\n",
-    "z_source = 2.0"
+    "z_src = 2.0"
    ]
   },
   {
@@ -487,23 +487,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nc_Sigmac = smd.sigma_critical(cosmo, z_source, z_cluster, z_cluster)\n",
+    "nc_Sigmac = smd.sigma_critical(cosmo, z_src, z_cluster, z_cluster)\n",
     "\n",
     "nc_rho_ein = np.array([nc_ein.eval_density(cosmo, r3d_i, z_cluster) for r3d_i in r3d])\n",
     "nc_Sigma_ein = np.array([smd.sigma(nc_ein, cosmo, r3d_i, z_cluster) for r3d_i in r3d])\n",
     "nc_SigmaMean_ein = np.array([smd.sigma_mean(nc_ein, cosmo, r3d_i, z_cluster) for r3d_i in r3d])\n",
     "nc_DeltaSigma_ein = nc_SigmaMean_ein - nc_Sigma_ein\n",
     "nc_gammat_ein = np.array(\n",
-    "    [smd.shear(nc_ein, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.shear(nc_ein, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_kappa_ein = np.array(\n",
-    "    [smd.convergence(nc_ein, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.convergence(nc_ein, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_gt_ein = np.array(\n",
-    "    [smd.reduced_shear(nc_ein, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.reduced_shear(nc_ein, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_mu_ein = np.array(\n",
-    "    [smd.magnification(nc_ein, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.magnification(nc_ein, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")"
    ]
   },
@@ -525,16 +525,16 @@
     "nc_SigmaMean_her = np.array([smd.sigma_mean(nc_her, cosmo, r3d_i, z_cluster) for r3d_i in r3d])\n",
     "nc_DeltaSigma_her = nc_SigmaMean_her - nc_Sigma_her\n",
     "nc_gammat_her = np.array(\n",
-    "    [smd.shear(nc_her, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.shear(nc_her, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_kappa_her = np.array(\n",
-    "    [smd.convergence(nc_her, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.convergence(nc_her, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_gt_her = np.array(\n",
-    "    [smd.reduced_shear(nc_her, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.reduced_shear(nc_her, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")\n",
     "nc_mu_her = np.array(\n",
-    "    [smd.magnification(nc_her, cosmo, r3d_i, z_source, z_cluster, z_cluster) for r3d_i in r3d]\n",
+    "    [smd.magnification(nc_her, cosmo, r3d_i, z_src, z_cluster, z_cluster) for r3d_i in r3d]\n",
     ")"
    ]
   },
@@ -601,15 +601,15 @@
    "source": [
     "config_dict_ein = {\n",
     "    \"dl\": dist.angular_diameter(cosmo, z_cluster) * cosmo.RH_Mpc(),\n",
-    "    \"z_source\": z_source,\n",
+    "    \"z_src\": z_src,\n",
     "    \"cluster_mass\": cluster_mass,\n",
     "    \"G[m3/km.s2]\": Ncm.C.G(),\n",
     "    \"aexp_cluster\": 1.0 / (1.0 + z_cluster),\n",
     "    \"density_profile_parametrization\": \"einasto\",\n",
-    "    \"ds\": dist.angular_diameter(cosmo, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"ds\": dist.angular_diameter(cosmo, z_src) * cosmo.RH_Mpc(),\n",
     "    \"cosmo_Ob0\": cosmo_ccl[\"Omega_b0\"],\n",
-    "    \"aexp_source\": 1.0 / (1.0 + z_source),\n",
-    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"aexp_source\": 1.0 / (1.0 + z_src),\n",
+    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_src) * cosmo.RH_Mpc(),\n",
     "    \"z_cluster\": z_cluster,\n",
     "    \"mass_Delta\": mass_Delta,\n",
     "    \"lightspeed[km/s]\": Ncm.C.c() / 1000.0,\n",
@@ -637,15 +637,15 @@
    "source": [
     "config_dict_her = {\n",
     "    \"dl\": dist.angular_diameter(cosmo, z_cluster) * cosmo.RH_Mpc(),\n",
-    "    \"z_source\": z_source,\n",
+    "    \"z_src\": z_src,\n",
     "    \"cluster_mass\": cluster_mass,\n",
     "    \"G[m3/km.s2]\": Ncm.C.G(),\n",
     "    \"aexp_cluster\": 1.0 / (1.0 + z_cluster),\n",
     "    \"density_profile_parametrization\": \"hernquist\",\n",
-    "    \"ds\": dist.angular_diameter(cosmo, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"ds\": dist.angular_diameter(cosmo, z_src) * cosmo.RH_Mpc(),\n",
     "    \"cosmo_Ob0\": cosmo_ccl[\"Omega_b0\"],\n",
-    "    \"aexp_source\": 1.0 / (1.0 + z_source),\n",
-    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_source) * cosmo.RH_Mpc(),\n",
+    "    \"aexp_source\": 1.0 / (1.0 + z_src),\n",
+    "    \"dsl\": dist.angular_diameter_z1_z2(cosmo, z_cluster, z_src) * cosmo.RH_Mpc(),\n",
     "    \"z_cluster\": z_cluster,\n",
     "    \"mass_Delta\": mass_Delta,\n",
     "    \"lightspeed[km/s]\": Ncm.C.c() / 1000.0,\n",
diff --git a/tests/test_clusterensemble.py b/tests/test_clusterensemble.py
index 7751d7bb6..63c3d3958 100644
--- a/tests/test_clusterensemble.py
+++ b/tests/test_clusterensemble.py
@@ -23,7 +23,7 @@ def test_cluster_ensemble():
     ra_source = [120.1, 119.9]
     dec_source = [41.9, 42.2]
     theta_source = [0.0025, 0.015]
-    z_source = [1.0, 2.0]
+    z_src = [1.0, 2.0]
     shear1 = [0.2, 0.4]
     shear2 = [0.3, 0.5]
     w_ls = [1.0e-30, 1.0e-31]
@@ -33,7 +33,7 @@ def test_cluster_ensemble():
     names = ("ra", "dec", "theta", "w_ls", "e1", "e2", "z")
 
     galcat = clmm.GCData(
-        [ra_source, dec_source, theta_source, w_ls, shear1, shear2, z_source], names=names
+        [ra_source, dec_source, theta_source, w_ls, shear1, shear2, z_src], names=names
     )
     # create cluster
     cluster = clmm.GalaxyCluster(
diff --git a/tests/test_cosmo_parent.py b/tests/test_cosmo_parent.py
index b65666e97..67de814db 100644
--- a/tests/test_cosmo_parent.py
+++ b/tests/test_cosmo_parent.py
@@ -269,7 +269,7 @@ def test_eval_sigma_crit(modeling_data):
     cosmo, testcase, _ = load_validation_config()
 
     assert_allclose(
-        cosmo.eval_sigma_crit(testcase["z_cluster"], testcase["z_source"]),
+        cosmo.eval_sigma_crit(testcase["z_cluster"], testcase["z_src"]),
         testcase["nc_Sigmac"],
         reltol,
     )
@@ -278,7 +278,7 @@ def test_eval_sigma_crit(modeling_data):
     assert_raises(ValueError, cosmo.eval_sigma_crit, 0.2, -0.3)
     # Check behaviour when sources are in front of the lens
     z_cluster = 0.3
-    z_source = 0.2
-    assert_allclose(cosmo.eval_sigma_crit(z_cluster, z_source), np.inf, 1.0e-10)
-    z_source = [0.2, 0.12, 0.25]
-    assert_allclose(cosmo.eval_sigma_crit(z_cluster, z_source), [np.inf, np.inf, np.inf], 1.0e-10)
+    z_src = 0.2
+    assert_allclose(cosmo.eval_sigma_crit(z_cluster, z_src), np.inf, 1.0e-10)
+    z_src = [0.2, 0.12, 0.25]
+    assert_allclose(cosmo.eval_sigma_crit(z_cluster, z_src), [np.inf, np.inf, np.inf], 1.0e-10)
diff --git a/tests/test_dataops.py b/tests/test_dataops.py
index aae203de5..f60ac795b 100644
--- a/tests/test_dataops.py
+++ b/tests/test_dataops.py
@@ -6,6 +6,7 @@
 from clmm import GCData
 from scipy.stats import multivariate_normal
 import clmm.dataops as da
+from clmm.theory import compute_critical_surface_density_eff
 
 TOLERANCE = {"rtol": 1.0e-7, "atol": 1.0e-7}
 
@@ -444,35 +445,9 @@ def test_compute_tangential_and_cross_components(modeling_data):
     # Check behaviour for the deltasigma option.
     cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.275, Omega_b0=0.025)
 
-    # test missing info for is_deltasigma=True
-    assert_raises(
-        TypeError,
-        da.compute_tangential_and_cross_components,
-        ra_lens=ra_lens,
-        dec_lens=dec_lens,
-        ra_source=gals["ra"],
-        dec_source=gals["dec"],
-        shear1=gals["e1"],
-        shear2=gals["e2"],
-        is_deltasigma=True,
-        cosmo=None,
-        z_lens=z_lens,
-        z_source=gals["z"],
-    )
-    assert_raises(
-        TypeError,
-        da.compute_tangential_and_cross_components,
-        ra_lens=ra_lens,
-        dec_lens=dec_lens,
-        ra_source=gals["ra"],
-        dec_source=gals["dec"],
-        shear1=gals["e1"],
-        shear2=gals["e2"],
-        is_deltasigma=True,
-        cosmo=cosmo,
-        z_lens=None,
-        z_source=gals["z"],
-    )
+    # check values for DeltaSigma
+    sigma_c = cosmo.eval_sigma_crit(z_lens, gals["z"])
+    # check validation between is_deltasigma and sigma_c
     assert_raises(
         TypeError,
         da.compute_tangential_and_cross_components,
@@ -482,13 +457,9 @@ def test_compute_tangential_and_cross_components(modeling_data):
         dec_source=gals["dec"],
         shear1=gals["e1"],
         shear2=gals["e2"],
-        is_deltasigma=True,
-        cosmo=cosmo,
-        z_lens=z_lens,
-        z_source=None,
+        is_deltasigma=False,
+        sigma_c=sigma_c,
     )
-
-    # test  missing info for use_pdz=True
     assert_raises(
         TypeError,
         da.compute_tangential_and_cross_components,
@@ -499,30 +470,9 @@ def test_compute_tangential_and_cross_components(modeling_data):
         shear1=gals["e1"],
         shear2=gals["e2"],
         is_deltasigma=True,
-        use_pdz=True,
-        cosmo=cosmo,
-        z_lens=z_lens,
-        z_source=None,
+        sigma_c=None,
     )
-
-    # Trying to got through line 173 of dataops/__init__.py
-    da.compute_tangential_and_cross_components(
-        ra_lens=ra_lens,
-        dec_lens=dec_lens,
-        ra_source=gals["ra"][0],
-        dec_source=gals["dec"][0],
-        shear1=gals["e1"][0],
-        shear2=gals["e2"][0],
-        is_deltasigma=True,
-        use_pdz=True,
-        cosmo=cosmo,
-        z_lens=z_lens,
-        z_source=None,
-        pzbins=[[0.55, 0.6, 0.65, 0.7, 0.75]],
-        pzpdf=[[0.01, 1, 0.01, 0.001, 0.0001]],
-    )
-
-    # check values for DeltaSigma
+    # test values
     for geometry, expected in geo_tests:
         angsep_DS, tDS, xDS = da.compute_tangential_and_cross_components(
             ra_lens=ra_lens,
@@ -532,9 +482,7 @@ def test_compute_tangential_and_cross_components(modeling_data):
             shear1=gals["e1"],
             shear2=gals["e2"],
             is_deltasigma=True,
-            cosmo=cosmo,
-            z_lens=z_lens,
-            z_source=gals["z"],
+            sigma_c=sigma_c,
             geometry=geometry,
         )
         assert_allclose(
@@ -583,15 +531,20 @@ def test_compute_tangential_and_cross_components(modeling_data):
             err_msg="Cross Shear not correct when using cluster method",
         )
 
+    # test basic weights functionality
+    cluster.compute_galaxy_weights()
+    expected = np.array([1.0, 1.0])
+    assert_allclose(cluster.galcat["w_ls"], expected, **TOLERANCE)
+
 
 def test_compute_background_probability():
     """test for compute background probability"""
     z_lens = 0.1
-    z_source = np.array([0.22, 0.35, 1.7])
+    z_src = np.array([0.22, 0.35, 1.7])
 
     # true redshift
     p_bkg = da.compute_background_probability(
-        z_lens, z_source=z_source, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True
+        z_lens, z_src=z_src, use_pdz=False, pzpdf=None, pzbins=None, validate_input=True
     )
     expected = np.array([1.0, 1.0, 1.0])
     assert_allclose(p_bkg, expected, **TOLERANCE)
@@ -599,7 +552,7 @@ def test_compute_background_probability():
         ValueError,
         da.compute_background_probability,
         z_lens,
-        z_source=None,
+        z_src=None,
         use_pdz=False,
         pzpdf=None,
         pzbins=None,
@@ -608,13 +561,13 @@ def test_compute_background_probability():
 
     # photoz + deltasigma
     pzbin = np.linspace(0.0001, 5, 100)
-    pzbins = [pzbin for i in range(z_source.size)]
-    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_source]
+    pzbins = [pzbin for i in range(z_src.size)]
+    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_src]
     assert_raises(
         ValueError,
         da.compute_background_probability,
         z_lens,
-        z_source=z_source,
+        z_src=z_src,
         use_pdz=True,
         pzpdf=None,
         pzbins=pzbins,
@@ -626,27 +579,23 @@ def test_compute_galaxy_weights():
     """test for compute galaxy weights"""
     cosmo = clmm.Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0)
     z_lens = 0.1
-    z_source = [0.22, 0.35, 1.7]
+    z_src = [0.22, 0.35, 1.7]
     shape_component1 = np.array([0.143, 0.063, -0.171])
     shape_component2 = np.array([-0.011, 0.012, -0.250])
     shape_component1_err = np.array([0.11, 0.01, 0.2])
     shape_component2_err = np.array([0.14, 0.16, 0.21])
 
     # true redshift + deltasigma
+    sigma_c = cosmo.eval_sigma_crit(z_lens, z_src)
     weights = da.compute_galaxy_weights(
-        z_lens,
-        cosmo,
-        z_source=z_source,
-        use_pdz=False,
-        pzpdf=None,
-        pzbins=None,
+        is_deltasigma=True,
+        sigma_c=sigma_c,
         use_shape_noise=False,
         shape_component1=shape_component1,
         shape_component2=shape_component2,
         use_shape_error=False,
         shape_component1_err=shape_component1_err,
         shape_component2_err=shape_component2_err,
-        is_deltasigma=True,
         validate_input=True,
     )
     expected = np.array([4.58644320e-31, 9.68145632e-31, 5.07260777e-31])
@@ -654,22 +603,23 @@ def test_compute_galaxy_weights():
 
     # photoz + deltasigma
     pzbin = np.linspace(0.0001, 5, 100)
-    pzbins = [pzbin for i in range(len(z_source))]
-    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_source]
-    weights = da.compute_galaxy_weights(
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
+    pzbins = [pzbin for i in range(len(z_src))]
+    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_src]
+    sigma_c_eff = compute_critical_surface_density_eff(
+        cosmo=cosmo,
+        z_cluster=z_lens,
         pzbins=pzbins,
+        pzpdf=pzpdf,
+    )
+    weights = da.compute_galaxy_weights(
+        is_deltasigma=True,
+        sigma_c=sigma_c_eff,
         use_shape_noise=False,
         shape_component1=shape_component1,
         shape_component2=shape_component2,
         use_shape_error=False,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=True,
         validate_input=True,
     )
 
@@ -679,21 +629,16 @@ def test_compute_galaxy_weights():
     # photoz + deltasigma - shared bins
     pzbin = np.linspace(0.0001, 5, 100)
     pzbins = pzbin
-    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_source]
+    pzpdf = [multivariate_normal.pdf(pzbin, mean=z, cov=0.3) for z in z_src]
     weights = da.compute_galaxy_weights(
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
+        is_deltasigma=True,
+        sigma_c=sigma_c_eff,
         use_shape_noise=False,
         shape_component1=shape_component1,
         shape_component2=shape_component2,
         use_shape_error=False,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=True,
         validate_input=True,
     )
 
@@ -702,19 +647,14 @@ def test_compute_galaxy_weights():
 
     # test with noise
     weights = da.compute_galaxy_weights(
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
+        is_deltasigma=True,
+        sigma_c=sigma_c_eff,
         use_shape_noise=True,
         shape_component1=shape_component1,
         shape_component2=shape_component2,
         use_shape_error=False,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=True,
         validate_input=True,
     )
 
@@ -723,19 +663,14 @@ def test_compute_galaxy_weights():
 
     # test with is_deltasigma=False and geometric weights only
     weights = da.compute_galaxy_weights(
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
+        is_deltasigma=False,
+        sigma_c=None,
         use_shape_noise=False,
         shape_component1=shape_component1,
         shape_component2=shape_component2,
         use_shape_error=False,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=False,
         validate_input=True,
     )
 
@@ -746,19 +681,14 @@ def test_compute_galaxy_weights():
     assert_raises(
         ValueError,
         da.compute_galaxy_weights,
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
+        is_deltasigma=True,
+        sigma_c=sigma_c_eff,
         use_shape_noise=True,
         shape_component1=None,
         shape_component2=None,
         use_shape_error=False,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=False,
         validate_input=True,
     )
 
@@ -766,52 +696,12 @@ def test_compute_galaxy_weights():
     assert_raises(
         ValueError,
         da.compute_galaxy_weights,
-        z_lens,
-        cosmo,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
-        use_shape_noise=False,
-        use_shape_error=True,
-        shape_component1_err=None,
-        shape_component2_err=None,
-        is_deltasigma=False,
-        validate_input=True,
-    )
-
-    assert_raises(
-        TypeError,
-        da.compute_galaxy_weights,
-        z_lens,
-        cosmo=None,
-        z_source=None,
-        use_pdz=False,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
-        use_shape_noise=True,
-        shape_component1=None,
-        shape_component2=None,
-        use_shape_error=False,
-        shape_component1_err=None,
-        shape_component2_err=None,
         is_deltasigma=True,
-        validate_input=True,
-    )
-    assert_raises(
-        TypeError,
-        da.compute_galaxy_weights,
-        z_lens,
-        cosmo=None,
-        z_source=None,
-        use_pdz=True,
-        pzpdf=pzpdf,
-        pzbins=pzbins,
+        sigma_c=sigma_c_eff,
         use_shape_noise=False,
-        use_shape_error=False,
+        use_shape_error=True,
         shape_component1_err=None,
         shape_component2_err=None,
-        is_deltasigma=True,
         validate_input=True,
     )
 
diff --git a/tests/test_galaxycluster.py b/tests/test_galaxycluster.py
index a3987255f..04d73c6dc 100644
--- a/tests/test_galaxycluster.py
+++ b/tests/test_galaxycluster.py
@@ -140,12 +140,12 @@ def test_print_gc():
 def test_integrity_of_lensfuncs():
     """test integrity of lensfuncs"""
     ra_source, dec_source = [120.1, 119.9, 119.9], [41.9, 42.2, 42.2]
-    id_source, z_source = [1, 2, 3], [1, 1, 1]
+    id_source, z_src = [1, 2, 3], [1, 1, 1]
     shape_component1 = np.array([0.143, 0.063, -0.171])
     shape_component2 = np.array([-0.011, 0.012, -0.250])
 
     galcat = GCData(
-        [ra_source, dec_source, z_source, id_source, shape_component1, shape_component2],
+        [ra_source, dec_source, z_src, id_source, shape_component1, shape_component2],
         names=("ra", "dec", "z", "id", "e1", "e2"),
     )
     galcat_noz = GCData([ra_source, dec_source, id_source], names=("ra", "dec", "id"))
@@ -168,8 +168,8 @@ def test_integrity_of_lensfuncs():
     pzbins = np.linspace(0.0001, 5, 100)
     cluster = clmm.GalaxyCluster(unique_id="1", ra=161.3, dec=34.0, z=0.3, galcat=galcat)
     cluster.galcat.pzpdf_info["zbins"] = pzbins
-    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_source))]
-    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_source]
+    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
+    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_src]
 
     for pztype in ("individual_bins", "shared_bins"):
         cluster.galcat.pzpdf_info["type"] = pztype
@@ -177,21 +177,20 @@ def test_integrity_of_lensfuncs():
         cluster.compute_tangential_and_cross_components(
             is_deltasigma=True, use_pdz=True, cosmo=cosmo, add=True
         )
-        assert_equal(cluster.galcat.meta["sigmac_type"], "effective")
+        for comp_name in ("et", "ex"):
+            assert_equal(cluster.galcat.meta[f"{comp_name}_sigmac_type"], "effective")
 
 
 def test_integrity_of_probfuncs():
     """test integrity of prob funcs"""
     ra_source, dec_source = [120.1, 119.9, 119.9], [41.9, 42.2, 42.2]
-    id_source, z_sources = [1, 2, 3], [1, 1, 1]
+    id_source, z_srcs = [1, 2, 3], [1, 1, 1]
     cluster = clmm.GalaxyCluster(
         unique_id="1",
         ra=161.3,
         dec=34.0,
         z=0.3,
-        galcat=GCData(
-            [ra_source, dec_source, z_sources, id_source], names=("ra", "dec", "z", "id")
-        ),
+        galcat=GCData([ra_source, dec_source, z_srcs, id_source], names=("ra", "dec", "z", "id")),
     )
     # true redshift
     cluster.compute_background_probability(use_pdz=False, p_background_name="p_bkg_true")
@@ -202,8 +201,8 @@ def test_integrity_of_probfuncs():
     assert_raises(TypeError, cluster.compute_background_probability, use_photoz=True)
     pzbins = np.linspace(0.0001, 5, 1000)
     cluster.galcat.pzpdf_info["zbins"] = pzbins
-    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_sources))]
-    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.01) for z in z_sources]
+    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_srcs))]
+    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.01) for z in z_srcs]
     for pztype in ("individual_bins", "shared_bins"):
         cluster.galcat.pzpdf_info["type"] = pztype
         cluster.compute_background_probability(use_pdz=True, p_background_name="p_bkg_pz")
@@ -214,7 +213,7 @@ def test_integrity_of_weightfuncs():
     """test integrity of weight funcs"""
     cosmo = clmm.Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0)
     z_lens = 0.1
-    z_source = [0.22, 0.35, 1.7]
+    z_src = [0.22, 0.35, 1.7]
     shape_component1 = np.array([0.143, 0.063, -0.171])
     shape_component2 = np.array([-0.011, 0.012, -0.250])
     shape_component1_err = np.array([0.11, 0.01, 0.2])
@@ -231,7 +230,7 @@ def test_integrity_of_weightfuncs():
                 shape_component2,
                 shape_component1_err,
                 shape_component2_err,
-                z_source,
+                z_src,
             ],
             names=("e1", "e2", "e1_err", "e2_err", "z"),
         ),
@@ -245,8 +244,8 @@ def test_integrity_of_weightfuncs():
     # photoz + deltasigma
     pzbins = np.linspace(0.0001, 5, 100)
     cluster.galcat.pzpdf_info["zbins"] = pzbins
-    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_source))]
-    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_source]
+    cluster.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
+    cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_src]
     for pztype in ("individual_bins", "shared_bins"):
         cluster.galcat.pzpdf_info["type"] = pztype
         cluster.compute_galaxy_weights(
@@ -271,11 +270,11 @@ def test_integrity_of_weightfuncs():
 def test_pzpdf_random_draw():
     """test draw_gal_z_from_pdz"""
     z_lens = 0.1
-    z_source = [0.22, 0.35, 1.7]
+    z_src = [0.22, 0.35, 1.7]
     shape_component1 = np.array([0.143, 0.063, -0.171])
     shape_component2 = np.array([-0.011, 0.012, -0.250])
     cluster_kwargs = dict(unique_id="1", ra=161.3, dec=34.0, z=z_lens)
-    gcat_args = [shape_component1, shape_component2, z_source]
+    gcat_args = [shape_component1, shape_component2, z_src]
     gcat_kwargs = {"names": ("e1", "e2", "z")}
 
     # set up photoz
@@ -288,19 +287,17 @@ def test_pzpdf_random_draw():
         assert_raises(TypeError, cluster.draw_gal_z_from_pdz)
 
         cluster.galcat.pzpdf_info["zbins"] = pzbins
-        cluster.galcat["pzbins"] = [pzbins for i in range(len(z_source))]
+        cluster.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
         assert_raises(TypeError, cluster.draw_gal_z_from_pdz)
 
         cluster.galcat.pzpdf_info.pop("zbins")
         cluster.galcat.remove_column("pzbins")
-        cluster.galcat["pzpdf"] = [
-            multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_source
-        ]
+        cluster.galcat["pzpdf"] = [multivariate_normal.pdf(pzbins, mean=z, cov=0.3) for z in z_src]
         assert_raises(TypeError, cluster.draw_gal_z_from_pdz)
 
         # add pzbins back to galcat
         cluster.galcat.pzpdf_info["zbins"] = pzbins
-        cluster.galcat["pzbins"] = [pzbins for i in range(len(z_source))]
+        cluster.galcat["pzbins"] = [pzbins for i in range(len(z_src))]
         # test raising TypeError when the name of the new column is already in cluster.galcat
         # also test default overwrite=False and zcol_out='z'
         assert_raises(TypeError, cluster.draw_gal_z_from_pdz)
@@ -349,7 +346,7 @@ def test_plot_profiles():
     ra_lens, dec_lens, z_lens = 120.0, 42.0, 0.5
     ra_source = [120.1, 119.9]
     dec_source = [41.9, 42.2]
-    z_source = [1.0, 2.0]
+    z_src = [1.0, 2.0]
     shear1 = [0.2, 0.4]
     shear2 = [0.3, 0.5]
     # Set up radial values
@@ -362,7 +359,7 @@ def test_plot_profiles():
         dec=dec_lens,
         z=z_lens,
         galcat=GCData(
-            [ra_source, dec_source, shear1, shear2, z_source], names=("ra", "dec", "e1", "e2", "z")
+            [ra_source, dec_source, shear1, shear2, z_src], names=("ra", "dec", "e1", "e2", "z")
         ),
     )
     cluster.compute_tangential_and_cross_components()
diff --git a/tests/test_theory.py b/tests/test_theory.py
index 9c482317c..ad8562b33 100644
--- a/tests/test_theory.py
+++ b/tests/test_theory.py
@@ -94,7 +94,7 @@ def load_validation_config(halo_profile_model=None):
         "mdelta": testcase["cluster_mass"],
         "cdelta": testcase["cluster_concentration"],
         "z_cluster": testcase["z_cluster"],
-        "z_source": testcase["z_source"],
+        "z_src": testcase["z_src"],
         "delta_mdef": testcase["mass_Delta"],
         "halo_profile_model": testcase["density_profile_parametrization"],
         "z_src_info": "discrete",
@@ -102,7 +102,7 @@ def load_validation_config(halo_profile_model=None):
 
     return {
         "TEST_CASE": testcase,
-        "z_source": testcase["z_source"],
+        "z_src": testcase["z_src"],
         "cosmo": cosmo,
         "cosmo_pars": {k.replace("cosmo_", ""): v for k, v in testcase.items() if "cosmo_" in k},
         "RHO_PARAMS": RHO_PARAMS,
@@ -473,7 +473,7 @@ def helper_physics_functions(func, additional_kwargs={}):
         "mdelta": 1.0e15,
         "cdelta": 4.0,
         "z_cluster": 0.2,
-        "z_source": 0.45,
+        "z_src": 0.45,
         "cosmo": theo.Cosmology(Omega_dm0=0.25, Omega_b0=0.05, H0=70.0),
     }
     kwargs.update(additional_kwargs)
@@ -547,19 +547,19 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         cfg_inf = load_validation_config()
 
         # compute some values
-        cfg_inf["GAMMA_PARAMS"]["z_source"] = 1000.0
+        cfg_inf["GAMMA_PARAMS"]["z_src"] = 1000.0
         beta_s_mean = compute_beta_s_mean(
-            cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_source"], cosmo
+            cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_src"], cosmo
         )
         beta_s_square_mean = compute_beta_s_square_mean(
-            cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_source"], cosmo
+            cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_src"], cosmo
         )
 
         gammat_inf = theo.compute_tangential_shear(cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"])
         kappa_inf = theo.compute_convergence(cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"])
 
         # test z_src = chang2013 distribution
-        cfg_inf["GAMMA_PARAMS"]["z_source"] = chang2013
+        cfg_inf["GAMMA_PARAMS"]["z_src"] = chang2013
         cfg_inf["GAMMA_PARAMS"]["z_src_info"] = "distribution"
 
         # store original values
@@ -762,7 +762,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
 
         # test z_src_info = 'beta'
         beta_s_mean, beta_s_square_mean = 0.9, 0.6
-        cfg_inf["GAMMA_PARAMS"]["z_source"] = (beta_s_mean, beta_s_square_mean)
+        cfg_inf["GAMMA_PARAMS"]["z_src"] = (beta_s_mean, beta_s_square_mean)
         cfg_inf["GAMMA_PARAMS"]["z_src_info"] = "beta"
         # tangential shear
         assert_allclose(
@@ -836,7 +836,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         # First, check for a array of radius and single source z
         radius = np.logspace(-2, 2, 10)
         z_cluster = 0.3
-        z_source = 0.2
+        z_src = 0.2
 
         assert_allclose(
             theo.compute_convergence(
@@ -844,7 +844,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
             np.zeros(len(radius)),
@@ -856,7 +856,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
             np.zeros(len(radius)),
@@ -868,7 +868,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
             np.zeros(len(radius)),
@@ -880,7 +880,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
             np.ones(len(radius)),
@@ -893,7 +893,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
             np.ones(len(radius)),
@@ -902,17 +902,17 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
 
         # Second, check a single radius and array of source z
         radius = 1.0
-        z_source = [0.25, 0.1, 0.14, 0.02]
+        z_src = [0.25, 0.1, 0.14, 0.02]
         assert_allclose(
             theo.compute_convergence(
                 radius,
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
-            np.zeros(len(z_source)),
+            np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
@@ -921,10 +921,10 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
-            np.zeros(len(z_source)),
+            np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
@@ -933,10 +933,10 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
-            np.zeros(len(z_source)),
+            np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
@@ -945,10 +945,10 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
-            np.ones(len(z_source)),
+            np.ones(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
@@ -958,10 +958,10 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
                 mdelta=1.0e15,
                 cdelta=4.0,
                 z_cluster=z_cluster,
-                z_source=z_source,
+                z_src=z_src,
                 cosmo=cosmo,
             ),
-            np.ones(len(z_source)),
+            np.ones(len(z_src)),
             1.0e-10,
         )
 
@@ -996,7 +996,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         profile_pars = [
             cfg["GAMMA_PARAMS"]["r_proj"],
             cfg["GAMMA_PARAMS"]["z_cluster"],
-            cfg["GAMMA_PARAMS"]["z_source"],
+            cfg["GAMMA_PARAMS"]["z_src"],
         ]
         # Validate tangential shear
         gammat = mod.eval_tangential_shear(*profile_pars)
@@ -1113,49 +1113,49 @@ def test_shear_convergence_unittests(modeling_data, profile_init):
         # First, check for a array of radius and single source z
         radius = np.logspace(-2, 2, 10)
         z_cluster = 0.3
-        z_source = 0.2
+        z_src = 0.2
 
         assert_allclose(
-            mod.eval_convergence(radius, z_cluster, z_source), np.zeros(len(radius)), 1.0e-10
+            mod.eval_convergence(radius, z_cluster, z_src), np.zeros(len(radius)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_tangential_shear(radius, z_cluster, z_source), np.zeros(len(radius)), 1.0e-10
+            mod.eval_tangential_shear(radius, z_cluster, z_src), np.zeros(len(radius)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_reduced_tangential_shear(radius, z_cluster, z_source),
+            mod.eval_reduced_tangential_shear(radius, z_cluster, z_src),
             np.zeros(len(radius)),
             1.0e-10,
         )
         assert_allclose(
-            mod.eval_magnification(radius, z_cluster, z_source), np.ones(len(radius)), 1.0e-10
+            mod.eval_magnification(radius, z_cluster, z_src), np.ones(len(radius)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_magnification_bias(radius, z_cluster, z_source, alpha),
+            mod.eval_magnification_bias(radius, z_cluster, z_src, alpha),
             np.ones(len(radius)),
             1.0e-10,
         )
 
         # Second, check a single radius and array of source z
         radius = 1.0
-        z_source = [0.25, 0.1, 0.14, 0.02]
+        z_src = [0.25, 0.1, 0.14, 0.02]
 
         assert_allclose(
-            mod.eval_convergence(radius, z_cluster, z_source), np.zeros(len(z_source)), 1.0e-10
+            mod.eval_convergence(radius, z_cluster, z_src), np.zeros(len(z_src)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_tangential_shear(radius, z_cluster, z_source), np.zeros(len(z_source)), 1.0e-10
+            mod.eval_tangential_shear(radius, z_cluster, z_src), np.zeros(len(z_src)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_reduced_tangential_shear(radius, z_cluster, z_source),
-            np.zeros(len(z_source)),
+            mod.eval_reduced_tangential_shear(radius, z_cluster, z_src),
+            np.zeros(len(z_src)),
             1.0e-10,
         )
         assert_allclose(
-            mod.eval_magnification(radius, z_cluster, z_source), np.ones(len(z_source)), 1.0e-10
+            mod.eval_magnification(radius, z_cluster, z_src), np.ones(len(z_src)), 1.0e-10
         )
         assert_allclose(
-            mod.eval_magnification_bias(radius, z_cluster, z_source, alpha),
-            np.ones(len(z_source)),
+            mod.eval_magnification_bias(radius, z_cluster, z_src, alpha),
+            np.ones(len(z_src)),
             1.0e-10,
         )