adding codes for planck source subtraction

Serena Giardiello · Serena Giardiello · commit 6ad2234caf15 · 2025-11-19T16:38:24.000Z
diff --git a/project/SO/pISO/python/dust/global_dust.dict b/project/SO/pISO/python/dust/global_dust.dict
@@ -20,6 +20,8 @@ sq_win_alms_dir = "/scratch/c.spxsg6/iso/dust/sq_win_alms"
 cov_dir = "/scratch/c.spxsg6/iso/dust/covariances"
 montecarlo_dir = "/scratch/c.spxsg6/iso/dust/montecarlo"
 noise_alms_dir = "/scratch/c.spxsg6/iso/dust/noise_alms"
+catalog_dir = "/scratch/c.spxsg6/iso/dust/catalogs"
+
 
 deconvolve_pixwin = True
 pixwin_Planck = {"pix": 'HEALPIX', "nside": 2048}
@@ -57,13 +59,15 @@ passband_dir_planck = release_dir + "pspipe/for_planck/passbands/"
 freq_info_Planck_f143 = {"freq_tag": 143, "passband": passband_dir_planck + "passband_npipe_f143.dat"}
 freq_info_Planck_f353 = {"freq_tag": 353, "passband": passband_dir_planck + "passband_npipe_f353.dat"}
 
-planck_fits_beam_path = release_dir + "pspipe/for_planck/beams/"
+planck_fits_beam_path = release_dir + "pspipe/for_planck/beams"
+
+beam_dir = "/scratch/c.spxsg6/iso/dust/beams"
 
-beam_T_Planck_f143 = "beams/npipe_DR6_AxB/bl_T_npipe_DR6_143Ax143B.dat" # Get from running pipeline_dust.yml
-beam_T_Planck_f353 = "beams/npipe_DR6_AxB/bl_T_npipe_DR6_353Ax353B.dat"
+beam_T_Planck_f143 = beam_dir + '/npipe_DR6_AxB/bl_T_npipe_DR6_AxB_143Ax143B.dat'  # Get from running pipeline_dust.yml
+beam_T_Planck_f353 = beam_dir + '/npipe_DR6_AxB/bl_T_npipe_DR6_AxB_353Ax353B.dat'
+beam_pol_Planck_f143 = beam_dir + '/npipe_DR6_AxB/bl_pol_npipe_DR6_AxB_143Ax143B.dat'
+beam_pol_Planck_f353 = beam_dir + '/npipe_DR6_AxB/bl_pol_npipe_DR6_AxB_353Ax353B.dat'
 
-beam_pol_Planck_f143 = "beams/npipe_DR6_AxB/bl_pol_npipe_DR6_143Ax143B.dat"
-beam_pol_Planck_f353 = "beams/npipe_DR6_AxB/bl_pol_npipe_DR6_353Ax353B.dat"
 
 include_beam_chromaticity_effect_in_best_fit = False
 
diff --git a/project/SO/pISO/python/planck/check_src_subtraction.py b/project/SO/pISO/python/planck/check_src_subtraction.py
@@ -0,0 +1,133 @@
+"""
+Iterates over the sources in the MF ACT catalog
+and plot a submap centered on individual sources
+of the residual map (i.e. the source subtracted map)
+"""
+from pspy import so_map, pspy_utils, so_map, so_mpi, so_dict
+import numpy as np
+import pandas as pd
+import os, sys
+import pspipe
+from pspipe_utils import log
+
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+map_dir = "planck_projected"
+
+planck_version = d["planck_version"]
+
+out_dir = d["plots_base_dir"] + "/source_sub_check"
+pspy_utils.create_directory(out_dir)
+
+# Read catalog
+cat_file = d["source_catalog"]
+cat = pd.read_table(cat_file, escapechar="#", sep="\s+")
+cat = cat.shift(1, axis=1)
+
+# Sort by flux at 150 GHz
+cat.sort_values(by="flux_T2", inplace=True, ascending=False)
+
+# Cuts
+snr_cut = 5
+flux_cut = 150 #mJy
+n_src_max = 100
+
+mask = (cat.snr_T >= snr_cut) & (cat.flux_T2 >= flux_cut)
+cat = cat[mask]
+
+ra_list = [ra for i, ra in enumerate(cat.ra) if i < n_src_max]
+dec_list = [dec for i, dec in enumerate(cat.dec) if i < n_src_max]
+
+n = len(ra_list)
+log.info(f"Display {n} sources with 150 GHz flux >= {flux_cut:.1f} mJy and SNR >= {snr_cut:.1f}")
+
+if planck_version == "legacy":
+    map_root = "HFI_SkyMap_2048_R3.01_halfmission-{}_f{}_map.fits"
+    map_splits = [1, 2]
+
+if planck_version == "npipe":
+    map_root = "npipe6v20{}_f{}_map.fits"
+    map_splits = ['A', 'B']
+
+map_freqs = [100, 143, 217]
+
+for freq in map_freqs:
+    for split in map_splits:
+
+        map_file = map_root.format(split, freq)
+        map_path = f"{map_dir}/{map_file}"
+
+        m = so_map.read_map(map_path)
+        m_srcfree = so_map.read_map(map_path.replace("map.fits", "map_srcfree.fits"))
+
+        so_mpi.init(True)
+        subtasks = so_mpi.taskrange(imin=0, imax=n-1)
+
+        for task in subtasks:
+
+            log.info(f"Task {task} of {n-1}")
+
+            ra, dec = ra_list[task], dec_list[task]
+
+            ra0, ra1 = ra - 1., ra + 1.
+            dec0, dec1 = dec - 1., dec + 1.
+
+            box = so_map.get_box(ra0, ra1, dec0, dec1)
+            sub = so_map.get_submap_car(m, box)
+            sub_srcfree = so_map.get_submap_car(m_srcfree, box)
+
+            sub.plot(file_name=f"{out_dir}/{planck_version}_f{freq}_{split}_{task:03d}",
+                     color_range=[250, 100, 100],
+                     ticks_spacing_car=0.6
+            )
+            sub_srcfree.plot(file_name=f"{out_dir}/{planck_version}_srcfree_f{freq}_{split}_{task:03d}",
+                     color_range=[250, 100, 100],
+                     ticks_spacing_car=0.6
+            )
+
+
+# Generate HTML file for visualization
+ids_src = np.arange(n)
+multistep_path = os.path.join(os.path.dirname(pspipe.__file__), "js")
+os.system(f"cp {multistep_path}/multistep2.js {out_dir}/")
+
+maps = [f"{freq}_{split}" for freq in map_freqs for split in map_splits]
+
+filename = f"{out_dir}/{planck_version}_sources.html"
+g = open(filename, mode='w')
+g.write('<html>\n')
+g.write('<head>\n')
+g.write(f'<title> Source subtraction for Planck {planck_version} </title>\n')
+g.write(f'<script src="multistep2.js"></script>\n')
+g.write('<script> add_step("src_id", ["c","v"]) </script> \n')
+g.write('<script> add_step("srcfree", ["j","k"]) </script> \n')
+g.write('<script> add_step("map", ["a","z"]) </script> \n')
+g.write('</head> \n')
+g.write('<body> \n')
+g.write(f'<h1> Source subtraction for Planck {planck_version} </h1>')
+g.write('<p> This webpage display the 100 brightest sources at 150 GHz, with a flux higher than 150mJy detected with a SNR > 5. </p> \n')
+g.write('<p> You can switch between sources (c/v), between the map and source free map (j/k) and between the maps (a/z). </p> \n')
+g.write('<div class=src_id> \n')
+for src_id in ids_src:
+
+    g.write('<div class=srcfree>\n')
+    for type in ["", "_srcfree"]:
+
+        prefix = f"{planck_version}{type}"
+        g.write('<div class=map>\n')
+        for map_name in maps:
+
+            freq, split = map_name.split("_")
+            file_name = f"{prefix}_f{freq}_{split}_{src_id:03d}_T.png"
+            g.write(f'<h2> {freq} GHz split {split} - {prefix} [Source no {src_id:03d}] </p> \n')
+            g.write('<img src="' + file_name + '" width="40%" /> \n')
+        g.write('</div>\n')
+
+    g.write('</div>\n')
+g.write('</div> \n')
+g.write('</body> \n')
+g.write('</html> \n')
+g.close()
diff --git a/project/SO/pISO/python/planck/extract_planck_beam.py b/project/SO/pISO/python/planck/extract_planck_beam.py
@@ -0,0 +1,144 @@
+"""
+Read the Planck beam files and save the beam as tables in .dat files.
+both the intensity, the polarisation and the leakage beam files
+"""
+
+import sys
+
+import numpy as np
+import pylab as plt
+import healpy as hp
+from astropy.io import fits
+from pspy import pspy_utils, so_dict
+from pspipe_utils import log
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+planck_fits_beam_path = d["planck_fits_beam_path"]
+beam_dir_d = d["beam_dir"]
+
+freqs = [100, 143, 217, 353]
+lmax = 3030
+lmax_for_plot = 2000
+releases = ["legacy", "npipe_DR6_AxB"]
+
+for release in releases:
+    
+    beam_dir = beam_dir_d + f"/{release}"
+    pspy_utils.create_directory(beam_dir)
+
+    for freq in freqs:
+        
+        if release == "npipe_DR6_AxB":
+            s1, s2 = "A", "B"
+            if freq != 353:
+                Wl = fits.open(f"{planck_fits_beam_path}/QP_dr6_pa6_f150/Wl_npipe6v20_{freq}{s1}x{freq}{s2}.fits")
+            else:
+                print(f"use full sky beam for {freq} GHz")
+                Wl = fits.open(f"{planck_fits_beam_path}/quickpol/Wl_npipe6v20_{freq}{s1}x{freq}{s2}.fits")
+
+        if release == "legacy":
+            s1, s2 = "hm1", "hm2"
+            if freq != 353:
+                Wl = fits.open(f"{planck_fits_beam_path}/BeamWf_HFI_R3.01/Wl_R3.01_plikmask_{freq}{s1}x{freq}{s2}.fits")
+            else:
+                print(f"use full sky beam for {freq} GHz")
+                Wl = fits.open(f"{planck_fits_beam_path}/BeamWf_HFI_R3.01/Wl_R3.01_fullsky_{freq}{s1}x{freq}{s2}.fits")
+
+        l = np.arange(lmax)
+        
+        Wl_TT_2_TT = Wl[1].data["TT_2_TT"][0, :lmax]
+        Wl_EE_2_EE = Wl[2].data["EE_2_EE"][0, :lmax]
+            
+        # extract beam and polarised beam
+        bl_T = np.sqrt(Wl_TT_2_TT)
+        bl_pol = np.sqrt(Wl_EE_2_EE)
+        
+        plt.figure(figsize=(12,8))
+        plt.subplot(2, 1, 1)
+        plt.plot(l[:lmax_for_plot], bl_T[:lmax_for_plot], label="temperature beam", color="lightblue")
+        plt.errorbar(l[:lmax_for_plot], bl_pol[:lmax_for_plot],  fmt="+", markevery=50, label="pol beam", color="red")
+        plt.ylabel(r"$ B_{\ell}$", fontsize=14)
+        plt.legend()
+        plt.subplot(2, 1, 2)
+        plt.plot(l[:lmax_for_plot], (bl_T[:lmax_for_plot] / bl_pol[:lmax_for_plot]) ** 2)
+        plt.ylabel(r"$ (B^{\rm T}_{\ell}/B^{\rm pol}_{\ell})^{2} $", fontsize=14)
+        plt.xlabel(r"$\ell$", fontsize=14)
+        plt.savefig(f"{beam_dir}/beam_{freq}.png")
+        plt.clf()
+        plt.close()
+        
+        np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_T]))
+        np.savetxt(f"{beam_dir}/bl_pol_{release}_{freq}{s1}x{freq}{s2}.dat", np.transpose([l, bl_pol]))
+
+        # extract leakage beam
+        Wl_TE_2_TE = Wl[4].data["TE_2_TE"][0, :lmax]
+        gamma_TE = Wl[1].data["TT_2_TE"][0, :lmax] / Wl_TE_2_TE
+        gamma_ET = Wl[1].data["TT_2_ET"][0, :lmax] / Wl_TE_2_TE
+
+        gamma_TB = Wl[1].data["TT_2_TB"][0, :lmax] / Wl_TE_2_TE
+        gamma_BT = Wl[1].data["TT_2_BT"][0, :lmax] / Wl_TE_2_TE
+
+        plt.figure(figsize=(12,8))
+        plt.title(f"{freq} GHz x {freq} GHz", fontsize=14)
+        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)
+        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)
+        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)
+
+        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)
+        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)
+        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)
+
+        plt.ylim(-0.8, 0.8)
+        plt.ylabel(r"$ \gamma_{\ell}$", fontsize=14)
+        plt.legend()
+        plt.xlabel(r"$\ell$", fontsize=14)
+        plt.savefig(f"{beam_dir}/beam_leakage_{freq}.png")
+        plt.clf()
+        plt.close()
+
+        zeros = np.zeros(len(l))
+    
+        np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s2}_t2e.dat", np.transpose([l, gamma_TE, zeros, zeros, zeros]))
+        np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s1}_t2e.dat", np.transpose([l, gamma_ET, zeros, zeros, zeros]))
+
+        np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s2}_t2b.dat", np.transpose([l, gamma_TB, zeros, zeros, zeros]))
+        np.savetxt(f"{beam_dir}/gamma_{release}_{freq}{s1}_t2b.dat", np.transpose([l, gamma_BT, zeros, zeros, zeros]))
+
+        gamma_mean_TE = (gamma_TE + gamma_ET) / 2
+        gamma_mean_TB = (gamma_TB + gamma_BT) / 2
+
+        np.savetxt(f"{beam_dir}/gamma_mean_{release}_{freq}{s1}{s2}_t2e.dat", np.transpose([l, gamma_mean_TE,  zeros, zeros,  zeros]))
+        np.savetxt(f"{beam_dir}/gamma_mean_{release}_{freq}{s1}{s2}_t2b.dat", np.transpose([l, gamma_mean_TB,  zeros, zeros,  zeros]))
+
+# these beams will be used for source-sub, they have high ell max
+
+release = "npipe"
+beam_dir = beam_dir_d + f"/{release}"
+pspy_utils.create_directory(beam_dir)
+
+
+pl = hp.pixwin(2048)
+
+for freq in freqs:
+
+    Wl_AA = fits.open(f"{planck_fits_beam_path}/quickpol/Wl_npipe6v20_{freq}Ax{freq}A.fits")
+    Wl_BB = fits.open(f"{planck_fits_beam_path}/quickpol/Wl_npipe6v20_{freq}Bx{freq}B.fits")
+
+    bl_T_A = np.sqrt(Wl_AA[1].data["TT_2_TT"][0, :])
+    bl_T_B = np.sqrt(Wl_BB[1].data["TT_2_TT"][0, :])
+    
+    bl_T_coadd = (bl_T_A + bl_T_B) / 2
+
+    bl_T_coadd = bl_T_coadd[:len(pl)]
+    bl_T_coadd_pixwin = bl_T_coadd * pl
+    
+    l = np.arange(len(bl_T_coadd))
+    
+
+    np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}_coadd.dat", np.transpose([l, bl_T_coadd]))
+    np.savetxt(f"{beam_dir}/bl_T_{release}_{freq}_coadd_pixwin.dat", np.transpose([l, bl_T_coadd_pixwin]))
+
+
diff --git a/project/SO/pISO/python/planck/reformat_source_catalog.py b/project/SO/pISO/python/planck/reformat_source_catalog.py
@@ -0,0 +1,72 @@
+"""
+This script is used to reformat the ACT multifrequency
+point source catalog `cat_skn_multifreq_20220526_nightonly`
+in three different monofrequency catalogs (at 90,150 and 220 GHz)
+"""
+from pixell import utils
+import pandas as pd
+import numpy as np
+from pspy import pspy_utils, so_dict
+from pspipe_utils import log
+import sys
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+cat_file = d["source_catalog"]
+out_dir = d["catalog_dir"]
+
+pspy_utils.create_directory(out_dir)
+
+input_catalog = pd.read_table(cat_file, escapechar="#", sep="\s+")
+input_catalog = input_catalog.shift(1, axis=1)
+
+flux_id = {90: 1, 150: 2, 220: 3}
+
+# https://phy-wiki.princeton.edu/polwiki/pmwiki.php?n=BeamAnalysis.BeamDistributionCenter
+beam_areas = {
+    90: [484.17, 487.03],
+    150: [228.44, 215.24, 221.88], #nsr
+    220: [107.34]
+}
+
+for freq, freq_id in flux_id.items():
+
+    flux = {
+        m: getattr(input_catalog, f"flux_{m}{freq_id}") for m in "TQU"
+    }
+    dflux = {
+        m: getattr(input_catalog, f"dflux_{m}{freq_id}") for m in "TQU"
+    }
+
+    # nsr to sr
+    beam_area = np.mean(beam_areas[freq])/1e9
+
+    amp = {}
+    damp = {}
+    for m in ["T", "Q", "U"]:
+        # utils.dplanck converts Jy to Kelvins
+        # flux is given in mJy and we want amps in mK
+        amp[m] = flux[m] / (utils.dplanck(freq*1e9) * beam_area)
+        damp[m] = dflux[m] / (utils.dplanck(freq*1e9) * beam_area)
+
+    ca = input_catalog.ca
+    status = np.where(ca < 3, 1, 2)
+
+    snr = amp["T"] / damp["T"]
+    # npix is not used anywhere
+    npix = np.zeros_like(snr)
+
+    output_catalog = [input_catalog.ra, input_catalog.dec, snr, amp["T"], damp["T"], amp["Q"], damp["Q"], amp["U"], damp["U"],
+                      flux["T"], dflux["T"], flux["Q"], dflux["Q"], flux["U"], dflux["U"],
+                      npix, status]
+
+    output_catalog = np.transpose(output_catalog)
+
+    # ordering by snr
+    output_catalog = output_catalog[output_catalog[:,2].argsort()[::-1]]
+
+    header = "ra dec SNR Tamp dTamp Qamp dQamp Uamp dUamp Tflux dTflux Qflux dQflux Uflux dUflux npix status"
+    fmt = "%11.6f %11.6f %8.3f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f %5.0f %2.0f"
+    np.savetxt(f"{out_dir}/cat_skn_{freq:03d}_20220526_nightonly_ordered.txt", output_catalog, header=header, fmt=fmt)
diff --git a/project/SO/pISO/python/planck/run_npipe_src_subtraction.sh b/project/SO/pISO/python/planck/run_npipe_src_subtraction.sh
