implement Lazar et al. (2021) correction to correlated structure cont…

pyHalo/Rendering/two_halo.py

@@ -2,6 +2,7 @@
 from copy import deepcopy
 from pyHalo.Rendering.line_of_sight import LineOfSightNoSheet
 from scipy.integrate import quad
+from pyHalo.concentration_models import ConcentrationDiemerJoyce
 
 
 class TwoHaloContribution(LineOfSightNoSheet):
@@ -34,25 +35,26 @@ def __init__(self, mass_function_model, kwargs_mass_function, spatial_distributi
         else:
             raise Exception('must specify the host halo mass through keyword argument host_m200 or log_m_host (base 10)'
                             'when adding the two-halo term!')
+        concentration_model = ConcentrationDiemerJoyce(lens_cosmo.cosmo.astropy, scatter=False)
+        c_host = concentration_model.nfw_concentration(host_m200, lens_cosmo.z_lens)
+        _, _, r200_host_kpc = lens_cosmo.NFW_params_physical(host_m200, c_host, lens_cosmo.z_lens)
+        r200_host_mpc = r200_host_kpc * 1e-3
         z_eval = lens_cosmo.z_lens
         idx = np.argmin(abs(np.array(lens_plane_redshifts) - z_eval))
         delta_z = delta_z_list[idx]
-        boost = two_halo_enhancement_factor(z_eval, delta_z, lens_cosmo, host_m200)
-        relative_enhancement = boost - 1.
+        boost = two_halo_enhancement_factor(z_eval, delta_z, lens_cosmo, host_m200, r200_host_mpc)
         kwargs_mass_function_scaled = deepcopy(kwargs_mass_function)
-        kwargs_mass_function_scaled['LOS_normalization'] *= relative_enhancement
+        kwargs_mass_function_scaled['LOS_normalization'] *= boost
         self._delta_z = delta_z
         super(TwoHaloContribution, self).__init__(mass_function_model, kwargs_mass_function_scaled, spatial_distribution_model,
                  geometry, lens_cosmo, lens_plane_redshifts, delta_z_list)
 
     def render(self):
-
         """
-        Generates halo masses and positions for objects along the line of sight
-        (except for halos from the two-halo contribution)
-        :return: mass (in Msun), x (arcsec), y (arcsec), r3d (kpc), redshift
+        Generates halo masses and positions for correlated structure near the main deflector.
+        These objects are placed at the lens redshift
+        :return: mass (in Msun), x (arcsec), y (arcsec), r3d (kpc), redshift, subhalo_flag (bool)
         """
-
         mfunc_model = self._get_mass_function_model(self._lens_cosmo.z_lens, self._delta_z)
         masses = mfunc_model.draw()
         nhalos = len(masses)
@@ -62,19 +64,43 @@ def render(self):
         r3d = np.array([None] * nhalos)
         return masses, x, y, r3d, redshifts, subhalo_flag
 
-def two_halo_enhancement_factor(z_lens, z_step, lens_cosmo, overdensity_m200):
+def modificationLazar2021(r, r200, b_e=0.1, s_e=4.0):
+    """
+    Modification to the two-halo contribution around the main deflector using Equation 3 from Lazar et al. (2021)
+    https://ui.adsabs.harvard.edu/abs/2021MNRAS.502.6064L/abstract
+    :param r: distance in Mpc
+    :param r200: host halo virial radius in Mpc
+    :param b_e: calibrated parameter
+    :param s_e: calibrated parameter
+    :return: the correction factor for the background density that gives a better match to simulations
+    """
+    return b_e * (r / r200 / 5) ** -s_e
+
+def _boost_integrand(r, m200_host, z, r200_host, lens_cosmo):
+    """
+    The local enhancement of the background density of the Universe from the two-halo term and a correction factor
+    see Lazar et al. 2021
+    :param r: distance from the host in Mpc
+    :param m200_host: host halo mass
+    :param z: host halo redshift
+    :param r200_host: host halo virial radius
+    :param lens_cosmo: an instance of LensCosmo
+    :return: the local enhancement to the background density
     """
+    return lens_cosmo.twohaloterm(r, m200_host, z) + modificationLazar2021(r, r200_host)
 
-    :param z_lens:
-    :param lens_plane_redshifts:
-    :param delta_z_list:
-    :param lens_cosmo:
-    :param overdensity_m200:
-    :return:
+def two_halo_enhancement_factor(z_lens, z_step, lens_cosmo, overdensity_m200, r200_host):
+    """
+    Calculates the enhancement of the background density due to correlated structure around the main deflector. Includes
+    a contribution from the two-halo term and a correction factor calibrated by Lazar et al. (2021)
+    :param z_lens: deflector redshift
+    :param z_step: redshift spacing
+    :param lens_cosmo: instance of LensCosmo
+    :param overdensity_m200: host halo mass
+    :return: the enhancement of the background density per unit length
     """
     rmax = lens_cosmo.cosmo.D_C_transverse(z_lens + z_step) - lens_cosmo.cosmo.D_C_transverse(z_lens)
     rmin = min(rmax, 0.5)
-    mean_boost = 2 * quad(lens_cosmo.twohaloterm, rmin, rmax, args=(overdensity_m200, z_lens))[0] / (rmax - rmin)
-    two_halo_boost = 1 + mean_boost
+    args = (overdensity_m200, z_lens, r200_host, lens_cosmo)
+    two_halo_boost = 2 * quad(_boost_integrand, rmin, rmax, args=args)[0] / (rmax - rmin)
     return two_halo_boost
-

tests/test_rendering/test_2halo.py

@@ -35,14 +35,9 @@ def test_boost(self):
 
         z_step = 0.02
         mhost = 10**13
-        boost = two_halo_enhancement_factor(self.zlens, z_step, self.lens_cosmo, mhost)
-        npt.assert_array_less(1.1, boost)
-
-        z_step = 0.2
-        mhost = 10 ** 13
-        boost = two_halo_enhancement_factor(self.zlens, z_step, self.lens_cosmo, mhost)
-        # as the redshift increment increases, the relative contribution from the two halo term diminishes
-        npt.assert_array_less(boost, 1.1)
+        r200_host = 0.4
+        boost = two_halo_enhancement_factor(self.zlens, z_step, self.lens_cosmo, mhost, r200_host)
+        npt.assert_equal(boost>0, True)
 
     def test_render(self):