Ctrl+K

Halo Occupation Distribution from extragalactic catalogs

Halo Occupation Distribution from extragalactic catalogs#

Notebook owner: Yao-Yuan Mao @yymao. Last run: Mar 8, 2024 by @patricialarsen

In this notebook we demostrate how to plot the halo occupation distribution of the cosmoDC2/skysim/roman_rubin galaxy catalogs.

Learning objectives#

  • Use GCRCatalogs to access the cosmoDC2, skysim or roman_rubin catalogs.

  • Access cosmology in the extragalactic catalogs.

  • Use CCL to predict Halo Mass Function.

import GCRCatalogs

import pyccl as ccl

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Catalog selection#

Uncomment the line corresponding to the catalog you want to inspect

#gc = GCRCatalogs.load_catalog('cosmoDC2_v1.1.4_small')
gc = GCRCatalogs.load_catalog('skysim5000_v1.1.2_small')
#gc = GCRCatalogs.load_catalog('roman_rubin_2023_v1.1.3_elais')

zmax = 0.25
mass_bins = np.logspace(10, 15, 21)
mass_center = np.sqrt(mass_bins[1:] * mass_bins[:-1])

Get the data#

Retrieve the data given the maximum redshift value determined above

data = gc.get_quantities(['halo_mass', 'Mag_true_r_lsst_z0', 'redshift','halo_id'], filters=['redshift < {}'.format(zmax)])

Cosmology#

Retrieve the cosmological parameters and use them to create a ccl cosmology object

cosmo = ccl.Cosmology(
    Omega_c=gc.cosmology.Om0-gc.cosmology.Ob0, 
    Omega_b=gc.cosmology.Ob0, 
    h=gc.cosmology.h, 
    sigma8=gc.cosmology.sigma8, 
    n_s=gc.cosmology.n_s, 
    transfer_function='bbks',
)

Get expected halo mass function from CCL#

We approximate the hmf by using the mean redshift, create a mass definition and decide which mass function to use

# approximate hmf using mean redshift
mean_scale_factor = 1.0/(1.0+data['redshift'].mean())

hmd_fof = ccl.halos.MassDefFof
mf = ccl.halos.MassFuncSheth99(mass_def=hmd_fof)
hmf_dn_dlogm =  mf(cosmo, mass_center, mean_scale_factor)

Determine expected number of halos#

Use the halo mass function and the volume (given by 4/3 pi * fsky * r^3) to get the total number of halos expected in the volume. We also select on halo ids to get the true number of halos in the volume.

d = gc.cosmology.comoving_distance(zmax).to('Mpc').value
volume = np.deg2rad(np.deg2rad(gc.sky_area)) * d**3 / 3.0

dlogm = np.ediff1d(np.log10(mass_bins))
nhalo_expected = hmf_dn_dlogm * volume * dlogm

hid,idx,count=np.unique(data['halo_id'],return_index=True,return_counts=True)

Plot halo occupation distribution#

for Mr_thres, color in zip((-21.5, -21, -20.5, -20), plt.cm.tab20c.colors):
    plt.loglog(
        mass_center, 
        np.histogram(data['halo_mass'][data['Mag_true_r_lsst_z0'] < Mr_thres], mass_bins)[0] / nhalo_expected,
        label=r'$CCL:\ M_r < {}$'.format(Mr_thres),
        c=color,
    );
    plt.loglog(
        mass_center, 
        np.histogram(data['halo_mass'][data['Mag_true_r_lsst_z0'] < Mr_thres], mass_bins)[0]/np.histogram(data['halo_mass'][idx], mass_bins)[0],
        label=r'$M_r < {}$'.format(Mr_thres),
        c=color,
        linestyle='--'
    );


plt.xlabel(r'${\rm M}_h \,/\, {\rm M}_\odot$');
plt.ylabel(r'$\langle N_{\rm gal} \,|\, {\rm M}_h \rangle$');
plt.title(r'HOD $(z < 0.25)$');
plt.ylim(0.01, None)
plt.axhline(1, lw=0.5, c='k');
plt.legend();
Contents