Skip to content

Commit

Permalink
SZ sugar (#784)
Browse files Browse the repository at this point in the history 
* SZ tracer

* SZ profile

* arnaud declare

* divide

* old notation

* misnomer

* final missing factor of h70

* small cleanup

* gnfw

* unit test from class_z

* doc fixes

* benchmarked

* explicit update

* minibug

* doc comoving, better test

* units doc
  • Loading branch information
@damonge
damonge committed Jun 23, 2020
1 parent 3cf3e68 commit 076d828
Show file tree
Hide file tree
Showing 12 changed files with 470 additions and 53 deletions.
62 changes: 62 additions & 0 deletions benchmarks/data/codes/class_sz_P13.ini
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
path_to_class = /home/alonso/Science/LSST/CCL/CCL/benchmarks/data/codes/class_sz
h = 0.67
tau_reio = 0.07
n_s = 0.9645
k_pivot = 0.05
Omega_b = 0.05
Omega_cdm = 0.25
A_s = 2.02e-9
N_ncdm = 0
N_ur = 3.046
deg_ncdm = 3
m_ncdm = 0.0
T_ncdm = 0.71611
output = tSZ_1h
root = ../sz_cl_P13_
write sz results to files = yes
component of tSZ power spectrum = total
B = 1.41
M1SZ = 100000000000.0
M2SZ = 1e+16
z1SZ = 1e-5
z2SZ = 3.
z_max_pk = 3.
x_outSZ = 6.
mass function = M500
pressure profile = Custom. GNFW
integration method (pressure profile) = gsl
HMF_prescription_NCDM = CDM
temperature mass relation = standard
units for tSZ spectrum = dimensionless
multipoles_sz = ell_mock
ell_max_mock = 15000.0
ell_min_mock = 2.0
dlogell = 0.4
mass_range = 0
has_completeness_for_cc = 0
max redshift for cluster counts = 10.
integration method (mass) = patterson
mass_epsrel = 1e-3
mass_epsabs = 1e-30
redshift_epsrel = 1e-3
redshift_epsabs = 1e-30
patterson_show_neval = 0
number of mass bins for cov(Y,N) = 50
experiment = 0
sz_verbose = 1
background_verbose = 1
headers = yes
write primordial = no
write background = no
write parameters = no
create_ref_trispectrum_for_cobaya = NO
f_sky = 0.47
ndimSZ = 100
n_arraySZ = 100
P0GNFW = 6.41
c500 = 1.81
gammaGNFW = 0.31
alphaGNFW = 1.33
betaGNFW = 4.13
pressure_profile_epsabs = 1e-09
pressure_profile_epsrel = 0.01
11 changes: 11 additions & 0 deletions benchmarks/data/codes/sz_cl_benchmark.sh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
#!/bin/bash

# This generates the SZ power spectrum benchmarks
git clone git@github.com:borisbolliet/class_sz.git
cd class_sz
make clean ; make -j8 class
cd ..

# Planck 2013
./class_sz/class class_sz_P13.ini
rm ../sz_cl_P13_pk_cb.dat ../sz_cl_P13_pk.dat ../sz_cl_P13_redshift_dependent_functions.txt ../sz_cl_P13_tk.dat
52 changes: 52 additions & 0 deletions benchmarks/data/sz_cl_P13_szpowerspectrum.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
# Input mass bias b = 2.907801e-01
# sigma8 = 7.719890e-01
# Omega_m = 3.000000e-01
# h = 6.700000e-01
# HMF: Tinker et al 2008 @ M500c
# Pressure Profile: Custom. GNFW
# P0GNFW = 6.410000e+00
# c500 = 1.810000e+00
# gammaGNFW = 3.100000e-01
# alphaGNFW = 1.330000e+00
# betaGNFW = 4.130000e+00
# Dimensions for y power: 'dimensionless'
# Omega_survey = 1.93917e+04 deg2
# f_sky = 0.470000

# Columns:
# 1:multipole
# 2:10^12*ell*(ell+1)/(2*pi)*C_l^tSZ (1-halo term)
# 3:Unbinned Gaussian sampling variance (sigma_g_C_l^2)
# 4:Diagonal elements of non-Gaussian sampling variance (trispectrum, T_ll/Omega_survey)
# 5:Binned Gaussian sampling variance (sigma_g_C_l^2_binned)
# 6:Binned total std dev (Gaussian + non-Gaussian)
# 7:2-halo term 10^12*ell*(ell+1)/(2*pi)*C_l^tSZ (2-halo term)
# 8:SZ temperature, Te [in keV]
# 9:b_kSZ_kSZ_gal_1halo [TBC]
# 10:ell*(ell+1)/(2*pi)*C_l^y-phi (1-halo term) [muK] [TBC]
# 11:ell*(ell+1)/(2*pi)*C_l^isw-phi [muK] [TBC]
# 11:ell*(ell+1)/(2*pi)*C_l^isw-y [muK] [TBC]

2.000000e+00 1.048539e-03 4.822658e-07 0.000000e+00 1.274111e-06 1.077891e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
2.983649e+00 1.492550e-03 1.787004e-07 0.000000e+00 3.164670e-07 1.064176e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
4.451082e+00 2.147762e-03 6.247919e-08 0.000000e+00 7.416868e-08 1.051668e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
6.640234e+00 3.121290e-03 2.092838e-08 0.000000e+00 1.665340e-08 1.041985e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
9.906065e+00 4.574784e-03 6.802614e-09 0.000000e+00 3.628485e-09 1.035743e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.477811e+01 6.750958e-03 2.166072e-09 0.000000e+00 7.744700e-10 1.032753e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
2.204635e+01 1.001408e-02 6.801854e-10 0.000000e+00 1.630201e-10 1.032475e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.288929e+01 1.491038e-02 2.115903e-10 0.000000e+00 3.399318e-11 1.034268e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
4.906506e+01 2.224960e-02 6.534532e-11 0.000000e+00 7.037086e-12 1.037107e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
7.319647e+01 3.321457e-02 2.003658e-11 0.000000e+00 1.446386e-12 1.039527e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.091963e+02 4.947596e-02 6.084241e-12 0.000000e+00 2.944077e-13 1.039129e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.629017e+02 7.326050e-02 1.818965e-12 0.000000e+00 5.899963e-14 1.032177e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
2.430208e+02 1.072496e-01 5.302416e-13 0.000000e+00 1.152872e-14 1.013400e-03 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.625445e+02 1.540801e-01 1.486125e-13 0.000000e+00 2.165932e-15 9.762509e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
5.408528e+02 2.151579e-01 3.930730e-14 0.000000e+00 3.840120e-16 9.140143e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
8.068576e+02 2.887880e-01 9.598174e-15 0.000000e+00 6.285538e-17 8.224756e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.203690e+03 3.678893e-01 2.110177e-15 0.000000e+00 9.263070e-18 7.024050e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.795695e+03 4.392708e-01 4.074343e-16 0.000000e+00 1.198880e-18 5.622310e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
2.678862e+03 4.850089e-01 6.725160e-17 0.000000e+00 1.326487e-19 4.161350e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
3.996392e+03 4.887228e-01 9.244278e-18 0.000000e+00 1.222238e-20 2.810883e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
5.961916e+03 4.438364e-01 1.032035e-18 0.000000e+00 9.146589e-22 1.711174e-04 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
8.894133e+03 3.598244e-01 9.181228e-20 0.000000e+00 5.454411e-23 9.299298e-05 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1.326849e+04 2.582760e-01 6.402343e-21 0.000000e+00 2.549577e-24 4.474349e-05 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
27 changes: 27 additions & 0 deletions benchmarks/test_sz.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
import numpy as np
import pyccl as ccl


def test_szcl():
COSMO = ccl.Cosmology(
Omega_b=0.05,
Omega_c=0.25,
h=0.67,
n_s=0.9645,
A_s=2.02E-9,
Neff=3.046)
bm = np.loadtxt("benchmarks/data/sz_cl_P13_szpowerspectrum.txt",
unpack=True)
l_bm = bm[0]
cl_bm = bm[1]
cl_bm *= (2*np.pi) / (1E12*l_bm*(l_bm+1))
mass_def = ccl.halos.MassDef(500, 'critical')
hmf = ccl.halos.MassFuncTinker08(COSMO, mass_def=mass_def)
hbf = ccl.halos.HaloBiasTinker10(COSMO, mass_def=mass_def)
hmc = ccl.halos.HMCalculator(COSMO, hmf, hbf, mass_def)
prf = ccl.halos.HaloProfilePressureGNFW(mass_bias=1./1.41)
pk = ccl.halos.halomod_Pk2D(COSMO, hmc, prf, get_2h=False)
tr = ccl.tSZTracer(COSMO, z_max=4.)
cl = ccl.angular_cl(COSMO, tr, tr, l_bm, p_of_k_a=pk)

assert np.all(np.fabs(cl/cl_bm-1) < 2E-2)
2 changes: 1 addition & 1 deletion pyccl/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@

# Cl's and tracers
from .tracers import Tracer, NumberCountsTracer, WeakLensingTracer, CMBLensingTracer, \
get_density_kernel, get_kappa_kernel, get_lensing_kernel
tSZTracer, get_density_kernel, get_kappa_kernel, get_lensing_kernel
from .cls import angular_cl

# Useful constants and unit conversions
Expand Down
3 changes: 2 additions & 1 deletion pyccl/halos/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,8 @@
from .profiles import ( # noqa
HaloProfile, HaloProfileGaussian,
HaloProfilePowerLaw, HaloProfileNFW,
HaloProfileEinasto, HaloProfileHernquist)
HaloProfileEinasto, HaloProfileHernquist,
HaloProfilePressureGNFW)

# Halo profile 2-point cumulants
from .profiles_2pt import ( # noqa
Expand Down
213 changes: 213 additions & 0 deletions pyccl/halos/profiles.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
from .. import ccllib as lib
from ..core import check
from ..background import h_over_h0
from ..power import sigmaM
from ..pyutils import resample_array, _fftlog_transform
from .concentration import Concentration
Expand Down Expand Up @@ -843,3 +844,215 @@ def _real(self, cosmo, r, M, a, mass_def):
if np.ndim(M) == 0:
prof = np.squeeze(prof, axis=0)
return prof


class HaloProfilePressureGNFW(HaloProfile):
""" Generalized NFW pressure profile by Arnaud et al.
(2010A&A...517A..92A).
The parametrization is:
.. math::
P_e(r) = C\\times P_0 h_{70}^E (c_{500} x)^{-\\gamma}
[1+(c_{500}x)^\\alpha]^{(\\gamma-\\beta)/\\alpha},
where
.. math::
C = 1.65\\,h_{70}^2\\left(\\frac{H(z)}{H_0}\\right)^{8/3}
\\left[\\frac{h_{70}\\tilde{M}_{500}}
{3\\times10^{14}\\,M_\\odot}\\right]^{2/3+0.12},
:math:`x = r/\\tilde{r}_{500}`, :math:`h_{70}=h/0.7`, and the
exponent :math:`E` is -1 for SZ-based profile normalizations
and -1.5 for X-ray-based normalizations. The biased mass
:math:`\\tilde{M}_{500}` is related to the true overdensity
mass :math:`M_{500}` via the mass bias parameter :math:`(1-b)`
as :math:`\\tilde{M}_{500}=(1-b)M_{500}`. :math:`\\tilde{r}_{500}`
is the overdensity halo radius associated with :math:`\\tilde{M}_{500}`
(note the intentional tilde!), and the profile is defined for
a halo overdensity :math:`\\Delta=500` with respect to the
critical density.
The default arguments (other than `mass_bias`), correspond to the
profile parameters used in the Planck 2013 (V) paper. The profile
is calculated in physical (non-comoving) units of eV/cm^3.
Args:
mass_bias (float): the mass bias parameter :math:`1-b`.
P0 (float): profile normalization.
c500 (float): concentration parameter.
alpha (float): profile shape parameter.
beta (float): profile shape parameter.
gamma (float): profile shape parameter.
alpha_P (float): additional mass dependence exponent
P0_hexp (float): power of `h` with which the normalization
parameter should scale (-1 for SZ-based normalizations,
-3/2 for X-ray-based ones).
qrange (tuple): limits of integration to be used when
precomputing the Fourier-space profile template, as
fractions of the virial radius.
nq (int): number of points over which the
Fourier-space profile template will be sampled.
"""
def __init__(self, mass_bias=0.8, P0=6.41,
c500=1.81, alpha=1.33, alpha_P=0.12,
beta=4.13, gamma=0.31, P0_hexp=-1.,
qrange=(1e-3, 1e3), nq=128):
self.qrange = qrange
self.nq = nq
self.mass_bias = mass_bias
self.P0 = P0
self.c500 = c500
self.alpha = alpha
self.alpha_P = alpha_P
self.beta = beta
self.gamma = gamma
self.P0_hexp = P0_hexp

# Interpolator for dimensionless Fourier-space profile
self._fourier_interp = None
super(HaloProfilePressureGNFW, self).__init__()

def update_parameters(self, mass_bias=None, P0=None,
c500=None, alpha=None, beta=None, gamma=None,
alpha_P=None, P0_hexp=None):
""" Update any of the parameters associated with
this profile. Any parameter set to `None` won't be updated.
.. note:: A change in `alpha`, `beta` or `gamma` will trigger
a recomputation of the Fourier-space template, which can be
slow.
Args:
mass_bias (float): the mass bias parameter :math:`1-b`.
P0 (float): profile normalization.
c500 (float): concentration parameter.
alpha (float): profile shape parameter.
beta (float): profile shape parameters.
gamma (float): profile shape parameters.
alpha_P (float): additional mass dependence exponent.
P0_hexp (float): power of `h` with which the normalization should \
scale (-1 for SZ-based normalizations, -3/2 for \
X-ray-based ones).
"""
if mass_bias is not None:
self.mass_bias = mass_bias
if c500 is not None:
self.c500 = c500
if alpha_P is not None:
self.alpha_P = alpha_P
if P0 is not None:
self.P0 = P0
if P0_hexp is not None:
self.P0_hexp = P0_hexp

# Check if we need to recompute the Fourier profile.
re_fourier = False
if alpha is not None:
if alpha != self.alpha:
re_fourier = True
self.alpha = alpha
if beta is not None:
if beta != self.beta:
re_fourier = True
self.beta = beta
if gamma is not None:
if gamma != self.gamma:
re_fourier = True
self.gamma = gamma

if re_fourier and (self._fourier_interp is not None):
self._fourier_interp = self._integ_interp()

def _form_factor(self, x):
# Scale-dependent factor of the GNFW profile.
f1 = (self.c500*x)**(-self.gamma)
exponent = -(self.beta-self.gamma)/self.alpha
f2 = (1+(self.c500*x)**self.alpha)**exponent
return f1*f2

def _integ_interp(self):
# Precomputes the Fourier transform of the profile in terms
# of the scaled radius x and creates a spline interpolator
# for it.
from scipy.interpolate import interp1d
from scipy.integrate import quad

def integrand(x):
return self._form_factor(x)*x

q_arr = np.geomspace(self.qrange[0], self.qrange[1], self.nq)
# We use the `weight` feature of quad to quickly estimate
# the Fourier transform. We could use the existing FFTLog
# framework, but this is a lot less of a kerfuffle.
f_arr = np.array([quad(integrand,
a=1e-4, b=np.inf, # limits of integration
weight="sin", # fourier sine weight
wvar=q)[0] / q
for q in q_arr])
Fq = interp1d(np.log(q_arr), f_arr,
fill_value="extrapolate",
bounds_error=False)
return Fq

def _norm(self, cosmo, M, a, mb):
# Computes the normalisation factor of the GNFW profile.
# Normalisation factor is given in units of eV/cm^3.
# (Bolliet et al. 2017).
h70 = cosmo["h"]/0.7
C0 = 1.65*h70**2
CM = (h70*M*mb/3E14)**(2/3+self.alpha_P) # M dependence
Cz = h_over_h0(cosmo, a)**(8/3) # z dependence
P0_corr = self.P0 * h70**self.P0_hexp # h-corrected P_0
return P0_corr * C0 * CM * Cz

def _real(self, cosmo, r, M, a, mass_def):
# Real-space profile.
# Output in units of eV/cm^3
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)

# Comoving virial radius
# (1-b)
mb = self.mass_bias
# R_Delta*(1+z)
R = mass_def.get_radius(cosmo, M_use * mb, a) / a

nn = self._norm(cosmo, M_use, a, mb)
prof = self._form_factor(r_use[None, :] / R[:, None])
prof *= nn[:, None]

if np.ndim(r) == 0:
prof = np.squeeze(prof, axis=-1)
if np.ndim(M) == 0:
prof = np.squeeze(prof, axis=0)
return prof

def _fourier(self, cosmo, k, M, a, mass_def):
# Fourier-space profile.
# Output in units of eV * Mpc^3 / cm^3.

# Tabulate if not done yet
if self._fourier_interp is None:
self._fourier_interp = self._integ_interp()

# Input handling
M_use = np.atleast_1d(M)
k_use = np.atleast_1d(k)

# hydrostatic bias
mb = self.mass_bias
# R_Delta*(1+z)
R = mass_def.get_radius(cosmo, M_use * mb, a) / a

ff = self._fourier_interp(np.log(k_use[None, :] * R[:, None]))
nn = self._norm(cosmo, M_use, a, mb)

prof = (4*np.pi*R**3 * nn)[:, None] * ff

if np.ndim(k) == 0:
prof = np.squeeze(prof, axis=-1)
if np.ndim(M) == 0:
prof = np.squeeze(prof, axis=0)
return prof

0 comments on commit 076d828

Please sign in to comment.