Skip to content

Commit 89c809c

Browse files
Louis ThibautLouis Thibaut
Louis Thibaut
authored and
Louis Thibaut
committed
use new convention for leakage file
1 parent 5cf8d23 commit 89c809c

1 file changed

Lines changed: 86 additions & 94 deletions

File tree

‎project/ACT_DR6/python/planck/extract_planck_beam.py‎

Lines changed: 86 additions & 94 deletions
Original file line numberDiff line numberDiff line change
@@ -2,118 +2,110 @@
22
Read the Planck beam files and save the beam as tables in .dat files.
33
both the intensity, the polarisation and the leakage beam files
44
"""
5+
6+
import sys
7+
58
import numpy as np
69
import pylab as plt
710
from astropy.io import fits
8-
from pspy import pspy_utils
11+
from pspy import pspy_utils, so_dict
912

10-
freqs = [100, 143, 217]
11-
lmax_for_plot = 2000
12-
beam_path = "/global/cfs/cdirs/act/data/tlouis/dr6v4/beams/planck_beams/"
13-
releases = ["npipe", "legacy", "npipe_DR6"]
1413

14+
d = so_dict.so_dict()
15+
d.read_from_file(sys.argv[1])
16+
log = log.get_logger(**d)
17+
18+
planck_fits_beam_path = d["planck_fits_beam_path"]
19+
20+
freqs = [100, 143, 217, 353]
21+
lmax = 3030
22+
lmax_for_plot = 2000
23+
releases = ["legacy", "npipe_DR6_AxB"]
1524

1625
for release in releases:
1726
beam_dir = f"beams/{release}"
1827
pspy_utils.create_directory(beam_dir)
28+
1929
for freq in freqs:
2030

21-
if release == "npipe":
22-
split_pairs = [("A", "A"), ("A", "B"), ("B", "B")]
23-
file_root = f"{beam_path}/quickpol/Wl_npipe6v20"
24-
25-
if release == "npipe_DR6":
26-
split_pairs = [("A", "A"), ("A", "B"), ("B", "B")]
27-
file_root = f"{beam_path}/QP_dr6_pa6_f150/Wl_npipe6v20"
31+
if release == "npipe_DR6_AxB":
32+
s1, s2 = "A", "B"
33+
Wl = fits.open(f"{planck_fits_beam_path}/QP_dr6_pa6_f150/Wl_npipe6v20_{freq}{s1}x{freq}{s2}.fits")
2834

2935
if release == "legacy":
30-
split_pairs = [("hm1", "hm1"), ("hm1", "hm2"), ("hm2", "hm2")]
31-
file_root = f"{beam_path}/BeamWf_HFI_R3.01/Wl_R3.01_plikmask"
36+
s1, s2 = "hm1", "hm2"
37+
if freq != 353:
38+
Wl = fits.open(f"{planck_fits_beam_path}/BeamWf_HFI_R3.01/Wl_R3.01_plikmask_{freq}{s1}x{freq}{s2}.fits")
39+
else:
40+
print(f"use full sky beam for {freq} GHz")
41+
Wl = fits.open(f"{planck_fits_beam_path}/BeamWf_HFI_R3.01/Wl_R3.01_fullsky_{freq}{s1}x{freq}{s2}.fits")
42+
43+
l = np.arange(lmax)
44+
45+
Wl_TT_2_TT = Wl[1].data["TT_2_TT"][0, :lmax]
46+
Wl_EE_2_EE = Wl[2].data["EE_2_EE"][0, :lmax]
47+
48+
# extract beam and polarised beam
49+
bl_T = np.sqrt(Wl_TT_2_TT)
50+
bl_pol = np.sqrt(Wl_EE_2_EE)
51+
52+
plt.figure(figsize=(12,8))
53+
plt.subplot(2, 1, 1)
54+
plt.plot(l[:lmax_for_plot], bl_T[:lmax_for_plot], label="temperature beam", color="lightblue")
55+
plt.errorbar(l[:lmax_for_plot], bl_pol[:lmax_for_plot], fmt="+", markevery=50, label="pol beam", color="red")
56+
plt.ylabel(r"$ B_{\ell}$", fontsize=14)
57+
plt.legend()
58+
plt.subplot(2, 1, 2)
59+
plt.plot(l[:lmax_for_plot], (bl_T[:lmax_for_plot] / bl_pol[:lmax_for_plot]) ** 2)
60+
plt.ylabel(r"$ (B^{\rm T}_{\ell}/B^{\rm pol}_{\ell})^{2} $", fontsize=14)
61+
plt.xlabel(r"$\ell$", fontsize=14)
62+
plt.savefig(f"{beam_dir}/beam_{freq}.png")
63+
plt.clf()
64+
plt.close()
65+
66+
np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_T]))
67+
np.savetxt(f"{beam_dir}/bl_pol_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_pol]))
68+
69+
# extract leakage beam
70+
Wl_TE_2_TE = Wl[4].data["TE_2_TE"][0, :lmax]
71+
gamma_TE = Wl[1].data["TT_2_TE"][0, :lmax] / Wl_TE_2_TE
72+
gamma_ET = Wl[1].data["TT_2_ET"][0, :lmax] / Wl_TE_2_TE
73+
74+
gamma_TB = Wl[1].data["TT_2_TB"][0, :lmax] / Wl_TE_2_TE
75+
gamma_BT = Wl[1].data["TT_2_BT"][0, :lmax] / Wl_TE_2_TE
76+
77+
plt.figure(figsize=(12,8))
78+
plt.title(f"{freq} GHz x {freq} GHz", fontsize=14)
79+
plt.errorbar(l[:lmax_for_plot], 100 * gamma_TE[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ E_{\rm %s}$" % (release, s1, s2), color="blue", fmt="-", markevery=50)
80+
plt.errorbar(l[:lmax_for_plot], 100 * gamma_ET[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ E_{\rm %s}$" % (release, s2, s1), color="navy", fmt="+", markevery=50)
81+
plt.errorbar(l[:lmax_for_plot], 100 * (gamma_TE[:lmax_for_plot] + gamma_ET[:lmax_for_plot]) / 2, label=r"%s $T \ x \ E$" % (release), color="black", fmt="--", markevery=50)
82+
83+
plt.errorbar(l[:lmax_for_plot], 100 * gamma_TB[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ B_{\rm %s}$" % (release, s1, s2), color="red", fmt="-.", markevery=50)
84+
plt.errorbar(l[:lmax_for_plot], 100 * gamma_BT[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ B_{\rm %s}$" % (release, s2, s1), color="orange", fmt="*", markevery=50)
85+
plt.errorbar(l[:lmax_for_plot], 100 * (gamma_TB[:lmax_for_plot] + gamma_BT[:lmax_for_plot]) / 2, label=r"%s $T \ x \ B$" % (release), color="gray", fmt="--", markevery=50)
86+
87+
plt.ylim(-0.8, 0.8)
88+
plt.ylabel(r"$ \gamma_{\ell}$", fontsize=14)
89+
plt.legend()
90+
plt.xlabel(r"$\ell$", fontsize=14)
91+
plt.savefig(f"{beam_dir}/beam_leakage_{freq}.png")
92+
plt.clf()
93+
plt.close()
94+
95+
zeros = np.zeros(len(l))
3296

33-
bl_dict = {}
34-
for s1, s2 in split_pairs:
97+
np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s2}_t2e.dat", np.transpose([l, gamma_TE, zeros, zeros, zeros]))
98+
np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s1}_t2e.dat", np.transpose([l, gamma_ET, zeros, zeros, zeros]))
3599

36-
Wl = fits.open(f"{file_root}_{freq}{s1}x{freq}{s2}.fits")
37-
Wl_TT_2_TT = Wl[1].data["TT_2_TT"][0, :]
38-
Wl_EE_2_EE = Wl[2].data["EE_2_EE"][0, :]
39-
40-
l = np.arange(len(Wl_TT_2_TT))
100+
np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s2}_t2b.dat", np.transpose([l, gamma_TB, zeros, zeros, zeros]))
101+
np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s1}_t2b.dat", np.transpose([l, gamma_BT, zeros, zeros, zeros]))
41102

42-
# extract beam and polarised beam
43-
bl_T = np.sqrt(Wl_TT_2_TT)
44-
bl_pol = np.sqrt(Wl_EE_2_EE)
45-
46-
if s1 != s2:
47-
plt.figure(figsize=(12,8))
48-
plt.subplot(2, 1, 1)
49-
plt.plot(l[:lmax_for_plot], bl_T[:lmax_for_plot], label="temperature beam", color="lightblue")
50-
plt.errorbar(l[:lmax_for_plot], bl_pol[:lmax_for_plot], fmt="+", markevery=50, label="pol beam", color="red")
51-
plt.ylabel(r"$ B_{\ell}$", fontsize=14)
52-
plt.legend()
53-
plt.subplot(2, 1, 2)
54-
plt.plot(l[:lmax_for_plot], (bl_T[:lmax_for_plot] / bl_pol[:lmax_for_plot]) ** 2)
55-
plt.ylabel(r"$ (B^{\rm T}_{\ell}/B^{\rm pol}_{\ell})^{2} $", fontsize=14)
56-
plt.xlabel(r"$\ell$", fontsize=14)
57-
plt.savefig(f"{beam_dir}/beam_{freq}.png")
58-
plt.clf()
59-
plt.close()
60-
61-
bl_dict[s1, s2, "T"] = bl_T
62-
bl_dict[s1, s2, "pol"] = bl_pol
63-
np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_T]))
64-
np.savetxt(f"{beam_dir}/bl_pol_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_pol]))
65-
66-
# extract leakage beam
67-
Wl_TE_2_TE = Wl[4].data["TE_2_TE"][0]
68-
gamma_TE = Wl[1].data["TT_2_TE"][0] / Wl_TE_2_TE
69-
gamma_ET = Wl[1].data["TT_2_ET"][0] / Wl_TE_2_TE
70-
71-
gamma_TB = Wl[1].data["TT_2_TB"][0] / Wl_TE_2_TE
72-
gamma_BT = Wl[1].data["TT_2_BT"][0] / Wl_TE_2_TE
73-
74-
if s1 != s2:
75-
plt.figure(figsize=(12,8))
76-
plt.title(f"{freq} GHz x {freq} GHz", fontsize=14)
77-
plt.errorbar(l[:lmax_for_plot], 100 * gamma_TE[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ E_{\rm %s}$" % (release, s1, s2), color="blue", fmt="-", markevery=50)
78-
plt.errorbar(l[:lmax_for_plot], 100 * gamma_ET[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ E_{\rm %s}$" % (release, s2, s1), color="navy", fmt="+", markevery=50)
79-
plt.errorbar(l[:lmax_for_plot], 100 * (gamma_TE[:lmax_for_plot] + gamma_ET[:lmax_for_plot]) / 2, label=r"%s $T \ x \ E$" % (release), color="black", fmt="--", markevery=50)
80-
81-
plt.errorbar(l[:lmax_for_plot], 100 * gamma_TB[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ B_{\rm %s}$" % (release, s1, s2), color="red", fmt="-.", markevery=50)
82-
plt.errorbar(l[:lmax_for_plot], 100 * gamma_BT[:lmax_for_plot], label=r"%s $T_{\rm %s} \ x \ B_{\rm %s}$" % (release, s2, s1), color="orange", fmt="*", markevery=50)
83-
plt.errorbar(l[:lmax_for_plot], 100 * (gamma_TB[:lmax_for_plot] + gamma_BT[:lmax_for_plot]) / 2, label=r"%s $T \ x \ B$" % (release), color="gray", fmt="--", markevery=50)
84-
85-
plt.ylim(-0.8, 0.8)
86-
plt.ylabel(r"$ \gamma_{\ell}$", fontsize=14)
87-
plt.legend()
88-
plt.xlabel(r"$\ell$", fontsize=14)
89-
plt.savefig(f"{beam_dir}/beam_leakage_{freq}.png")
90-
plt.clf()
91-
plt.close()
92-
93-
94-
zeros = np.zeros(len(l))
95-
96-
np.savetxt(f"{beam_dir}/gamma_TP_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, gamma_TE, gamma_TB, zeros, zeros]))
97-
np.savetxt(f"{beam_dir}/gamma_PT_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, gamma_ET, gamma_BT, zeros, zeros]))
98-
99-
gamma_mean_TE = (gamma_TE + gamma_ET) / 2
100-
gamma_mean_TB = (gamma_TB + gamma_BT) / 2
103+
gamma_mean_TE = (gamma_TE + gamma_ET) / 2
104+
gamma_mean_TB = (gamma_TB + gamma_BT) / 2
101105

102-
np.savetxt(f"{beam_dir}/gamma_mean_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, gamma_mean_TE, gamma_mean_TB, zeros, zeros]))
106+
np.savetxt(f"{beam_dir}/gamma_mean_{release}_{freq}{s1}{s2}_t2e.dat", np.transpose([l, gamma_mean_TE, zeros, zeros, zeros]))
107+
np.savetxt(f"{beam_dir}/gamma_mean_{release}_{freq}{s1}{s2}_t2b.dat", np.transpose([l, gamma_mean_TB, zeros, zeros, zeros]))
103108

104109

105-
np.savetxt(f"{beam_dir}/error_modes_gamma_TP_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, zeros, zeros, zeros, zeros, zeros, zeros]))
106-
np.savetxt(f"{beam_dir}/error_modes_gamma_PT_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, zeros, zeros, zeros, zeros, zeros, zeros]))
107-
108-
np.savetxt(f"{beam_dir}/error_modes_gamma_mean_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, zeros, zeros, zeros, zeros, zeros, zeros]))
109110

110-
if "npipe" in release:
111-
bl_T_coadd = (bl_dict["A", "A", "T"] + bl_dict["B", "B", "T"]) / 2.
112-
bl_pol_coadd = (bl_dict["A", "A", "pol"] + bl_dict["B", "B", "pol"]) / 2.
113-
114-
elif release == "legacy":
115-
bl_T_coadd = (bl_dict["hm1", "hm1", "T"] + bl_dict["hm2", "hm2", "T"]) / 2.
116-
bl_pol_coadd = (bl_dict["hm1", "hm1", "pol"] + bl_dict["hm2", "hm2", "pol"]) / 2.
117111

118-
np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}_coadd.dat", np.transpose([l, bl_T_coadd]))
119-
np.savetxt(f"{beam_dir}/bl_pol_{release}_{freq}_coadd.dat", np.transpose([l, bl_pol_coadd]))

0 commit comments

Comments
 (0)