correct bandpass construction in the presence of shifts

mflike/MFLike.yaml

@@ -71,6 +71,7 @@ data:
 # (if bandpass in sacc is a single freq, no band integration)
 # the default is to read bandpasses from file, to build top-hat uncomment the
 # parameters of the block!
+# Bandpass has to be in RJ units
 top_hat_band:
 #  nsteps: 1
 #  bandwidth: 0

mflike/tests/test_mflike.py

@@ -158,8 +158,8 @@ def get_model(nsteps, bandwidth):
 
         # chi2 reference results for the different models and different bandshifts
         chi2s = {
-            "model1": [2265.168, 3742.528, 43753.810],
-            "model2": [2345.508, 4224.293, 41764.366],
+            "model1": [2265.168, 3742.528, 43753.809],
+            "model2": [2345.508, 4458.148, 47533.477],
         }
 
         for model, chi2 in chi2s.items():

mflike/theoryforge.py

@@ -5,7 +5,10 @@
 from cobaya.log import LoggedError
 
 
-# Converts from cmb units to brightness. Numerical factors not included,
+# Converts from cmb units to brightness.
+# There is an additional nu**2 factor to convert the bandpasses
+# from RJ to brightness units.
+# Numerical factors not included,
 # it needs proper normalization when used.
 def _cmb2bb(nu):
     # NB: numerical factors not included
@@ -78,6 +81,7 @@ def __init__(self, mflike=None):
 
     # Takes care of the bandpass construction. It returns a list of nu-transmittance
     # for each frequency or an array with the effective freqs.
+    # bandpasses saved in the sacc file have to be in RJ units
     def _bandpass_construction(self, **params):
         data_are_monofreq = False
         self.bandint_freqs = []
@@ -93,15 +97,17 @@ def _bandpass_construction(self, **params):
                 if self.bandint_nsteps > 1:
                     bandlow = fr * (1 - self.bandint_width[iexp] * 0.5)
                     bandhigh = fr * (1 + self.bandint_width[iexp] * 0.5)
-                    nubtrue = np.linspace(bandlow, bandhigh, self.bandint_nsteps, dtype=float)
+                    #nubtrue = np.linspace(bandlow, bandhigh, self.bandint_nsteps, dtype=float)
                     nub = np.linspace(
                         bandlow + params[bandpar],
                         bandhigh + params[bandpar],
                         self.bandint_nsteps,
                         dtype=float,
                     )
                     tranb = _cmb2bb(nub)
-                    tranb_norm = np.trapz(_cmb2bb(nubtrue), nubtrue)
+                    # normalization integral to be evaluated at the shifted freqs 
+                    # in order to have cmb component calibrated to 1
+                    tranb_norm = np.trapz(_cmb2bb(nub), nub)
                     self.bandint_freqs.append([nub, tranb / tranb_norm])
                 # in case we don't want to do band integration, e.g. when we have multifreq bandpass in sacc file
                 if self.bandint_nsteps == 1:
@@ -116,14 +122,15 @@ def _bandpass_construction(self, **params):
                     data_are_monofreq = True
                     self.bandint_freqs.append(nub[0])
                 else:
-                    trans_norm = np.trapz(bp * _cmb2bb(nu_ghz), nu_ghz)
+                    trans_norm = np.trapz(bp * _cmb2bb(nub), nub)
                     trans = bp / trans_norm * _cmb2bb(nub)
                     self.bandint_freqs.append([nub, trans])
 
         # fgspectra can't mix monofrequency with [nu, bp]. If one channel is mono-frequency then we
         # assume all of them and pass to fgspectra an array (not list!!) of frequencies
         if data_are_monofreq:
             self.bandint_freqs = np.asarray(self.bandint_freqs)
+            self.log.info("bandpass is delta function, no band integration performed")
 
     def get_modified_theory(self, Dls, **params):
         fg_params = {k: params[k] for k in self.expected_params_fg}