LSSTDESC / CCL

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

Cross-correlations of two disjoint galaxy survey bins #1157

Closed AndreaRubiola97 closed 4 months ago

AndreaRubiola97 commented 4 months ago

Dear all,

I am implementing a multitracer code with a redshift survey divided in multiple, non-overlapping redshift bins.

As a consistency check, I asked pyccl to calculate the cross-correlation between the two bins, expecting it to be zero. Instead, the angular power spectrum is nonzero, although the product of the window, as expected, is. The issue is appearing also for a generic, overly simplified window, e.g. a constant. I am attaching below a minimal example with a simple window and a linear power spectrum as Pk2d input. Thank you very much, Best Andrea

`cosmo_test = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.7, n_s=0.95, sigma8=0.8, transfer_function='bbks',matter_power_spectrum='linear')

z_test = np.linspace(0,2,100) a_test = 1/(1+z_for_test) a_test = np.flip(a_test)

print(a_test)

chi_test = np.flip(cosmo_test.comoving_radial_distance(a_test))

print(chi_test)

k_test = np.logspace(-3,2,100)

Create window

window_1 = [1 if chis <= chi_test[int( len(chi_test)/2)] else 0. for chis in chi_test ] window_1 = np.array(window_1) window_2 = [1 if chis > chi_test[int( len(chi_test)/2)] else 0. for chis in chi_test ] window_2 = np.array(window_2)

plt.plot(chi_test,window_1) plt.plot(chi_test,window_2) plt.show() plt.plot(chi_test,window_1*window_2) # test window product is identically zero plt.show()

Implement simple linear pk

pk_test = np.array([[cosmo_test.linear_matter_power(k_val,a_val) for k_val in k_test ] for a_val in a_test ])

print(pk_test)

pk_test_val = ccl.pk2d.Pk2D(a_arr=a_test, lk_arr=np.log(k_test), pk_arr = pk_test, is_logp = False, extrap_order_lok=2, extrap_order_hik=2, cosmo=cosmo_test, empty=False)

Implement tracer

window_1_tracer = ccl.Tracer() window_1_tracer.add_tracer(cosmo_test,(chi_test,window_1))

window_2_tracer = ccl.Tracer() window_2_tracer.add_tracer(cosmo_test,(chi_test,window_2))

# l_arr = np.unique((np.geomspace(2,10000, 64)).astype(int)) #definire gli l plt.loglog(l_arr,ccl.angular_cl(cosmo_test,window_1_tracer,window_2_tracer,l_arr,pk_test_val)) plt.show()`

tilmantroester commented 4 months ago

My suspicion is that the internal spline representation of the windows is non-zero between the chi samples where the window is 0 and the one where it is 1. Does this still give a non-zero cross-correlation if the windows are separated by >=2 \Delta\chi?

AndreaRubiola97 commented 4 months ago

Dear Tillman,

thank you very much for the prompt answer. Indeed I forgot to specify exactly what you are saying: in an earlier test, the very first I did, my kernels were not contiguous as in this case, but divided by some decimals in redshift and the correlation was correctly null. The issue rose exactly when the kernels became contiguous, as I need in reality. So a possible solution might be to introduce some "epsilon" additional displacement between the last chi of the first bin and the first chi of the second one?

AndreaRubiola97 commented 4 months ago

Indeed, if I edit the implementation of the second window as follows

[1 if chis > chi_test[int( len(chi_test)/2)+1] else 0. for chis in chi_test ]

I correctly get zero in the power spectrum. Of course it makes sense as long as the chi_test array is dense enough that the distance between two points is very small.

Also notice that if I use the window.get_kernel() method, the product of the two windows is zero also in the original implementations, the problem appears only when the c_l calculator method is called.

damonge commented 4 months ago

probably. Alternatively, sample the kernels more finely to reduce spline ringing around the edges.

Since we seem to know what's causing this, I'll close the issue, but please reopen it if that's not the case.

AndreaRubiola97 commented 4 months ago

Dear David, The issue appears also when my chis are pushed up to 1000 or 10000 points, unfortunately.

tilmantroester commented 4 months ago

Fundamentally, the window is only defined at the sample points. Any finite sampling will lead to some bleed because the window isn't defined between the sample points.