Added calculation of power and trispectrum for dusty sources

pspipe_utils/radio_sources.py → pspipe_utils/poisson_sources.py

@@ -7,6 +7,8 @@
 from pixell import enmap, utils
 from scipy.interpolate import InterpolatedUnivariateSpline
 from matplotlib import rcParams
+
+#alex hacking this
 from . import get_data_path
 
 
@@ -16,8 +18,60 @@
 rcParams["axes.titlesize"] = 16
 
 ref_freq_radio_GHz = 148
+ref_freq_dusty_GHz = 217
+
+
 tucci_file_name = f"{get_data_path()}/radio_source/ns_148GHz_modC2Ex.dat"
 
+bethermin_file_name = f"{get_data_path()}/radio_source/Bethermin2012model_grids.save"
+
+def read_bethermin_source_distrib(lam=1380.0, plot_fname=None):
+    """
+    Read in the source distribution from Bethermin et al 2012
+    The model is described in this paper https://arxiv.org/abs/1208.6512 and the data products come from https://irfu.cea.fr/dap/Phocea/Page/index.php?id=537
+    
+    Parameters
+    ----------
+    lam : float
+        The wavelenth, in um, at which we would like the counts.  Must be one of 
+        160., 250., 350., 500., 550., 850., 1100., 1380., 2097., 3000., 210000.    
+        Default is 1380 um, i.e. 217GHz.
+    plot_fname : str
+        if not None save the plot in plot_fname
+    """
+    from  scipy.io import readsav
+
+    bethermin = readsav(bethermin_file_name)
+
+    #Figure out the index of the Bethermin model grid of the wavelength that we are interested in.  Note that 1100 um = 220 GHz
+
+    lam_idx = np.where(bethermin['lambda'   ] == lam)[0][0]
+
+    #bethermin['dndsnudz'] has size [n_lambda, n_z, n_S] or namely [16, 210, 200].
+    #Units are gal/Jy/sr, so there is actually a 'dOmega' in the denominator
+    #Let's sum over redshift
+    dz = np.gradient(bethermin['z'])
+    #Note the broadcasting of dz here, from shape=(210) to shape=(210,1)
+    dNdSdOmega = np.nansum(bethermin['dndsnudz'][lam_idx, :, :] * dz[:,None], axis = 0)
+
+
+    S = bethermin['snu']
+
+
+    if plot_fname is not None:
+        plt.figure(figsize=(15, 10))
+        plt.loglog()
+        plt.ylim(1e-3, 200)
+        plt.xlim(0.003, 1.5)
+        plt.plot(S, dNdSdOmega * S ** (5/2), "--", color="red")
+        plt.xlabel("S (Flux [Jy])", fontsize=16)
+        plt.ylabel(r"$S^{5/2}\frac{dN}{dS d\Omega} [Jy^{3/2} sr^{-1}]$", fontsize=22)
+        plt.savefig(plot_fname, bbox_inches="tight")
+        plt.clf()
+        plt.close()
+
+    return S, dNdSdOmega
+
 def read_tucci_source_distrib(plot_fname=None):
     """
     Read the source distribution from tucci et al (https://arxiv.org/pdf/1103.5707.pdf) with flux (S) at ref frequency 148 GHz
@@ -46,7 +100,8 @@ def read_tucci_source_distrib(plot_fname=None):
         plt.close()
 
     return S, dNdSdOmega
-
+
+
 def convert_Jy_per_str_to_muK_cmb(freq_GHz):
     """
     Convert from Jy/str to muK CMB at the observation frequency in GHz
@@ -60,7 +115,7 @@ def convert_Jy_per_str_to_muK_cmb(freq_GHz):
 
     return 10 ** 6 / K_to_Jy_per_str
 
-def get_poisson_power(S, dNdSdOmega, plot_fname=None):
+def get_poisson_power(S, dNdSdOmega, plot_fname=None, ref_freq_GHz=ref_freq_radio_GHz):
 
     """
     Get the expected poisson power as a function of Smax: the maximal flux to consider in the integral
@@ -75,9 +130,11 @@ def get_poisson_power(S, dNdSdOmega, plot_fname=None):
         source distribution dN/(dS dOmega) corresponding to S
     plot_fname : str
         if not None save the plot in plot_fname
+    ref_freq_GHz : float
+        Reference frequency, defaults to ref_freq_ratio_GHz defined above
     """
 
-    Jy_per_str_to_muK = convert_Jy_per_str_to_muK_cmb(ref_freq_radio_GHz)
+    Jy_per_str_to_muK = convert_Jy_per_str_to_muK_cmb(ref_freq_GHz)
     dS = np.gradient(S)
     poisson_power = (Jy_per_str_to_muK) ** 2 * np.cumsum(dS * S ** 2 * dNdSdOmega)
 
@@ -91,14 +148,14 @@ def get_poisson_power(S, dNdSdOmega, plot_fname=None):
         plt.ylim(10 ** -1, 5 * 10 ** 2)
         plt.plot(S * 10 ** 3, poisson_power * fac0)
         plt.xlabel(r"$S_{max}$ (Flux cut [mJy])", fontsize=16)
-        plt.ylabel(r"$a_{s} [ {%d} GHz] $" % ref_freq_radio_GHz, fontsize=22)
+        plt.ylabel(r"$a_{s} [ {%d} GHz] $" % ref_freq_GHz, fontsize=22)
         plt.savefig(plot_fname, bbox_inches="tight")
         plt.clf()
         plt.close()
 
     return poisson_power
 
-def get_trispectrum(S, dNdSdOmega):
+def get_trispectrum(S, dNdSdOmega, ref_freq_GHz = ref_freq_radio_GHz):
 
     """
     Get the expected trispectrum  as a function of Smax: the maximal flux to consider in the integral
@@ -109,9 +166,10 @@ def get_trispectrum(S, dNdSdOmega):
         source flux at the ref frequency in Jy
     dNdSdOmega : 1d array
         source distribution dN/(dS dOmega) corresponding to S
+    ref_freq_GHz : Reference frequency at which to convert from Jy/Sr to uK.  If not explicitly given, this is assumed to be the value of ref_freq_radio_GHz above.
     """
 
-    Jy_per_str_to_muK = convert_Jy_per_str_to_muK_cmb(ref_freq_radio_GHz)
+    Jy_per_str_to_muK = convert_Jy_per_str_to_muK_cmb(ref_freq_GHz)
     dS = np.gradient(S)
     trispectrum = (Jy_per_str_to_muK) ** 4 * np.cumsum(dS * S ** 4 * dNdSdOmega)
 

setup.py

@@ -16,7 +16,7 @@
     python_requires=">=3.5",
     install_requires=[
         "pspy>=1.5.3",
-        "mflike>=0.8.2",
+        # "mflike>=0.8.2",
     ],
     package_data={"": ["data/**"]},
 )

tutorials/tuto_dusty.py

@@ -0,0 +1,53 @@
+"""
+this is a small simulation script to check that the  radio power and trispectrum of simualtion match the theoretical expectation
+"""
+
+import numpy as np
+import pylab as plt
+import scipy
+from pixell import enmap, curvedsky
+from pspy import so_map, so_window, so_mcm, sph_tools, so_spectra, pspy_utils, so_cov
+from pspipe_utils import poisson_sources, get_data_path
+
+test_dir = "test_poisson"
+pspy_utils.create_directory(test_dir)
+ref_freq_dusty_GHz = 217
+ref_lam_dusty_um = 1380.0
+#check that this 
+if not np.isclose((ref_freq_dusty_GHz  * 1e9 )* (ref_lam_dusty_um * 1e-6), 3e8, rtol = .01):
+    raise "Check values of ref freq and lambda"
+
+
+S, dNdSdOmega = poisson_sources.read_bethermin_source_distrib(lam = ref_lam_dusty_um, plot_fname=f"{test_dir}/source_distrib_dusty.png")
+
+#Scale the 15 mJy cut at 148 GHz to 217 GHz.  Asssume we are in the RJ regime and beta = 1.6
+Smax_148GHz = 0.015
+Smax_dusty_217GHz = (0.015) * (217 / 148)**(2 + 1.6)
+
+poisson_power = poisson_sources.get_poisson_power(S, dNdSdOmega, plot_fname=f"{test_dir}/as_dusty.png", ref_freq_GHz=ref_freq_dusty_GHz)
+trispectrum = poisson_sources.get_trispectrum(S, dNdSdOmega, ref_freq_GHz=ref_freq_dusty_GHz)
+
+poisson_power_cut, trispectrum_cut = poisson_sources.get_power_and_trispectrum_at_Smax(S, poisson_power, trispectrum, Smax=Smax_dusty_217GHz)
+print(f"poisson_power_cut: {poisson_power_cut}", f"trispectrum_cut: {trispectrum_cut}")
+
+dunkley_dusty_data = {"S_cut_Jy" : 15e-3, "D3000" : 90, "sigmaD3000" : 10} # extracted from Dunkley et al. 
+
+l0 = 3000 # pivot scale for the fg amplitude
+fac0 = (l0 * (l0 + 1)) / (2 * np.pi)
+
+plt.figure(figsize=(12, 10))
+plt.loglog(S * 1e3, poisson_power * fac0) # a_s in Dunkley is computed at l=3000, in Dl unit
+plt.ylabel("$D^{217GHz}_{\ell = 3000}$  [$\mu K^{2}$]", fontsize=22)
+plt.xlabel("$S_{max}$ (mJy)", fontsize=22)
+plt.xlim([1e0, 1e3])
+plt.ylim(0.1, 1000)
+plt.errorbar(dunkley_dusty_data["S_cut_Jy"] * 1e3, dunkley_dusty_data["D3000"], dunkley_dusty_data["sigmaD3000"], fmt=".", label = "ACT 2013 dusty (Dunkley+)")
+plt.errorbar(dunkley_dusty_data["S_cut_Jy"] * 1e3, poisson_power_cut * fac0, fmt=".", label = "code prediction at 15 mJY")
+plt.legend(fontsize=22)
+plt.savefig(f"{test_dir}/check_dunkley_dusty.png")
+plt.clf()
+plt.close()
+
+#Note, the rest of the tuto_radio file contains code to make Poisson sims and check these predictions.  
+#These are not included here - there are far more dusty sources per unit area than radio, 
+#so this will take too long and would ultimately show the same thing.

tutorials/tuto_radio.py

@@ -7,18 +7,18 @@
 import scipy
 from pixell import enmap, curvedsky
 from pspy import so_map, so_window, so_mcm, sph_tools, so_spectra, pspy_utils, so_cov
-from pspipe_utils import radio_sources, get_data_path
+from pspipe_utils import poisson_sources, get_data_path
 
 test_dir = "test_poisson"
 pspy_utils.create_directory(test_dir)
 
-S, dNdSdOmega = radio_sources.read_tucci_source_distrib(plot_fname=f"{test_dir}/source_distrib.png")
+S, dNdSdOmega = poisson_sources.read_tucci_source_distrib(plot_fname=f"{test_dir}/source_distrib_radio.png")
 
 
-poisson_power = radio_sources.get_poisson_power(S, dNdSdOmega, plot_fname=f"{test_dir}/as.png")
-trispectrum = radio_sources.get_trispectrum(S, dNdSdOmega)
+poisson_power = poisson_sources.get_poisson_power(S, dNdSdOmega, plot_fname=f"{test_dir}/as_radio.png")
+trispectrum = poisson_sources.get_trispectrum(S, dNdSdOmega)
 
-poisson_power_15mJy, trispectrum_15mJy = radio_sources.get_power_and_trispectrum_at_Smax(S, poisson_power, trispectrum, Smax=0.015)
+poisson_power_15mJy, trispectrum_15mJy = poisson_sources.get_power_and_trispectrum_at_Smax(S, poisson_power, trispectrum, Smax=0.015)
 print(f"poisson_power_15mJy: {poisson_power_15mJy}", f"trispectrum_15mJy: {trispectrum_15mJy}")
 
 dunkley_radio_data = {"S_cut_Jy" : 15e-3, "D3000" : 3.1, "sigmaD3000" : .4} # extracted from Dunkley et al
@@ -35,7 +35,7 @@
 plt.errorbar(dunkley_radio_data["S_cut_Jy"] * 1e3, dunkley_radio_data["D3000"], dunkley_radio_data["sigmaD3000"], fmt=".", label = "ACT 2013 radio(Dunkley+)")
 plt.errorbar(dunkley_radio_data["S_cut_Jy"] * 1e3, poisson_power_15mJy * fac0, fmt=".", label = "code prediction at 15 mJY")
 plt.legend(fontsize=22)
-plt.savefig(f"{test_dir}/check_dunkley.png")
+plt.savefig(f"{test_dir}/check_dunkley_radio.png")
 plt.clf()
 plt.close()
 
@@ -48,7 +48,7 @@
 binning_file = f"{test_dir}/binning.dat"
 n_splits = 2
 ref_freq = 148
-Jy_per_str_to_muK = radio_sources.convert_Jy_per_str_to_muK_cmb(ref_freq)
+Jy_per_str_to_muK = poisson_sources.convert_Jy_per_str_to_muK_cmb(ref_freq)
 
 
 # try to emulate the actual DR6 window (although at 8 times lower res)
@@ -95,7 +95,7 @@
         Clth_dict[id1 + id2] = Clth_dict[id1 + id2][:-2]
 
 
-source_mean_numbers = radio_sources.get_mean_number_of_source(template_car, S, dNdSdOmega, plot_fname=f"{test_dir}/N_source.png")
+source_mean_numbers = poisson_sources.get_mean_number_of_source(template_car, S, dNdSdOmega, plot_fname=f"{test_dir}/N_source.png")
 
 mbb_inv, Bbl = so_mcm.mcm_and_bbl_spin0(window, binning_file, lmax=lmax, type=ps_type, niter=niter)