From 5ace6f6d1eaa6a35f4fbff145ee927d695b97446 Mon Sep 17 00:00:00 2001 From: Louis Thibaut Date: Tue, 17 Sep 2024 11:27:50 +0200 Subject: [PATCH] all comparison scripts between ACT and PLanck --- project/data_analysis/dr6xplanck.rst | 5 +- .../python/planck/AxP_comparison.py | 65 ++--- .../data_analysis/python/planck/AxP_plots.py | 100 +++++++ .../python/planck/AxP_residuals.py | 273 ++++++++++++++++++ .../data_analysis/python/planck/AxP_utils.py | 8 +- 5 files changed, 403 insertions(+), 48 deletions(-) create mode 100644 project/data_analysis/python/planck/AxP_plots.py create mode 100644 project/data_analysis/python/planck/AxP_residuals.py diff --git a/project/data_analysis/dr6xplanck.rst b/project/data_analysis/dr6xplanck.rst index b87654e0..ff1d9191 100644 --- a/project/data_analysis/dr6xplanck.rst +++ b/project/data_analysis/dr6xplanck.rst @@ -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. diff --git a/project/data_analysis/python/planck/AxP_comparison.py b/project/data_analysis/python/planck/AxP_comparison.py index aaeb2dee..c2242c4e 100644 --- a/project/data_analysis/python/planck/AxP_comparison.py +++ b/project/data_analysis/python/planck/AxP_comparison.py @@ -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 @@ -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 = {} @@ -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"]] @@ -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"]] @@ -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) @@ -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() - - - - diff --git a/project/data_analysis/python/planck/AxP_plots.py b/project/data_analysis/python/planck/AxP_plots.py new file mode 100644 index 00000000..cf06b086 --- /dev/null +++ b/project/data_analysis/python/planck/AxP_plots.py @@ -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() + + + + + + + diff --git a/project/data_analysis/python/planck/AxP_residuals.py b/project/data_analysis/python/planck/AxP_residuals.py new file mode 100644 index 00000000..2e6f44aa --- /dev/null +++ b/project/data_analysis/python/planck/AxP_residuals.py @@ -0,0 +1,273 @@ +""" +Plot a suit of residuals for EE and TE between ACT and Planck +""" +from pspy import so_dict, pspy_utils, so_cov, so_spectra +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 +import sys +import AxP_utils +import matplotlib + +d = so_dict.so_dict() +d.read_from_file(sys.argv[1]) + +spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] +remove_first_bin = True +part = "EE" +matplotlib.rcParams["font.family"] = "serif" +matplotlib.rcParams["font.size"] = "20" + + +legacy_dir = "dr6xlegacy/" +npipe_dir = "dr6xnpipe/" + +cov_dir_run_a = f"{legacy_dir}/covariances" +cov_dir_run_b = f"{npipe_dir}/covariances" +bestfit_dir = "best_fits" +spec_dir = "spectra_leak_corr_planck_bias_corr" +cov_dir = "covariances" + +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}/{bestfit_dir}" +bf_dir_run_b = f"{npipe_dir}/{bestfit_dir}" + +run_a_name = "legacy" +run_b_name = "NPIPE" + +plot_dir = "AxP_plots" + +pspy_utils.create_directory(plot_dir) + + +null_list = {} +null_list["EE"] = [] +null_list["EE"] += [["EE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f100", "Planck_f100"]] +null_list["EE"] += [["EE", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f100", "Planck_f100"]] +null_list["EE"] += [["EE", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f100", "Planck_f100"]] +null_list["EE"] += [["EE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f100", "Planck_f100"]] + +null_list["EE"] += [["EE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f143", "Planck_f143"]] +null_list["EE"] += [["EE", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f143", "Planck_f143"]] +null_list["EE"] += [["EE", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f143", "Planck_f143"]] +null_list["EE"] += [["EE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f143", "Planck_f143"]] + +null_list["EE"] += [["EE", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f217", "Planck_f217"]] +null_list["EE"] += [["EE", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f217", "Planck_f217"]] +null_list["EE"] += [["EE", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f217", "Planck_f217"]] +null_list["EE"] += [["EE", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f217", "Planck_f217"]] + + +null_list["ET"] = [] + +null_list["ET"] += [["ET", "dr6_pa5_f090", "Planck_f100", "Planck_f100", "Planck_f100"]] +null_list["ET"] += [["ET", "dr6_pa5_f150", "Planck_f100", "Planck_f100", "Planck_f100"]] +null_list["ET"] += [["ET", "dr6_pa6_f090", "Planck_f100", "Planck_f100", "Planck_f100"]] +null_list["ET"] += [["ET", "dr6_pa6_f150", "Planck_f100", "Planck_f100", "Planck_f100"]] + +null_list["ET"] += [["ET", "dr6_pa5_f090", "Planck_f143", "Planck_f143", "Planck_f143"]] +null_list["ET"] += [["ET", "dr6_pa5_f150", "Planck_f143", "Planck_f143", "Planck_f143"]] +null_list["ET"] += [["ET", "dr6_pa6_f090", "Planck_f143", "Planck_f143", "Planck_f143"]] +null_list["ET"] += [["ET", "dr6_pa6_f150", "Planck_f143", "Planck_f143", "Planck_f143"]] + +null_list["ET"] += [["ET", "dr6_pa5_f090", "Planck_f217", "Planck_f217", "Planck_f217"]] +null_list["ET"] += [["ET", "dr6_pa5_f150", "Planck_f217", "Planck_f217", "Planck_f217"]] +null_list["ET"] += [["ET", "dr6_pa6_f090", "Planck_f217", "Planck_f217", "Planck_f217"]] +null_list["ET"] += [["ET", "dr6_pa6_f150", "Planck_f217", "Planck_f217", "Planck_f217"]] + +#null_list["ET"] += [["ET", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f100", "Planck_f100"]] +#null_list["ET"] += [["ET", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f100", "Planck_f100"]] +#null_list["ET"] += [["ET", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f100", "Planck_f100"]] +#null_list["ET"] += [["ET", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f100", "Planck_f100"]] + +#null_list["ET"] += [["ET", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f143", "Planck_f143"]] +#null_list["ET"] += [["ET", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f143", "Planck_f143"]] +#null_list["ET"] += [["ET", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f143", "Planck_f143"]] +#null_list["ET"] += [["ET", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f143", "Planck_f143"]] + +#null_list["ET"] += [["ET", "dr6_pa5_f090", "dr6_pa5_f090", "Planck_f217", "Planck_f217"]] +#null_list["ET"] += [["ET", "dr6_pa6_f090", "dr6_pa6_f090", "Planck_f217", "Planck_f217"]] +#null_list["ET"] += [["ET", "dr6_pa5_f150", "dr6_pa5_f150", "Planck_f217", "Planck_f217"]] +#null_list["ET"] += [["ET", "dr6_pa6_f150", "dr6_pa6_f150", "Planck_f217", "Planck_f217"]] + + + + +cov_type_list = ["analytic_cov", "mc_cov", "leakage_cov"] + +map_set_list = pspipe_list.get_map_set_list(d) + +lb, all_ps_run_a, all_cov_run_a = AxP_utils.read_data(map_set_list, spec_dir_run_a, cov_dir_run_a, cov_type_list, spectra) +lb, all_ps_run_b, all_cov_run_b = AxP_utils.read_data(map_set_list, spec_dir_run_b, cov_dir_run_b, cov_type_list, spectra) + +# Load foreground best fits +fg_file_name = bf_dir_run_a + "/fg_{}x{}.dat" +l_fg, fg_dict = best_fits.fg_dict_from_files(fg_file_name, map_set_list, d["lmax"], spectra=spectra) + +operations = {"diff": "ab-cd"} + +multipole_range, l_pows, y_lims = AxP_utils.get_plot_params() +y_lims_spec = {} +y_lims_spec["TT"] = [0, 2*10**9] +y_lims_spec["TE"] = [-150, 150] +y_lims_spec["ET"] = [-150, 150] +y_lims_spec["EE"] = [-10, 70] + + +print(f"we will do {len(null_list)} null tests") + + +plt.figure(figsize=(20, 25)) + +count = 0 +for null in null_list[part]: + + mode, ms1, ms2, ms3, ms4 = null + res_cov_dict_run_a, res_cov_dict_run_b = {}, {} + + l, Db1_run_a, sigma1_run_a = AxP_utils.read_ps_and_sigma(spec_dir_run_a, cov_dir_run_a, ms1, ms2, mode, cov_type_list) + l, Db2_run_a, sigma2_run_a = AxP_utils.read_ps_and_sigma(spec_dir_run_a, cov_dir_run_a, ms3, ms4, mode, cov_type_list) + + l, Db1_run_b, sigma1_run_b = AxP_utils.read_ps_and_sigma(spec_dir_run_b, cov_dir_run_b, ms1, ms2, mode, cov_type_list) + l, Db2_run_b, sigma2_run_b = AxP_utils.read_ps_and_sigma(spec_dir_run_b, cov_dir_run_b, ms3, ms4, mode, cov_type_list) + + if remove_first_bin: + Db1_run_a, sigma1_run_a, Db2_run_a, sigma2_run_a = Db1_run_a[1:], sigma1_run_a[1:], Db2_run_a[1:], sigma2_run_a[1:] + Db1_run_b, sigma1_run_b, Db2_run_b, sigma2_run_b = Db1_run_b[1:], sigma1_run_b[1:], Db2_run_b[1:], sigma2_run_b[1:] + + for cov in cov_type_list: + lb, res_ps_run_a, res_cov_dict_run_a[cov] = consistency.compare_spectra([ms1, ms2, ms3, ms4], + "ab-cd", + all_ps_run_a, + all_cov_run_a[cov], + mode = mode, + return_chi2=False) + + lb, res_ps_run_b, res_cov_dict_run_b[cov] = consistency.compare_spectra([ms1, ms2, ms3, ms4], + "ab-cd", + all_ps_run_b, + all_cov_run_b[cov], + mode = mode, + return_chi2=False) + + if remove_first_bin: + lb, res_cov_dict_run_a[cov] = lb[1:], res_cov_dict_run_a[cov][1:,1:] + res_ps_run_a = res_ps_run_a[1:] + res_cov_dict_run_b[cov] = res_cov_dict_run_b[cov][1:,1:] + res_ps_run_b = res_ps_run_b[1:] + + + res_fg = fg_dict[ms1, ms2][mode] - fg_dict[ms3, ms4][mode] + lb_fg, res_fg_b = pspy_utils.naive_binning(l_fg, res_fg, d["binning_file"], d["lmax"]) + + if remove_first_bin: + lb_fg, res_fg_b = lb_fg[1:], res_fg_b[1:] + + + name = "analytical" + if "mc_cov" in cov_type_list: + r_cov_run_a = covariance.correct_analytical_cov(res_cov_dict_run_a["analytic_cov"], + res_cov_dict_run_a["mc_cov"], + only_diag_corrections=True) + + r_cov_run_b = covariance.correct_analytical_cov(res_cov_dict_run_b["analytic_cov"], + res_cov_dict_run_b["mc_cov"], + only_diag_corrections=True) + + name += "+mc_corr" + if "beam_cov" in cov_type_list: + r_cov_run_a += res_cov_dict_run_a["beam_cov"] + r_cov_run_b += res_cov_dict_run_b["beam_cov"] + + name += "+beam" + + if "leakage_cov" in cov_type_list: + r_cov_run_a += res_cov_dict_run_a["leakage_cov"] + r_cov_run_b += res_cov_dict_run_b["leakage_cov"] + + name += "+leakage" + + corr_run_a = so_cov.cov2corr(r_cov_run_a) + + lmin, lmax = AxP_utils.get_lmin_lmax(null, multipole_range) + + # Plot residual and get chi2 + lrange = np.where((lb >= lmin) & (lb <= lmax))[0] + fname = f"diff_{mode}_{ms1}x{ms2}_{ms3}x{ms4}" + l_pow = l_pows[mode] + + expected_res = 0. + ylims = y_lims[mode] + if "pa4" in fname: + ylims = (y_lims[mode][0] * 5, y_lims[mode][1] * 5) + + + chi2_run_a = (res_ps_run_a[lrange] - res_fg_b[lrange]) @ np.linalg.inv(r_cov_run_a[np.ix_(lrange, lrange)]) @ (res_ps_run_a[lrange] - res_fg_b[lrange]) + chi2_run_b = (res_ps_run_b[lrange] - res_fg_b[lrange]) @ np.linalg.inv(r_cov_run_b[np.ix_(lrange, lrange)]) @ (res_ps_run_b[lrange] - res_fg_b[lrange]) + + ndof = len(lb[lrange]) + ndof -= 1 + + pte_run_a = 1 - ss.chi2(ndof).cdf(chi2_run_a) + pte_run_b = 1 - ss.chi2(ndof).cdf(chi2_run_b) + + xleft, xright = lb[lrange][0], lb[lrange][-1] + + if mode == "TT": + l_pow_plain = 2 + else: + l_pow_plain = 0 + + plt.subplot(6, 2, count+1) + plot_title = f"{ms3}x{ms4} - {ms1}x{ms2}" + plot_title = plot_title.replace("dr6_", "") + plot_title = plot_title.replace("-"," vs ") + plot_title = plot_title.replace("_", "-") + + plt.title(plot_title, fontsize=20) + + plt.errorbar(lb - 8, -(res_ps_run_a - res_fg_b) * lb ** l_pow, + yerr=np.sqrt(r_cov_run_a.diagonal()) * lb ** l_pow, + ls="None", marker = ".", + label=f"{run_a_name} p= {pte_run_a:.3f}", color="lightseagreen") + + plt.errorbar(lb + 8, -(res_ps_run_b - res_fg_b) * lb ** l_pow, + yerr=np.sqrt(r_cov_run_b.diagonal()) * lb ** l_pow, + ls="None", marker = ".", + label=f"{run_b_name} p= {pte_run_b:.3f}", color="blue", alpha=0.7) + + plt.plot(lb_fg, lb_fg*0, linestyle="--", color="black") + plt.legend(fontsize=15, loc="lower right") + plt.axvspan(xmin=0, xmax=xleft, color="gray", alpha=0.7) + if xright != lb[-1]: + plt.axvspan(xmin=xright, xmax=lb[-1], color="gray", alpha=0.7) + + # plt.title(plot_title.replace("dr6_", "")) + plt.xlim(300, 2000) + if ylims is not None: + plt.ylim(*ylims) + plt.xlabel(r"$\ell$", fontsize=25) + + if l_pow != 0: + plt.ylabel(r"$\ell^{%d} \Delta D_\ell^\mathrm{%s}$" % (l_pow, mode), fontsize=25) + else: + plt.ylabel(r"$\Delta D_\ell^\mathrm{%s}$" % (mode), fontsize=25) + + print(f"{ms3}x{ms4} - {ms1}x{ms2}, {run_a_name} p= {pte_run_a:.3f}") + print(f"{ms3}x{ms4} - {ms1}x{ms2}, {run_b_name} p= {pte_run_b:.3f}") + + count += 1 + +plt.tight_layout() +plt.savefig(f"{plot_dir}/{part}.png", bbox_inches='tight') +plt.clf() +plt.close() + + + + + diff --git a/project/data_analysis/python/planck/AxP_utils.py b/project/data_analysis/python/planck/AxP_utils.py index 20da346e..ff1887fe 100644 --- a/project/data_analysis/python/planck/AxP_utils.py +++ b/project/data_analysis/python/planck/AxP_utils.py @@ -99,17 +99,17 @@ def get_plot_params(): }, "Planck_f100": { "T": [350, 1500], - "E": [350, 1200], - "B": [350, 1200] + "E": [350, 1500], + "B": [350, 1500] }, "Planck_f143": { "T": [350, 2200], - "E": [350, 2000], + "E": [350, 2100], "B": [350, 2000], }, "Planck_f217": { "T": [350, 2200], - "E": [350, 2000], + "E": [350, 2100], "B": [350, 2000], }