DX12 SMICA release¶
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'font.size': 14})
matplotlib.rcParams.update({'legend.labelspacing':0.01})
matplotlib.rcParams.update({'figure.subplot.bottom': 0.1, 'figure.subplot.top': 0.9,
'figure.subplot.right': 0.9, 'figure.subplot.left': 0.125 }) # the bottom of the subplots of the figure
import spherelib as sp
import pipetools as pt
import smicatools as st
import plotutils as pu
reload(pu)
import healpy as hp
rdir = "/workdir/lejeune/planck/outputs/dx12"
Weight derived plots¶
In [228]:
def plotweight(case='T'):
""" Paper plot of SMICA weights.
"""
bins, lmean = pu.get_bins_plot(20)
plt.figure(figsize=(15,10))
for iplot,icase,zname in pu.getcases(case):
ax = plt.subplot(1,2,iplot)
names, colors = pu.get_names(icase, zone=zname)
W, A = pu.extractWl(rdir, icase, zname)
nmap,lmax = np.shape(W)
for m in range(nmap):
fq = sp.bin_cl(W[m,:]*1e-3/A[m], bins)
plt.plot(np.sqrt(lmean), fq, linewidth=2, label=names[m], color=colors[m])#, marker='+')
xarr, xlab = pu.get_xticks(ax)
plt.title(icase+"_"+zname)
plt.axis([np.sqrt(3),np.sqrt(4001),-2,3 ])
leg = plt.legend(loc="lower center",bbox_to_anchor=(0.5, -0.021))
leg,ax = pu.format_legend_and_axes(leg, ax)
plt.xlabel("$\ell$",fontsize=16)
plt.ylabel(r"Weight $\left[ \mu \mathrm{K} / \mu \mathrm{K^{RJ}} \right]$",fontsize=20)
plotweight("T")
plotweight("E")
In [230]:
def plotnoise(case='T'):
""" Make a plot of the detector noise contribution as a function of the multipole.
"""
bins_plot,lmean = pu.get_bins_plot()
plt.figure(figsize=(15,10))
for iplot,icase,zname in pu.getcases(case):
names, colors = pu.get_names(icase, zone=zname)
W,A = pu.extractWlraw(rdir, icase, zname)
N = pu.extractN(rdir, icase, zname)
nmap,lmaxp1 = N.shape
sumn = np.zeros((lmaxp1))
ax = plt.subplot(1,2,iplot)
plt.title(icase+"_"+zname)
for m in range(nmap):
data = W[m,:]**2*N[m,:]
sumn += data
plt.semilogy(np.sqrt(lmean), sp.bin_cl(data, bins_plot), color=colors[m], label=names[m])
plt.semilogy(np.sqrt(lmean), sp.bin_cl(sumn, bins_plot), lw=2, color='k', label='sum')
xarr, xlab = pu.get_xticks(ax)
plt.title(icase+"_"+zname)
leg= plt.legend(loc="lower left")
leg,ax = pu.format_legend_and_axes(leg, ax)
plt.xlabel("$\ell$",fontsize=16)
plt.ylabel(r"$\mathbf{w_i}^2 N_i \/ \left[ \mu \mathrm{K}^2 \right]$",fontsize=20)
plt.axis([None,None,1e-12,None])
plotnoise("T")
plotnoise("E")
In [20]:
def plotcmb(case='T', lmax=2500):
""" Make a plot of the detector contribution to the cmb map as a function of the multipole.
"""
bins_plot,lmean = pu.get_bins_plot(lmax=lmax)
plt.figure(figsize=(15,8))
ells = np.arange(lmax+1)
for iplot,icase,zname in pu.getcases(case):
names, colors = pu.get_names(icase, zone=zname)
W, A = pu.extractWlraw(rdir, icase, zname)
cmbbf = pt.get_reference_cmb_powspec(icase+icase)[0:lmax+1]
ax = plt.subplot(1,2,iplot)
plt.title(icase+"_"+zname)
for m in range(len(names)):
G = np.array([np.sum(W[m,ell]*np.squeeze(A[m])) for ell in ells])
data = sp.bin_cl(ells*(ells+1)*cmbbf*G**2/(2*np.pi), bins_plot)
plt.semilogy(np.sqrt(lmean), data, color=colors[m], label=names[m])
plt.axis([None,np.sqrt(lmax+10),np.max(ells**2*cmbbf/6.)*1e-6,None])
xarr, xlab = pu.get_xticks(ax)
plt.title(icase+"_"+zname)
leg= plt.legend(loc="upper right")
leg = pu.format_legend_and_axes(leg,ax)
plt.xlabel("$\ell$",fontsize=16)
plt.ylabel(r"$D_\ell G_\ell^2 \/ \left[ \mu \mathrm{K}^2 \right]$",
fontsize=20)
#plotcmb("T")
#plotcmb("E")
TT CMB power spectrum compared to plik on the 100GHz plik mask¶
In [125]:
def plotdlplik(rdir, maskfn, mllfn, iref, lmax=2508):
"""
"""
mll = np.load(mllfn)
fsky = mll[0,:].sum()
clhm1xhm2 = hp.read_cl(pu.getclfn(rdir, maskfn, "hm1xhm2"))
clfull = hp.read_cl(pu.getclfn(rdir, maskfn, "full")) - hp.read_cl(pu.getclfn(rdir, maskfn, "hmhd"))
dlplik, dlcmb = pu.load_plik(iref, iref)
ells = np.arange(lmax+1, dtype=float)
clplik = dlplik*(2*np.pi)/(ells*(ells+1))
clplik[np.isnan(clplik)] = 0
clplik = np.dot(mll, clplik/fsky)
bins = sp.uniformbin(30,1500,30)
lmean = np.squeeze(sp.get_lmean(bins))
W,A = pu.extractWlraw(rdir, "T", "X2")
N = pu.extractN(rdir, "T", "X2")
Nc = np.array([np.sum(W[:,ell]**2*N[:,ell]) for ell in ells])
plt.figure(figsize=(15,10))
for subp in range(2):
ax = plt.subplot(2,1,subp+1)
plt.title("mask plik %sGHz"%iref)
for cl,label in zip([clhm1xhm2, clhm1xhm2*1.002, clfull],
["smica hm1xhm2","1.002x smica hm1xhm2", "[smica full- smica hmhd]"]):
if subp==0:
plt.plot(lmean, sp.bin_cl(ells*(ells+1)*(clplik-cl)/(2*np.pi), bins),label="plik - "+label)
else:
plt.plot(lmean, sp.bin_cl(clplik/cl, bins),label="plik / "+label)
if subp==0:
plt.plot(lmean, sp.bin_cl(ells*(ells+1)*(Nc)/(2*np.pi), bins), label="W^T N W")
plt.axis([2,1500,-30,30])
plt.ylabel(r"$D_\ell \/ \left[ \mu \mathrm{K}^2 \right]$",fontsize=20)
else:
plt.axis([2,1500,0.99,1.01])
plt.ylabel(r"$C_\ell \/ \left[ \mu \mathrm{K}^2 \right]$",fontsize=20)
plt.grid()
leg= plt.legend(loc="lower center")
leg = pu.format_legend_and_axes(leg,ax)
plt.xlabel("$\ell$",fontsize=16)
for iref in ["100", "143", "217"]:
maskfn = "/workdir/lejeune/planck/data/benabed/plik_common/temperature_%s.fits"%iref
mllfn = "/workdir/lejeune/planck/data_script/mll%s_bin.txt.npy"%iref
lmax = 2508
#pu.compute_cls(rdir, maskfn, lmax=lmax)
plotdlplik(rdir, maskfn, mllfn, iref, lmax=lmax)
Note: On the top panel we plot $\ell (\ell+1) [C_\ell^{plik} - \alpha C_\ell^{smica}]/(2\pi)$ with $\alpha=$1 (blue) or 1.002 (green). When one makes the difference between the 2 by eye, it goes to $\ell(\ell+1)C_\ell^{smica}(1.002-1)/(2pi) \rightarrow 0$ as $C_\ell$ does
TT inter calibration¶
In [42]:
def plot_calib(lstcase, a_ratios, errors, names):
""" """
nmap = len(names)
plt.figure(figsize=(8,5))
ax = plt.subplot(111)
plt.plot(np.arange(nmap), np.ones((nmap)),color="k")
colors = ["b","r"]
lines = list()
labs = list()
for case,a,e,j in zip(lstcase, a_ratios, errors, range(len(lstcase))):
strtoprint = "\n".join(["%s %.5f pm %.6f"%(names[m], a[m], e[m]) for m in range(nmap)])
plt.errorbar(np.arange(nmap)+0.2*j, a, yerr=e, color=colors[j], fmt='o', label=case)
plt.xticks(np.arange(nmap), tuple(names),rotation=45)
plt.axis([-0.5, nmap-1+0.5-1, 0.98, 1.05])
plt.grid(True)
plt.axis([-0.5, nmap-0.5, 0.999, 1.01])
leg = plt.legend(loc='upper center')
pu.format_legend_and_axes(leg,ax)
dcalib = dict({"dx11":([1.00524,1.00079,1,1.0029,1.008,1.01740],[0.00042,0.00016,0,0.00026,0.0016,0.0162]),
"dx12":([1,1.00017,1,1.00166,1.00256,1.01807],[0,0.00015,0,0.00023,0.0017,0.0183])})
names = ["070", "100", "143", "217", "353", "545"]
plot_calib(dcalib.keys(), [v[0] for v in dcalib.values()], [v[1] for v in dcalib.values()],names)
Polarized foreground SED and angular power spectra¶
In [54]:
freqs = [30,44,70,100,143,217,353]
xrange = np.arange(len(freqs))
xfit, dref = plotutils.computeSED(rdir, freqs)
xref, info = plotutils.getrefSED(freqs, dref=dref)
colors = dict({"sync":"r", "dust":"b"})
lws = dict({"E":1.5, "B":"0.5"})
plt.figure(figsize=(12,5))
for subp in range(2):
ax = plt.subplot(1,2,subp+1)
for iis in range(0,6,3):
data = xref["f0"][iis+1:iis+3]
comp = xref["name"][iis][0:4]
data /= xref["f0"][iis] if subp==1 else 1.
plt.fill_between(xrange, data[0], data[1], alpha=0.25, color=colors[comp], label=info[comp])
for iis in [0,2,5,7]:
data = xfit["f0"][iis:iis+2] ; data[data==0] = None
comp = xfit["name"][iis][0:4]
rsed = xref["f0"][xref["name"]==comp][0]
if subp==1:
if xfit[iis][1][-1]=="B": # rescale B in second plot
data *= rsed[dref[comp][0]]/data[0][dref[comp][0]]
data /= rsed
plt.errorbar(xrange+iis*2e-2, data[0], data[1], marker='^', color=colors[comp], lw=lws[xfit["name"][iis][-1]])
if subp==0:
leg = plt.legend(loc="upper left")
plt.axis([-0.25, 6.25, 1e-7, 1e-3])
ax.set_yscale("log")
plt.ylabel(r"SED $\left[ \mu \mathrm{K} \right]$",fontsize=20)
else:
leg = plt.legend([plt.Line2D([0],[0],lw=lw) for lw in lws.values()], lws.keys(),loc="lower left")
plt.axis([-0.25, 6.25, 0.5, 1.15])
plt.ylabel(r"SED ratio $\left[ \mu \mathrm{K} \right]$",fontsize=20)
plt.xticks(xrange, tuple([str(freq) for freq in freqs]),rotation=45)
pu.format_legend_and_axes(leg,ax)
In [119]:
import plotutils
reload(plotutils)
xfit, lmean, dref = pu.computeAPS(rdir, freqs)
print "Dust power spectra is computed at %dGHz"%freqs[dref["dust"]]
print "Sync power spectra is computed at %dGHz"%freqs[dref["sync"]]
colors = dict({"sync":"r", "dust":"b", "$\ell^{-2.42}$":"k"})
lws = dict({"E":1.5, "B":"0.5"})
fac = dict({"sync":1., "dust":1.})
powerlaw = lmean**-2.42
l50 = np.where(lmean>80)[0][0]
Ad = "Dust BB amplitude at ell=80 is %e uK^2\n"%(xfit['f0'][xfit['name']=="dustB"][0][l50])
print Ad
xref = dict({"sync":powerlaw*xfit['f0'][xfit['name']=="syncE"][0][l50]/powerlaw[l50],
"dust":powerlaw*xfit['f0'][xfit['name']=="dustE"][0][l50]/powerlaw[l50]})
plt.figure(figsize=(12,5))
for subp in range(2):
ax = plt.subplot(1,2,subp+1)
for iis in range(0,8,2):
comp = xfit['name'][iis][0:4]
data = xfit['f0'][iis:iis+2].copy()
data /= xref[comp] if subp==1 else 1/fac[comp]
data[0] += np.float(fac.keys().index(comp)) if subp==1 else 0
plt.errorbar(lmean, data[0], data[1], color=colors[comp], lw=lws[xfit['name'][iis][-1]])
if subp==0:
ax.set_yscale("log")
ax.set_xscale("log")
for c in ['sync', 'dust']:
for f,lw in zip([1.,0.5], [1.5, 0.5]):
plt.plot(lmean, xref[c]*fac[c]*f, color='k', lw=lw)
leg = plt.legend([plt.Line2D([0],[0],color=color) for color in colors.values()], colors.keys(),loc="upper right")
plt.axis([3,lmean[-1]+50, 2e-4, 2e4])
plt.ylabel(r"$C_\ell$ $\left[ \mu \mathrm{K^2} \right]$",fontsize=20)
else:
for f,lw in zip([0.5,1.,1.5,2.], [0.5,1.5,0.5,1.5]):
plt.plot(lmean, f*np.ones((len(lmean))), color='k', lw=lw)
leg = plt.legend([plt.Line2D([0],[0],lw=lw) for lw in lws.values()], lws.keys(),loc="upper left")
plt.ylabel(r"$C_\ell$ ratio",fontsize=20)
pu.format_legend_and_axes(leg,ax)
plt.xlabel(r"$\ell$",fontsize=20)
Difference maps between dx11 and dx12¶
In [5]:
m2015 = hp.read_map("/workdir/lejeune/planck/outputs/dx11d/dx11d_smica_harmonic_cmb_I.fits")
m2015d = hp.read_map("/workdir/lejeune/planck/outputs/new/dx11d/dx11d_smica_harmonic_int_cmb_raw_full.fits")
m2017 = hp.read_map("/workdir/lejeune/planck/outputs/dx12/dx12_smica_harmonic_int_cmb_raw_full.fits")
mrc4 = hp.read_map("/workdir/lejeune/planck/outputs/rd12_rc3/rd12_rc3_smica_harmonic_int_cmb_raw_full.fits")
pmask15 = hp.read_map("/workdir/lejeune/planck/outputs/dx11d/dx11d_smica_harmonic_mask_I.fits")
pmask15d = hp.read_map("/workdir/lejeune/planck/outputs/new/dx11d/dx11d_smica_harmonic_int_mask.fits.gz")
pmask17 = hp.read_map("/workdir/lejeune/planck/outputs/dx12/dx12_smica_harmonic_int_mask.fits.gz")
dmap15 = hp.smoothing(m2015, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
dmap15d = hp.smoothing(m2015d, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
dmap17 = hp.smoothing(m2017, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
dmaprc = hp.smoothing(mrc4, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
In [6]:
fignum = 0
plt.figure(fignum,figsize=(10,8))
hp.mollview(pmask15*pmask17*(m2015-m2017), sub=(3,2,1), min=-10, max=10, title="2015-2017")
hp.mollview(pmask15*pmask17*(dmap15-dmap17), sub=(3,2,2), min=-10, max=10, title="2015-2017 @80arcmin")
hp.mollview(pmask15d*pmask17*(m2015d-m2017), sub=(3,2,3), min=-10, max=10, title="2015d-2017")
hp.mollview(pmask15d*pmask17*(dmap15d-dmap17), sub=(3,2,4), min=-10, max=10, title="2015d-2017 @80arcmin")
hp.mollview(pmask15d*pmask17*(mrc4-m2017), sub=(3,2,5), min=-10, max=10, title="rc3-2017")
hp.mollview(pmask15d*pmask17*(dmaprc-dmap17), sub=(3,2,6), min=-10, max=10, title="rc3-2017 @80arcmin")
The 2015d map is the one which has been built after the 2nd release. The pipeline is exactly the same as in 2015 except for the recalibration of the 44GHz and 70GHz. To be more precise, the SMICA fit was kept fixed, only the CMB mixing column was corrected in the ILC filter computation.
The main effect causing the difference between 2015 and 2017 is this LFI recalibration bug found in 2016. Correcting the calibration factor in the ILC filter suppress most of the difference we see (compare 2015-2017 with 2015d-2017). The remaining signal in the 2015d-2017 difference can also be attributed to this bug, as it correspond to feature seen in the 44GHz (known free-free spot near the galactic plane and ophiucus region). One may argue that a wrong CMB mixing column has also impacted the SMICA fit of the covariance matrices which enter in the ILC filter computation.
In any case, this power vanishes in the rc3-2017 difference where a new SMICA fit is done on the same LFI maps (DX11) but with RC3 HFI channels.