Adds prescription of Costanzi et al. 2013 (JCAP12(2013)012) for evalu…

…ating the halo mass function in nuCDM cosmologies
LSSTDESC · Jun 7, 2024 · d963d16 · d963d16
1 parent d762d99
commit d963d16
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 3 deletions.
diff --git a/pyccl/boltzmann.py b/pyccl/boltzmann.py
@@ -35,6 +35,15 @@ def get_camb_pk_lin(cosmo, *, nonlin=False):
     except (KeyError, TypeError):
         pass
 
+    # Get extra parameters to evaluate
+    # the halo mass function in cosmology with massive neutrinos
+    # using prescription by Costanzi et al. 2013(JCAP12(2013)012)
+    use_costanzi13 = False
+    try:
+        use_costanzi13 = cosmo["extra_parameters"]["use_costanzi13"]
+    except (KeyError, TypeError):
+        pass
+
     # z sampling from CCL parameters
     na = lib.get_pk_spline_na(cosmo.cosmo)
     status = 0
@@ -151,8 +160,13 @@ def get_camb_pk_lin(cosmo, *, nonlin=False):
     camb_res = camb.get_transfer_functions(cp)
 
     def construct_Pk2D(camb_res, nonlin=False):
-        k, z, pk = camb_res.get_linear_matter_power_spectrum(
-            hubble_units=True, nonlinear=nonlin)
+        if not use_costanzi13:
+            k, z, pk = camb_res.get_linear_matter_power_spectrum(
+                hubble_units=True, nonlinear=nonlin)
+        else:
+            k, z, pk = camb_res.get_linear_matter_power_spectrum(
+                var1='delta_nonu', var2='delta_nonu',
+                hubble_units=True, nonlinear=nonlin)
 
         # convert to non-h inverse units
         k *= cosmo['h']

diff --git a/pyccl/halos/halo_model_base.py b/pyccl/halos/halo_model_base.py
@@ -193,7 +193,23 @@ def __call__(self, cosmo, M, a):
         logM, sigM, dlns_dlogM = self._get_logM_sigM(
             cosmo, M_use, a, return_dlns=True)
 
-        rho = (const.RHO_CRITICAL * cosmo['Omega_m'] * cosmo['h']**2)
+        # prescription by Costanzi et al. 2013(JCAP12(2013)012)
+        # to evaluate the halo mass function in nuCDM cosmology
+
+        use_costanzi13 = False
+
+        try:
+            use_costanzi13 = cosmo["extra_parameters"]["use_costanzi13"]
+        except (KeyError, TypeError):
+            pass
+
+        if not use_costanzi13:
+            rho = (const.RHO_CRITICAL * cosmo['Omega_m'] * cosmo['h']**2)
+        else:
+            rho = (const.RHO_CRITICAL
+                   * (cosmo['Omega_c'] + cosmo['Omega_b'])
+                   * cosmo['h']**2)
+
         f = self._get_fsigma(cosmo, sigM, a, 2.302585092994046 * logM)
         mf = f * rho * dlns_dlogM / M_use
         if np.ndim(M) == 0: