implement host halo concentration dependence

pyHalo/Halos/HaloModels/TNFW.py

@@ -114,17 +114,10 @@ def bound_mass(self):
         Computes the mass inside the virial radius (with truncation effects included)
         :return: the mass inside r = c * r_s
         """
-        def _integrand(r, rs, rhos, rt):
-            x = r/rs
-            tau = rt/rs
-            return 4 * tau **2 * np.pi * r ** 2 * rhos / x / (1+x)**2 / (x**2 + tau**2)
-
-        params_physical = self.params_physical
-        rhos = params_physical['rhos']
-        rs = params_physical['rs']
-        rt = params_physical['r_trunc_kpc']
-        r200 = rs * self.c
-            #tau = params_physical['r_trunc_kpc'] / params_physical['rs']
-        #f = tnfw_mass_fraction(tau, self.c) * self._rescale_norm
-        bound_mass = quad(_integrand, 0, r200, args=(rs, rhos, rt))[0]
-        return bound_mass
+        if hasattr(self, '_kwargs_lenstronomy'):
+            tau = self._kwargs_lenstronomy[0]['r_trunc'] / self._kwargs_lenstronomy[0]['Rs']
+        else:
+            params_physical = self.params_physical
+            tau = params_physical['r_trunc_kpc'] / params_physical['rs']
+        f = tnfw_mass_fraction(tau, self.c)
+        return f * self.mass

pyHalo/Halos/galacticus_truncation/interp_mass_loss.py

@@ -80,28 +80,33 @@ class InterpGalacticus(object):
     def __init__(self):
         from pyHalo.Halos.galacticus_truncation.johnsonSUparams import a_fit, \
             b_fit
-        log10c_values = np.linspace(np.log10(2.0), np.log10(384), 25)
-        t_inf_values = np.linspace(0.0, 8.1, 25)
-        a_values = np.array(a_fit).reshape(25, 25)
-        b_values = np.array(b_fit).reshape(25, 25)
-        _points = (t_inf_values, log10c_values)
+        log10c_values = np.linspace(np.log10(2.0), np.log10(384), 10)
+        t_inf_values = np.linspace(0.0, 8.1, 10)
+        chost_values = np.linspace(3.0, 9.0, 10)
+        a_values = np.array(a_fit).reshape(10, 10, 10)
+        b_values = np.array(b_fit).reshape(10, 10, 10)
+        _points = (t_inf_values, log10c_values, chost_values)
         self._a_interp = RegularGridInterpolator(_points, a_values, bounds_error=False,
                                                fill_value=None)
         self._b_interp = RegularGridInterpolator(_points, b_values, bounds_error=False,
                                            fill_value=None)
 
-    def __call__(self, log10_concentration_infall, time_since_infall):
+    def __call__(self, log10_concentration_infall, time_since_infall, chost):
         """
         Evaluates the prediction from galacticus for subhalo bound mass
         :param log10_concentration_infall: log10(c) where c is the halo concentration at infall
         :param time_since_infall: the time ellapsed since infall and the deflector redshift
+        :param chost: host halo concentration at z=0.5
         :return: the log10(bound mass divided by the infall mass), plus scatter
         """
         log10_concentration_infall = max(np.log10(2), log10_concentration_infall)
         log10_concentration_infall = min(np.log10(384), log10_concentration_infall)
         time_since_infall = max(0.0, time_since_infall)
         time_since_infall = min(time_since_infall, 8.1)
-        p = (time_since_infall, log10_concentration_infall)
+        #chost = max(3.0, chost)
+        chost = max(6.0, chost) # we're getting some unphysical behavior at chost = 3.0
+        chost = min(9.0, chost)
+        p = (time_since_infall, log10_concentration_infall, chost)
         a, b = self._a_interp(p), self._b_interp(p)
         output = float(johnsonsu.rvs(a, b))
         return output

pyHalo/Halos/tidal_truncation.py

@@ -321,12 +321,13 @@ def truncation_radius(self, halo_mass, infall_concentration,
 
 class TruncationGalacticus(object):
 
-    def __init__(self, lens_cosmo):
+    def __init__(self, lens_cosmo, c_host):
         """
 
-        :param lens_cosmo:
+        :param lens_cosmo: an instance of LensCosmo
+        :param c_host: host halo concentration at z=0.5
         """
-
+        self._chost = c_host
         self._lens_cosmo = lens_cosmo
         self._mass_loss_interp = InterpGalacticus()
         self._tau_mf_interpolation = tau_mf_interpolation()
@@ -354,9 +355,8 @@ def truncation_radius_halo(self, halo):
         halo_mass = halo.mass
         infall_concentration = halo.c
         time_since_infall = halo.time_since_infall
-
         log10c = np.log10(infall_concentration)
-        log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall)
+        log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall, self._chost)
         m_bound = halo_mass * 10 ** log10mbound_over_minfall
         _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass,
                                                            infall_concentration,
@@ -377,7 +377,7 @@ def truncation_radius(self, halo_mass, infall_concentration,
         :return:
         """
         log10c = np.log10(infall_concentration)
-        log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall)
+        log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall, self._chost)
         m_bound = halo_mass * 10 ** log10mbound_over_minfall
         _, rs, r200 = self._lens_cosmo.NFW_params_physical(halo_mass,
                                                         infall_concentration,

pyHalo/PresetModels/cdm.py

@@ -17,7 +17,7 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig
         shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3,  r_tidal=0.25,
         LOS_normalization=1.0, two_halo_contribution=True, delta_power_law_index=0.0,
         geometry_type='DOUBLE_CONE', kwargs_cosmo=None, host_scaling_factor=0.5,
-        redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, draw_poisson=True):
+        redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, draw_poisson=True, c_host=6.0):
     """
     This class generates realizations of dark matter structure in Cold Dark Matter
 
@@ -57,6 +57,7 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig
     :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos
     :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the
     main deflector presented by Lazar et al. (2021)
+    :param c_host: manually set host halo concentration
     :return: a realization of dark matter halos
     """
 
@@ -90,15 +91,16 @@ def CDM(z_lens, z_source, sigma_sub=0.025, log_mlow=6., log_mhigh=10., log10_sig
     model_fieldhalos, kwargs_mc_field = preset_concentration_models(concentration_model_fieldhalos,
                                                                     kwargs_concentration_model_fieldhalos)
     concentration_model_fieldhalos = model_fieldhalos(**kwargs_mc_field)
-    c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens)
+    if c_host is None:
+        c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens)
 
     # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS
     kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo
     kwargs_truncation_model_fieldhalos['lens_cosmo'] = pyhalo.lens_cosmo
 
     model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos)
     kwargs_trunc_subs.update(kwargs_truncation_model_subhalos)
-    if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24':
+    if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']:
         kwargs_trunc_subs['c_host'] = c_host
     truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs)
 

pyHalo/PresetModels/wdm.py

@@ -18,7 +18,7 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l
         LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None,
         mdef_subhalos='TNFW', mdef_field_halos='TNFW', kwargs_density_profile={},
         host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True,
-        draw_poisson=True):
+        draw_poisson=True, c_host=None):
 
     """
     This class generates realizations of dark matter structure in Warm Dark Matter
@@ -63,6 +63,7 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l
     :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos
     :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the
     main deflector presented by Lazar et al. (2021)
+    :param c_host: manually set host halo concentration
     :return: a realization of dark matter halos
     """
     # FIRST WE CREATE AN INSTANCE OF PYHALO, WHICH SETS THE COSMOLOGY
@@ -104,14 +105,15 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l
     concentration_model_CDM = preset_concentration_models('DIEMERJOYCE19')[0]
     kwargs_mc_field['concentration_cdm_class'] = concentration_model_CDM
     concentration_model_fieldhalos = model_fieldhalos(**kwargs_mc_field)
-    c_host = concentration_model_fieldhalos.nfw_concentration(10 ** log_m_host, z_lens)
+    if c_host is None:
+        c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens)
 
     # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS
     kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo
     kwargs_truncation_model_fieldhalos['lens_cosmo'] = pyhalo.lens_cosmo
     model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos)
     kwargs_trunc_subs.update(kwargs_truncation_model_subhalos)
-    if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24':
+    if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']:
         kwargs_trunc_subs['c_host'] = c_host
     kwargs_trunc_subs['lens_cosmo'] = pyhalo.lens_cosmo
     truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs)
@@ -247,7 +249,7 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow=
         shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3, r_tidal=0.25,
         LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None,
         mdef_subhalos='TNFW', mdef_field_halos='TNFW', kwargs_density_profile={},
-               host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True):
+        host_scaling_factor=0.5, redshift_scaling_factor=0.3, two_halo_Lazar_correction=True, c_host=None):
 
     """
     This preset model implements a generalized treatment of warm dark matter, or any theory that produces a cutoff in
@@ -290,6 +292,7 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow=
     :param redshift_scaling_factor: the scaling with (1+z) of the projected number density of subhalos
     :param two_halo_Lazar_correction: bool; if True, adds the correction to the two-halo contribution from around the
     main deflector presented by Lazar et al. (2021)
+    :param c_host: manually fix the host halo concentration
     :return:
     """
     # FIRST WE CREATE AN INSTANCE OF PYHALO, WHICH SETS THE COSMOLOGY
@@ -320,15 +323,15 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow=
 
     concentration_model_fieldhalos = model_fieldhalos(**kwargs_concentration_model_fieldhalos)
     concentration_model_CDM = preset_concentration_models('DIEMERJOYCE19')[0](pyhalo.astropy_cosmo)
-    c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens)
+    if c_host is None:
+        c_host = concentration_model_CDM.nfw_concentration(10 ** log_m_host, z_lens)
     # SET THE TRUNCATION RADIUS FOR SUBHALOS AND FIELD HALOS
     kwargs_truncation_model_subhalos['lens_cosmo'] = pyhalo.lens_cosmo
     kwargs_truncation_model_fieldhalos['lens_cosmo'] = pyhalo.lens_cosmo
     model_subhalos, kwargs_trunc_subs = truncation_models(truncation_model_subhalos)
     kwargs_trunc_subs.update(kwargs_truncation_model_subhalos)
-    if truncation_model_subhalos == 'TRUNCATION_GALACTICUS_KEELEY24':
+    if truncation_model_subhalos in ['TRUNCATION_GALACTICUS_KEELEY24', 'TRUNCATION_GALACTICUS']:
         kwargs_trunc_subs['c_host'] = c_host
-
     truncation_model_subhalos = model_subhalos(**kwargs_trunc_subs)
     model_fieldhalos, kwargs_trunc_field = truncation_models(truncation_model_fieldhalos)
     kwargs_trunc_field.update(kwargs_truncation_model_fieldhalos)

tests/test_halos/test_galacticus_truncation.py

@@ -12,7 +12,7 @@ def setup_method(self):
 
         self.lens_cosmo_instance = LensCosmo(0.5, 1.5)
         chost = 5.0
-        self.tg = TruncationGalacticus(self.lens_cosmo_instance)
+        self.tg = TruncationGalacticus(self.lens_cosmo_instance, chost)
         self.tgk24 = TruncationGalacticusKeeley24(self.lens_cosmo_instance, chost)
 
     def test_interp_k24(self):
@@ -113,7 +113,7 @@ def rescale_normalization(self, rescale):
 
     def test_tg(self):
         lens_cosmo = LensCosmo(0.5, 2.0,None)
-        truncation_class = TruncationGalacticus(lens_cosmo)
+        truncation_class = TruncationGalacticus(lens_cosmo, c_host=4.0)
         concentration_class = ConcentrationDiemerJoyce(lens_cosmo.cosmo, scatter=False)
         mass = 10 ** 8
         x = 0

tests/test_halos/test_tnfw_halo.py

@@ -92,10 +92,10 @@ def test_density_mass(self):
         mtheory = 4 * np.pi * rho_s * rs ** 3 * (np.log(1 + c) - c/(1+c))
         npt.assert_almost_equal(mtheory, tnfw_fieldhalo.mass)
         rmax = c * rs
+
         m_calculated = tnfw_fieldhalo.mass_3d(rmax)
         npt.assert_almost_equal(mtheory/m_calculated, 1.0)
 
-
 if __name__ == '__main__':
     pytest.main()
 

tests/test_halos/test_truncation.py

@@ -24,7 +24,7 @@ def test_load_models(self):
                            'TRUNCATION_GALACTICUS',
                            'TRUNCATION_GALACTICUS_KEELEY24']
         kwargs_model_list = [{}, {'LOS_truncation_factor': 50.}, {'RocheNorm': 1.0, 'm_power': 1./3, 'RocheNu': 2.0/3.0}, {},
-                             {}, {}, {'rt_arcsec': 1.0}, {}, {'c_host': 5.0}]
+                             {}, {}, {'rt_arcsec': 1.0}, {'c_host': 5.0}, {'c_host': 9.0}]
         for model,kwargs in zip(model_name_list, kwargs_model_list):
             mod, kw = truncation_models(model)
             kwargs.update(kw)