From deb67c1b80c858723440a952b92ce4da86820b30 Mon Sep 17 00:00:00 2001 From: Louis Thibaut Date: Wed, 11 Dec 2024 07:25:53 +0100 Subject: [PATCH] write component to disk and add a plot of TT vs different components --- .../python/get_best_fit_mflike.py | 7 +- .../results_multifrequency_components_TT.py | 173 ++++++++++++++++++ 2 files changed, 179 insertions(+), 1 deletion(-) create mode 100644 project/data_analysis/python/paper_plots/results_multifrequency_components_TT.py diff --git a/project/data_analysis/python/get_best_fit_mflike.py b/project/data_analysis/python/get_best_fit_mflike.py index 0195f32e..7bc19aba 100644 --- a/project/data_analysis/python/get_best_fit_mflike.py +++ b/project/data_analysis/python/get_best_fit_mflike.py @@ -28,10 +28,14 @@ bestfit_dir = f"best_fits{tag}" plot_dir = f"plots/best_fits{tag}/" -log.info(f"save best fits in {bestfit_dir} folder") +components_dir = f"{bestfit_dir}/components" + + +log.info(f"save best fits in {bestfit_dir} folder") pspy_utils.create_directory(bestfit_dir) +pspy_utils.create_directory(components_dir) pspy_utils.create_directory(plot_dir) cosmo_params = d["cosmo_params"] @@ -147,6 +151,7 @@ for comp in fg_components[mode]: ax.plot(l_th, fg_dict[mode, comp, name1, name2]) + np.savetxt(f"{components_dir}/{mode}_{comp}_{cross}.dat", np.transpose([l_th, fg_dict[mode, comp, name1, name2]])) ax.plot(l_th, fg_dict[mode, "all", name1, name2], color="k") ax.plot(l_th, ps_dict[mode.upper()], color="gray") ax.set_title(cross) diff --git a/project/data_analysis/python/paper_plots/results_multifrequency_components_TT.py b/project/data_analysis/python/paper_plots/results_multifrequency_components_TT.py new file mode 100644 index 00000000..bc5fe6b8 --- /dev/null +++ b/project/data_analysis/python/paper_plots/results_multifrequency_components_TT.py @@ -0,0 +1,173 @@ +""" +Plot the TT power spectra with their respective foreground components +""" +import matplotlib +import sys + +import numpy as np +import pylab as plt +from pspipe_utils import best_fits, log, pspipe_list, beam_chromaticity, covariance +from pspy import pspy_utils, so_dict, so_spectra +from matplotlib import rcParams + + +rcParams["font.family"] = "serif" +rcParams["font.size"] = "12" +rcParams["xtick.labelsize"] = 24 +rcParams["ytick.labelsize"] = 24 +rcParams["axes.labelsize"] = 20 +rcParams["axes.titlesize"] = 12 + + +d = so_dict.so_dict() +d.read_from_file(sys.argv[1]) +log = log.get_logger(**d) + + +tag = d["best_fit_tag"] +fg_components = d["fg_components"] +binning_file = d["binning_file"] +lmax = d["lmax"] +spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"] + +bestfit_dir = f"best_fits{tag}" +components_dir = f"{bestfit_dir}/components" +cov_dir = "covariances" +spec_dir = f"spectra_leak_corr_ab_corr_cal{tag}" +result_dir = "plots/TT_components" + +pspy_utils.create_directory(result_dir) + + +bin_lo, bin_hi, lb, bin_size = pspy_utils.read_binning_file(binning_file, lmax) + + + +spec_name_list = pspipe_list.get_spec_name_list(d, delimiter = "_") + + + +## Read the data and extract TT +cov_xar = np.load(f"{cov_dir}/x_ar_final_cov_data.npy") +vec_xar = covariance.read_x_ar_spectra_vec(spec_dir, + spec_name_list, + "cross", + spectra_order = spectra, + type="Dl") + +spectra_cuts = {"dr6_pa4_f220": dict(T=[975, lmax], P=[lmax, lmax]), + "dr6_pa5_f150": dict(T=[775, lmax], P=[775, lmax]), + "dr6_pa6_f150": dict(T=[575, lmax], P=[575, lmax]), + "dr6_pa5_f090": dict(T=[975, lmax], P=[975, lmax]), + "dr6_pa6_f090": dict(T=[975, lmax], P=[975, lmax])} + +bin_out_dict, indices = covariance.get_indices(bin_lo, + bin_hi, + lb, + spec_name_list, + spectra_cuts=spectra_cuts, + spectra_order=spectra, + selected_spectra=["TT"]) + +cov_TT = cov_xar[np.ix_(indices, indices)] +vec_TT = vec_xar[indices] + +lb_dict, ps_TT, sigma_TT = {}, {}, {} +for spec_select in bin_out_dict.keys(): + my_spec = spec_select[0] + my_id, lb_dict[my_spec] = bin_out_dict[spec_select] + ps_TT[my_spec] = vec_TT[my_id] + sigma_TT[my_spec] = np.sqrt(cov_TT[np.ix_(my_id, my_id)].diagonal()) + +# Done + + + + +#### remove pa6 for clarity +d["arrays_dr6"] = ["pa4_f220", "pa5_f090", "pa5_f150"] +#spectra_list = pspipe_list.get_spec_name_list(d, delimiter = "_") + +narrays, _, _ = pspipe_list.get_arrays_list(d) + +l_th, ps_dict = so_spectra.read_ps(f"{bestfit_dir}/cmb.dat", spectra=spectra) + + +fg_components["tt"].remove("tSZ_and_CIB") +for comp in ["tSZ", "cibc", "tSZxCIB"]: + fg_components["tt"].append(comp) + + + +spectra_list = ["dr6_pa5_f090xdr6_pa5_f090", "dr6_pa5_f090xdr6_pa5_f150", "dr6_pa4_f220xdr6_pa5_f090", "dr6_pa5_f150xdr6_pa5_f150", "dr6_pa4_f220xdr6_pa5_f150", "dr6_pa4_f220xdr6_pa4_f220"] + +comp_color = ["darkorange", "green", "red", "blue", "brown", "darkmagenta", "magenta"] +fig, axes = plt.subplots(narrays, narrays, sharex=True, sharey=True, figsize=(24, 18)) +axes = np.atleast_2d(axes) +indices = np.triu_indices(narrays)[::-1] +for i, cross in enumerate(spectra_list): + name1, name2 = cross.split("x") + idx = (indices[0][i], indices[1][i]) + ax = axes[idx] + + if i == 0: + ax.errorbar(lb_dict[cross], ps_TT[cross], sigma_TT[cross], fmt=".", label="data", markersize=12) + else: + ax.errorbar(lb_dict[cross], ps_TT[cross], sigma_TT[cross], fmt=".", markersize=12) + + l_th, fg_all = so_spectra.read_ps(f"{bestfit_dir}/fg_{cross}.dat", spectra=spectra) + + if i==0: + ax.plot(l_th, ps_dict["TT"], color="gray", linestyle="--", label="CMB", linewidth=2) + ax.plot(l_th, ps_dict["TT"] + fg_all["TT"], color="gray", label="CMB + fg", linewidth=2) + else: + ax.plot(l_th, ps_dict["TT"], color="gray", linestyle="--", linewidth=2) + ax.plot(l_th, ps_dict["TT"] + fg_all["TT"], color="gray", linewidth=2) + + + for comp, col in zip(fg_components["tt"], comp_color): + l_th, fg_comp = np.loadtxt(f"{components_dir}/tt_{comp}_{cross}.dat", unpack=True) + if comp == "tSZxCIB": + fg_comp = np.abs(fg_comp) + + + if i==0: + if comp == "cibp": + label = "CIB-Poisson" + elif comp == "cibc": + label = "CIB-Clustered" + elif comp == "tSZxCIB": + label = "|tSZxCIB|" + else: + label=comp + + ax.plot(l_th, fg_comp, label=label, linewidth=2, color=col) + else: + ax.plot(l_th, fg_comp, linewidth=2, color=col) + + + title_ax = cross.replace("dr6_", "") + title_ax = title_ax.replace("_", "-") + if "pa4-f220" in title_ax: + a, b = title_ax.split("x") + title_ax = f"{b}x{a}" + + ax.set_title(title_ax, fontsize=30) + ax.set_yscale("log") + ax.set_ylim(1, 1e4) + ax.set_xlim(200, 7800) + + +for idx in zip(*np.triu_indices(narrays, k=1)): + ax = axes[idx] + fig.delaxes(ax) + +for i in range(narrays): + axes[-1, i].set_xlabel(r"$\ell$", fontsize=35) + axes[i, 0].set_ylabel(r"$D_\ell$", fontsize=35) +fig.legend(bbox_to_anchor=(0.94,1), fontsize=30) +plt.tight_layout() +plt.savefig(f"{result_dir}/TT_per_components.png") +plt.clf() +plt.close() +