
%load_ext autoreload
%autoreload 2


#load modules
import numpy as np, glob, sys

show_plot = True
if show_plot:
    from pylab import *
    #plot(); show()
    rcParams['figure.dpi'] = 150
    rcParams['font.family'] = 'serif'
    rcParams['figure.facecolor'] = 'white'
    #plot(); show()


fname = 'publish/bandpowers_and_cov/spt_100d_dl_yy_bandpowers_and_cov_all_estimators.npy'
spt_100d_yy_dict = np.load(fname, allow_pickle = True).item()
ilc_key_arr = [('ymv', 'ymv'), ('ycibfree', 'ycibfree'), ('ycibfree', 'ymv')]
el_eff = spt_100d_yy_dict['el_eff']
el_min = spt_100d_yy_dict['el_min']
el_max = spt_100d_yy_dict['el_max']
bpwf = spt_100d_yy_dict['bpwf']
print( '\n100d SPT yy results.\n\nEstimators are: ', list(spt_100d_yy_dict.keys()) )

#data
for m1m2 in ilc_key_arr:
    m1, m2 = m1m2
    print('\n\tEstimator = %sx%s' %(m1, m2))
    dl_yy = spt_100d_yy_dict[m1m2]['bandpower'] #Dl_yy [1e12]
    dl_yy_err = spt_100d_yy_dict[m1m2]['bandpower_error'] #Dl_yy_err [1e12] 
    dl_yy_cov = spt_100d_yy_dict[m1m2]['cov'] #Dl_yy_cov [1e12]
    print('\t\t', 'el_eff =', el_eff)
    print('\t\t', 'dl_yy =', dl_yy)
    print('\t\t', 'dl_yy_err =', dl_yy_err)

100d SPT yy results.

Estimators are:  ['el_eff', 'bpwf', 'el_min', 'el_max', 'delta_el', ('ymv', 'ymv'), ('ycibfree', 'ycibfree'), ('ycibfree', 'ymv')]

	Estimator = ymvxymv
		 el_eff = [ 750. 1250. 1750. 2250. 2750. 3250. 3750. 4250. 4750.]
		 dl_yy = [0.3874191  0.41697756 0.40722179 0.57300277 0.67155125 0.64040909
 0.59087548 0.654273   0.63533072]
		 dl_yy_err = [0.20979787 0.14851846 0.1269094  0.11088423 0.09341208 0.08276539
 0.08479397 0.08692268 0.09399785]

	Estimator = ycibfreexycibfree
		 el_eff = [ 750. 1250. 1750. 2250. 2750. 3250. 3750. 4250. 4750.]
		 dl_yy = [0.56661888 0.50512548 0.46756056 0.54718716 0.66731365 0.62197424
 0.58014743 0.65117204 0.62898178]
		 dl_yy_err = [0.21750301 0.14848823 0.13705187 0.11486615 0.09548166 0.08249556
 0.08623394 0.08811555 0.09478625]

	Estimator = ycibfreexymv
		 el_eff = [ 750. 1250. 1750. 2250. 2750. 3250. 3750. 4250. 4750.]
		 dl_yy = [0.48262144 0.46112836 0.42705692 0.54636736 0.66819478 0.6208671
 0.58355304 0.65307235 0.6313517 ]
		 dl_yy_err = [0.21164065 0.1456658  0.12998236 0.1109142  0.09293964 0.08177848
 0.08486734 0.08702453 0.09408163]


el_eff = spt_100d_yy_dict['el_eff']
el_min = spt_100d_yy_dict['el_min']
el_max = spt_100d_yy_dict['el_max']
bpwf = spt_100d_yy_dict['bpwf']

el_unbinned = np.arange( el_min, el_max + 1)
clf()
fsval = 14
color_arr = [cm.RdYlBu_r(int(d)) for d in np.linspace(0, 255, len(bpwf))]
for b in range(len(bpwf)):
    plot( el_unbinned, bpwf[b], color = color_arr[b])
xlabel(r'Multipole $\ell$', fontsize = fsval)
ylabel(r'BPWF $w_{\ell}$', fontsize = fsval)
show();

# Compare with FLAMINGO simulations
#flamingo_fname = 'data/flamingo_LS8_fgas_8sigma.txt'
#flamingo_el, flamingo_dl = np.loadtxt(flamingo_fname, unpack = True, delimiter = ',')
import pickle as pk
flamingo_fname = 'data/FLAMINGO_tSZ_power_spectra/LS8_fgas-8sig.pickle'
flamingo_data = np.asarray( pk.load( open(flamingo_fname,'rb') ) ).T
flamingo_el, flamingo_cl = flamingo_data[0], flamingo_data[1] * 1e12

#interpolate to get the power in the unbinned \ells.
flamingo_cl_unbinned = np.interp( el_unbinned, flamingo_el, flamingo_cl) 

#bin theory using BPWF
flamingo_cl_binned = np.dot( bpwf, flamingo_cl_unbinned )


clf()
reqd_m1m2 = ('ycibfree', 'ycibfree')
dl_fac = el_eff * (el_eff+1)/2/np.pi
ax = subplot(111)
fsval = 14
#data
data_dl_yy = spt_100d_yy_dict[reqd_m1m2]['bandpower'] #Dl_yy [1e12]
data_dl_yy_err = spt_100d_yy_dict[reqd_m1m2]['bandpower_error'] #Dl_yy_err [1e12] 
errorbar(el_eff, data_dl_yy, 
         yerr = data_dl_yy_err, 
         color = 'tab:green', 
         marker = 'o', 
         ms = 4., 
         capsize = 0,
         label = r'Data', 
         ls = 'None', 
         zorder = 1000,
        )

#FLAMINGO
dl_fac_unbinned = el_unbinned * (el_unbinned+1)/2/np.pi
flamingo_name = r'LS8$+f_{\rm gas}-8\sigma$'
plot(el_unbinned, dl_fac_unbinned * flamingo_cl_unbinned, color = 'silver', label = r'FLAMINGO (%s): Unbinned' %(flamingo_name), ls = '-', lw = 2.)
plot(el_eff, dl_fac * flamingo_cl_binned, color = 'black', label = r'FLAMINGO (%s): Binned' %(flamingo_name), ls = (4, (4, 0.5)), lw = 1.)
xlim(10., 5000.); ylim(0.1, 1.2)
legend(loc = 1, fontsize = fsval - 3)
xlabel(r'Multipole $\ell$', fontsize = fsval)
ylabel(r'$\ell(\ell+1)/2\pi C_{\ell}^{yy}$ [$10^{12}$]', fontsize = fsval)
show()

/var/folders/08/mkcy0rls72j2q09krs7qx2jw0000gn/T/ipykernel_11179/574348054.py:21: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  flamingo_data = np.asarray( pk.load( open(flamingo_fname,'rb') ) ).T


map_folder = 'publish/maps_and_masks/'
spt_100d_y_mask_dict_fname = '%s/masks.npy' %(map_folder)
spt_100d_y_mask_dict = np.load( spt_100d_y_mask_dict_fname, allow_pickle = True ).item()

spt_100d_y_boundary_mask = spt_100d_y_mask_dict['apod_mask']
spt_100d_y_boundary_plus_source_mask = spt_100d_y_mask_dict['mask']

if show_plot:
    imshow(spt_100d_y_boundary_mask); colorbar(); 
    title(r'SPT + {\it Herschel}-SPIRE: Boundary mask', fontsize = 10 ); 
    show()
    
    imshow(spt_100d_y_boundary_plus_source_mask); colorbar(); 
    title(r'Boundary + Point source ($S_{150} \ge $ 2 mJy) mask ', fontsize = 10 ); 
    show()


from astropy.io import fits
map_folder = 'publish/maps_and_masks/'
ilc_name_dict = {'mv': r'MV', 'cibfree': r'CIB-min'}
data_splits_name_dict = {'full': r'Full', 'half1': r'Half1', 'half2': r'Half2'}

reqd_ilc_name_arr = ['mv', 'cibfree']

for reqd_ilc_name in reqd_ilc_name_arr:
    reqd_full_half_iden = 'full'

    assert reqd_ilc_name in ilc_name_dict #['mv', 'cibfree']
    assert reqd_full_half_iden in data_splits_name_dict #['full', 'half1', 'half2']

    #read map
    spt_100d_y_map_fname = '%s/y%s_%s_proj0.fits' %(map_folder, reqd_ilc_name, reqd_full_half_iden)
    spt_100d_y_map_rec = fits.open( spt_100d_y_map_fname )
    spt_100d_y_map = spt_100d_y_map_rec[1].data

    #mapspecs
    ra0, dec0 = spt_100d_y_map_rec[1].header['ALPHA0'], spt_100d_y_map_rec[1].header['DELTA0']
    map_resol_degrees = spt_100d_y_map_rec[1].header['CDELT1']
    #map_resol_arcmins = map_resol_degrees * 60.
    print('Map centres: (ra, dec) = (%g, %g) degrees' %(ra0, dec0))
    ny, nx = spt_100d_y_map.shape
    boxsize_degrees = ny * map_resol_degrees

    #get masks
    spt_100d_y_mask_dict_fname = '%s/masks.npy' %(map_folder)
    spt_100d_y_mask_dict = np.load( spt_100d_y_mask_dict_fname, allow_pickle = True ).item()
    spt_100d_y_boundary_mask = spt_100d_y_mask_dict['apod_mask']


    if show_plot:
        clf()
        fsval = 14
        x1, x2 = ra0-boxsize_degrees/2, ra0+boxsize_degrees/2
        y1, y2 = dec0-boxsize_degrees/2, dec0+boxsize_degrees/2
        extent = [x1, x2, y1, y2]
        spt_100d_y_map[spt_100d_y_boundary_mask == 0] = None
        cmap = cm.RdYlBu_r
        vmin, vmax = -10., 15.
        scalefac = 1e6
        imshow(spt_100d_y_map * scalefac, cmap = cmap, extent = extent, vmin = vmin, vmax = vmax); 
        cbar = colorbar(pad = 0.01); 
        cbar.set_label(r'Compton-$y$: %s: %s [$10^{%g}$]' 
                       %(ilc_name_dict[reqd_ilc_name], data_splits_name_dict[reqd_full_half_iden], np.log10(scalefac)),
                        fontsize = fsval-1, )
        xlabel(r'R.A. [degrees]', fontsize = fsval)
        ylabel(r'Declination [degrees]', fontsize = fsval)
        show()

Map centres: (ra, dec) = (352.5, -55) degrees

Map centres: (ra, dec) = (352.5, -55) degrees


def get_lxly(flatskymapparams):

    """
    returns lx, ly based on the flatskymap parameters
    input:
    flatskymyapparams = [nx, ny, dx] where ny, nx = flatskymap.shape; and reso_am is the pixel resolution in arcminutes.
    for example: [100, 100, 0.5] is a 50' x 50' flatskymap that has dimensions 100 x 100 with dx = dy = 0.5 arcminutes.

    output:
    lx, ly
    """

    nx, ny, reso_am = flatskymapparams
    reso_rad = np.radians(reso_am/60.)

    lx, ly = np.meshgrid( np.fft.fftfreq( nx, reso_rad ), np.fft.fftfreq( ny, reso_rad ) )
    lx *= 2* np.pi
    ly *= 2* np.pi

    return lx, ly

#read beam and filters
beam_filter_fname = 'publish/bandpowers_and_cov/beam_filter.npy'
beam_filter_dict = np.load( beam_filter_fname, allow_pickle=True ).item()
bl_eff_2D = beam_filter_dict['beam_2D_effective']
filter_2D = beam_filter_dict['filter_2D']

#get lx, ly#, and ell_grid
ny, nx = bl_eff_2D.shape
reso_am = 0.5
mapparams = [ny, nx, reso_am]
lx, ly = get_lxly(mapparams)
##ell_grid = np.str( lx**2 + ly**2 )
print(beam_filter_dict.keys())

#plot
axlim = 7000
extent_val = [np.min(lx), np.max(lx), np.min(ly), np.max(ly)]
fsval = 14
clf()
imshow( np.fft.fftshift(bl_eff_2D), extent = extent_val); title(r'Effective beam'); 
colorbar()
xlabel(r'$\ell_{x}$', fontsize = fsval); ylabel(r'$\ell_{y}$', fontsize = fsval)
xlim(-axlim, axlim); ylim(-axlim, axlim)
show()

clf()
imshow( np.fft.fftshift(filter_2D), extent = extent_val); title(r'Filter'); 
colorbar()
xlabel(r'$\ell_{x}$', fontsize = fsval); ylabel(r'$\ell_{y}$', fontsize = fsval)
xlim(-axlim, axlim); ylim(-axlim, axlim)
show()

dict_keys(['beam_2D_effective', 'filter_2D'])

Read SPT 100d $yy$ data products¶

¶

Bandpowers, covariance matrices, and Bandpower window functions (BPWF)¶

¶

¶

Let us now plot the BPWF¶

The dimension of BPWF is length(el_unbinned) x total_el_bins¶

In our case: $\ell_{\rm min}$ = 500; $\ell_{\rm max}$ = 5000; $\Delta_\ell$ = 500.¶

Thus, el_unbinned = np.arange(lmin, lmax+1); total_el_bins = 9; BPWF is 5001 x 9 matrix¶

The bandpowers are binned with equal weights in $C_{\ell}$ space.¶

As a result, BPWF is simply 1/$\Delta_\ell$ (=1/500 = 0.002) in each bin.¶

Comparing theory spectrum with bandpowers¶

In case you have some theory spectrum and would like to compare that to the measured bandpowers, you must use the BPWF.¶

You can do this by¶

Below, we will one show such example using a theory curve from FLAMINGO simulations. Note that this is not a fit! Just an example.¶

The FLAMINGO bandpowers come from Fig. 6 of McCarthy et al. 2024. Thanks to Ian McCarthy for providing these.¶

¶

¶

Masks¶

We provide both the boundary mask and the source mask.¶

Source mask is used to mask all sources with flux $S_{150} \ge 2\ {\rm mJy}$.¶

(More details can be found in the paper.)¶

¶

¶

Maps (Sanson-Flamsteed projection).¶

We release two sets of LC maps: (A) MV and (B) CIB-minimized.¶

If you need CMB-free and/or Radio-minimized, please write Srini. These maps exist but have not been fully quantified, and hence not released publicly at the moment.¶

For each LC type, we release three sets of maps:¶

Full coadd, that includes all the observations.¶

Coadd of half1: Includes 50 per cent of the observations.¶

Coadd of half2: Includdes the other 50 per cent of the observations.¶

Cross-spectrum of B and C will help to remove the noise bias.¶

Note that in the paper, we perform crosses of 100 such bundles and just half1 and half2.¶

Important: Half1 x Half2 will result in slighlty different high-$\ell$ spectra compared to the 100 bundle-crosses due to detector time constants that are not deconvolved.¶

¶

Beams and filters¶

We release the effective beam for the LC maps and the filters. The below code will allow you to read/plot them.¶

If you want to calculate the power spectrum of one of the maps, the beam and the filter must be deconvolved to get an unbiased estimate.¶

(Please note again that the bandpowers reported in the paper are not from these full coadd or half-split maps. They are from the cross-spectra of 100 different bundles.)¶