```
import healpy as hp
import numpy as np
import astropy.units as u
hp.disable_warnings()
```

If you are analyzing a map from an instrument with a specific beam width, you can correct the power spectrum by the smoothing factor caused by that beam and obtain a better approximation of the power spectrum of the orignal sky.

```
= hp.read_map(
m, h "https://portal.nersc.gov/project/cmb/so_pysm_models_data/equatorial/dust_T_ns512.fits", h=True
)
```

`min=0, max=1000, title="Dust map", unit="uK_RJ") hp.mollview(m, `

`= hp.anafast(m) cl `

In this case we assume that the dust map from PySM is the true sky, then we apply a smoothing caused by the beam

`= 30 * u.arcmin beam `

`= hp.smoothing(m, fwhm=beam.to_value(u.radian)) m_smoothed `

`= hp.anafast(m_smoothed) cl_smoothed `

We can get the transfer function of the beam, generally referred as \(B_\ell\):

`= hp.gauss_beam(fwhm=beam.to_value(u.radian), lmax=len(cl)-1) bl `

`import matplotlib.pyplot as plt`

```
plt.loglog(bl)"Beam window function")
plt.title("$\ell$"); plt.xlabel(
```

`min=0, max=1000) hp.mollview(m_smoothed, `

We can recover the input \(C_\ell\) as $C_^{input} = $:

```
"seaborn-talk")
plt.style.use(="cl")
plt.loglog(cl, label="cl smoothed")
plt.plot(cl_smoothed, label/bl**2, label="cl smoothed corrected")
plt.plot(cl_smoothed1, len(cl)+100)
plt.xlim(1100, color="gray", ls="--", label="$\ell=1100$");
plt.axvline(
plt.legend(); plt.grid()
```

However, once the smoothed \(C_\ell\) reach machine precision, there is no more signal left, the beam deconvolution then causes extra noise. We need to limit our analysis to a range of \(\ell\) before the numerical error dominates, in this case for example, \(\ell=1100\).