CAMB vs CCL angular power spectrum

camb_ccl.py

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

import camb
import pyccl as ccl
from pyccl import ccllib as lib


#redshift of sources

z0=0.311
zs = np.linspace(0.5, 0.7, 50)
W = np.exp(-zs/z0)*(zs/z0)**2/2/z0
bias = 1 + 0.84*zs

lmax = 500


#define cosmology

#ccl cosmo
cosmo = ccl.Cosmology(Omega_c = 0.25, Omega_b = 0.05, h = 0.7, sigma8 = 0.8, n_s = 0.96, transfer_function='boltzmann_camb', baryons_power_spectrum='nobaryons', matter_power_spectrum='linear')

#camb cosmo

# =====  BELOW IS THE SETTING OF CAMB COSMOLOGY TO THE INPUT CCL COSMOLOGY ===== #
#see https://ccl.readthedocs.io/en/latest/_modules/pyccl/boltzmann.html#get_camb_pk_lin

nonlin = False
kmax = 10
extra_camb_params = {}
if np.isfinite(cosmo["A_s"]):
    A_s_fid = cosmo["A_s"]
elif np.isfinite(cosmo["sigma8"]):
    # in this case, CCL will internally normalize for us when we init
    # the linear power spectrum - so we just get close
    #TODO what?
    A_s_fid = 2.43e-9 * (cosmo["sigma8"] / 0.87659)**2
else:
    raise CCLError(
        "Could not normalize the linear power spectrum! "
        "A_s = %f, sigma8 = %f" % (
            cosmo['A_s'], cosmo['sigma8']))

# init camb params
cp = camb.model.CAMBparams()

# turn some stuff off
cp.WantCls = True
cp.DoLensing = False
cp.Want_CMB = False
cp.Want_CMB_lensing = False
cp.Want_cl_2D_array = True
cp.WantTransfer = True

# basic background stuff
h2 = cosmo['h']**2
cp.H0 = cosmo['h'] * 100
cp.ombh2 = cosmo['Omega_b'] * h2
cp.omch2 = cosmo['Omega_c'] * h2
cp.omk = cosmo['Omega_k']

# "constants"
cp.TCMB = cosmo['T_CMB']


# neutrinos
# We maually setup the CAMB neutrinos to match the adjustments CLASS
# makes to their temperatures.
cp.share_delta_neff = False
cp.omnuh2 = cosmo['Omega_nu_mass'] * h2
cp.num_nu_massless = cosmo['N_nu_rel']
cp.num_nu_massive = int(cosmo['N_nu_mass'])
cp.nu_mass_eigenstates = int(cosmo['N_nu_mass'])

delta_neff = cosmo['Neff'] - 3.046  # used for BBN YHe comps

# CAMB defines a neutrino degeneracy factor as T_i = g^(1/4)*T_nu
# where T_nu is the standard neutrino temperature from first order
# computations
# CLASS defines the temperature of each neutrino species to be
# T_i_eff = TNCDM * T_cmb where TNCDM is a fudge factor to get the
# total mass in terms of eV to match second-order computations of the
# relationship between m_nu and Omega_nu.
# We are trying to get both codes to use the same neutrino temperature.
# thus we set T_i_eff = T_i = g^(1/4) * T_nu and solve for the right
# value of g for CAMB. We get g = (TNCDM / (11/4)^(-1/3))^4
g = np.power(
    lib.cvar.constants.TNCDM / np.power(11.0/4.0, -1.0/3.0),
    4.0)

if cosmo['N_nu_mass'] > 0:
    nu_mass_fracs = cosmo['m_nu'][:cosmo['N_nu_mass']]
    nu_mass_fracs = nu_mass_fracs / np.sum(nu_mass_fracs)

    cp.nu_mass_numbers = np.ones(cosmo['N_nu_mass'], dtype=np.int)
    cp.nu_mass_fractions = nu_mass_fracs
    cp.nu_mass_degeneracies = np.ones(int(cosmo['N_nu_mass'])) * g
else:
    cp.nu_mass_numbers = []
    cp.nu_mass_fractions = []
    cp.nu_mass_degeneracies = []

# get YHe from BBN
cp.bbn_predictor = camb.bbn.get_predictor()
cp.YHe = cp.bbn_predictor.Y_He(
    cp.ombh2 * (camb.constants.COBE_CMBTemp / cp.TCMB) ** 3,
    delta_neff)

camb_de_models = ['DarkEnergyPPF', 'ppf', 'DarkEnergyFluid', 'fluid']
camb_de_model = extra_camb_params.get('dark_energy_model', 'fluid')
if camb_de_model not in camb_de_models:
    raise ValueError("The only dark energy models CCL supports with"
                    " camb are fluid and ppf.")
cp.set_classes(
    dark_energy_model=camb_de_model
)

if camb_de_model not in camb_de_models[:2] and cosmo['wa'] and \
        (cosmo['w0'] < -1 - 1e-6 or
            1 + cosmo['w0'] + cosmo['wa'] < - 1e-6):
    raise ValueError("If you want to use w crossing -1,"
                    " then please set the dark_energy_model to ppf.")
cp.DarkEnergy.set_params(
    w=cosmo['w0'],
    wa=cosmo['wa']
)

nonlin = False
cp.set_matter_power(
    redshifts=[_z for _z in zs],
    kmax=extra_camb_params.get("kmax", kmax),
    nonlinear=nonlin)
assert cp.NonLinear == camb.model.NonLinear_none

cp.InitPower.set_params(
        As=A_s_fid,
        ns=cosmo['n_s'])


#calculate and compare cls

#camb
cp.SourceWindows = [camb.sources.SplinedSourceWindow(dlog10Ndm=0, z=zs, W=W, bias_z = bias)]
cp.set_for_lmax(lmax)
cp.SourceTerms.limber_phi_lmin = 500
#cp.SourceTerms.SourceLimberBoost = limber_l
cp.SourceTerms.counts_redshift = False
results = camb.get_results(cp)


ell = np.arange(2,lmax+1)
cls = results.get_source_cls_dict(raw_cl=True)
spectrum = 'W1xW1'
cl_camb = cls[spectrum][2:lmax+1]

#ccl


clu1 = ccl.NumberCountsTracer(cosmo, has_rsd=False, dndz=(zs,W), bias=(zs,bias))
cl_ccl = ccl.angular_cl(cosmo, clu1, clu1, ell = ell)



fig,  [ax,ax_rat] =  plt.subplots(2,  figsize = (8,8), sharex=True)
ax.loglog(ell, cl_camb, 'g--', lw = 3, alpha = 0.7, label = 'CAMB: Full treatment')
ax.loglog(ell, cl_ccl, 'r--', lw = 1, alpha = 0.8,  label = 'CCL: Limber approx')



ax_rat.plot(ell, 100*(cl_camb/cl_ccl-1), 'g--', lw = 3, alpha = 0.7, label = 'CAMB/CCL')
ax_rat.set_ylabel(r'$\Delta C_\ell/C_\ell$')
ax_rat.grid()
ax_rat.set_ylim(-10,10)
ax.legend()

ax_rat.set_xlabel('$\ell$')
ax.set_xlabel('$\ell$')
ax.set_ylabel('Cl')