Power spectrum discrepancy between CAMB and CLASS

[Gabriel-Ducrocq]

Hi,

I have noticed some discrepancy between the CAMB computed power spectrum and the Class - with option "format=camb" - computed power spectrum. I can't find where it comes from.
Would you have any idea ?

Here is the code I am using and the graphics showing the difference.

from classy import Class
import matplotlib.pyplot as plt
import numpy as np
import camb
import time


cosmo = Class()

COSMO_PARAMS_NAMES = ["n_s", "omega_b", "omega_cdm", "H0", "ln10^{10}A_s", "tau_reio"]
COSMO_PARAMS_CLASS = np.array([0.9665, 0.02242, 0.11933, 68.2, 3.047, 0.0561])
COSMO_PARAMS_CAMB = np.array([0.9665, 0.02242, 0.11933, 68.2, 3.047e-9, 0.0561])

L_MAX_SCALARS=2500
LENSING = 'yes'
OUTPUT_CLASS = 'tCl pCl lCl'

def generate_cls_class(theta):
    params = {'output': OUTPUT_CLASS,
              'l_max_scalars': L_MAX_SCALARS,
              'lensing': LENSING,
              'format':'camb'}
    d = {name:val for name, val in zip(COSMO_PARAMS_NAMES, theta)}
    params.update(d)
    cosmo.set(params)
    start = time.time()
    cosmo.compute()
    end = time.time()
    print("Computation time CLASS:")
    print(end-start)
    cls = cosmo.lensed_cl(L_MAX_SCALARS)
    #10^12 parce que les cls sont exprimés en kelvin carré, du coup ça donne une stdd en 10^6
    cls["tt"] *= 1e12
    cosmo.struct_cleanup()
    cosmo.empty()
    return cls["tt"]


def generate_cls_camb(theta):
    ns, omega_b, omega_cdm, H0, As, tau_reio = theta
    pars = camb.CAMBparams()
    pars.set_cosmology(H0=H0, ombh2=omega_b, omch2=omega_cdm, tau=tau_reio)
    pars.InitPower.set_params(As=As, ns=ns)
    pars.set_for_lmax(config.L_MAX_SCALARS)

    print("Computation time CAMB:")
    start = time.time()
    results = camb.get_results(pars)
    pow_spec = results.get_cmb_power_spectra(pars, CMB_unit='muK')
    end = time.time()
    print(end-start)
    pow_spec = pow_spec["total"][: L_MAX_SCALARS+1, 0]
    return pow_spec


cls_class = generate_cls_class(COSMO_PARAMS_CLASS)
cls_camb = generate_cls_camb(COSMO_PARAMS_CAMB)
scale = np.array([l*(l+1)/(2*np.pi) for l in range(L_MAX_SCALARS+1)])

plt.plot(cls_class*scale, label="CLASS")
plt.plot(cls_camb, label="CAMB")
plt.legend(loc="upper right")
plt.title("C_ell")
plt.show()

Thank you,
Gabriel.

[Stefan-Heimersheim]

Hey,

thanks for that very detailed issue, that's super nice to reproduce!

At the point where you multiply the C_ell with 1e12 you should multiply with the CMB temperature (2.7K) if you want these units cls["tt"] *= T_cmb**2. In uK this is 2.7255e6**2 and gives you the following result:

Now the CLASS curve is correctly showing D_ell = l*(l+1)/(2pi)*C_ell, compare e.g. to the Planck paper here:
https://arxiv.org/pdf/1807.06209.pdf (Figure 1)

Edit: I don't know how CAMB is supposed to be used to reproduce the Planck plot, there is a factor 1.5 missing ...
And CLASS's format=camb option apparently only affects the output files (when running class from the terminal) and not the python output.

Here's your edited code:

from classy import Class
import matplotlib.pyplot as plt
import numpy as np
import camb
import time


cosmo = Class()

COSMO_PARAMS_NAMES = ["n_s", "omega_b", "omega_cdm", "H0", "ln10^{10}A_s", "tau_reio"]
COSMO_PARAMS_CLASS = np.array([0.9665, 0.02242, 0.11933, 68.2, 3.047, 0.0561])
COSMO_PARAMS_CAMB = np.array([0.9665, 0.02242, 0.11933, 68.2, 3.047e-9, 0.0561])

L_MAX_SCALARS=2500
LENSING = 'yes'
OUTPUT_CLASS = 'tCl pCl lCl'

def generate_cls_class(theta):
    params = {'output': OUTPUT_CLASS,
              'l_max_scalars': L_MAX_SCALARS,
              'lensing': LENSING,
              'format':'camb'}
    d = {name:val for name, val in zip(COSMO_PARAMS_NAMES, theta)}
    params.update(d)
    cosmo.set(params)
    start = time.time()
    cosmo.compute()
    end = time.time()
    print("Computation time CLASS:")
    print(end-start)
    cls = cosmo.lensed_cl(L_MAX_SCALARS)
    #10^12 parce que les cls sont exprimés en kelvin carré, du coup ça donne une stdd en 10^6
    T_cmb = 2.7255e6 #uK
    cls["tt"] *= T_cmb**2
    cosmo.struct_cleanup()
    cosmo.empty()
    return cls["tt"]


def generate_cls_camb(theta):
    ns, omega_b, omega_cdm, H0, As, tau_reio = theta
    pars = camb.CAMBparams()
    pars.set_cosmology(H0=H0, ombh2=omega_b, omch2=omega_cdm, tau=tau_reio)
    pars.InitPower.set_params(As=As, ns=ns)
    pars.set_for_lmax(L_MAX_SCALARS)

    print("Computation time CAMB:")
    start = time.time()
    results = camb.get_results(pars)
    pow_spec = results.get_cmb_power_spectra(pars, CMB_unit='muK')
    end = time.time()
    print(end-start)
    pow_spec = pow_spec["total"][: L_MAX_SCALARS+1, 0]
    return pow_spec


cls_class = generate_cls_class(COSMO_PARAMS_CLASS)
cls_camb = generate_cls_camb(COSMO_PARAMS_CAMB)
scale = np.array([l*(l+1)/(2*np.pi) for l in range(L_MAX_SCALARS+1)])

plt.plot(cls_class*scale, label="CLASS")
plt.plot(cls_camb, label="CAMB")
plt.legend(loc="upper right")
plt.title("C_ell")
plt.show()

PS: Also changed config.L_MAX_SCALARS to L_MAX_SCALARS

[Gabriel-Ducrocq]

Hi,

Thanks for your quick reply !

At the point where you multiply the C_ell with 1e12 you should multiply with the CMB temperature (2.7K) if you want these units cls["tt"] *= T_cmb2. In uK this is 2.7255e62 and gives you the following result:

Should I understand that the units of Cls are usually expressed in 2.7^2 K^2 and not K^2 as I expected ?

And CLASS's format=camb option apparently only affects the output files (when running class from the terminal) and not the python output.

Do you have any idea why ?

I will create an issue on the CAMB repo to understand where the 1.5 factors come from.

Thanks again !

[lesgourg]

Dear Gabriel,

Let me add details to Stefan's answer. The C_l can be defined in two ways:

1. as the power spectrum of the map of CMB relative anisotropies, deltaT/T, in which case the a_lm's and the C_l's are both dimensionless; this convention is more diffuse in the theoretical cosmology community;

2. as the power spectrum of the map of CMB absolute anisotropies deltaT, in which case the a_lm can be in K, or muK, etc., and the C__l in K^2, or muK^2, etc; this convention is more diffuse in the observational cosmology community.

At the beginning, both CLASS and CAMB were giving their output in files, not through python wrappers. They were both writing D_l=l(l+1)C_l/(2pi) in their output files, but CLASS had the convention 1) and CAMB the convention 2) with muK^2. Thus the two spectra were different by a factor T_cmb^2 in [muK^2], i.e., by 2.726e6^2. To allow for a quick comparison, I had coded the option "camb_format", such that the CLASS results were multiplied by 2.726e6^2 before being written in the output file, and thus, the two output files were identical. In addition, when one uses CLASS with output files, the conventions and the units are written in crystal-clear way in the headers at the top fo the file, there are no ambiguities.

Now we all use python wrappers. This reshuffled things a bit. I don't think that we propagated the option "camb_format" in the functions of the wrapper that deals with C_l's output. Thus I don't think that it makes a difference to pass 'format':'camb' in your script. Of course we could implement this option at the level of the wrapper. But it should be anyway straightforward to mutiply the wrapper's output by 2.726e6^2 if you want your output in [muK^2]. Then it should agree with CAMB if the cosmological parameters are the same. I don't know where your additional 1.5 factor comes from, did you check that you passed the correct A_s or log[10^10A_s] and tau_reio or z_reio to CAMB?

[Gabriel-Ducrocq]

Hi Lesgourg !

Okay thank you, that's very clear :)

Yes, I finally found where the factor came from: I made a mistake when passing A_s to CAMB. Instead of passing np.exp(3.047)1e-10 I passed 3.047e-9.
Modifying this, everything works as expected !

Thank you !

Gabriel.