So I've spent the last few days trying to work out how to use the code supplied in plc_2.0 (COM_Likelihood_Code...) in the Planck 2015 Likelihood code with the plik data to construct a data covariance matrix. This is including the off-diagonal components due to multipole mixing as well.

Does anyone have any idea where to start? I've been looking at cliklike.F90 in the CosmoMC code (most recent version) as well, but there doesn't seem to be any clear links. I've only seen covariance codes set up for the foreground routines, in the cmb_only component of the code and as I say I don't really understand how to use it.

To be clear: trying to create a covariance matrix between the Cls including binned-ls and the mixing between different multipoles.

I've been looking at the clik header files and the computational aspects included in clik_compute but I cannot find the routines where the covariance is calculated??

Not sure this is still useful but the Cl-Cl' covariance can be accessed by reading in the criterion_gauss_mat file that lives in the plik/plik_dx11dr2_HM_v18_TT.clik/clik/lkl_0/ folder. In python you can do this for example by using astropy.io.fits.getdata.

This will give you an array of floats that can be reshaped into a 765x765 covariance matrix - on each side the first 136 entries are the bandpower bins for 100x100 GHz, the next 199 are 143x143, then 215 for 143x217, and finally another 215 for 217x217 GHz.

This matrix has the signal (cosmic variance) and noise contributions included as described in the Planck 2015 likelihood paper. Note that the dependence of the CMB cosmic variance contribution on the actual cosmological parameters is ignored in the Planck likelihood, meaning the covariance matrix is fixed rather than having to be updated as the values of cosmological params change.

What exactly are you trying to do? Planck made a CMB-only spectrum and covariance marginalized over foregrounds and other nuisance parameters that they called 'plik_lite'. Are you after something like that?