I have recently been working on the cross-correlation of CIB intensity and CMB lensing from the Planck data. I've had a few hiccups with the format of the lensing data and thus decided to document this properly in a notebook, just in case others ran into similar issues. I use the PR1 and PR2 data release, each of which has a slightly different format.

The goal is to understand the data format, units, and conventions to use this map for CIB-lensing cross correlations. Furthermore, I demonstrate how to deal with maps, $a_{lm}$, and $C_l$ in healpy while also dealing with sky fractions and window functions.

A big thanks goes to Duncan Hanson, author of the Planck 2015 lensing paper, for help with the data format and for clarifying various issues that I ran into.

In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In [2]:
import numpy as np
import healpy as hp
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['figure.figsize'] = (12, 6)
rcParams['text.usetex'] = True

from joblib import Memory
mem = Memory(cachedir='.joblib')

import warnings
warnings.filterwarnings('ignore')

We start by defining some parameters for this notebook. Update the datapath if you want to download the data yourself and re-run the notebook.

In [3]:
NSIDE = 2048
LMAX = 2048  # The lensing data are limited to 10 < l < 2048
l = np.arange(LMAX)
datapath = '/users/dlenz/projects/planck/data/cmb_lensing/'
pixfunc = hp.pixwin(NSIDE)

The pixel window function can be generated with healpy and is used to correct the $C_l$ for the finite pixelsize, which we demonstrate later.

In [4]:
plt.loglog(l, pixfunc[:LMAX])
plt.xlabel('l')
plt.ylabel('Pixel window function')
Out[4]:
<matplotlib.text.Text at 0x10aba5240>

We furthermore create a helper to smooth the mask, also caching the result.

In [5]:
@mem.cache
def smooth_mask(mask):
    smoothing_scale = 2.  # in degree
    mask_apod = hp.smoothing(mask.astype(float), fwhm=np.radians(smoothing_scale))
    mask_apod = np.clip(mask_apod, a_min=0., a_max=None)
    return mask_apod

PR1 from 2013

Load data

The PR1 data are taken from here and described in Planck collaboration (2013) XVII. The data product is also accessible from the Planck Legacy Archive (PLA).

First, load the map and mask, convert to ring ordering, and apodize the mask.

In [6]:
pr1_ext1 = fits.getdata(datapath + 'pr1/COM_CompMap_Lensing_2048_R1.10.fits', ext=1)
phibar = hp.reorder(pr1_ext1['PHIBAR'], n2r=True)
mask_pr1 = hp.reorder(pr1_ext1['MASK'], n2r=True)
In [7]:
# Note that smoothing will take a while (~10 minutes) with this high-resolution mask and the large kernel
mask_pr1_apod = smooth_mask(mask_pr1)
skyfrac_pr1 = mask_pr1_apod.sum()/mask_pr1_apod.size
WARNING:root:[MemorizedFunc(func=<function smooth_mask at 0x10ac4cae8>, cachedir='.joblib/joblib')]: Clearing cache .joblib/joblib/__main__--Users-dlenz-projects-cibxphi-notebooks-__ipython-input__/smooth_mask
________________________________________________________________________________
[Memory] Calling __main__--Users-dlenz-projects-cibxphi-notebooks-__ipython-input__.smooth_mask...
smooth_mask(array([1, ..., 1], dtype=uint8))
Sigma is 50.959308 arcmin (0.014823 rad) 
-> fwhm is 120.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
____________________________________________________smooth_mask - 525.7s, 8.8min

Inspect the lensing potential:

In [8]:
hp.mollview(phibar * mask_pr1_apod, title='', unit=r'$\bar\Phi$')

Second, also load the response function $R_l^{\Phi\Phi}$ and the noise term $N_l^{\Phi\Phi}$. These are later used to correct the derived angular power spectrum of the lensing potential.

In [9]:
pr1_ext2 = fits.getdata(datapath + 'pr1/COM_CompMap_Lensing_2048_R1.10.fits', ext=2)
rlpp = pr1_ext2['RLPP']
nlpp = pr1_ext2['NLPP']
In [10]:
plt.plot(rlpp, label=r'$R_l^{\Phi\Phi}$')
plt.plot(nlpp, label=r'$N_l^{\Phi\Phi}$')
plt.loglog()

# labels & legend
plt.legend()
plt.xlabel(r'$l')
plt.ylabel(r'$X_l$');

Estimate $C_l^{\Phi\Phi}$

We simply use hp.anafast() here to compute the $C_l$ for a given map. Ideally, you would also want to correct for the mode-mode coupling through e.g. the MASTER method, which is done in polspice or in David Alonso's NaMaster code (github).

Below, we compute $C_l^{\bar{\Phi}\bar{\Phi}}$ and correct its amplitude for the sky fraction that is covered by the mask.

In [11]:
Cl_pr1 = hp.anafast(phibar * mask_pr1_apod, lmax=LMAX-1) / skyfrac_pr1

To correct this power spectrum for the noise spectrum $N_l$ and the response function $R_l$, we compute $C_l^{\Phi\Phi} = \left(C_l^{\bar{\Phi}\bar{\Phi}} - N_l^{\Phi\Phi}\right) / (R_l^{\Phi\Phi})^2$.

In [12]:
Cl_pr1_pp = (Cl_pr1[:LMAX] - nlpp[:LMAX])/(rlpp[:LMAX])**2

If we correct for this and apply these two changes, we get the following, which looks like the paper result (modulo the binning). Note that most commonly, the plots show $\sim l^4\,C_l$.

In [13]:
plt.plot(l, 1.e7 *  (l *(l+1.))**2 / 2. / np.pi * Cl_pr1_pp)
plt.semilogx()

# limits
# The PR1 data are difficult to interpret at the largest scale, hence we focus only on l > 15 
plt.xlim([20, LMAX])
plt.ylim([-3, 5])

# labes
plt.xlabel(r'$l$')
plt.ylabel(r'$[l(l+1)]^2/2\pi\, C_l\ [\times 10^7]$');

PR2 from 2014

For the PR2, the data comes in two different formats:

  • $C_l$ and $C_l + N_l$ for the lensing convergence $\kappa$
  • $a_{lm}$ and a mask for the lensing convergence $\kappa$

They are taken from here, and described in Planck collaboration (2015) XV.

The lensing convergence $\kappa$ is related to the lensing potential $\Phi_{l}$ through:

$\kappa_{l} = 0.5 l(l+1) \Phi_{l}$

Load data

We begin by loading all the available data

In [14]:
alm_pr2 = hp.read_alm(datapath + 'pr2/dat_klm.fits')
mask_pr2 = hp.read_map(datapath + 'pr2/mask.fits.gz')
NSIDE = 2048
ORDERING = RING in fits file
INDXSCHM = IMPLICIT
In [15]:
# Note that smoothing will take a while (~10 minutes) with this high-resolution mask and the large kernel
mask_pr2_apod = smooth_mask(mask_pr2)
________________________________________________________________________________
[Memory] Calling __main__--Users-dlenz-projects-cibxphi-notebooks-__ipython-input__.smooth_mask...
smooth_mask(array([ 1., ...,  1.]))
Sigma is 50.959308 arcmin (0.014823 rad) 
-> fwhm is 120.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
____________________________________________________smooth_mask - 451.1s, 7.5min
In [16]:
hp.mollview(mask_pr2_apod)

Powerspectra

After the initial inspection, we compute the angular power spectra $C_l^{\Phi\Phi}$ for the lensing map. Since we want to understand the transformations, this concistency check is well-suited.

There are two different formats available The $\hat{\kappa}_{lm}$ directly contain the lensing convergence, and the noise and signal+noise power spectrum of $\hat\kappa$ is also available.