DX12 SMICA release: impact of the 217GHz recalibration¶
- X1 refers to low ell, LFI+HFI zone
- X2 refers to high ell HFI only plik zone
- dx12c217 refers to a pipeline where the 100GHz and the 217GHz a recalibrated with SMICA (only the 217GHz moves significantly ie 2 per mil)
- dx12no217 is obtained with fixed, official calibration
In [2]:
%matplotlib inline
import os
os.environ["RELEASE"] = "dx12c217"
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
rdir1 = "/workdir/lejeune/planck/outputs/dx12c217"
rdir2 = "/workdir/lejeune/planck/outputs/dx12no217"
Weight derived plots¶
In [11]:
def plotweight(r1, r2, case='T'):
""" Paper plot of SMICA weights.
"""
bins, lmean = pu.get_bins_plot(2)
for iplot,icase,zname in pu.getcases(case):
plt.figure(figsize=(15,10))
for rdir in [r1,r2]:
os.environ["RELEASE"] = rdir.split("/")[-1]
reload(pu)
ax = plt.subplot(1,2,[r1,r2].index(rdir)+1)
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+"_"+rdir.split("/")[-1])
plt.axis([np.sqrt(3),np.sqrt(4001),-3,3.5 ])
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(rdir1, rdir2, "T")
In [14]:
def plotnoise(r1,r2,case='T'):
""" Make a plot of the detector noise contribution as a function of the multipole.
"""
bins_plot,lmean = pu.get_bins_plot()
for iplot,icase,zname in pu.getcases(case):
plt.figure(figsize=(15,10))
for rdir in [r1,r2]:
os.environ["RELEASE"] = rdir.split("/")[-1]
reload(pu)
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,[r1,r2].index(rdir)+1)
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+"_"+rdir.split("/")[-1])
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(rdir1, rdir2, "T")
In [20]:
TT CMB power spectrum compared to plik on the 100GHz plik mask¶
In [36]:
def plotdlplik(rdir1, rdir2, maskfn, mllfn, iref, lmax=2508):
"""
"""
mll = np.load(mllfn)
fsky = mll[0,:].sum()
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
clplik0 = clplik.copy()
clplik = np.dot(mll, clplik/fsky)
bins = sp.uniformbin(30,1500,30)
lmean = np.squeeze(sp.get_lmean(bins))
plt.figure(figsize=(15,10))
ax = plt.subplot(111)
for rdir,co in zip([rdir1, rdir2], ["b", "r"]):
os.environ["RELEASE"] = rdir.split("/")[-1]
reload(pu)
cfit,cefit,lmean_fit, bins_fit = pu.extractCMB(rdir, "T", z="X2")
cqplik = sp.bin_cl(clplik0, bins_fit)
cl = hp.read_cl(pu.getclfn(rdir, maskfn, "full")) - hp.read_cl(pu.getclfn(rdir, maskfn, "hmhd"))
plt.title("mask plik %sGHz"%iref)
label = "[smica full- smica hmhd]"
plt.plot(lmean, sp.bin_cl(ells*(ells+1)*(clplik-cl)/(2*np.pi), bins),label=rdir.split("/")[-1],color=co)
#ax.errorbar(lmean_fit, lmean_fit*(lmean_fit+1)*(cqplik-cfit)/(2*np.pi), yerr=0*lmean_fit*(lmean_fit+1)*cefit/(2*np.pi),color=co, fmt='+')
plt.plot(lmean, np.zeros((len(lmean))), color="k")
#plt.plot(lmean, sp.bin_cl(clplik/cl, bins),label="plik / "+label)
plt.axis([2,1500,-15,15])
plt.ylabel(r"$D_\ell^{plik}-D_\ell^{ilc} \/ \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"]:
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
#for rdir in [rdir1, rdir2]:
# os.environ["RELEASE"] = rdir.split("/")[-1]
# reload(pu)
# pu.compute_cls(rdir, maskfn, lmax=lmax)
plotdlplik(rdir1, rdir2, 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)/(2\pi) \rightarrow 0$ as $C_\ell$ does
TT inter calibration¶
In [133]:
def draw_calib(ax, lstcase, a_ratios, errors, names, axis):
nmap = len(names)
plt.plot(np.arange(nmap), np.ones((nmap)),color="k")
colors = ["b","r","c","g","m","y","k","grey","brown","orange","navy"]
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)])
ax.errorbar(np.arange(nmap)+0.06*j, a, yerr=e, color=colors[j], fmt='o', label=case)
plt.xticks(np.arange(nmap), tuple(names),rotation=45)
ax.axis()
plt.grid(True)
ax.axis(axis)
leg = ax.legend(loc='upper center')
pu.format_legend_and_axes(leg,ax)
def plot_calib(lstcase, a_ratios, errors, names, filename=None):
""" """
nmap = len(names)
plt.figure(figsize=(16,6))
ax = plt.subplot(121)
draw_calib(ax, lstcase, a_ratios, errors, names, [0.5, nmap-1+0.5, 0.95, 1.05])
ax2 = plt.subplot(122)
draw_calib(ax2, lstcase, a_ratios, errors, names, [0.5, nmap-1+0.5-1, 0.995, 1.006])
if filename is not None:
plt.savefig(filename)
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)
Difference maps between dx11 and dx12¶
In [43]:
%%capture
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")
m12c217 = hp.read_map("/workdir/lejeune/planck/outputs/dx12c217/dx12c217_smica_harmonic_int_cmb_raw_full.fits")
m12no217 = hp.read_map("/workdir/lejeune/planck/outputs/dx12no217/dx12no217_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)
dmapc217 = hp.smoothing(m12c217, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
dmapno217 = hp.smoothing(m12no217, fwhm=80*np.pi/(180*60.), lmax=2000, iter=0)
In [38]:
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.
Impact of the 217GHz recalibration on CMB maps¶
In [42]:
fignum = 0
plt.figure(fignum,figsize=(10,8))
hp.mollview(pmask15d*pmask17*(m12no217-m12c217), sub=(1,2,1), min=-10, max=10, title="no217-c217")
hp.mollview(pmask15d*pmask17*(dmapno217-dmapc217), sub=(1,2,2), min=-10, max=10, title="no217-c217 @80arcmin")
