import healpy
import spherelib
import spherelib.astro as astro
from pipetools import make_psmask
import planckdata as data
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as pl
[docs]def get_mask(nside, fsky=0, radius=60):
""" Return
- point source mask -> to be applied to I and fastfilled
- apodized galactic mask -> to be applied to T, E and B maps
"""
# ps mask
psmask = make_psmask(get_data_fn("", seg="ps_fit")+"../", nside)
plot_map_in_pipe(psmask, "ps_common_mask", str(nside), norm=None)
fn = get_data_fn("psmask_common_%d.fits"%nside)
healpy.write_map(fn , psmask)
logged_subprocess(["gzip", fn])
# galactic mask
if fsky>0:
gmask = healpy.read_map (data.get_galmask(fsky=fsky, radius=radius))
else:
gmask = healpy.read_map (data.get_ps2galmask(radius=radius))
if nside!=2048:
gmask[gmask!=0] = 1
gmask = healpy.ud_grade (gmask, nside)
gmask[gmask!=1] = 0
spherelib.apodize_mask(gmask, radius)
plot_map_in_pipe(gmask, "gal_common_mask", str(nside), norm=None)
fn = get_data_fn("gmask_apodized_%s.fits"%nside)
healpy.write_map(fn , gmask)
logged_subprocess(["gzip", fn])
return psmask, gmask
[docs]def do_ps_subtraction(mapF, name, dset):
""" Perform point source subtraction as fitted on full data set
"""
cat_fn = glob_parent("cat_fit*.txt")[0]
fit_cat = np.loadtxt(cat_fn)
load_param("ps_fit", globals(), ["latcut", "maxxi2"])
beam = data.get_fwhm_beam(name)
freq = data.get_freq(name)
logger.info("subtracting ps %s : beam %e , latcut %d , maxxi2 %d"%(cat_fn, beam,latcut,maxxi2[freq]))
spherelib.ps_subtract (mapF, fit_cat, beam, nsigma=1, latcut=latcut, maxxi2=maxxi2[freq])
dp = dict({})
load_products(glob_parent("fit_cat_postproc.pkl")[0], dp)
psmask = healpy.read_map(glob_parent("psmask*.fits.gz")[0])
plot_ps_subtraction(mapF*psmask, name, dp['fit_cat'], dp['idx_bad'], dp['idx_good'], beam, dset=dset)
return mapF
[docs]def main():
"""
Apply source subtraction and masking
Compute alm from map
"""
name = get_input("ps_fit")
dset = get_input("dset")
freq = data.get_freq(name)
inmap = data.get_mapfile(name, dset)
nside = int(spherelib.read_key(inmap, "NSIDE"))
unit_to = 'K_RJ'
unit_fr = data.parse_unit(inmap)
logger.info ("original unit is %s"%unit_fr)
cv_unit = data.get_unit_factor(name, fr=unit_fr, to=unit_to)
logger.info ("using conversion factor %e"%cv_unit)
hook("config", globals()) # fsky [= 0 for full sky map , 40 for calibration] , do_polar [False for calibration]
polar = data.is_polarized(name) and do_polar
fields = ['I','Q','U'] if polar else ['I']
fieldsht = ['T','E','B'] if polar else ['T']
lmax = 4000
logger.info("fsky = %d for galactic mask "%fsky_target)
psmask, gmask = get_mask(nside, fsky=fsky_target)
# compute alm of full sky map
maps = list()
alms = list()
for field in fields:
mapF = healpy.read_map(inmap, field=fields.index(field))
mapF[mapF<-1e30] = 0 # replace undef with 0
if field=='I':
mn = np.mean(np.double(mapF[mapF!=0]))
mapF[mapF!=0]-= mn
if freq<545:
mapF = do_ps_subtraction(mapF, name, dset)
mapF*= psmask
spherelib.fastfill(mapF, 1e-12, 1000, 256)# fill point sources holes
mapF*= cv_unit #convert in new unit
plot_map_in_pipe(mapF, "map_full_sky", "%s-%s-%s"%(dset,name,field))
maps.append(mapF) # full sky map
alms = healpy.map2alm(maps, lmax, iter=1) # alm of full sky map
alms = [alms] if len(alms)!=3 else list(alms)
write_alm(alms, fieldsht, "alm_full_sky_%s_%s"%(dset,name))
# compute alm of masked map
for field, i_f in zip(fieldsht, range(len(fieldsht))):
if field=='T':
mapF = maps[i_f]
spherelib.regress (mapF, gmask, 2)
else:
mapF = healpy.alm2map (alms[i_f], nside)
mapF*= gmask
plot_map_in_pipe(mapF, "map_masked", "%s-%s-%s"%(dset,name,field))
alms[i_f] = healpy.map2alm (mapF, lmax, iter=1)
write_alm(alms, fieldsht, "alm_masked_%s_%s"%(dset,name))
# compute transfer function
bl = data.get_bl(name)[0:lmax+1]
pwf = healpy.pixwin(nside)[0:lmax+1]
ref_beam = 5
max_inv = 1000
bl_ref = spherelib.gaussbeam(ref_beam, lmax)
f = np.squeeze(bl_ref)/(np.squeeze(bl))
f[f>max_inv] = max_inv
f /= pwf
healpy.write_cl(get_data_fn("ibl%s.fits"%name), f)
# plot spectra
fsky = np.mean(gmask**2)
plt.figure()
plt.subplot(211)
for field, j in zip(fieldsht, range (len(fieldsht))):
cl = healpy.alm2cl(alms[j])/fsky
plt.semilogy(cl[2:], label="masked %s"%field)
plt.semilogy(cl[2:]*f[2:]**2, label="unbeamed %s"%field)
plt.subplot(212)
plt.semilogy(f)
plt.savefig(get_data_fn("cls.png"))
set_output([(dset,name)])
save_param(["ref_beam","lmax", "unit_to", "fsky"])
fsky_str = str(fsky)[0:5]
globals().update(locals())
expose(["lmax", "fsky_str"])
if __name__=="__builtin__":
main()