From d963d16f97dc01804d9e9313ed330d5aababccd6 Mon Sep 17 00:00:00 2001
From: Delon Shen <delon@stanford.edu>
Date: Fri, 7 Jun 2024 15:48:38 -0700
Subject: [PATCH] Adds prescription of Costanzi et al. 2013 (JCAP12(2013)012)
 for evaluating the halo mass function in nuCDM cosmologies

---
 pyccl/boltzmann.py             | 18 ++++++++++++++++--
 pyccl/halos/halo_model_base.py | 18 +++++++++++++++++-
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/pyccl/boltzmann.py b/pyccl/boltzmann.py
index 65e085675..947d97c5a 100644
--- 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
index 8c9b5103c..96569bec9 100644
--- 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: