simonsobs / LAT_MFLike

Multifrequency Likelihood for SO Large Aperture Telescope
https://lat-mflike.readthedocs.io
Other
4 stars 12 forks source link

Implementing beam chromaticity #76

Closed sgiardie closed 1 month ago

sgiardie commented 5 months ago

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 commented 5 months ago

@sgiardie Can you try this

sudo ln -s /opt/homebrew/bin/gfortran-11 /usr/local/bin/gfortran
codecov-commenter commented 5 months ago

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 :warning:
Additional details and impacted files [![Impacted file tree graph](https://app.codecov.io/gh/simonsobs/LAT_MFLike/pull/76/graphs/tree.svg?width=650&height=150&src=pr&token=qrrVcbNCs5&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simonsobs)](https://app.codecov.io/gh/simonsobs/LAT_MFLike/pull/76?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simonsobs) ```diff @@ 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](https://app.codecov.io/gh/simonsobs/LAT_MFLike/pull/76?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simonsobs) | Coverage Δ | | |---|---|---| | [mflike/mflike.py](https://app.codecov.io/gh/simonsobs/LAT_MFLike/pull/76?src=pr&el=tree&filepath=mflike%2Fmflike.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simonsobs#diff-bWZsaWtlL21mbGlrZS5weQ==) | `83.08% <100.00%> (+0.32%)` | :arrow_up: | | [mflike/foreground.py](https://app.codecov.io/gh/simonsobs/LAT_MFLike/pull/76?src=pr&el=tree&filepath=mflike%2Fforeground.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simonsobs#diff-bWZsaWtlL2ZvcmVncm91bmQucHk=) | `88.08% <87.50%> (-1.17%)` | :arrow_down: |
sgiardie commented 5 months ago

@sgiardie Can you try this

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

That worked, thanks!

sgiardie commented 5 months ago

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.

cmbant commented 4 months ago

healpy is loaded dynamically, so I think can be an optional dependency (then can work on Windows as long as those functions not used).

sgiardie commented 1 month ago

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.