Skip to content

Commit

Permalink
all comparison scripts between ACT and PLanck
Browse files Browse the repository at this point in the history
  • Loading branch information
Louis Thibaut committed Sep 17, 2024
1 parent 565e53a commit 5ace6f6
Show file tree
Hide file tree
Showing 5 changed files with 403 additions and 48 deletions.
5 changes: 3 additions & 2 deletions project/data_analysis/dr6xplanck.rst
Original file line number Diff line number Diff line change
Expand Up @@ -157,11 +157,12 @@ The first code is similar to the standard simulation spectra code, but it's resi
Comparison of ACT and Planck
*******************

In order to compare ACT and Planck power spectrum, once you have computed both ACTxNpipe and ACTxlegacy, grab the script in the planck folder and run
In order to compare ACT and Planck power spectrum, once you have computed all products for legacy and NPIPE, create two folders "dr6xlegacy" and "dr6xnpipe" containing the folders "spectra_leak_corr_planck_bias_corr", "covariances" and "best_fits" and run

.. code:: shell
salloc -N 1 -C cpu -q interactive -t 00:30:00
OMP_NUM_THREADS=256 srun -n 1 -c 256 --cpu_bind=cores python AxP_plots.py global_dr6v4xlegacy.dict
OMP_NUM_THREADS=256 srun -n 1 -c 256 --cpu_bind=cores python AxP_residuals.py global_dr6v4xlegacy.dict
OMP_NUM_THREADS=256 srun -n 1 -c 256 --cpu_bind=cores python AxP_comparison.py global_dr6v4xlegacy.dict
note that you have to specify in this script the location of the spectra and covariances of the npipe and legacy run.
65 changes: 23 additions & 42 deletions project/data_analysis/python/planck/AxP_comparison.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@

from pspy import so_dict, pspy_utils, so_cov, so_spectra
"""
Compare ACT DR6 spectra with Planck NPIPE and legacy spectra
"""
from pspy import so_dict, pspy_utils, so_cov
from pspipe_utils import consistency, best_fits, covariance, pspipe_list
import pickle
import scipy.stats as ss
import matplotlib.pyplot as plt
import numpy as np
Expand All @@ -13,31 +14,30 @@

spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
remove_first_bin = True
part = "part_c"
part = "part_b"


legacy_dir = "dr6xlegacy/"
npipe_dir = "dr6xnpipe/"
bestfit_dir = "best_fits"
spec_dir = "spectra_leak_corr_planck_bias_corr"
cov_dir = "covariances"

cov_dir_run_a = f"{legacy_dir}/covariances"
cov_dir_run_b = f"{npipe_dir}/covariances"
cov_dir_run_a = f"{legacy_dir}/{cov_dir}"
cov_dir_run_b = f"{npipe_dir}/{cov_dir}"

spec_dir_run_a = f"{legacy_dir}/spectra_leak_corr_planck_bias_corr"
spec_dir_run_b = f"{npipe_dir}/spectra_leak_corr_planck_bias_corr"
spec_dir_run_a = f"{legacy_dir}/{spec_dir}"
spec_dir_run_b = f"{npipe_dir}/{spec_dir}"

bf_dir_run_a = f"{legacy_dir}/best_fits"
bf_dir_run_b = f"{npipe_dir}/best_fits"
bf_dir_run_a = f"{legacy_dir}/{bestfit_dir}"
bf_dir_run_b = f"{npipe_dir}/{bestfit_dir}"

run_a_name = "legacy [new]"
run_b_name = "npipe [new]"
run_a_name = "legacy"
run_b_name = "npipe"

plot_dir = "AxP_plots"
pspy_utils.create_directory(f"{plot_dir}/{part}")

null_test_dir = "null_test"
plot_dir = "plots/planck_nulls"
bestfit_dir = "best_fits"

pspy_utils.create_directory(plot_dir)
pspy_utils.create_directory(null_test_dir)


null_list = {}
Expand All @@ -56,6 +56,8 @@
null_list["part_b"] += [["EE", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f143", "Planck_f143"]]
null_list["part_b"] += [["EE", "dr6_pa6_f150", "Planck_f143", "Planck_f143", "Planck_f143"]]
null_list["part_b"] += [["EE", "dr6_pa5_f150", "Planck_f143", "Planck_f143", "Planck_f143"]]
null_list["part_b"] += [["EE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f217", "Planck_f217"]]
null_list["part_b"] += [["EE", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f217", "Planck_f217"]]

null_list["part_b"] += [["ET", "dr6_pa6_f090", "Planck_f100", "Planck_f100", "Planck_f100"]]
null_list["part_b"] += [["ET", "dr6_pa5_f090", "Planck_f100", "Planck_f100", "Planck_f100"]]
Expand All @@ -68,18 +70,8 @@
null_list["part_b"] += [["TE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f143", "Planck_f143"]]
null_list["part_b"] += [["TE", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f143", "Planck_f143"]]
null_list["part_b"] += [["TE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f143", "Planck_f143"]]

null_list["part_c"] = []
null_list["part_c"] += [["TT", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["TT", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["TT", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["TT", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["TT", "dr6_pa4_f220", "Planck_f217", "Planck_f217", "Planck_f217"]]

null_list["part_c"] += [["EE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["EE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f100", "Planck_f100"]]
null_list["part_c"] += [["ET", "dr6_pa6_f150", "Planck_f143", "Planck_f143", "Planck_f143"]]
null_list["part_c"] += [["ET", "dr6_pa5_f090", "Planck_f100", "Planck_f100", "Planck_f100"]]
null_list["part_b"] += [["TE", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f217", "Planck_f217"]]
null_list["part_b"] += [["TE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f217", "Planck_f217"]]



Expand Down Expand Up @@ -242,7 +234,6 @@
ls="None", marker = ".",
label=f"[{run_b_name} run] {name} [$\chi^2 = {{{chi2_run_b:.1f}}}/{{{ndof}}}$ (${{{pte_run_b:.4f}}}$)]", color="darkorange", alpha=0.7)

# plt.plot(lb_fg, res_fg_b * lb_fg ** l_pow, color="gray")
plt.plot(lb_fg, lb_fg*0, color="black")
plt.legend()
plt.axvspan(xmin=0, xmax=xleft, color="gray", alpha=0.7)
Expand All @@ -256,17 +247,7 @@
plt.xlabel(r"$\ell$", fontsize=18)
plt.ylabel(r"$\ell^{%d} \Delta D_\ell^\mathrm{%s}$" % (l_pow, mode), fontsize=18)
plt.tight_layout()
# plt.subplot(2,2,3)
# plt.imshow(corr, interpolation="nearest", cmap="coolwarm")
# plt.xticks(np.arange(len(lb)), lb, rotation=90, fontsize=6)
# plt.yticks(np.arange(len(lb)), lb, fontsize=6)
# plt.colorbar()
# plt.show()
plt.savefig(f"{plot_dir}/{part}_{fname}.png", bbox_inches='tight')
plt.savefig(f"{plot_dir}/{part}/{fname}.png", bbox_inches='tight')
plt.clf()
plt.close()





100 changes: 100 additions & 0 deletions project/data_analysis/python/planck/AxP_plots.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
"""
plot the Planck and DR6 EE power and TE power spectra
"""
from pspy import so_dict, pspy_utils, so_spectra
from pspipe_utils import best_fits
import matplotlib.pyplot as plt
import numpy as np
import sys
import AxP_utils
import matplotlib


matplotlib.rcParams["font.family"] = "serif"
matplotlib.rcParams["font.size"] = "40"

d = so_dict.so_dict()
d.read_from_file(sys.argv[1])

spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]

legacy_dir = "dr6xlegacy/"
npipe_dir = "dr6xnpipe/"
bestfit_dir = "best_fits"
spec_dir = "spectra_leak_corr_planck_bias_corr"
cov_dir = "covariances"

plot_dir = "AxP_plots"
pspy_utils.create_directory(plot_dir)

cov_dir_legacy = f"{legacy_dir}/{cov_dir}"
cov_dir_npipe = f"{npipe_dir}/{cov_dir}"

spec_dir_legacy = f"{legacy_dir}/{spec_dir}"
spec_dir_npipe = f"{npipe_dir}/{spec_dir}"

bf_dir_legacy = f"{legacy_dir}/{bestfit_dir}"


cov_type_list = ["analytic_cov", "mc_cov", "leakage_cov"]

lth, psth = so_spectra.read_ps(f"{bf_dir_legacy}/cmb.dat", spectra=spectra)


my_ylim = {}
my_ylim["EE"] = (0,45)
my_ylim["TE"] = (-150,130)



fig = plt.figure(figsize=(40,40))
ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
ax_list = [ax1, ax2]
count = 0
for ax, mode in zip(ax_list, ["EE", "TE"]):

l, Db_legacy, sigma_legacy = AxP_utils.read_ps_and_sigma(spec_dir_legacy, cov_dir_legacy, "Planck_f143", "Planck_f143", mode, cov_type_list)
l, Db_npipe, sigma_npipe = AxP_utils.read_ps_and_sigma(spec_dir_npipe, cov_dir_npipe, "Planck_f143", "Planck_f143", mode, cov_type_list)

l, Db_pa6_f150, sigma_pa6_f150 = AxP_utils.read_ps_and_sigma(spec_dir_legacy, cov_dir_legacy, "dr6_pa6_f150", "dr6_pa6_f150", mode, cov_type_list)
l, Db_pa5_f150, sigma_pa5_f150 = AxP_utils.read_ps_and_sigma(spec_dir_npipe, cov_dir_npipe, "dr6_pa5_f150", "dr6_pa5_f150", mode, cov_type_list)
l, Db_pa6_f090, sigma_pa6_f090 = AxP_utils.read_ps_and_sigma(spec_dir_legacy, cov_dir_legacy, "dr6_pa6_f090", "dr6_pa6_f090", mode, cov_type_list)
l, Db_pa5_f090, sigma_pa5_f090 = AxP_utils.read_ps_and_sigma(spec_dir_npipe, cov_dir_npipe, "dr6_pa5_f090", "dr6_pa5_f090", mode, cov_type_list)
#plt.subplot(2,1,1 +count)
ax.errorbar(l-15, Db_legacy, sigma_legacy, fmt=".", color="black", alpha=0.7, label="Planck* legacy 143 GHz")
ax.errorbar(l-10, Db_pa6_f150, sigma_pa6_f150, fmt=".", color="purple", label="dr6 pa6-f150")
ax.errorbar(l-5, Db_pa5_f150, sigma_pa5_f150, fmt=".", color="blue", label="dr6 pa5-f150")
ax.errorbar(l, Db_pa6_f090, sigma_pa6_f090, fmt=".", color="darkorange", label="dr6 pa6-f090")
ax.errorbar(l+5, Db_pa5_f090, sigma_pa5_f090, fmt=".", color="red", label="dr6 pa5-f090")
ax.errorbar(l+10, Db_npipe, sigma_npipe, fmt=".", color="dimgrey", alpha=0.7, label="Planck* NPIPE 143 GHz")

ax.plot(lth, psth[mode], color="gray", linestyle="--")
ax.set_ylim(my_ylim[mode])
ax.set_xlim(0, 2000)
if count == 0:
ax.legend(fontsize=40)
ax.set_xticks([])
yticks = ax.yaxis.get_major_ticks()
yticks[0].label1.set_visible(False)

ax.set_ylabel(r"$D^{%s}_\ell$" % mode, fontsize=50)

if count == 1:
ax.set_xlabel(r"$\ell$", fontsize=50)
yticks = ax.yaxis.get_major_ticks()
yticks[-1].label1.set_visible(False)

count+=1
plt.subplots_adjust(wspace=0, hspace=0)

plt.savefig(f"{plot_dir}/ACT_and_planck.png",bbox_inches='tight')
plt.clf()
plt.close()







Loading

0 comments on commit 5ace6f6

Please sign in to comment.