TE/ET fix (#71)

* Quick fix for TE/ET cross-frequencies

* Quick fix for TE/ET cross-frequencies (Hidde's solution)

* fixing typos

* Adding comments to the code

* fixing typo in comment

* Neater condition in TE/ET for mflike._get_power_spectra()

---------

Co-authored-by: sgiardie <grdsrn1@unife.it>

mflike/mflike.py

@@ -286,7 +286,9 @@ def get_sacc_names(pol, exp_1, exp_2):
                     # will all be sampled at the same ells.
                     self.l_bpws = ws.values
 
-                # Symmetrize if needed.
+                # Symmetrize if needed. If symmetrize = True, the "ET" polarization
+                # is eliminated by the polarization list and the TE spectrum becomes
+                # (TE + ET)/2. The associated spec_meta dict will have "hasYX_xsp": False
                 if (pol in ["TE", "ET"]) and symm:
                     pol2 = pol[::-1]
                     pols.remove(pol2)
@@ -315,8 +317,8 @@ def get_sacc_names(pol, exp_1, exp_2):
                     {
                         "ids": (index_sofar + np.arange(cls.size, dtype=int)),
                         "pol": ppol_dict[pol],
-                        "hasYX_xsp": pol
-                        in ["ET", "BE", "BT"],  # This is necessary for handling symmetrization
+                        "hasYX_xsp": pol        # this flag is true for pol = ET, BE, BT
+                        in ["ET", "BE", "BT"],  
                         "t1": exp_1,
                         "t2": exp_2,
                         "leff": ls,  #
@@ -357,7 +359,13 @@ def _get_power_spectra(self, cl, **params_values_nocosmo):
             p = m["pol"]
             i = m["ids"]
             w = m["bpw"].weight.T
-            clt = w @ DlsObs[p, m["t1"], m["t2"]]
+            # If symmetrize = False, the (ET, exp1, exp2) spectrum 
+            # will have the flag m["hasYX_xsp"] = True. 
+            # In this case, the power spectrum
+            # is computed as DlsObs["te", m["t2"], m["t1"]], to associate
+            # T --> exp2, E --> exp1
+            dls_obs = DlsObs[p, m["t2"], m["t1"]] if m["hasYX_xsp"] else DlsObs[p, m["t1"], m["t2"]]
+            clt = w @ dls_obs
             ps_vec[i] = clt
 
         return ps_vec

mflike/theoryforge.py

@@ -176,12 +176,16 @@ def get_modified_theory(self, Dls, **params):
             if p in ["tt", "ee", "bb"]:
                 dls_dict[p, m["t1"], m["t2"]] = cmbfg_dict[p, m["t1"], m["t2"]]
             else:  # ['te','tb','eb']
-                if m["hasYX_xsp"]:  # not symmetrizing
-                    dls_dict[p, m["t1"], m["t2"]] = cmbfg_dict[p, m["t2"], m["t1"]]
-                else:
+                if m["hasYX_xsp"]:  # case with symmetrize = False and ET/BT/BE spectra
+                    dls_dict[p, m["t2"], m["t1"]] = cmbfg_dict[p, m["t2"], m["t1"]]
+                else: # case of TE/TB/EB spectra, or symmetrize = True
                     dls_dict[p, m["t1"], m["t2"]] = cmbfg_dict[p, m["t1"], m["t2"]]
 
-                if self.defaults_cuts["symmetrize"]:  # we average TE and ET (as we do for data)
+                # if symmetrize = True, dls_dict has already been set 
+                # equal to cmbfg_dict[p, m["t1"], m["t2"]
+                # now we add cmbfg_dict[p, m["t2"], m["t1"] and we average them
+                # as we do for our data
+                if self.defaults_cuts["symmetrize"]:  
                     dls_dict[p, m["t1"], m["t2"]] += cmbfg_dict[p, m["t2"], m["t1"]]
                     dls_dict[p, m["t1"], m["t2"]] *= 0.5