In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches
import astropy.io.fits
import healpy as hp
import scipy.special
import scipy.interpolate
import os
import collections
pi = np.pi
In [ ]:
def create_tSZ_catalog(fields, y_map_file, y_map_mask_file, shear_path, pad, output_path, survey):
y_map = hp.read_map(y_map_file, field=0)
y_map_mask = hp.read_map(y_map_mask_file, field=4)
for field in fields:
shear_filename = shear_path + "{}.fits".format(field)
field_info = np.loadtxt(shear_filename + "_sizeinfo.dat")
print("Field: " + field)
bbox = np.array([[field_info[0,0]-pad, field_info[1,0]-pad],
[field_info[0,1]+pad, field_info[1,0]-pad],
[field_info[0,1]+pad, field_info[1,1]+pad],
[field_info[0,0]-pad, field_info[1,1]+pad]])
tSZ_output = output_path + "{}.fits".format(field)
tpcf_formats.create_standard_catalog_from_healpix_map(y_map, bbox, tSZ_output, y_map_mask, use_fits=True, field=field, survey=survey)
del y_map
del y_map_mask
In [ ]:
pad_RCSLenS = 4.0
y_map_file = "/shared_data/Planck_data/COM_Compton_SZMap_R2.00/milca_ymaps.fits"
mask_file = "/shared_data/Planck_data/COM_Compton_SZMap_R2.00/masks.fits"
RCSLenS_fields = ["CDE0047", "CDE0133", "CDE0310", "CDE0357", "CDE1040", "CDE1111", "CDE1303", "CDE1514", "CDE1613", "CDE1645", "CDE2143", "CDE2329", "CDE2338", "CSP0320"]
shear_path = "/data/aha25/shear_catalogs/RCSLenS/mag_18_30/"
output_path = "/data/aha25/tSZ/tSZ_data/RCSLenS/pad{}/".format(pad_RCSLenS)
create_tSZ_catalog(RCSLenS_fields, y_map_file, mask_file, shear_path, pad_RCSLenS, output_path, survey="RCSLenS")
In [ ]:
RCSLenS_fields = ["CDE0047", "CDE0133", "CDE0310", "CDE0357", "CDE1040", "CDE1111", "CDE1303", "CDE1514", "CDE1613", "CDE1645", "CDE2143", "CDE2329", "CDE2338", "CSP0320"]
for field in RCSLenS_fields:
lenses_filename = "/home/aha25/y_maps/{}_y_planck_milca_2048_big_new.fits".format(field)
tSZ_map_output = "/data/aha25/tSZ/tSZ_data/RCSLenS_grid/pad2/{}.fits".format(field)
tpcf_formats.create_catalog_from_grids(mode="lenses",
filename=lenses_filename, output_filename=tSZ_map_output,
pixel_size=1.0/60.0, weight_filename=None, remove_mean=True,
survey="RCSLenS", field=field)
In [ ]:
tpcf_path = "/data/aha25/tpcf/calculate_shear_2pcf"
tpcf_path_debug = "/data/aha25/tpcf/calculate_shear_2pcf_debug"
bootstrap_configs = [#{"n_resample" : 10000, "type" : "field-equal-weight-blocks", "n_x_block" : None, "n_y_block" : None, "n_block_total" : 40, "supersampling" : None, "n_run" : 1},
{"n_resample" : 10000, "type" : "field-equal-weight-blocks", "n_x_block" : None, "n_y_block" : None, "n_block_total" : 80, "supersampling" : None, "n_run" : 1},
{"n_resample" : 10000, "type" : "field-equal-weight-blocks", "n_x_block" : None, "n_y_block" : None, "n_block_total" : 160, "supersampling" : None, "n_run" : 1},
{"n_resample" : 10000, "type" : "field-rectangular-blocks", "n_x_block" : 5, "n_y_block" : 5, "n_block_total" : None, "supersampling" : None, "n_run" : 1},
{"n_resample" : 10000, "type" : "field-rectangular-blocks", "n_x_block" : 10, "n_y_block" : 10, "n_block_total" : None, "supersampling" : None, "n_run" : 1}]
n_threads = 10
In [ ]:
marks_file = "/data/aha25/tpcf_marks/CFHTLenS/CFHTLenS_kappa_scale_{}_theta_{}_{}_nbin_{}_{}_pad{}_marks.fits".format(scale_CFHTLenS, theta_min_CFHTLenS*60.0, theta_max_CFHTLenS*60.0, n_bin_CFHTLenS, bin_type_CFHTLenS, pad_CFHTLenS)
foreground_path = "/data/aha25/tSZ/tSZ_data/CFHTLenS/pad{}/".format(pad_CFHTLenS)
tpcf_output_path = "/data/aha25/tSZ/tpcf/CFHTLenS_kappa/scale{}/".format(scale_CFHTLenS)
tpcf_formats.mkdirs(tpcf_output_path)
tpcf_driver.run_marks_crosscorrelation("kappa", foreground_path=foreground_path, marks_file=marks_file, tpcf_output_path=tpcf_output_path,
fields=CFHTLenS_fields,
n_bin=n_bin_CFHTLenS, theta_min=theta_min_CFHTLenS, theta_max=theta_max_CFHTLenS, logspaced=logspaced_CFHTLenS,
bootstrap_configs=bootstrap_configs,
n_thread=n_threads, tpcf_path=tpcf_path, verbose=True)
In [ ]:
tSZ_path = os.path.join("/data/aha25/tSZ/tSZ_data/sims/", cosmology_sims, feedback_model_sims)
shear_path = os.path.join("/data/aha25/shear_catalogs/cosmo-OWLS/", n_z_sims, cosmology_sims, feedback_model_sims)
kappa_path = os.path.join("/data/aha25/massmaps/cosmo-OWLS/", n_z_sims, cosmology_sims, feedback_model_sims)
for i in range(n_LOS):
tpcf_output_base_path = "/data/aha25/tSZ/tpcf/cosmo-OWLS/{}/".format(i)
if not os.access(tpcf_output_base_path, os.F_OK):
os.mkdir(tpcf_output_base_path)
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims))
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims))
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims))
tpcf_output_path = os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims)
LOS = ["cone_{}".format(i)]
tpcf_driver.run_standard_crosscorrelation(mode="tangential-shear",
foreground_path=tSZ_path, background_path=shear_path, tpcf_output_path=tpcf_output_path,
fields=LOS,
n_bin=n_bin_sims, theta_min=theta_min_sims, theta_max=theta_max_sims, logspaced=logspaced_sims,
bootstrap_configs=None,
spherical_coordinates=False, left_handed_coordinates=True, kdtree=False,
calculate_tpcf=True, n_thread=20)
for i in range(n_LOS):
tpcf_output_base_path = "/data/aha25/tSZ/tpcf/cosmo-OWLS_kappa/{}/".format(i)
if not os.access(tpcf_output_base_path, os.F_OK):
os.mkdir(tpcf_output_base_path)
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims))
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims))
if not os.access(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims), os.F_OK):
os.mkdir(os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims))
tpcf_output_path = os.path.join(tpcf_output_base_path, n_z_sims, cosmology_sims, feedback_model_sims)
LOS = ["cone_{}".format(i)]
tpcf_driver.run_standard_crosscorrelation(mode="kappa",
foreground_path=tSZ_path, background_path=kappa_path, tpcf_output_path=tpcf_output_path,
fields=LOS,
n_bin=n_bin_sims, theta_min=theta_min_sims, theta_max=theta_max_sims, logspaced=logspaced_sims,
bootstrap_configs=None,
spherical_coordinates=False, left_handed_coordinates=False, kdtree=False,
calculate_tpcf=True, n_thread=20)
In [ ]:
C_ell_Planck_AGN80_bin2_raw = np.loadtxt("/home/aha25/RCSLenS_new/Ian_sims/FT/AGN_Planck_Cl_kappa_bin2_y_smooth.dat")
C_ell_Planck_AGN80_bin2_raw = C_ell_Planck_AGN80_bin2_raw[1:70]
C_ell_WMAP7_AGN80_bin2_raw = np.loadtxt("/home/aha25/RCSLenS_new/Ian_sims/FT/AGN_WMAP7_Cl_kappa_bin2_y_smooth.dat")
C_ell_WMAP7_AGN80_bin2_raw = C_ell_WMAP7_AGN80_bin2_raw[1:70]
C_ell_WMAP7_AGN80_CFHTLenS_raw = np.loadtxt("/data/aha25/tSZ/Cl_kappa_y_simulation/AGN_WMAP7_Cl_kappa_y.dat")
C_ell_WMAP7_AGN80_CFHTLenS_raw = C_ell_WMAP7_AGN80_CFHTLenS_raw[1:70]
In [ ]:
plt.loglog(C_ell_WMAP7_AGN80_bin2_raw[:,0], C_ell_WMAP7_AGN80_bin2_raw[:,0]*C_ell_WMAP7_AGN80_bin2_raw[:,1])
plt.loglog(C_ell_WMAP7_AGN80_CFHTLenS_raw[:,0], C_ell_WMAP7_AGN80_CFHTLenS_raw[:,0]*C_ell_WMAP7_AGN80_CFHTLenS_raw[:,1])
In [ ]:
ell_log_C_ell_Planck_AGN80_bin2_intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(C_ell_Planck_AGN80_bin2_raw[:,0]), np.log(C_ell_Planck_AGN80_bin2_raw[:,0]*C_ell_Planck_AGN80_bin2_raw[:,1]), ext="const")
ell_log_C_ell_WMAP7_AGN80_bin2_intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(C_ell_WMAP7_AGN80_bin2_raw[:,0]), np.log(C_ell_WMAP7_AGN80_bin2_raw[:,0]*C_ell_WMAP7_AGN80_bin2_raw[:,1]), ext="const")
ell_log_C_ell_WMAP7_AGN80_CFHTLenS_intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(C_ell_WMAP7_AGN80_CFHTLenS_raw[:,0]), np.log(C_ell_WMAP7_AGN80_CFHTLenS_raw[:,0]*C_ell_WMAP7_AGN80_CFHTLenS_raw[:,1]), ext="const")
C_ell_Planck_AGN80_bin2 = lambda ell: np.exp(ell_log_C_ell_Planck_AGN80_bin2_intp(np.log(ell)))/ell
C_ell_WMAP7_AGN80_bin2 = lambda ell: np.exp(ell_log_C_ell_WMAP7_AGN80_bin2_intp(np.log(ell)))/ell
C_ell_WMAP7_AGN80_CFHTLenS = lambda ell: np.exp(ell_log_C_ell_WMAP7_AGN80_CFHTLenS_intp(np.log(ell)))/ell
ell = np.arange(1, C_ell_Planck_AGN80_bin2_raw[-1,0]+1)
theta_th = np.linspace(0.1, 181.0, 80.0)
xi_th_kappa_Planck_AGN80_bin2 = np.zeros_like(theta_th)
xi_th_g_t_Planck_AGN80_bin2 = np.zeros_like(theta_th)
xi_th_kappa_WMAP7_AGN80_bin2 = np.zeros_like(theta_th)
xi_th_g_t_WMAP7_AGN80_bin2 = np.zeros_like(theta_th)
xi_th_kappa_WMAP7_AGN80_CFHTLenS = np.zeros_like(theta_th)
xi_th_g_t_WMAP7_AGN80_CFHTLenS = np.zeros_like(theta_th)
sigma_Planck = 9.5/60.0/180.0*pi/(2.0*np.sqrt(2.0* np.log(2.0)))
Planck_smoothing = np.exp(-sigma_Planck**2*ell**2/2.0)
sigma_RCSLenS = 2.0/60.0/180.0*pi
RCSLenS_smoothing = np.exp(-sigma_RCSLenS**2*ell**2/2.0)
for i in range(len(theta_th)):
xi_th_kappa_Planck_AGN80_bin2[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_Planck_AGN80_bin2(ell) * ell * scipy.special.jv(0, ell*theta_th[i]/60.0/180.0*pi))
xi_th_g_t_Planck_AGN80_bin2[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_Planck_AGN80_bin2(ell) * ell * scipy.special.jv(2, ell*theta_th[i]/60.0/180.0*pi))
xi_th_kappa_WMAP7_AGN80_bin2[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_bin2(ell) * ell * scipy.special.jv(0, ell*theta_th[i]/60.0/180.0*pi))
xi_th_g_t_WMAP7_AGN80_bin2[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_bin2(ell) * ell * scipy.special.jv(2, ell*theta_th[i]/60.0/180.0*pi))
xi_th_kappa_WMAP7_AGN80_CFHTLenS[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_CFHTLenS(ell)*Planck_smoothing * ell * scipy.special.jv(0, ell*theta_th[i]/60.0/180.0*pi))
xi_th_g_t_WMAP7_AGN80_CFHTLenS[i] = 1.0e-11*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_CFHTLenS(ell)*Planck_smoothing * ell * scipy.special.jv(2, ell*theta_th[i]/60.0/180.0*pi))
In [ ]:
ell = C_ell_WMAP7_AGN80_CFHTLenS_raw[:,0]
sigma_Planck = 9.5/60.0/180.0*pi/(2.0*np.sqrt(2.0* np.log(2.0)))
Planck_smoothing = np.exp(-sigma_Planck**2*ell**2/2.0)
sigma_RCSLenS = 2.0/60.0/180.0*pi
RCSLenS_smoothing = np.exp(-sigma_RCSLenS**2*ell**2/2.0)
for i in range(len(theta_th)):
xi_th_kappa_WMAP7_AGN80_CFHTLenS[i] = 1.0e-9*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_CFHTLenS_raw[:,1]*Planck_smoothing*RCSLenS_smoothing * ell * scipy.special.jv(0, ell*theta_th[i]/60.0/180.0*pi))
xi_th_g_t_WMAP7_AGN80_CFHTLenS[i] = 1.0e-9*1.0/(2.0*pi) * np.sum(C_ell_WMAP7_AGN80_CFHTLenS_raw[:,1]*Planck_smoothing*RCSLenS_smoothing * ell * scipy.special.jv(2, ell*theta_th[i]/60.0/180.0*pi))
In [ ]:
def load_xi(base_path, covariance, n_bin):
xi = np.loadtxt(base_path + "xi.dat")
xi_E_cov = np.loadtxt(base_path + "{}".format(covariance))[:n_bin,:]
xi_B_cov = np.loadtxt(base_path + "{}".format(covariance))[n_bin:,:]
return {"xi" : xi, "xi_E_cov" : xi_E_cov, "xi_B_cov" : xi_B_cov}
def stack_simulations(base_path, N_sims, filter_type, n_bin):
sims_xi = np.zeros((N_sims, n_bin, 3))
xi = np.zeros((n_bin, 3))
for i in range(N_sims):
tmp = np.loadtxt(base_path + "{}/{}/xi.dat".format(i, filter_type))
sims_xi[i] = tmp[:,1:4]
xi[:,0] = tmp[:,0]
xi[:,1] = np.sum(sims_xi[:,:,0]*sims_xi[:,:,2], axis=0)/np.sum(sims_xi[:,:,2], axis=0)
xi[:,2] = np.sum(sims_xi[:,:,1]*sims_xi[:,:,2], axis=0)/np.sum(sims_xi[:,:,2], axis=0)
xi_E_cov = np.cov(sims_xi[:,:,0].T, ddof=1)
xi_B_cov = np.cov(sims_xi[:,:,1].T, ddof=1)
return {"xi" : xi, "xi_E_cov" : xi_E_cov, "xi_B_cov" : xi_B_cov}
In [ ]:
xi_kappa = collections.OrderedDict()
xi_shear = collections.OrderedDict()