Implementing beam chromaticity

[sgiardie]

This PR introduces the beam chromaticity effect in mflike. This means considering the frequency dependence of beams (both in temperature and polarization) $b^{T/P}_{\ell}(\nu)$ and computing the foreground SEDs by taking that into account: $ \frac{\int d\nu \frac{\partial B_{\nu}(\nu)}{\partial T} \tau(\nu) b^{T/P}_{\ell}(\nu) f^{T/P}_{SED}(\nu)}{\int d\nu \frac{\partial B_{\nu}(\nu)}{\partial T} \tau(\nu) b^{T/P}_{\ell}(\nu) }$ (where $\tau(\nu)$ is the bandpass and $\frac{\partial B_{\nu}(\nu)}{\partial T} $ is the derivative of the black body emission). Since this integral is computed in fgspectra, we modified also fgspectra/frequency.py to allow the transmission factor $\frac{\frac{\partial B_{\nu}(\nu)}{\partial T} \tau(\nu) b^{T/P}_{\ell}(\nu) }{\int d\nu \frac{\partial B_{\nu}(\nu)}{\partial T} \tau(\nu) b^{T/P}_{\ell}(\nu) }$ to have dimension (freqs, ells).

Since we need to differentiate between temperature and polarization beam, we now need to distinguish between the transmission factors in T/P. This also means having to introduce a new class in fgspectra/cross.py, FactorizedCrossSpectrumTE, to separately pass the different transmissions for T and E for the same cross spectrum.

There are three ways to provide the beams for each frequency array/experiment: computing them in the code as simple Gaussian beams, reading them from external file or reading them from the data sacc file (currently this would require modifications to sacc, since sacc does not allow an array(freqs, ells) to be passed as a beam tracer tracer.beam. If we cannot modify that, we should drop the option to read from sacc, or find another way to do that). This choice is done by selecting the appropriate options in the MFLike.py block beam_profile.

We also introduced a correction factor to apply to the beam profiles in the presence of bandpass shift: $b^{T/P}_{\ell}(\nu+\Delta \nu) \simeq b^{T/P}_{\ell}(\nu) b^{T/P, Gauss \, corr}_{\ell}(\nu, \Delta \nu)$. For simplicity, this is the correction factor expected for a Gaussian beam.

[xgarrido]

@sgiardie Can you try this

sudo ln -s /opt/homebrew/bin/gfortran-11 /usr/local/bin/gfortran

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 88.18182% with 13 lines in your changes missing coverage. Please review.

Project coverage is 85.45%. Comparing base (5afa861) to head (8867e5e).
Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
mflike/foreground.py	87.50%	13 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #76      +/-   ##
==========================================
+ Coverage   85.24%   85.45%   +0.21%     
==========================================
  Files           3        3              
  Lines         454      550      +96     
==========================================
+ Hits          387      470      +83     
- Misses         67       80      +13     

Files with missing lines	Coverage Δ
mflike/mflike.py	83.08% <100.00%> (+0.32%)	⬆️
mflike/foreground.py	88.08% <87.50%> (-1.17%)	⬇️

[sgiardie]

@sgiardie Can you try this

sudo ln -s /opt/homebrew/bin/gfortran-11 /usr/local/bin/gfortran

That worked, thanks!

[sgiardie]

We changed a bit the approach since the beginning of this PR, due to discussions with @AdriJD.
Now, when Gaussian beams are computed within the code, we don't automatically assume that the FWHM is that of a diffraction limited experiment anymore. The FWHM is computed as $FWHM(\nu) = FWHM(\nu_0) \left( \frac{\nu}{\nu_0} \right)^{-\alpha/2}$, so the user needs to provide the parameters $\nu_0$, $FWHM(\nu_0) $ and $\alpha$ through a dictionary in MFLike.yaml:

beam_profile:
Gaussian_beam:
LAT_93_s0:
FWHM_0: ...
nu_0: ...
alpha: ...
LAT_93_s2:
FWHM_0: ...
nu_0: ...
alpha: ...
LAT_145_s0:
FWHM_0: ...
nu_0: ...
alpha: ...

There is still the option to read precomputed chromatic beams from an external file, which is activated in MFLike.yaml in this way:

beam_profile:
Gaussian_beam: null/False
beam_from_file: filename

Finally, we changed the way we deal with bandpass shifts. Before, we were applying the Gaussian correction expected for $\Delta \nu$ to $b_{\ell}(\nu)$. Now we assume that $b_{\ell}(\nu + \Delta \nu) = b_{\ell (\nu/\nu_0)^{-\alpha/2}}(\nu_0 + \Delta \nu)$. The monochromatic beam $b_{\ell}(\nu_0 + \Delta \nu)$ is computed from planet beams measurements assuming the bandpass shift $\Delta \nu$. Then the scaling $b_{\ell}(\nu_0 + \Delta \nu) \rightarrow b_{\ell (\nu/\nu_0)^{-\alpha/2}}(\nu_0 + \Delta \nu)$ is applied within the code. Thus the user needs to provide an external yaml file with the following structure:

LAT_93_s0:
beams: {..., '-2.0': b_ell(nu_0 -2),
'-1.0': b_ell(nu_0 -1),
...
'5.0': b_ell(nu_0 + 5),
...}
nu_0: ...
alpha: ...
LAT_93_s2:
beams: {'-10.0': b_ell(nu_0 - 10), ...}
nu_0: ...
alpha: ...
LAT_145_s0:
beams: ...
nu_0: ...
alpha: ...
...

This yaml file should be saved in the cobaya data directory and will be read by theoryforge by providing its filename in MFLike.yaml:

beam_profile:
Bandpass_shifted_beams: bandpass_shifted_beams
Gaussian_beam: dict/False/null
beam_from_file: filename/False/null

In the absence of measurements of these $b_{\ell}(\nu_0 + \Delta \nu)$, I have added a piece of code in mflike_utils/utils.py to compute them as simple Gaussian beams for a diffraction limited experiment and $\Delta \nu \in [-20, 20]$. This code also saves the corresponding yaml file.
I have also added more docs to explain all of this.

[sgiardie]

Recent changes to this branch: I have aligned it to the current mflike master and allowed the choice of leaving the beams unchanged in case of bandpass shifts (keeping $b_{\ell}(\nu)$ instead of computing $b_{\ell}(\nu + \Delta \nu)$). This is done for simplicity, since it should be a second order effect.

[xgarrido]

I have made some kozmetix comments but since the current unit test still pass (and you have added a new test for beam chromaticity), I guess we are safe. I'm only a bit concerned by the mflike_utils/utils.py script (see comment below)

[xgarrido]

I thought the use of abs function was not going to work because it had to deal with numpy array and a scalar float but it actually works and it is as fast as np.abs

[sgiardie]

If you prefer I can change it to np.abs!

[xgarrido]

Looks good to me