<jupyter_start><jupyter_text>Halo model power spectraThis example illustrates how to use CCL to calculate power spectra using the halo model.<jupyter_code>import numpy as np
import pyccl as ccl
import pylab as plt
import matplotlib.cm as cm
from scipy.special import erf
%matplotlib inline<jupyter_output><empty_output><jupyter_text>PreliminariesGenerate a cosmology object and a couple of $k$ (comoving wavenumber) and $a$ (scale factor) arrays.<jupyter_code># Cosmology
cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.83, n_s=0.96)

# Wavenumbers and scale factors
k_arr = np.geomspace(1E-4,1E1,256)
a_arr = np.linspace(0.1,1,32)<jupyter_output><empty_output><jupyter_text>Let us also initialize a series of halo model quantities. We will need:- A halo mass definition.- A concentration-mass relation.- A halo mass function parametrization.- A halo bias parametrization.- Halo profiles corresponding to the quantities we want to compute the power spectrum of.<jupyter_code># We will use a mass definition with Delta = 200 times the matter density
hmd_200m = ccl.halos.MassDef200m()

# The Duffy 2008 concentration-mass relation
cM = ccl.halos.ConcentrationDuffy08(hmd_200m)

# The Tinker 2008 mass function
nM = ccl.halos.MassFuncTinker08(cosmo, mass_def=hmd_200m)

# The Tinker 2010 halo bias
bM = ccl.halos.HaloBiasTinker10(cosmo, mass_def=hmd_200m)

# The NFW profile to characterize the matter density around halos
pM = ccl.halos.profiles.HaloProfileNFW(cM)<jupyter_output><empty_output><jupyter_text>Additionally, in order to showcase the flexibility of the current halo model implementation, let us create our own halo profile characterizing the galaxy overdensity using a Halo Occupation Distribution model. We will use the same parametrization used in Nicola et al. https://arxiv.org/abs/1912.08209.<jupyter_code>class HaloProfileHOD(ccl.halos.HaloProfileNFW):
    def __init__(self, c_M_relation,
                 lMmin=12.02, lMminp=-1.34,
                 lM0=6.6, lM0p=-1.43,
                 lM1=13.27, lM1p=-0.323):
        self.lMmin=lMmin
        self.lMminp=lMminp
        self.lM0=lM0
        self.lM0p=lM0p
        self.lM1=lM1
        self.lM1p=lM1p
        self.a0 = 1./(1+0.65)
        self.sigmaLogM = 0.4
        self.alpha = 1.
        super(HaloProfileHOD, self).__init__(c_M_relation)
        self._fourier = self._fourier_analytic_hod

    def _lMmin(self, a):
        return self.lMmin + self.lMminp * (a - self.a0)

    def _lM0(self, a):
        return self.lM0 + self.lM0p * (a - self.a0)

    def _lM1(self, a):
        return self.lM1 + self.lM1p * (a - self.a0)

    def _Nc(self, M, a):
        # Number of centrals
        Mmin = 10.**self._lMmin(a)
        return 0.5 * (1 + erf(np.log(M / Mmin) / self.sigmaLogM))

    def _Ns(self, M, a):
        # Number of satellites
        M0 = 10.**self._lM0(a)
        M1 = 10.**self._lM1(a)
        return np.heaviside(M-M0,1) * ((M - M0) / M1)**self.alpha

    def _fourier_analytic_hod(self, cosmo, k, M, a, mass_def):
        M_use = np.atleast_1d(M)
        k_use = np.atleast_1d(k)

        Nc = self._Nc(M_use, a)
        Ns = self._Ns(M_use, a)
        # NFW profile
        uk = self._fourier_analytic(cosmo, k_use, M_use, a, mass_def) / M_use[:, None]

        prof = Nc[:, None] * (1 + Ns[:, None] * uk)

        if np.ndim(k) == 0:
            prof = np.squeeze(prof, axis=-1)
        if np.ndim(M) == 0:
            prof = np.squeeze(prof, axis=0)
        return prof

    def _fourier_variance(self, cosmo, k, M, a, mass_def):
        # Fourier-space variance of the HOD profile
        M_use = np.atleast_1d(M)
        k_use = np.atleast_1d(k)

        Nc = self._Nc(M_use, a)
        Ns = self._Ns(M_use, a)
        # NFW profile
        uk = self._fourier_analytic(cosmo, k_use, M_use, a, mass_def) / M_use[:, None]

        prof = Ns[:, None] * uk
        prof = Nc[:, None] * (2 * prof + prof**2)

        if np.ndim(k) == 0:
            prof = np.squeeze(prof, axis=-1)
        if np.ndim(M) == 0:
            prof = np.squeeze(prof, axis=0)
        return prof

pg = HaloProfileHOD(cM)<jupyter_output><empty_output><jupyter_text>Halo model power spectraThe halo-model power spectrum for two quantities defined by their halo profiles $u$ and $v$ is given by\begin{equation}   P_{uv}(k) = \int dM\,n(M)\,\langle u(k|M)\,v(k|M)\rangle + \left[\int dM\,n(M)\,b(M)\,\langle u(k|M)\rangle\right]\,\left[\int dM\,n(M)\,b(M)\,\langle v(k|M)\rangle\right]\,P_{\rm lin}(k).\end{equation}The first part, the so-called "one-halo" term, involves the Fourier-space two-point cumulant of the two profiles ($\langle u(k|M)\,v(k|M)\rangle$). This is often approximated by the product of the profile means ($\langle u(k|M)\rangle\,\langle v(k|M)\rangle$), but this is not correct in general (e.g. for the HOD profile). For this reason, the computation of the halo model power spectrum in CCL requires three main ingredients:- The two profiles to be correlated.- The linear matter power spectrum. Together with the profile means it controls the 2-halo term.- The two-point cumulant of the two profiles, controlling the 1-halo term. By default we will set this to the product of the profile means, but we provide ways to define your own two-points.Let's start by computing the different auto- and cross-correlations between matter and galaxy overdensities at $z=0$.The first step is creating a `HMCalculator` object that essentially takes care of the mass integrals.<jupyter_code>hmc = ccl.halos.HMCalculator(cosmo, nM, bM, hmd_200m)<jupyter_output><empty_output><jupyter_text>Now let's start computing the different power spectra. Note that, because we are interested in overdensities, we normalize by the means of both quantities. This is done by setting the `normprof` parameters of all method calls to `True`. 1. Matter-matter<jupyter_code>pk_MM = ccl.halos.halomod_power_spectrum(cosmo, hmc, k_arr, 1.,
                                         pM, normprof1=True)<jupyter_output><empty_output><jupyter_text>2. Galaxy-matter<jupyter_code>pk_gM = ccl.halos.halomod_power_spectrum(cosmo, hmc, k_arr, 1.,
                                         pg, prof2=pM,
                                         normprof1=True, normprof2=True)<jupyter_output><empty_output><jupyter_text>3. Galaxy-galaxyIn this case we have a non-trivial two-point cumulant (see Eq. 3.9 of https://arxiv.org/abs/1909.09102). To deal with this we first create a `ProfileCovar` object that implements this.<jupyter_code>class Profile2ptHOD(ccl.halos.Profile2pt):
    def fourier_2pt(self, prof, cosmo, k, M, a,
                      prof2=None, mass_def=None):
        return prof._fourier_variance(cosmo, k, M ,a, mass_def)
HOD2pt = Profile2ptHOD()
    
pk_gg = ccl.halos.halomod_power_spectrum(cosmo, hmc, k_arr, 1.,
                                         pg, prof_2pt=HOD2pt,
                                         normprof1=True)<jupyter_output><empty_output><jupyter_text>OK, great, now let's plot the power spectra!<jupyter_code>plt.figure()
plt.plot(k_arr, pk_MM, 'r-', label='$P_{MM}(k)$')
plt.plot(k_arr, pk_gM, 'b-', label='$P_{gM}(k)$')
plt.plot(k_arr, pk_gg, 'y-', label='$P_{gg}(k)$')
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='lower left', fontsize=15)
plt.ylabel(r'$P(k)\,\,[{\rm Mpc}^3]$', fontsize=15)
plt.xlabel(r'$k\,\,[{\rm Mpc}^{-1}]$', fontsize=15);<jupyter_output><empty_output><jupyter_text>2D power spectra and $C_\ell$sThe `HMCalculator` objects can also provide `Pk2D` power spectrum objects that can then be used to perform Limber integrals and compute angular power spectra. Let's do so for the 3D cross correlations we just looked at.First, let's create two CCL tracers, a `NumberCountsTracer` that we will associate with the galaxy overdensity, and a `WeakLensing` tracer that we will associate with the matter overdensity.<jupyter_code>z_arr = np.linspace(0,0.5,256)
nz_arr = np.exp(-((z_arr - 0.25) / 0.05)**2)

# We give it unit bias because the HOD P(k) is already biased
t_g = ccl.NumberCountsTracer(cosmo, False, (z_arr, nz_arr),
                             bias=(z_arr, np.ones_like(z_arr)))
t_M = ccl.WeakLensingTracer(cosmo, (z_arr, nz_arr))<jupyter_output><empty_output><jupyter_text>Now let's create `Pk2D` objects for the three cross correlations<jupyter_code>pk_MMf = ccl.halos.halomod_Pk2D(cosmo, hmc, pM,
                                normprof1=True,
                                lk_arr=np.log(k_arr), a_arr=a_arr)
pk_gMf = ccl.halos.halomod_Pk2D(cosmo, hmc, pg, prof2=pM,
                                normprof1=True, normprof2=True,
                                lk_arr=np.log(k_arr), a_arr=a_arr)
pk_ggf = ccl.halos.halomod_Pk2D(cosmo, hmc, pg, prof_2pt=HOD2pt,
                                normprof1=True,
                                lk_arr=np.log(k_arr), a_arr=a_arr)<jupyter_output><empty_output><jupyter_text>Now let's use these to compute $C_\ell$s<jupyter_code>l_arr = np.unique(np.geomspace(2, 1000, 64).astype(int))

cl_MM = ccl.angular_cl(cosmo, t_M, t_M, l_arr, p_of_k_a=pk_MMf)
cl_gM = ccl.angular_cl(cosmo, t_g, t_M, l_arr, p_of_k_a=pk_gMf)
cl_gg = ccl.angular_cl(cosmo, t_g, t_g, l_arr, p_of_k_a=pk_ggf)<jupyter_output><empty_output><jupyter_text>Alright, plotting time!<jupyter_code>plt.figure()
lfac = l_arr * (l_arr + 1) / (2 * np.pi)
plt.plot(l_arr, lfac * cl_MM, 'r-', label='$\gamma\gamma$')
plt.plot(l_arr, lfac * cl_gM, 'b-', label='$g\gamma$')
plt.plot(l_arr, lfac * cl_gg, 'y-', label='$gg$')
plt.xscale('log')
plt.yscale('log')
plt.legend(loc='upper left', fontsize=15, ncol=2, frameon=False)
plt.ylabel(r'$\ell(\ell+1)C_\ell/(2\pi)$', fontsize=15)
plt.xlabel(r'$\ell$', fontsize=15);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12547724693014137, 'fraction_numerical': 0.025281188731813023, 'mean_word_length': 3.3037300177619895, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '662530', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3806, 'n_tokens_neox': 3538, 'n_words': 892, 'path': 'Halo-model-Pk.ipynb', 'repo_name': 'nikosarcevic/CCLX', 'size': '9691'}