# Investigate broken implementation of bandpass unit conversion in PySM 3

cosmology
Published

August 24, 2020

import pysm3
import numpy as np
import healpy as hp
from astropy.tests.helper import assert_quantity_allclose
import pysm3.units as u

Investigation of the broken implementation of the bandpass unit conversion function in PySM 3 and check of its impact on existing Simons Observatory and CMB-S4 simulated datasets.

See the related issue and the pull request with the fix. This affects all PySM 3 versions <3.2.2, it will be fixed in 3.3.0.

$= g() I()d$

If we consider emission of the CMB in $$K_{CMB}$$ units, $$I_{CMB}$$ is not a function of $$\nu$$:

${CMB}[K_{CMB}] = g() I{CMB}[K_{CMB}] d= I_{CMB}[K_{CMB}] g() d= I_{CMB}[K_{CMB}]$

${CMB}[K_{RJ}] = I{CMB}[K_{CMB}] g() C_{K_{CMB}}^{K_{RJ}}() d$

However we assume that the bandpass of the instrument $$g(\nu)$$ is given in power units, i.e. $$Jy/sr$$, the python normalize_weights functions does:

$g [K_{RJ}/K_{RJ}] =$

CMB bandpass integrated converted to $$K_{RJ}$$.

${CMB}[K_{RJ}] = I{CMB}[K_{CMB}] g()[K_{RJ}/K_{RJ}] C_{K_{CMB}}^{K_{RJ}}() d= I_{CMB}[K_{CMB}] C_{K_{CMB}}^{K_{RJ}}() d$

You can think the last equation as getting a value in $$K_{CMB}$$, turn it into $$K_{RJ}$$ and then to $$Jy~sr^{-1}$$, and do the relative weighting then.

However we assume that the bandpass of the instrument $$g(\nu)$$ is given in power units, i.e. $$Jy/sr$$, the python normalize_weights functions does:

$g [K_{RJ}/K_{RJ}] =$

## From $$K_{CMB}$$ to $$K_{RJ}$$:

If the output is requested in $$K_{CMB}$$, we basically have to undo that steps, so:

${CMB}[K_{RJ}] = I{CMB}[K_{CMB}] C_{K_{CMB}}^{K_{RJ}}() d$

### fixed_bandpass_unit_conversion

${CMB}[K_{CMB}] = {CMB}[K_{RJ}] { C_{K_{CMB}}{Jy~sr{-1}}() g() d}$

### broken_bandpass_unit_conversion

${CMB}[K_{CMB}] = {CMB}[K_{RJ}] C_{K_{RJ}}^{K_{CMB}}() g_{RJ}() d() = {CMB}[K_{RJ}] C{K_{RJ}}^{K_{CMB}}() d()$

from pysm3.utils import normalize_weights
def broken_bandpass_unit_conversion(freqs, weights, output_unit):
"""Unit conversion from uK_RJ to output unit given a bandpass

Parameters
----------
freqs : astropy.units.Quantity
Frequency array in a unit compatible with GHz
"""
freqs = check_freq_input(freqs)
factors = (np.ones(len(freqs), dtype=np.float) * u.uK_RJ).to_value(
output_unit, equivalencies=u.cmb_equivalencies(freqs * u.GHz)
)
if len(freqs) > 1:
w = normalize_weights(freqs, weights)
factor = np.trapz(factors * w, freqs)
else:
factor = factors[0]
return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)
    nside = 32
freqs = np.array([250, 300, 350]) * u.GHz
weights = np.ones(len(freqs))
sky = pysm3.Sky(nside=nside, preset_strings=["c2"])
CMB_rj_int = sky.get_emission(freqs, weights)
CMB_thermo_int = CMB_rj_int*fixed_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
"pysm_2/lensed_cmb.fits", field=(0, 1), nside=nside, unit=u.uK_CMB
)
for pol in [0, 1]:
#assert_quantity_allclose(expected_map[pol], CMB_thermo_int[pol], rtol=1e-4)
print("Error: {:.2f} %".format(100-np.median((CMB_thermo_int[pol]/expected_map[pol]).value)*100))
Error: -0.00 %
Error: -0.00 %
from pysm3.utils import check_freq_input
def fixed_bandpass_unit_conversion(freqs, weights, output_unit):
"""Unit conversion from uK_RJ to output unit given a bandpass
Parameters
----------
freqs : astropy.units.Quantity
Frequency array in a unit compatible with GHz
"""
freqs = check_freq_input(freqs)
weights_to_rj = (weights * u.uK_RJ).to_value(
(u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
)
weights_to_out = (weights * output_unit).to_value(
(u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz)
)
if len(freqs) > 1:
factor = np.trapz(weights_to_rj, freqs)/np.trapz(weights_to_out, freqs)
else:
factor = (1. * u.uK_RJ).to_value(
output_unit, equivalencies=u.cmb_equivalencies(freqs * u.GHz)
)
return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)

## Unit conversion error for 20% tophat bandpasses

perc_error = {}
for center_freq in np.array([10, 20, 40, 70, 100, 150, 200, 250, 270, 300, 350, 400])*u.GHz:
freqs = np.linspace(center_freq*.8, center_freq*1.2, 10)
weights = np.ones(10)
rj_to_cmb = broken_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
print("{: <3.0f}, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(center_freq, rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))
10  GHz, broken: 1.003 uK_CMB / uK_RJ, fixed: 1.003 uK_CMB / uK_RJ, error 0.00%
20  GHz, broken: 1.011 uK_CMB / uK_RJ, fixed: 1.011 uK_CMB / uK_RJ, error 0.00%
40  GHz, broken: 1.045 uK_CMB / uK_RJ, fixed: 1.045 uK_CMB / uK_RJ, error 0.01%
70  GHz, broken: 1.143 uK_CMB / uK_RJ, fixed: 1.142 uK_CMB / uK_RJ, error 0.08%
100 GHz, broken: 1.310 uK_CMB / uK_RJ, fixed: 1.306 uK_CMB / uK_RJ, error 0.32%
150 GHz, broken: 1.808 uK_CMB / uK_RJ, fixed: 1.782 uK_CMB / uK_RJ, error 1.47%
200 GHz, broken: 2.770 uK_CMB / uK_RJ, fixed: 2.661 uK_CMB / uK_RJ, error 4.06%
250 GHz, broken: 4.627 uK_CMB / uK_RJ, fixed: 4.260 uK_CMB / uK_RJ, error 8.62%
270 GHz, broken: 5.800 uK_CMB / uK_RJ, fixed: 5.221 uK_CMB / uK_RJ, error 11.10%
300 GHz, broken: 8.297 uK_CMB / uK_RJ, fixed: 7.178 uK_CMB / uK_RJ, error 15.59%
350 GHz, broken: 15.749 uK_CMB / uK_RJ, fixed: 12.560 uK_CMB / uK_RJ, error 25.39%
400 GHz, broken: 31.306 uK_CMB / uK_RJ, fixed: 22.588 uK_CMB / uK_RJ, error 38.59%

## Unit conversion error for Simons Observatory channels

import h5py
SO_chs = h5py.File("../mapsims/mapsims/data/simonsobs_instrument_parameters_2020.06.h5", mode='r')
SO_chs["LT0_UHF1"].keys()
<KeysViewHDF5 ['bandpass_frequency_GHz', 'bandpass_weight']>
SO_chs["LT0_UHF1"].attrs.keys()
<KeysViewHDF5 ['band', 'band_id', 'center_frequency_GHz', 'fwhm_arcmin', 'noise_band_index', 'telescope', 'tube', 'tube_id']>
perc_error = {}
for ch in SO_chs.values():
freqs = ch["bandpass_frequency_GHz"] * u.GHz
weights = ch["bandpass_weight"]
rj_to_cmb = broken_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
print("{} {:4s}, {:6.2f} GHz, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(ch.attrs["telescope"], ch.attrs["band"], ch.attrs["center_frequency_GHz"], rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))
LA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
LA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
LA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
LA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
LA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
LA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
LA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
LA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
LA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
LA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
LA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
LA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
LA LF1 ,  25.70 GHz, broken: 1.018 uK_CMB / uK_RJ, fixed: 1.018 uK_CMB / uK_RJ, error 0.00%
LA LF2 ,  38.90 GHz, broken: 1.043 uK_CMB / uK_RJ, fixed: 1.043 uK_CMB / uK_RJ, error 0.01%
SA UHF1, 225.70 GHz, broken: 3.377 uK_CMB / uK_RJ, fixed: 3.293 uK_CMB / uK_RJ, error 2.54%
SA UHF2, 285.40 GHz, broken: 6.166 uK_CMB / uK_RJ, fixed: 5.990 uK_CMB / uK_RJ, error 2.93%
SA MFF1,  92.00 GHz, broken: 1.248 uK_CMB / uK_RJ, fixed: 1.247 uK_CMB / uK_RJ, error 0.12%
SA MFF2, 147.50 GHz, broken: 1.729 uK_CMB / uK_RJ, fixed: 1.721 uK_CMB / uK_RJ, error 0.49%
SA MFS1,  88.60 GHz, broken: 1.229 uK_CMB / uK_RJ, fixed: 1.228 uK_CMB / uK_RJ, error 0.11%
SA MFS2, 146.50 GHz, broken: 1.720 uK_CMB / uK_RJ, fixed: 1.711 uK_CMB / uK_RJ, error 0.54%
SA LF1 ,  25.70 GHz, broken: 1.018 uK_CMB / uK_RJ, fixed: 1.018 uK_CMB / uK_RJ, error 0.00%
SA LF2 ,  38.90 GHz, broken: 1.043 uK_CMB / uK_RJ, fixed: 1.043 uK_CMB / uK_RJ, error 0.01%

## Unit conversion error for CMB-S4 channels

S4_chs = h5py.File("../s4mapbasedsims/202006_foregrounds_extragalactic_cmb_tophat/cmbs4_tophat.h5", mode='r')
perc_error = {}
for ch in S4_chs.values():
freqs = ch["bandpass_frequency_GHz"] * u.GHz
weights = ch["bandpass_weight"]
rj_to_cmb = broken_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
fixed_rj_to_cmb = fixed_bandpass_unit_conversion(
freqs, weights, u.uK_CMB
)
perc_error[center_freq] = (rj_to_cmb - fixed_rj_to_cmb) / fixed_rj_to_cmb * 100
print("{} {:5s}, {:6.2f} GHz, broken: {:.3f}, fixed: {:.3f}, error {:.2f}%".format(
ch.attrs["telescope"], ch.attrs["band"], ch.attrs["center_frequency_GHz"], rj_to_cmb, fixed_rj_to_cmb, perc_error[center_freq]))
LAT HFL1 , 225.00 GHz, broken: 3.364 uK_CMB / uK_RJ, fixed: 3.275 uK_CMB / uK_RJ, error 2.71%
LAT HFL2 , 278.00 GHz, broken: 5.632 uK_CMB / uK_RJ, fixed: 5.523 uK_CMB / uK_RJ, error 1.98%
SAT HFS1 , 220.00 GHz, broken: 3.163 uK_CMB / uK_RJ, fixed: 3.109 uK_CMB / uK_RJ, error 1.71%
SAT HFS2 , 270.00 GHz, broken: 5.269 uK_CMB / uK_RJ, fixed: 5.099 uK_CMB / uK_RJ, error 3.34%
LAT LFL1 ,  27.00 GHz, broken: 1.019 uK_CMB / uK_RJ, fixed: 1.019 uK_CMB / uK_RJ, error 0.00%
LAT LFL2 ,  39.00 GHz, broken: 1.044 uK_CMB / uK_RJ, fixed: 1.044 uK_CMB / uK_RJ, error 0.01%
SAT LFS1 ,  30.00 GHz, broken: 1.024 uK_CMB / uK_RJ, fixed: 1.024 uK_CMB / uK_RJ, error 0.00%
SAT LFS2 ,  40.00 GHz, broken: 1.044 uK_CMB / uK_RJ, fixed: 1.044 uK_CMB / uK_RJ, error 0.01%
SAT MFHS1,  95.00 GHz, broken: 1.263 uK_CMB / uK_RJ, fixed: 1.262 uK_CMB / uK_RJ, error 0.10%
SAT MFHS2, 155.10 GHz, broken: 1.822 uK_CMB / uK_RJ, fixed: 1.813 uK_CMB / uK_RJ, error 0.51%
LAT MFL1 ,  93.00 GHz, broken: 1.262 uK_CMB / uK_RJ, fixed: 1.259 uK_CMB / uK_RJ, error 0.22%
LAT MFL2 , 145.00 GHz, broken: 1.707 uK_CMB / uK_RJ, fixed: 1.697 uK_CMB / uK_RJ, error 0.62%
SAT MFLS1,  85.00 GHz, broken: 1.207 uK_CMB / uK_RJ, fixed: 1.206 uK_CMB / uK_RJ, error 0.06%
SAT MFLS2, 145.10 GHz, broken: 1.696 uK_CMB / uK_RJ, fixed: 1.690 uK_CMB / uK_RJ, error 0.40%
LAT ULFL1,  20.00 GHz, broken: 1.011 uK_CMB / uK_RJ, fixed: 1.011 uK_CMB / uK_RJ, error 0.00%