DX12 SMICA release: impact of masking the co in X2 analysis¶
- X1 refers to low ell, LFI+HFI zone
- X2 refers to high ell HFI only plik zone
- dx12c217no353 refers to the baseline where the co is masked
- dx12nomaskco is obtained with a new analysis mask where co is not masked
In [46]:
%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/dx12Tc217no353"
rdir2 = "/workdir/lejeune/planck/outputs/dx12Tnomaskco"
Weight derived plots¶
In [47]:
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 [56]:
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, ls in zip([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(111)#,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], ls=ls)
plt.semilogy(np.sqrt(lmean), sp.bin_cl(sumn, bins_plot), lw=2, color='k', label='sum '+rdir.split("/")[-1], ls=ls)
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 [50]:
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 [53]:
%%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/dx12Tc217no353/dx12Tc217no353_smica_harmonic_int_cmb_raw_full.fits")
m12no217 = hp.read_map("/workdir/lejeune/planck/outputs/dx12Tnomaskco/dx12Tnomaskco_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 [ ]:
Impact of the 217GHz recalibration on CMB maps¶
In [54]:
fignum = 0
plt.figure(fignum,figsize=(10,8))
hp.mollview(pmask15d*pmask17*(m12no217-m12c217), sub=(1,2,1), min=-10, max=10, title="wo - w mask co")
hp.mollview(pmask15d*pmask17*(dmapno217-dmapc217), sub=(1,2,2), min=-10, max=10, title="wo - w mask co@80arcmin")
