Add SZ likelihood to SOLikeT

setup.cfg

@@ -27,7 +27,7 @@ install_requires =
     syslibrary
 
 [options.package_data]
-soliket = *.yaml,*.bibtex,clusters/data/*,clusters/data/selFn_equD56/*,lensing/data/*.txt,lensing/data/*.fits,mflike/*.yaml,tests/*.yaml,data/xcorr_simulated/*.txt,data/CosmoPower/CP_paper/CMB/*.pkl
+soliket = *.yaml,*.bibtex,clusters/data/*,clusters/data/selFn_equD56/*,lensing/data/*.txt,lensing/data/*.fits,mflike/*.yaml,tests/*.yaml,data/xcorr_simulated/*.txt,data/CosmoPower/CP_paper/CMB/*.pkl,szlike/*.yaml
 testpaths = "soliket"
 text_file_format = rst
 

soliket/__init__.py

@@ -5,6 +5,7 @@
 from .mflike import MFLike  # noqa: F401
 from .mflike import TheoryForge_MFLike
 from .xcorr import XcorrLikelihood  # noqa: F401
+from .szlike import KSZLikelihood, TSZLikelihood
 from .foreground import Foreground
 from .bandpass import BandPass
 from .cosmopower import CosmoPower, CosmoPowerDerived

soliket/constants.py

@@ -7,10 +7,13 @@
 k_Boltzmann = constants.k
 elementary_charge = constants.e
 electron_mass_kg = constants.m_e
+proton_mass_kg = constants.m_p
 T_CMB = 2.72548 #in K
 MSUN_CGS = 1.98840987e+33
 G_CGS = constants.G * 1e3
 MPC2CM = constants.parsec * 1e8
+ST_CGS = constants.value('Thomson cross section') * 1.e4
+hydrogen_fraction = 0.76
 
 full_sky_area_sqdeg = 41252.9612
 SO_sky_area_sqdeg = 987.5

soliket/szlike/KSZLikelihood.yaml

@@ -0,0 +1,56 @@
+sz_data_file: gnfw_test_projections.txt
+beam_file: beam_f150_daynight.txt 
+transform_type: Hankel #FFT or Hankel
+#optional beam response file
+beam_response: False #False, or give file
+cov_ksz_file: cov_diskring_ksz_varweight_bootstrap.txt
+#optional 2-halo file
+twohalo_term: twohalo_cmass_average.txt #False, or give file. Need to flesh out this option more, but the idea is that if False it will use two_halo.py to calculate a 2h term
+
+redshift: 0.55 
+mass_halo_mean_Msol: 1.e13
+frequency_GHz: 150.
+#input_model: "gnfw"
+
+params:
+  Omega_m:
+    value: 0.25
+  Omega_b:
+    value: 0.044
+  hh:
+    value: 0.7
+  Omega_L:
+    value: 0.75
+  rhoc_0:
+    value: 2.77525e2
+  C_OVER_HUBBLE:
+    value: 2997.9
+  XH:
+    value: 0.76 #hydrogen fraction
+  v_rms:
+    value: 1.06e-3 #v_rms/c
+
+  gnfw_rho0:
+    prior:
+      min: 1.0
+      max: 6.0
+    ref: 3.1864410705359125
+    latex: \log_{10}\rho_0
+  gnfw_bt_ksz:
+    prior:
+      min: 1.0
+      max: 10.0
+    ref: 3.45493977635
+    latex: \beta_k
+  gnfw_A2h_ksz:
+    prior:
+      min: 0.01
+      max: 5.0
+    ref: 1.0
+    latex: A_{k2h}
+
+
+#prior: #al<bt for ksz
+#   ksz_al_less_than_bt: 'lambda gnfw_al_ksz,gnfw_bt_ksz:
+#                        stats.uniform.logpdf(gnfw_bt_ksz-gnfw_al_ksz,0,10.)'
+

soliket/szlike/OBBLikelihood.yaml

@@ -0,0 +1,64 @@
+sz_data_file: gnfw_test_projections.txt
+beam_file: beam_f150_daynight.txt 
+#optional beam response file
+beam_response: act_planck_s08_s18_cmb_f150_daynight_response_tsz.txt #False, or give file
+cov_ksz_file: cov_diskring_ksz_varweight_bootstrap.txt
+cov_tsz_file: cov_diskring_tsz_varweight_bootstrap.txt
+
+redshift: 0.55 
+mass_halo_mean_Msol: 10**13.
+frequency_GHz: 150.
+#input_model: "obb" #"gnfw" or "obb"
+
+params:
+  Omega_m:
+    value: 0.25
+  Omega_b:
+    value: 0.044
+  hh:
+    value: 0.7
+  Omega_L:
+    value: 0.75
+  rhoc_0:
+    value: 2.77525e2
+  C_OVER_HUBBLE:
+    value: 2997.9
+  XH:
+    value: 0.76 #hydrogen fraction
+  v_rms:
+    value: 1.06e-3 #v_rms/c
+
+  obb_Gamma: #how do we do gauss prior with min and max? it's truncated. maybe proposal for width?
+    prior:
+      min: 1.
+      max: 1.667 #5/3
+    ref:
+      dist: norm
+      loc: 1.2
+      scale: 0.2
+    latex: \Gamma
+  obb_alpha_Nth:
+    prior:
+      min: 0.01
+      max: 0.8
+    ref: 0.13
+    latex: \alpha_{Nth}
+  obb_logE:
+    prior:
+      min: -4.8
+      max: -4.0
+    latex: \log_{10}\epsilon
+  obb_Ak2h:
+    prior:
+      min: 0.01
+      max: 5.0
+    ref: 1.0
+    latex: A_{k2h}
+  obb_At2h:
+    prior:
+      min: 0.01
+      max: 5.0
+    ref: 1.0
+    latex: A_{t2h}
+
+

soliket/szlike/TSZLikelihood.yaml

@@ -0,0 +1,52 @@
+sz_data_file: gnfw_test_projections.txt
+beam_file: beam_f150_daynight.txt
+transform_type: Hankel #FFT or Hankel
+#optional beam response file
+beam_response: act_planck_s08_s18_cmb_f150_daynight_response_tsz.txt #False, or give file
+cov_tsz_file: cov_diskring_tsz_varweight_bootstrap.txt
+#optional 2-halo file
+twohalo_term: twohalo_cmass_average.txt #False, or give file. Need to flesh out this option more, but the idea is that if False it will use two_halo.py to calculate a 2h term
+
+redshift: 0.55 
+mass_halo_mean_Msol: 1.e13
+frequency_GHz: 150.
+#input_model: "gnfw" #"gnfw" or "obb"
+
+params:
+  Omega_m:
+    value: 0.25
+  Omega_b:
+    value: 0.044
+  hh:
+    value: 0.7
+  Omega_L:
+    value: 0.75
+  rhoc_0:
+    value: 2.77525e2
+  C_OVER_HUBBLE:
+    value: 2997.9
+  XH:
+    value: 0.76 #hydrogen fraction
+  v_rms:
+    value: 1.06e-3 #v_rms/c
+
+  gnfw_P0:
+    prior:
+      min: 0.1
+      max: 30.0
+    ref: 10
+    latex: P_0
+  gnfw_bt_tsz:
+    prior:
+      min: 1.0
+      max: 10.0
+    ref: 4.8
+    latex: \beta_t
+  gnfw_A2h_tsz:
+    prior:
+      min: 0.01
+      max: 5.0
+    ref: 1.0
+    latex: A_{t2h}
+
+

soliket/szlike/__init__.py

@@ -0,0 +1 @@
+from .szlike import KSZLikelihood, TSZLikelihood, OBBLikelihood

soliket/szlike/act_planck_s08_s18_cmb_f150_daynight_response_tsz.txt

@@ -0,0 +1,50 @@
+  1.0000000e+02  -2.7857403e+00   1.5997740e-04  -2.7856139e+00   4.8228018e-03  -2.7856313e+00   4.4830121e-03
+  1.1117845e+02  -2.7856534e+00   1.8469731e-04  -2.7854840e+00   5.5797143e-03  -2.7855074e+00   5.1870487e-03
+  1.2360647e+02  -2.7855449e+00   4.9142881e-04  -2.7842978e+00   1.5100645e-02  -2.7844523e+00   1.4155778e-02
+  1.3742375e+02  -2.7854108e+00   7.0369581e-04  -2.7828468e+00   2.1571018e-02  -2.7831586e+00   2.0256515e-02
+  1.5278559e+02  -2.7852452e+00   8.7837626e-04  -2.7812459e+00   2.6825564e-02  -2.7817249e+00   2.5229086e-02
+  1.6986465e+02  -2.7850344e+00   1.5038144e-03  -2.7732693e+00   4.4995450e-02  -2.7745783e+00   4.2622745e-02
+  1.8885287e+02  -2.7847723e+00   2.0925479e-03  -2.7628954e+00   5.9495480e-02  -2.7651742e+00   5.6765063e-02
+  2.0996369e+02  -2.7844454e+00   2.7415170e-03  -2.7488881e+00   7.2538906e-02  -2.7523072e+00   6.9826619e-02
+  2.3343437e+02  -2.7840349e+00   3.5484271e-03  -2.7298937e+00   8.3643641e-02  -2.7344838e+00   8.1512705e-02
+  2.5952870e+02  -2.7835323e+00   4.2377871e-03  -2.7145305e+00   8.8726171e-02  -2.7195742e+00   8.7379491e-02
+  2.8853998e+02  -2.7829162e+00   4.8612234e-03  -2.7019044e+00   9.0776623e-02  -2.7069713e+00   9.0141534e-02
+  3.2079427e+02  -2.7821455e+00   5.6812608e-03  -2.6884775e+00   9.1122398e-02  -2.6932317e+00   9.1207736e-02
+  3.5665408e+02  -2.7811880e+00   6.6176387e-03  -2.6764074e+00   8.9794280e-02  -2.6807110e+00   9.0419077e-02
+  3.9652247e+02  -2.7799884e+00   7.7989074e-03  -2.6651680e+00   8.7124256e-02  -2.6688968e+00   8.8126028e-02
+  4.4084752e+02  -2.7784696e+00   9.3872995e-03  -2.6544746e+00   8.3207176e-02  -2.6575164e+00   8.4420009e-02
+  4.9012742e+02  -2.7765397e+00   1.1432576e-02  -2.6452076e+00   7.8581491e-02  -2.6476749e+00   7.9883650e-02
+  5.4491605e+02  -2.7740662e+00   1.4060951e-02  -2.6375292e+00   7.3724573e-02  -2.6394098e+00   7.4960280e-02
+  6.0582920e+02  -2.7708470e+00   1.7493990e-02  -2.6309561e+00   6.8717002e-02  -2.6323211e+00   6.9804580e-02
+  6.7355149e+02  -2.7665546e+00   2.2070163e-02  -2.6257433e+00   6.4048074e-02  -2.6267047e+00   6.4961521e-02
+  7.4884408e+02  -2.7609028e+00   2.7463291e-02  -2.6214865e+00   5.9645292e-02  -2.6221452e+00   6.0387637e-02
+  8.3255321e+02  -2.7530970e+00   3.4386623e-02  -2.6174964e+00   5.5036680e-02  -2.6179442e+00   5.5634007e-02
+  9.2561972e+02  -2.7425625e+00   4.2091245e-02  -2.6142223e+00   5.0761329e-02  -2.6145549e+00   5.1270044e-02
+  1.0290896e+03  -2.7283183e+00   5.0157340e-02  -2.6111689e+00   4.6397002e-02  -2.6113931e+00   4.6813562e-02
+  1.1441258e+03  -2.7098290e+00   5.7049597e-02  -2.6082694e+00   4.1921900e-02  -2.6083909e+00   4.2233796e-02
+  1.2720213e+03  -2.6868474e+00   6.0953810e-02  -2.6054717e+00   3.7395426e-02  -2.6055174e+00   3.7615805e-02
+  1.4142136e+03  -2.6612665e+00   5.9414688e-02  -2.6028982e+00   3.3161759e-02  -2.6028859e+00   3.3294560e-02
+  1.5723007e+03  -2.6363129e+00   5.1621845e-02  -2.6006178e+00   2.9602854e-02  -2.6005606e+00   2.9644302e-02
+  1.7480594e+03  -2.6160545e+00   3.9732505e-02  -2.5986774e+00   2.7085557e-02  -2.5985817e+00   2.7047832e-02
+  1.9434653e+03  -2.6033194e+00   2.9486757e-02  -2.5970886e+00   2.5711662e-02  -2.5969452e+00   2.5617107e-02
+  2.1607146e+03  -2.5972109e+00   2.5379678e-02  -2.5957285e+00   2.5185161e-02  -2.5955600e+00   2.5068263e-02
+  2.4022489e+03  -2.5946020e+00   2.5081178e-02  -2.5944084e+00   2.5050874e-02  -2.5942559e+00   2.4938609e-02
+  2.6707830e+03  -2.5930062e+00   2.5208764e-02  -2.5929729e+00   2.5002656e-02  -2.5928193e+00   2.4892133e-02
+  2.9693350e+03  -2.5913793e+00   2.4983309e-02  -2.5913205e+00   2.4881939e-02  -2.5911674e+00   2.4778917e-02
+  3.3012605e+03  -2.5894704e+00   2.4660563e-02  -2.5893646e+00   2.4712072e-02  -2.5892230e+00   2.4616674e-02
+  3.6702901e+03  -2.5871192e+00   2.4389182e-02  -2.5870145e+00   2.4548655e-02  -2.5868742e+00   2.4446433e-02
+  4.0805715e+03  -2.5842129e+00   2.4137968e-02  -2.5841673e+00   2.4351556e-02  -2.5840452e+00   2.4266834e-02
+  4.5367160e+03  -2.5808629e+00   2.3879828e-02  -2.5809112e+00   2.4123850e-02  -2.5807834e+00   2.4048535e-02
+  5.0438504e+03  -2.5766085e+00   2.3814471e-02  -2.5766883e+00   2.4033056e-02  -2.5765489e+00   2.3969942e-02
+  5.6076745e+03  -2.5717533e+00   2.3787097e-02  -2.5717726e+00   2.3995997e-02  -2.5716337e+00   2.3935283e-02
+  6.2345254e+03  -2.5658345e+00   2.3816159e-02  -2.5659259e+00   2.4136946e-02  -2.5657792e+00   2.4076110e-02
+  6.9314484e+03  -2.5584699e+00   2.4062583e-02  -2.5586784e+00   2.4534829e-02  -2.5585448e+00   2.4483332e-02
+  7.7062767e+03  -2.5494999e+00   2.4614330e-02  -2.5496323e+00   2.5156997e-02  -2.5494980e+00   2.5107706e-02
+  8.5677186e+03  -2.5382574e+00   2.5989145e-02  -2.5380947e+00   2.6344473e-02  -2.5379297e+00   2.6295683e-02
+  9.5254565e+03  -2.5234453e+00   2.8452862e-02  -2.5225354e+00   2.8365270e-02  -2.5223466e+00   2.8324558e-02
+  1.0590254e+04  -2.5029798e+00   3.1972716e-02  -2.5008038e+00   3.1573322e-02  -2.5005982e+00   3.1531378e-02
+  1.1774080e+04  -2.4719823e+00   3.7183609e-02  -2.4687746e+00   3.6567599e-02  -2.4685400e+00   3.6524351e-02
+  1.3090240e+04  -2.4224620e+00   4.4503654e-02  -2.4202826e+00   4.3557684e-02  -2.4199965e+00   4.3511240e-02
+  1.4553525e+04  -2.3508047e+00   5.2823585e-02  -2.3498315e+00   5.1695877e-02  -2.3494733e+00   5.1608435e-02
+  1.6180383e+04  -2.2729733e+00   6.5550722e-02  -2.2736769e+00   6.5473622e-02  -2.2732954e+00   6.5303953e-02
+  1.7989098e+04  -2.2698761e+00   9.8036711e-02  -2.2712074e+00   9.7963778e-02  -2.2705603e+00   9.7797528e-02

soliket/szlike/beam.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+
+def create_synthetic_beam(fwhm=1.4):
+    """Returns beam of appropriate shape to use in mopc"""
+    ...
+
+
+"""from pixell.utils.beam_transform_to_profile"""
+
+
+def beam_transform_to_profile(bl, theta, normalize=False):
+    """Given the transform b(l) of a beam, evaluate its real space
+    angular profile
+    at the given radii theta."""
+    bl = np.asarray(bl)
+    ell = np.arange(bl.size)
+    x = np.cos(theta)
+    a = bl * (2 * ell + 1) / (4 * np.pi)
+    profile = np.polynomial.legendre.legval(x, a)
+    if normalize:
+        profile /= np.sum(a)
+    return profile
+
+
+def read_beam(filename, normalize=True):
+    """Reads the beam multipole info from file, returns real-space beam data"""
+    ell, b_ell = np.genfromtxt(filename, unpack=True)
+    theta_arcmin = np.linspace(0, 20, 100)  # arcmin
+    theta_rad = theta_arcmin / 60.0 * (np.pi / 180)  # radians
+    b = beam_transform_to_profile(b_ell, theta_rad)
+
+    if normalize:
+        integ = np.trapz(b * 2 * np.pi * theta_rad, x=theta_rad)
+        b /= integ
+
+    return b
+
+
+def f_beam(tht, b):
+    """Interpolates beam transform, b, from read_beam for a given theta value"""
+
+    theta_arcmin = np.linspace(0, 20, 100)
+    theta_rad = theta_arcmin / 60.0 * (np.pi / 180)
+
+    tht_in = theta_rad
+
+    return np.interp(tht, tht_in, b, period=np.pi)
+
+
+def f_beam_fft(beam_txt, ell):
+    """Reads in beam from txt file in ell, interpolates
+    Specific to ACT/CMASS
+    """
+    beam_file = np.genfromtxt(beam_txt)
+    ell_data = beam_file[:, 0]
+    beam_data = beam_file[:, 1]
+
+    return np.interp(ell, ell_data, beam_data)
+
+
+def f_response(beam_response_txt, ell):
+    """Reads in response from txt file in ell, interpolates
+    Specific to ACT/CMASS
+    """
+    response_file = np.genfromtxt(beam_response_txt)
+    ell_data = response_file[:, 0]
+    resp_data = response_file[:, 1]
+
+    return np.interp(ell, ell_data, resp_data)