script for combined null test

project/data_analysis/python/montecarlo/mc_analyze_combined_spectra.py

@@ -13,17 +13,25 @@
 d.read_from_file(sys.argv[1])
 log = log.get_logger(**d)
 
+include_syst = True
+if include_syst:
+    print(f"Use sim spectra with beam and leakage uncertainties")
+    add_str = "_syst"
+else:
+    add_str = ""
+
 bestfit_dir = f"best_fits"
-combined_spec_dir = "combined_sim_spectra"
-plot_dir = f"plots/combined_sim_spectra/"
+combined_spec_dir = f"combined_sim_spectra{add_str}"
+plot_dir = f"plots/combined_sim_spectra{add_str}/"
 
 pspy_utils.create_directory(plot_dir)
 
 binning_file = d["binning_file"]
 lmax = d["lmax"]
 type = d["type"]
 iStart = d["iStart"]
-iStop = d["iStop"]
+iStop =  d["iStop"]
+
 name = "all"
 
 spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
@@ -74,7 +82,6 @@
     plt.clf()
     plt.close()
 
-
     Dl_mean = np.mean(Dl_list, axis=0)
     Dl_fg_sub_mean = np.mean(Dl_list_fg_sub, axis=0)
 
@@ -105,7 +112,6 @@
     plt.clf()
     plt.close()
 
-
     plt.figure(figsize=(12,8))
     plt.subplot(2,1,1)
     plt.semilogy()

project/data_analysis/python/montecarlo/mc_get_combined_spectra.py

@@ -72,7 +72,6 @@ def combine_and_save_spectra(lb_ml, P_mat, cov_ml, i_sub_cov, sub_vec,  sub_vec_
                                                       P_mat,
                                                       sub_vec_fg_sub)
 
-    np.savetxt(f"{combined_spec_dir}/{type}_{name}.dat", np.transpose([lb_ml, vec_ml, np.sqrt(cov_ml.diagonal())]))
     np.savetxt(f"{combined_spec_dir}/{type}_{name}_cmb_only.dat", np.transpose([lb_ml, vec_ml_fg_sub, np.sqrt(cov_ml.diagonal())]))
 
 def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dict):
@@ -103,10 +102,17 @@ def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dic
 d.read_from_file(sys.argv[1])
 log = log.get_logger(**d)
 
-sim_spec_dir = f"sim_spectra"
+include_syst = True
+if include_syst:
+    print(f"Use sim spectra with beam and leakage uncertainties")
+    add_str = "_syst"
+else:
+    add_str = ""
+
+sim_spec_dir = f"sim_spectra{add_str}"
 bestfit_dir = f"best_fits"
 mcm_dir = "mcms"
-combined_spec_dir = f"combined_sim_spectra"
+combined_spec_dir = f"combined_sim_spectra{add_str}"
 plot_dir = f"plots/combined_sim_spectra/"
 
 pspy_utils.create_directory(combined_spec_dir)
@@ -124,6 +130,10 @@ def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dic
 
 cov = np.load("covariances/x_ar_final_cov_sim_gp.npy")
 
+if include_syst:
+    cov_beam = np.load("covariances/x_ar_beam_cov.npy")
+    cov_leakage = np.load("covariances/x_ar_leakage_cov.npy")
+    cov += cov_beam + cov_leakage
 
 vec_xar_th = covariance.read_x_ar_theory_vec(bestfit_dir,
                                              mcm_dir,
@@ -138,9 +148,6 @@ def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dic
                                                 spectra_order=spectra,
                                                 foreground_only=True)
 
-
-
-
 ########################################################################################
 spectra_cuts = {
     "dr6_pa4_f220": dict(T=[975, lmax], P=[lmax, lmax]),
@@ -243,9 +250,8 @@ def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dic
                                                        type=type)
 
             sub_vec = vec_xar[all_indices]
-            combine_and_save_spectra(lb_ml, P_mat, cov_ml, i_sub_cov, sub_vec,  sub_vec_fg_th, name)
+            combine_and_save_spectra(lb_ml, P_mat, cov_ml, i_sub_cov, sub_vec,  sub_vec_fg_th, sim_name)
 
-
 #### Now do the cross-freq spectra combination
 
 cross_freq_pairs = ["90x220", "90x150", "150x220"]
@@ -284,15 +290,13 @@ def ml_helper(cov, vec_xar_th, vec_xar_fg_th, bin_mean, all_indices, bin_out_dic
         for iii in range(iStart, iStop + 1):
             sim_name = name + f"_{iii:05d}"
             print(sim_name)
-
             vec_xar = covariance.read_x_ar_spectra_vec(sim_spec_dir,
                                                        spec_name_list,
                                                        f"cross_{iii:05d}",
                                                        spectra_order=spectra,
                                                        type=type)
 
-
             sub_vec = vec_xar[all_indices]
-            combine_and_save_spectra(lb_ml, P_mat, cov_ml, i_sub_cov, sub_vec,  sub_vec_fg_th, name)
+            combine_and_save_spectra(lb_ml, P_mat, cov_ml, i_sub_cov, sub_vec,  sub_vec_fg_th, sim_name)
 
 

project/data_analysis/python/paper_plots/results_get_combined_null.py

@@ -0,0 +1,104 @@
+"""
+This script test the bias and covariance of combined spectra
+"""
+
+from pspy import so_dict, pspy_utils, so_spectra, so_cov
+from pspipe_utils import log
+import numpy as np
+import pylab as plt
+import sys, os
+import scipy.stats as ss
+from matplotlib import rcParams
+
+
+
+rcParams["font.family"] = "serif"
+rcParams["font.size"] = "20"
+rcParams["xtick.labelsize"] = 20
+rcParams["ytick.labelsize"] = 20
+rcParams["axes.labelsize"] = 20
+rcParams["axes.titlesize"] = 20
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+tag = d["best_fit_tag"]
+combined_spec_dir = f"combined_spectra{tag}"
+
+combined_sim_spec_dir = f"combined_sim_spectra_syst"
+plot_dir = f"plots/combined_null_test/"
+
+pspy_utils.create_directory(plot_dir)
+
+binning_file = d["binning_file"]
+lmax = d["lmax"]
+type = d["type"]
+iStart = 0
+iStop =  1639
+
+freq_pairs = ["90x90", "90x150", "150x150"]
+
+ylim = {}
+ylim["TT"] = (-10,10)
+ylim["TE"] = (-3,3)
+ylim["TB"] = (-3,3)
+ylim["EE"] = (-2,2)
+ylim["EB"] = (-1,1)
+ylim["BB"] = (-1,1)
+
+#plt.figure(figsize=(18,20))
+
+combined_spectra = ["TT", "TE", "TB", "EE", "EB", "BB"]
+for ispec, spectrum in enumerate(combined_spectra):
+
+    plt.figure(figsize=(18,6))
+
+    count=0
+
+    for fp1 in freq_pairs:
+        for fp2 in freq_pairs:
+            if fp1 <= fp2: continue
+            l1, Dl1_fg_sub, _ = np.loadtxt(f"{combined_spec_dir}/Dl_{fp1}_{spectrum}_cmb_only.dat", unpack=True)
+            l2, Dl2_fg_sub, _ = np.loadtxt(f"{combined_spec_dir}/Dl_{fp2}_{spectrum}_cmb_only.dat", unpack=True)
+
+            min_l_null = np.maximum(l1[0], l2[0])
+            id1 = np.where(l1 >= min_l_null)
+            id2 = np.where(l2 >= min_l_null)
+            l1 = l1[id1]
+            l2 = l2[id2]
+            diff = (Dl1_fg_sub[id1] - Dl2_fg_sub[id2])
+
+            mc_diff_list = []
+            for iii in range(iStart, iStop + 1):
+                l1_, Dl1_fg_sub_sim, _ = np.loadtxt(f"{combined_sim_spec_dir}/Dl_{fp1}_{spectrum}_{iii:05d}_cmb_only.dat", unpack=True)
+                l2_, Dl2_fg_sub_sim, _ = np.loadtxt(f"{combined_sim_spec_dir}/Dl_{fp2}_{spectrum}_{iii:05d}_cmb_only.dat", unpack=True)
+                diff_sim = (Dl1_fg_sub_sim[id1] - Dl2_fg_sub_sim[id2])
+                mc_diff_list  += [diff_sim]
+
+            cov = np.cov(mc_diff_list, rowvar=False)
+            corr = so_cov.cov2corr(cov, remove_diag=True)
+
+            chi2 = diff @ np.linalg.inv(cov) @ diff
+            ndof = len(diff)
+            pte = 1 - ss.chi2(ndof).cdf(chi2)
+
+            std = np.sqrt(cov.diagonal())
+
+            fa1, fb1 = fp1.split("x")
+            fa2, fb2 = fp2.split("x")
+
+            plt.errorbar(l1 - 8 + 8*count, diff, std, fmt="o", label=f"{fa1} GHz x {fb1} GHz - {fa2} GHz x {fb2} GHz, PTE: {pte*100:0.2f} %")
+            count += 1
+
+    lth = np.linspace(0, 8500)
+    plt.plot(lth, lth*0, color="black", linestyle="--", alpha=0.5)
+    plt.ylim(ylim[spectrum])
+    plt.xlim(800, 4000)
+    plt.legend(fontsize=16)
+    plt.xlabel(r"$\ell$", fontsize=30)
+    plt.ylabel(r"$D^{%s}_{\ell}$" % spectrum, fontsize=30)
+
+    plt.savefig(f"{plot_dir}/null_{spectrum}.png", bbox_inches="tight")
+    plt.clf()
+    plt.close()