#depend .planck_release
import os
import healpy
import spherelib
import planckdata as data
import numpy as np
from matplotlib import pylab as pl
import matplotlib.pyplot as plt
[docs]def fit_postproc(fit_cat, maxxi2, high_snr):
""" Post processing of the fit catalog includes snr sorting and chi2 cut
"""
idx1 = pl.find(fit_cat[:,2]/fit_cat[:,5]>high_snr) # high snr
idx2 = pl.find(fit_cat[:,-1]<maxxi2)
fit_cat[np.intersect1d(idx1, idx2), -1] = np.inf # reject high snr good fit
idxsort = np.argsort(fit_cat[:,2]/fit_cat[:,5])
fit_cat = fit_cat[idxsort, :]
idxsort = len(idxsort)-1-np.arange(len(idxsort))
fit_cat = fit_cat[idxsort, :] # sort by snr
idx_good = pl.find(fit_cat[:,-1]<maxxi2)
idx_bad = pl.find(fit_cat[:,-1]>=maxxi2)
idx_cut = np.intersect1d(pl.find(fit_cat[:,-1]>=maxxi2), pl.find(np.isfinite(fit_cat[:,-1]))) # count number of sources in the galactic cut
logger.info ("After fitting : %d good vs %d bad sources (including %d in the latitude cut)"%(len(idx_good), len(idx_bad), len(idx_cut)))
save_products(get_data_fn("fit_cat_postproc.pkl"), locals(), ["idx_good", "idx_bad", "fit_cat"])
fig = plt.figure()
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(fit_cat[pl.find(fit_cat[:,-1]<3*maxxi2),-1], 100, facecolor='green')
plt.title ("Fitted point sources: maxxi2=%d"%int(maxxi2))
plt.savefig (get_data_fn("histogram-fitted.png"))
return idx_good, idx_bad, fit_cat
[docs]def compute_variance (name):
""" Compute variance map from half ring
"""
Vmap = healpy.read_map(data.get_mapfile(name, 'hr1'))
Vmap[Vmap<-1e30] = 0 # replace undef with 0
Vmap *= -1.
Vmap += healpy.read_map(data.get_mapfile(name, 'hr2'))
Vmap[Vmap<-1e30] = 0 # replace undef with 0
Vmap /= 2. # noise is (first-last)/2.
Vmap = Vmap**2
nside = healpy.get_nside(Vmap)
Vmap = healpy.ud_grade(Vmap, nside/8) # compute variance over big healpix pixel
Vmap = healpy.ud_grade(Vmap, nside)
return Vmap
[docs]def main():
"""
Perform point source subtraction, and compute corresponding point source mask
For freq>=545, use a point source mask only.
"""
name = get_input()
freq = data.get_freq(name)
beam = data.get_fwhm_beam(name)
# Read input maps
Imap = healpy.read_map(data.get_mapfile(name, 'full'))
Imap[Imap<-1e30] = 0 # replace undef with 0
Imap[Imap!=0] -= np.mean(np.double(Imap[Imap!=0])) # remove monopole
plot_map_in_pipe(Imap, "input_map", name)
# Sources selection
snr_thresh = 5 if freq<=353 else 7.5
cat,theta,phi,snr = data.load_catalog(name, snr_thresh)
write_cat(cat, "det%d"%snr_thresh, name)
# Gaussian fit
latcut = 5 # galactic latitude cut above which sources are fitted (degree)
maxxi2 = dict({30:2,44:2,70:2,100:3,143:3,217:6,353:6}) # chi2 maximum value of the gaussian fit for subtraction
high_snr = 80 # snr thresh above which the fit can't be trusted
if freq<545:
Vmap = compute_variance(name)
plot_map_in_pipe(Vmap, "variance_map", name)
fit_cat = spherelib.ps_gaussfit (Imap, Vmap, cat, beam, nsigma=1, latcut=latcut, maxxi2=maxxi2[freq])
write_cat(fit_cat, "fit%d"%snr_thresh, name)
idx_good, idx_bad, fit_cat = fit_postproc(fit_cat, maxxi2[freq], high_snr)
# Make point source mask
beam = data.get_fwhm_beam(name)
sigma_beam = 1.5 # point source radius cut used for masking
mask_snr_thresh = 8 # snr thresh above which the sources are masked
if freq<545:
snr = fit_cat[idx_bad,2]/fit_cat[idx_bad,5]
idx = pl.find(snr>=mask_snr_thresh)
theta, phi = tuple(fit_cat[idx_bad[idx],0:2].T)
snr = snr[idx]
else:
phi = phi[snr>mask_snr_thresh]
theta = theta[snr>mask_snr_thresh]
snr = snr[snr>mask_snr_thresh]
psmask = np.zeros((len(Imap)))
spherelib.cat2mask(psmask, theta, phi, sigma_beam*beam, snr=snr)
fn = get_data_fn("psmask-%s.fits"%(name))
healpy.write_map(fn, psmask)
logged_subprocess(["gzip", fn])
# Make some plot
plot_map_in_pipe(psmask, "psmask", name)
if freq<545:
plot_ps_subtraction(Imap*psmask, name, fit_cat, idx_bad, idx_good, beam)
set_output(get_input())
save_param(["snr_thresh","sigma_beam", "maxxi2" ,"latcut", "high_snr", "mask_snr_thresh"])
release = os.environ["RELEASE"]
expose(["release"])
globals().update(locals())
if __name__=="__builtin__":
main()