Original examples

The examples below allow to reproduce results from the original Lacasa & Grain article (arXiv:1809.05437).

They use the “turbo mode” Sij function which is the original one but is limited to top-hat redshift kernels.

Users are now encouraged to use the general function which allows arbitrary probe kernels. See the corresponding examples in the other notebooks.

Import modules

import math ; pi=math.pi
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
import time

# Import PySSC module
import PySSC

Compute the Sij matrix

Define boundaries (stakes) of the redshift bins

#zstakes = np.linspace(0.2,1.5,num=14)
zstakes = np.array([0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5]) #Explicitely
zmin = np.min(zstakes) ; zmax = np.max(zstakes)

Compute the matrix

t0 = time.process_time()
Sij = PySSC.turboSij(zstakes=zstakes)  #Uses the default cosmology of the article (arXiv:1809.05437) : Planck 2013 LCDM
t1 = time.process_time()
print(t1-t0)

7.993676833999999

# If you want to change cosmology, specify the parameters with a dictionnary in the format of CLASS :
params = {'omega_b':0.022,'omega_cdm':0.12,'H0':67.,'n_s':0.96,'sigma8':0.81}
Sijbis = PySSC.turboSij(zstakes=zstakes,cosmo_params=params)

Applications

Example 1: Plot the Sij matrix

# Sij can be negative (anti-correlation between bins), and varies by some order of magnitude due to redshift evolution.
# Let's first plot ln|Sij|
fig = plt.figure(figsize=(5.5,5))
P = plt.imshow(np.log(abs(Sij)),interpolation='none',cmap='bwr',extent=[zmin,zmax,zmax,zmin])
plt.xticks([]) ; plt.yticks([])
ax1 = fig.add_axes([0.89, 0.1, 0.035, 0.8])
cbar = plt.colorbar(P,ax1)
cbar.ax.tick_params(labelsize=15)
plt.show()

# Let's now plot the correlation matrix : Sij/sqrt(Sii*Sjj)
#Compute the correlation matrix
nzbins = len(zstakes) - 1
correl = np.zeros((nzbins,nzbins))
for i in range(nzbins):
    for j in range(nzbins):
        correl[i,j] = Sij[i,j] / np.sqrt(Sij[i,i]*Sij[j,j])
#Plot it
fig = plt.figure(figsize=(5.5,5))
P = plt.imshow(correl,interpolation='none',cmap='bwr',vmin=-1,vmax=1,extent=[zmin,zmax,zmax,zmin])
plt.xticks([]) ; plt.yticks([])
ax1 = fig.add_axes([0.89, 0.1, 0.035, 0.8])
cbar = plt.colorbar(P,ax1)
cbar.ax.tick_params(labelsize=15)
plt.show()

Example 2: Find the characteristic multipole ell_SSC depending on redshift

# The multipole where the SSC decreases the S/N by a factor 2 compared to the Gaussian cosmic variance-limited case
# Defined in Eq.46 of the article (arXiv:1809.05437)
#Compute it
fsky    = 1.
Nprobes = 1.
Resp    = 5.
ell_SSC = np.zeros(nzbins)
for i in range(nzbins):
    ell_SSC[i] = np.sqrt(2./(Nprobes*Resp**2*fsky*Sij[i,i]))

# Plot it as a function of redshift
zcenter = (zstakes[1:]+zstakes[:-1])/2.
plt.plot(zcenter,ell_SSC,lw=3)
plt.xlabel('z',fontsize=20) ; plt.ylabel(r'$\ell_\mathrm{SSC}$',fontsize=20)
plt.show()

so SSC starts to dominate for multipoles above a few hundred, i.e. sub-degree scale in real space

Example 3: Find the critical density of clusters where SSC/sample variance surpasses shot-noise/Poisson

# If Ncl is the cluster count (per steradian), then the shot-noise/Poisson variance is Cov_shot(Ncl)=Ncl/4pi
# and the SSC/sample variance is Cov_SSC = b^2 Ncl^2 * Sij
# the latter dominates shot-noise when Ncl >= Ncrit = 1./(4pi*b^2*Sij)

#Compute Ncrit
bcl = 5.                                  # Average cluster bias
Ncrit = np.zeros(nzbins)                  # Number density per sterad
for i in range(nzbins):
    Ncrit[i] = 1./(4.*pi*bcl**2*Sij[i,i])
Ncrit_deg2 = Ncrit * (pi/180.)**2         # Number density per deg^2

# Plot it as a function of redshift
plt.plot(zcenter,Ncrit_deg2,lw=3)
plt.xlabel('z',fontsize=20) ; plt.ylabel(r'$N_\mathrm{crit}$ [deg$^{-2}$]',fontsize=20)
#plt.savefig("Ncrit-halos-vs-z.png",bbox_inches='tight')
plt.show()

so SSC dominates when we detect more than a few clusters per square degree