Fix massdef in miscentering calculations and virial delta_mdef

clmm/theory/ccl.py

@@ -147,6 +147,12 @@ def _get_einasto_alpha(self, z_cl=None):
             a_cl = self.cosmo.get_a_from_z(z_cl)
         return self.hdpm._get_alpha(self.cosmo.be_cosmo, self.__mdelta_cor, a_cl)
 
+    def _get_delta_mdef_virial(self, z_cl):
+        "Gets the virial overdensity"
+        a_cl = self.cosmo.get_a_from_z(z_cl)
+        # This fuction can only return the virial overdensity
+        return self.mdef.get_Delta_vir(self.cosmo.be_cosmo, a_cl)
+
     def _eval_3d_density(self, r3d, z_cl):
         """eval 3d density"""
         return self._call_ccl_profile_lens(self.hdpm.real, r3d, z_cl, ndim=3)

clmm/theory/numcosmo.py

@@ -116,6 +116,12 @@ def _get_einasto_alpha(self, z_cl=None):
         # Note that z_cl is needed for CCL<2.6 only
         return self.hdpm.props.alpha
 
+    def _get_delta_mdef_virial(self, z_cl):
+        "Gets the virial overdensity"
+        # This fuction can return the overdensity for all mass definitions,
+        # but in CLMM, we only use it in the case of massdef='virial'.
+        return self.hdpm.Delta(self.cosmo.be_cosmo, z_cl)
+
     def _eval_reduced_tangential_shear_core(self, r_proj, z_cl, z_src):
         """eval reduced tangential shear considering a single redshift plane
         for background sources"""

clmm/theory/parent_class.py

@@ -186,6 +186,10 @@ def _get_einasto_alpha(self, z_cl=None):
         if defined"""
         raise NotImplementedError
 
+    def _get_delta_mdef_virial(self, z_cl):
+        "Gets the overdensity delta value for massdef='virial' from the backends"
+        raise NotImplementedError
+
     def _set_projected_quad(self, use_projected_quad):
         """Implemented for the CCL backend only"""
         raise NotImplementedError
@@ -269,18 +273,20 @@ def _eval_excess_surface_density_2h(
         )
 
     def _eval_rdelta(self, z_cl):
-        return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, self.delta_mdef)
+        delta_mdef = self._get_delta_mdef(z_cl)
+        return compute_rdelta(self.mdelta, z_cl, self.cosmo, self.massdef, delta_mdef)
 
     def _eval_mass_in_radius(self, r3d, z_cl):
         alpha = self._get_einasto_alpha(z_cl) if self.halo_profile_model == "einasto" else None
+        delta_mdef = self._get_delta_mdef(z_cl)
         return compute_profile_mass_in_radius(
             r3d,
             z_cl,
             self.cosmo,
             self.mdelta,
             self.cdelta,
             self.massdef,
-            self.delta_mdef,
+            delta_mdef,
             self.halo_profile_model,
             alpha,
         )
@@ -289,13 +295,14 @@ def _convert_mass_concentration(
         self, z_cl, massdef=None, delta_mdef=None, halo_profile_model=None, alpha=None
     ):
         alpha1 = self._get_einasto_alpha(z_cl) if self.halo_profile_model == "einasto" else None
+        delta_mdef1 = self._get_delta_mdef(z_cl)
         return convert_profile_mass_concentration(
             self.mdelta,
             self.cdelta,
             z_cl,
             self.cosmo,
             massdef=self.massdef,
-            delta_mdef=self.delta_mdef,
+            delta_mdef=delta_mdef1,
             halo_profile_model=self.halo_profile_model,
             alpha=alpha1,
             massdef2=massdef,
@@ -304,6 +311,11 @@ def _convert_mass_concentration(
             alpha2=alpha,
         )
 
+    def _get_delta_mdef(self, z_cl):
+        if self.massdef == "virial":
+            return int(self._get_delta_mdef_virial(z_cl))
+        return self.delta_mdef
+
     # 3.1. Miscentering functions
 
     def _eval_surface_density_miscentered(self, r_proj, z_cl, r_mis, mis_from_backend):
@@ -363,12 +375,16 @@ def _miscentering_params(self, z_cl, mis_from_backend):
             params = 1, (z_cl,)
 
         else:
-            rho_def = self.cosmo.get_rho_m(z_cl)
+            rho_def = (
+                self.cosmo.get_rho_m(z_cl) if self.massdef == "mean" else self.cosmo.get_rho_c(z_cl)
+            )
             r_s = self.eval_rdelta(z_cl) / self.cdelta
 
+            delta_mdef = self._get_delta_mdef(z_cl)
+
             if self.halo_profile_model == "nfw":
                 rho_s = (
-                    self.delta_mdef
+                    delta_mdef
                     / 3.0
                     * self.cdelta**3.0
                     * rho_def
@@ -379,7 +395,7 @@ def _miscentering_params(self, z_cl, mis_from_backend):
             elif self.halo_profile_model == "einasto":
                 alpha_ein = self._get_einasto_alpha(z_cl)
                 rho_s = (
-                    self.delta_mdef
+                    delta_mdef
                     / 3.0
                     * self.cdelta**3.0
                     * rho_def
@@ -395,7 +411,7 @@ def _miscentering_params(self, z_cl, mis_from_backend):
 
             elif self.halo_profile_model == "hernquist":
                 rho_s = (
-                    self.delta_mdef
+                    delta_mdef
                     / 3.0
                     * self.cdelta**3.0
                     * rho_def
@@ -493,7 +509,7 @@ def set_halo_density_profile(self, halo_profile_model="nfw", massdef="mean", del
             Mass definition, supported options are 'mean', 'critical', 'virial'
             (letter case independent)
         delta_mdef: int
-            Overdensity number
+            Overdensity number. No effect if massdef='virial'.
         """
         # make case independent
         validate_argument(locals(), "massdef", str)
@@ -1608,6 +1624,7 @@ def eval_rdelta(self, z_cl):
         """
         if self.validate_input:
             validate_argument(locals(), "z_cl", float, argmin=0)
+
         return self._eval_rdelta(z_cl)
 
     def eval_mass_in_radius(self, r3d, z_cl, verbose=False):

tests/test_theory.py

@@ -1276,3 +1276,12 @@ def test_mass_conversion(modeling_data, profile_init):
         if halo_profile_model == "einasto":
             profile._get_einasto_alpha = lambda z_cl: None
             assert_raises(ValueError, profile.convert_mass_concentration, z_cl)
+
+
+def test_delta_mdef_virial(modeling_data):
+    if modeling_data["nick"] in ["nc", "ccl"]:
+        cfg = load_validation_config()
+        cosmo = cfg['cosmo']
+        mod = theo.Modeling(massdef="virial")
+        mod.set_cosmo(cosmo)
+        assert_equal(mod._get_delta_mdef(0.1), 111)

tests/test_theory_parent.py

@@ -21,6 +21,7 @@ def test_unimplemented(modeling_data):
     assert_raises(NotImplementedError, mod._update_halo_density_profile)
     assert_raises(NotImplementedError, mod._set_einasto_alpha, 0.5)
     assert_raises(NotImplementedError, mod._get_einasto_alpha)
+    assert_raises(NotImplementedError, mod._get_delta_mdef_virial, 0.1)
     assert_raises(NotImplementedError, mod._set_projected_quad, True)
     assert_raises(NotImplementedError, mod.eval_3d_density, [0.3], 0.3)
     assert_raises(NotImplementedError, mod.eval_surface_density, [0.3], 0.3)