LSSTDESC / CCL

DESC Core Cosmology Library: cosmology routines with validated numerical accuracy
BSD 3-Clause "New" or "Revised" License
145 stars 68 forks source link

angular 2pcf from angular power spectrum #916

Closed SergeiDBykov closed 2 years ago

SergeiDBykov commented 2 years ago

Hello!

I have been playing with the pyccl.correlations and tried to reproduce pyccl.correlations.correlation 2pcf for some input angular spectrum. I found difference from the result of CCL and from my implementation of summation over Legendre polynomials. For my implementation see the code below. It uses pyhtools package for calculating Legendre polynomials.

Then I tried to check the output of from CAMB's cl2corr function and it matched with my implementation, see the plot below. a3e6b4c1-463e-4083-b605-eb5d672afaaf

Why is that? I tried different method argument for CCL function and it has no effect on the 2pcf. The formula I used is eq. (38) from CCL paper.

The code is below: (you need pyhtools package to run my calculation, but comparison with CAMB is sufficient).

import numpy as np
import matplotlib.pyplot as plt
import camb.correlations
import pyccl as ccl
import pyshtools as pysh
%matplotlib inline 
fiducial_params = {'Omega_c': 0.222, 'Omega_b': 0.0449,
                 'h': 0.71, 'sigma8': 0.801, 'n_s': 0.963,
                 'transfer_function': 'boltzmann_camb',
                   'baryons_power_spectrum': 'nobaryons',
                   'matter_power_spectrum': 'linear'}
cosmo = ccl.Cosmology(**fiducial_params)

zarr = np.linspace(0.5, 1.5, 100)
dndz_bin = np.exp(-0.5*(zarr-1)**2/0.1**2)

tracer = ccl.NumberCountsTracer(cosmo, has_rsd = True, dndz = (zarr, dndz_bin), bias = (zarr, 2*np.ones_like(zarr)))
ell = np.arange(1000)
cell = ccl.cls.angular_cl(cosmo, tracer, tracer, ell)

theta_test = np.linspace(1e-2, 1, 50)
ccl_res = ccl.correlation(cosmo, ell, cell,theta = theta_test, type = 'NN', method='Legendre', ) 
# Legendre or fftlog or Bessel

wtheta_camb = camb.correlations.cl2corr(np.repeat(np.atleast_2d(cell*(ell*(ell+1)/(2*np.pi))).T, 4, axis = 1), np.cos(np.deg2rad(theta_test)))[:,0]

Pl_i = [pysh.legendre.legendre(ell[-1], np.cos(np.deg2rad(x)), normalization = 'schmidt', csphase = 1)[:,0] for x in theta_test]
wtheta_i = [np.sum(x*cell*(2*ell+1)/(4*np.pi)) for x in Pl_i]

plt.figure(figsize=(10,10))
plt.plot(theta_test, wtheta_i, '.',label = 'pyhtools')
plt.plot(theta_test, ccl_res, '--',label = 'ccl', lw = 3)
plt.plot(theta_test, wtheta_camb, '--',label = 'camb', lw = 3)
plt.legend()
plt.xlabel(r'$\theta$')
plt.ylabel(r'$w(\theta)$')
SergeiDBykov commented 2 years ago

Can it be due to the CCL interpolation of an input angular power spectrum to large ell values (sect 3.4 in the ccl paper)? I tried ell up to 5000 and all the calculations seem to return the same results within 1%!

damonge commented 2 years ago

Hi @SergeiDBykov . This is likely the case. CCL extrapolates the C_ell to a large ell in order to give a converged correlation function. This is likely not a bad thing to do, since power spectra shouldn't magically go to zero after a specific ell, and you can see that it actually ends up mattering at surprisingly large angles.

Have you tried supplying both CCL and your code with a C_ell that is sampled to a very large ell (e.g. 10,000 or 20,000)?

SergeiDBykov commented 2 years ago

Hi @damonge, thanks! Indeed, for a very large ell values the codes agree very well. I conclude that this is due to the interpolation. May be it is worth noting in the docs that the input Cell are interpolated when you use pyccl.correlations. correlation? In the CCL paper there are just a few sentences about that step, easy to miss. Anyway, thank you very much for your help.