Skip to content

Commit ba9cd63

Browse files
author
Louis Thibaut
committed
add simulation scripts and helper for modifying alm_mask
1 parent 60d26a6 commit ba9cd63

6 files changed

Lines changed: 236 additions & 1 deletion

File tree

‎project/SPT_dev/SPT.rst‎

Lines changed: 22 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,28 @@ Running the main pipeline
2121
The pipeline take roughly 30 minutes, all time is spent doing map2alms since spt has 30 splits per frequency.
2222
An alternative dictfile using half-missions is provided that allows for much faster investigations.
2323

24-
You can then compare with the spt released spectra by running spt_plot_spectra.py, note that you need to install candl and spt_candl_data first.
24+
You can then compare with the spt released spectra by running spt_plot_spectra.py and spt_plot_noise_spectra.py, note that you need to install candl and spt_candl_data first.
25+
Another useful routine to modify spt alm mask (and reduce E-> B leakage) is provided in spt_modify_alm_mask.py.
26+
27+
Simulations
28+
-------------------------------------------------------
29+
30+
SPT provides 500 unfiltered simulations and 500 filtered simulations, allowing us to estimate the transfer function.
31+
We have written script to do both 2d transfer function analysis
32+
33+
.. code:: shell
34+
35+
salloc --nodes 4 --qos interactive --time 02:00:00 --constraint cpu
36+
OMP_NUM_THREADS=16 srun -n 16 -c 64 --cpu-bind=cores python spt_mc_get_tf2d.py global_spt.dict
37+
38+
and one dimensional transfer function analysis
39+
40+
.. code:: shell
41+
42+
salloc --nodes 4 --qos interactive --time 02:00:00 --constraint cpu
43+
OMP_NUM_THREADS=16 srun -n 16 -c 64 --cpu-bind=cores python spt_mc_get_tf_spectra.py global_spt.dict
44+
45+
2546
2647
.. code:: shell
2748

‎project/SPT_dev/global_spt.dict‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,3 +79,6 @@ window_pol_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
7979
cosmo_params = {"cosmomc_theta":0.0104059, "logA": 3.053, "ombh2": 0.02258, "omch2": 0.1242, "ns": 0.9666, "Alens": 1.0, "tau": 0.0567}
8080

8181

82+
iStart = 0
83+
iStop = 15
84+

‎project/SPT_dev/global_spt_half.dict‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,3 +78,6 @@ window_pol_spt_220 = window_dir + "pixel_mask_apodized_borders_only.fits"
7878
cosmo_params = {"cosmomc_theta":0.0104059, "logA": 3.053, "ombh2": 0.02258, "omch2": 0.1242, "ns": 0.9666, "Alens": 1.0, "tau": 0.0567}
7979

8080

81+
iStart = 0
82+
iStop = 15
83+

‎project/SPT_dev/spt_mc_get_tf2d.py‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
"""
2+
This script uses spt simulations to compute 2d (l,m) transfer functions
3+
"""
4+
5+
from pspy import pspy_utils, so_dict, so_map, sph_tools, so_mpi
6+
from pspipe_utils import log
7+
import numpy as np
8+
import healpy as hp
9+
import sys
10+
import time
11+
12+
d = so_dict.so_dict()
13+
d.read_from_file(sys.argv[1])
14+
log = log.get_logger(**d)
15+
16+
survey = "spt"
17+
lmax = d["lmax"]
18+
niter = d["niter"]
19+
type = d["type"]
20+
release_dir = d["release_dir"]
21+
alm_conv = d[f"alm_conv_{survey}"]
22+
23+
24+
tf_dir = "tf2d"
25+
26+
pspy_utils.create_directory(tf_dir)
27+
28+
spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
29+
spin_pairs = ["spin0xspin0", "spin0xspin2", "spin2xspin0", "spin2xspin2"]
30+
arrays_spt = d["arrays_spt"]
31+
32+
33+
so_mpi.init(True)
34+
subtasks = so_mpi.taskrange(imin=d["iStart"], imax=d["iStop"])
35+
36+
for iii in subtasks:
37+
log.info(f"Simulation n° {iii:05d}/{d['iStop']:05d}")
38+
log.info(f"-------------------------")
39+
40+
master_alms = {}
41+
42+
for ar in arrays_spt:
43+
44+
sim_in = so_map.read_map(f"{release_dir}/simulated_maps/input_maps/input_maps_realization{iii:03d}_{ar}ghz.fits")
45+
sim_out = so_map.read_map(f"{release_dir}/simulated_maps/output_maps/masking_yes/output_maps_masking_yes_realization{iii:03d}_{ar}ghz.fits")
46+
47+
win_T = so_map.read_map(d[f"window_T_{survey}_{ar}"])
48+
49+
fsky = np.sum(win_T.data) / (12 * win_T.nside ** 2)
50+
51+
win_pol = so_map.read_map(d[f"window_pol_{survey}_{ar}"])
52+
53+
window_tuple = (win_T, win_pol)
54+
55+
master_alms["nofilter"] = sph_tools.get_alms(sim_in, window_tuple, niter, lmax, alm_conv=alm_conv)
56+
master_alms["filter"] = sph_tools.get_alms(sim_out, window_tuple, niter, lmax, alm_conv=alm_conv)
57+
58+
l = np.arange(lmax)
59+
for i, comp in enumerate(["TT", "EE", "BB"]):
60+
for filt in ["nofilter", "filter"]:
61+
62+
master_alms[filt][i] = hp.sphtfunc.almxfl(master_alms[filt][i], l)
63+
ps_2d = np.abs(master_alms[filt][i] * np.conjugate(master_alms[filt][i]))/(2 * np.pi * fsky)
64+
65+
np.save(f"{tf_dir}/tf2d_{comp}_{ar}_{filt}_{iii:05d}.npy", ps_2d)
66+

‎project/SPT_dev/spt_mc_get_tf_spectra.py‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
"""
2+
This script compute the one dimensional power spectra of spt simulations
3+
"""
4+
5+
from pspy import pspy_utils, so_dict, so_map, sph_tools, so_mcm, so_spectra, so_mpi
6+
from pspipe_utils import pspipe_list, log
7+
import numpy as np
8+
import healpy as hp
9+
import sys
10+
import time
11+
12+
d = so_dict.so_dict()
13+
d.read_from_file(sys.argv[1])
14+
log = log.get_logger(**d)
15+
16+
survey = "spt"
17+
lmax = d["lmax"]
18+
niter = d["niter"]
19+
type = d["type"]
20+
binning_file = d["binning_file"]
21+
binned_mcm = d["binned_mcm"]
22+
release_dir = d["release_dir"]
23+
alm_conv = d[f"alm_conv_{survey}"]
24+
25+
26+
mcm_dir = "mcms"
27+
tf_dir = "sim_spectra_for_tf"
28+
29+
pspy_utils.create_directory(tf_dir)
30+
31+
spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
32+
spin_pairs = ["spin0xspin0", "spin0xspin2", "spin2xspin0", "spin2xspin2"]
33+
arrays_spt = d["arrays_spt"]
34+
35+
36+
so_mpi.init(True)
37+
subtasks = so_mpi.taskrange(imin=d["iStart"], imax=d["iStop"])
38+
39+
for iii in subtasks:
40+
log.info(f"Simulation n° {iii:05d}/{d['iStop']:05d}")
41+
log.info(f"-------------------------")
42+
43+
t0 = time.time()
44+
45+
master_alms = {}
46+
47+
for ar in arrays_spt:
48+
49+
sim_in = so_map.read_map(f"{release_dir}/simulated_maps/input_maps/input_maps_realization{iii:03d}_{ar}ghz.fits")
50+
sim_out = so_map.read_map(f"{release_dir}/simulated_maps/output_maps/masking_yes/output_maps_masking_yes_realization{iii:03d}_{ar}ghz.fits")
51+
52+
win_T = so_map.read_map(d[f"window_T_{survey}_{ar}"])
53+
win_pol = so_map.read_map(d[f"window_pol_{survey}_{ar}"])
54+
55+
window_tuple = (win_T, win_pol)
56+
57+
master_alms[survey, ar, "nofilter"] = sph_tools.get_alms(sim_in, window_tuple, niter, lmax, alm_conv=alm_conv)
58+
master_alms[survey, ar, "filter"] = sph_tools.get_alms(sim_out, window_tuple, niter, lmax, alm_conv=alm_conv)
59+
60+
alm_mask = hp.read_alm(d[f"alm_mask_{survey}_{ar}"], hdu=1)
61+
alm_mask = hp.sphtfunc.resize_alm(alm_mask, d["lmax_mask"], d["lmax_mask"], lmax, lmax)
62+
63+
master_alms[survey, ar, "nofilter_mask"] = master_alms[survey, ar, "nofilter"] * alm_mask
64+
master_alms[survey, ar, "filter_mask"] = master_alms[survey, ar, "filter"] * alm_mask
65+
66+
67+
_, _, lb, _ = pspy_utils.read_binning_file(binning_file, lmax)
68+
69+
n_spec, sv1_list, ar1_list, sv2_list, ar2_list = pspipe_list.get_spectra_list(d)
70+
71+
for i_spec in range(n_spec):
72+
sv1, ar1, sv2, ar2 = sv1_list[i_spec], ar1_list[i_spec], sv2_list[i_spec], ar2_list[i_spec]
73+
spec_name = f"{type}_{sv1}_{ar1}x{sv2}_{ar2}"
74+
75+
mbb_inv, Bbl = so_mcm.read_coupling(prefix=f"{mcm_dir}/{sv1}_{ar1}x{sv2}_{ar2}", spin_pairs=spin_pairs)
76+
77+
for filt in ["filter", "filter_mask", "nofilter", "nofilter_mask"]:
78+
79+
l, ps_master = so_spectra.get_spectra_pixell(master_alms[sv1, ar1, filt], master_alms[sv2, ar2, filt], spectra=spectra)
80+
81+
lb, ps = so_spectra.bin_spectra(l,
82+
ps_master,
83+
binning_file,
84+
lmax,
85+
type=type,
86+
mbb_inv=mbb_inv,
87+
spectra=spectra,
88+
binned_mcm=binned_mcm)
89+
90+
so_spectra.write_ps(tf_dir + f"/{spec_name}_{filt}_{iii:05d}.dat", lb, ps, type, spectra=spectra)
91+
92+
log.info(f"[{iii}] Simulation n° {iii:05d} done in {time.time()-t0:.02f} s")

‎project/SPT_dev/spt_modify_alm_mask.py‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
"""
2+
script to modify spt3g alm mask
3+
"""
4+
import sys
5+
6+
import numpy as np
7+
import healpy as hp
8+
from pspipe_utils import log
9+
from pspy import pspy_utils, so_dict, sph_tools
10+
11+
12+
13+
def create_mask_lm(l, m, l_cutoff, m_cutoff, width=20):
14+
m_mask = np.where(m <= (m_cutoff - width), 1.0,
15+
np.where(m >= m_cutoff, 0.0,
16+
0.5 * (1 + np.cos(np.pi * (m - (m_cutoff - width)) / width))))
17+
l_mask = np.where(l <= l_cutoff, 0.0,
18+
np.where(l >= (l_cutoff + width), 1.0,
19+
0.5 * (1 - np.cos(np.pi * (l - l_cutoff) / width))))
20+
21+
return l_mask * m_mask
22+
23+
24+
25+
d = so_dict.so_dict()
26+
d.read_from_file(sys.argv[1])
27+
log = log.get_logger(**d)
28+
29+
30+
lmax_mask = d["lmax_mask"]
31+
surveys = d["surveys"]
32+
release_dir = d["release_dir"]
33+
plot_dir = "plots"
34+
pspy_utils.create_directory(plot_dir)
35+
36+
37+
for sv in surveys:
38+
arrays = d[f"arrays_{sv}"]
39+
for ar in arrays:
40+
alm_mask = hp.read_alm(release_dir + f"ancillary_products/specific_to_c25/alm_mask_{ar}ghz.fits", hdu=(1,2,3))
41+
sph_tools.show_alm_triangle(alm_mask, lmax_mask, 0, 1, xlims=[0,2000], ylims=[0,2000], fig_file=f"{plot_dir}/old_alm_mask")
42+
43+
l, m = hp.Alm.getlm(lmax_mask)
44+
45+
soft_mask = create_mask_lm(l, m, l_cutoff=400, m_cutoff=270, width=50)
46+
alm_mask *= (1.0 - soft_mask)
47+
48+
sph_tools.show_alm_triangle(alm_mask, lmax_mask, 0, 1, xlims=[0,2000], ylims=[0,2000], fig_file=f"{plot_dir}/new_alm_mask")
49+
50+
hp.fitsfunc.write_alm(f"alm_mask_{ar}ghz_modif.fits", alm_mask, overwrite=True)

0 commit comments

Comments
 (0)