lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
230 stars 285 forks source link

Temperature-Density from CLASS and CLASSY #401

Closed ShazAlvi closed 3 years ago

ShazAlvi commented 3 years ago

I calculated the Temperature-Temperature (TT) and Temperature-Density (TD)power spectra between CLASSY and CLASS. The CLASS param file is as follows,

ln10^{10}A_s = 3.054
l_max_lss = 500
tau_reio = 0.0590
omega_cdm = 0.1201
omega_b = 0.022237
lensing = no
n_s = 0.9651
100*theta_s = 1.04090
z_max_pk = 2.5
non linear = halofit
output = dCl tCl pCl mPk
P_k_max_1/Mpc = 15.705
selection = gaussian
selection_mean = 0.49
selection_width = 0.5

The code that I use to work with CLASSY is the following,

from classy import Class
import numpy as np
import matplotlib.pyplot as plt
MyClass = Class()

params = {'output': 'tCl dCl pCl mPk', 'P_k_max_1/Mpc': 15.705,
            'z_max_pk': 2.5, 'non linear':  'halofit',
            'lensing': 'no', 'l_max_lss': 500, 'omega_b': 0.022237,
            'omega_cdm': 0.1201, 'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
            '100*theta_s': 1.04090, 'tau_reio': 0.0590,
            'selection': 'gaussian', 'selection_mean': '0.49',
            'selection_width': 0.5}

MyClass.set(params)
MyClass.compute()
Cl_DD_TD = MyClass.density_cl()
Cl_TT_TE = MyClass.raw_cl()
CLASS = np.loadtxt('../class_public/output/MP_Params_LCDM_cl.dat')
# plt.loglog(CLASS[0:499, 0], CLASS[0:499, 6])
plt.plot(CLASS[0:499, 0], (CLASS[0:499, 0]*(CLASS[0:499, 0] + 1)*(2*np.pi)**-1*Cl_TT_TE['tt'][2:501])/CLASS[0:499, 1])
plt.xlabel(r"$\ell$, Multipole")
plt.ylabel(r"CLASSY/CLASS")
plt.title("Temperature-Temperature")
plt.figure()
plt.plot(CLASS[0:499, 0], (CLASS[0:499, 0]*(CLASS[0:499, 0] + 1)*(2*np.pi)**-1*Cl_DD_TD['td'][2:501])/CLASS[0:499, 6])
plt.xlabel(r"$\ell$, Multipole")
plt.ylabel(r"CLASSY/CLASS")
plt.title("Temperature-Density")
plt.show()

On plotting the ratio of the CLASSY and CLASS power spectra for TT and TD I get these plots. image image The individual power spectra look like this, image image While the TT spectra look fine to me, the TD spectra differ significantly (for my application). The two files show that there is no mismatch between the parameters.

It will be very helpful if someone can point out what I am missing here.

ShazAlvi commented 3 years ago

I did some exploration into this tension with respect to the precision settings in CLASS (precision.h) and the culprit, primarily, appears to be the parameter transfer_neglect_delta_k_S_t0. This parameter is by default set to 0.15. When I increase this parameter to 30. I get the following plots for the ratio of the power spectra from CLASSY and CLASS, image This result can be further improved by reducing the parameter perturb_sampling_stepsize to 0.01 which gives the following plots (from now on transfer_neglect_delta_k_S_t0 is always set to 30 unless mentioned), image However, upon further reducing this quantity (thus increasing the resolution of the conformal time integral of source functions) the tension between the plots seems to increase as is apparent in the comparison below, image

I checked that the CLASS precision settings (even when directly hardcoded in precision.h file) do not seem to affect the output from CLASSY and furthermore, there does not seem to be a way to pass the precision settings in the dictionary when I use the same keywords as I use in the CLASS ini file. A little clarification on this matter will help my work with CLASSY and perhaps even MontePython.

The comparison between the CLASS output for the two precision settings is shown in the following plot, image The similarity between the last two plots indicates that the default precision setting for CLASSY for the conformal time step and transfer_neglect_delta_k_S_t0 is 0.01 and 30, respectively (?) which is why the tension between the two outputs seems to be minimum for this precision setting in CLASS. If so, I wonder how and where this is parameter is passed to CLASS as I could not find a mention of relevant parameter keywords in MontePython or CLASS files (CLASS files which are relevant to MontePython). Some clarification on this matter will be very useful for my future work (with thanks).

I keep this issue open to wait for one of the developers or people who are more experienced on these tools like @ThomasTram @baudren to comment on these plots and my questions.