Estimating RMS maps by smoothing and differencing

[lukashergt]

Premise

A bunch of our frequency maps are fairly old and don't come with their own RMS maps. Typically the paper will quote a zero-level-offset error (i.e. the error in the monopole) and a scale error (i.e. the error in the gain). These errors should ideally not go into the RMS map used with Commander3, since Commander3 can handle monopole and gain estimation on its own. Also the confusion limit should not go into the RMS map, instead it should be considered to be part of the signal.

Of course we cannot vary all the monopoles and gains simultaneously, we need to anchor them somehow. In BeyondPlanck the Haslam map was chosen for that purpose, i.e. its monopole and gain were fixed. In this project, we will try and use the lowest absolutely calibrated ARCADE2 band as our anchor.

However, we are left with the question of how to generate an appropriate RMS map that we can supply to Commander3, given that the only errors quoted in older papers are not the values that we are looking for.

Compute RMS from smoothed map

Here, I am presenting one option of estimating the RMS value:

• smooth the original map with a large beam
• create a galaxy mask by discarding the upper quartile of the map
• create a difference map between smoothed and original and apply the galaxy mask
• compute the root mean square of the resulting masked difference map

This is the above recipe turned into a python function that I applied to various of our synch maps:

def get_rms(m, mask, fwhm_deg, smoothing_order=3):
    """Compute RMS map.
    
    Parameters
    ----------
        m: array
            original map
        mask: array
            mask of original map
        fwhm_deg: float
            the beam full width at half maximum (fwhm) in degrees of the original map
        smoothing_order: float
            factor to be multiplied with the fwhm of the original map to create the beam used for smoothing
            
    Returns
    -------
        rms: array
            RMS map matching the shape of m and mask
    """
    s = np.where(mask, m, np.nanmedian(m))                       # inpaint masked area with median value before smoothing to reduce edge effects
    s = hp.smoothing(s, fwhm=deg2rad(fwhm_deg*smoothing_order))  # smooth map with some multiple of the beam fwhm
    s = np.where(mask, s, np.nan)                                # reintroduce the masked area
    d = np.abs(m-s)                                              # take the difference between smoothed map and original
    mask_galaxy = np.where(s<np.nanquantile(s, 0.75), 1, 0)      # create a galaxy mask
    d = np.where(mask_galaxy, d, np.nan)                         # apply galaxy mask to difference map
    rms = np.sqrt(np.nanmean(d**2))                              # compute RMS outside the galaxy
    return rms

Drawbacks

The estimation via the smoothing approach is very dependent on the choice of smoothing scale. I went with a smoothing order of 3, i.e. I applied a smoothing beam three times the beam size of the corresponding experiment. This appeared to produce values in an expectant ballpark, however, it will be good to revisit these values once we have reliable residual maps.

Note, that in my approach I was focusing on finding something that can ideally be applied to all experiments, rather than optimising for one specific experiment. Hence, it would certainly be good if everyone could critically look into how this applies to their respective maps.