Merge pull request #20 from mathias77515/Tom-dev2

Tom dev2

src/CMM/preset/preset_acquisition.py

@@ -3,6 +3,7 @@
 
 from lib.Qacquisition import *
 from lib.Qnoise import *
+import fgb.component_model as c
 
 
 class PresetAcquisition:
@@ -247,6 +248,27 @@ def get_preconditioner(self, A_qubic, A_ext, precond=True, thr=0):
         else:
             return None
 
+    def _get_scalar_acquisition_operator(self):
+        """
+        Function that will compute "scalar acquisition operatord" by applying the acquisition operators to a vector full of ones.
+        These scalar operators will be used to compute the resolutions in the case where we do not add convolutions during reconstruction.
+        """
+        ### Import the acquisition operators
+        acquisition_operators = self.preset_qubic.joint_out.qubic.H
+
+        ### Create the vector full of ones which will be used to compute the scalar operators
+        vector_ones = np.ones(acquisition_operators[0].shapein)
+
+        ### Apply each sub_operator on the vector
+        scalar_acquisition_operators = np.empty(
+            len(self.preset_qubic.joint_out.qubic.allnus)
+        )
+        for freq in range(len(self.preset_qubic.joint_out.qubic.allnus)):
+            scalar_acquisition_operators[freq] = np.mean(
+                acquisition_operators[freq](vector_ones)
+            )
+        return scalar_acquisition_operators
+
     def get_convolution(self):
         """Convolutions.
 
@@ -279,7 +301,7 @@ def get_convolution(self):
 
         # Check if convolution_out is True
         if self.preset_qubic.params_qubic["convolution_out"]:
-            swhm_mapmaking = np.sqrt(
+            fwhm_mapmaking = np.sqrt(
                 self.preset_qubic.joint_in.qubic.allfwhm**2
                 - np.min(self.preset_qubic.joint_in.qubic.allfwhm) ** 2
             )
@@ -299,7 +321,53 @@ def get_convolution(self):
             not self.preset_qubic.params_qubic["convolution_in"]
             and not self.preset_qubic.params_qubic["convolution_out"]
         ):
-            fwhm_rec = 0
+            scalar_acquisition_operators = self._get_scalar_acquisition_operator()
+            self.fwhm_reconstructed = np.zeros(len(self.preset_fg.components_model_out))
+            for comp in range(len(self.preset_fg.components_model_out)):
+                if self.preset_fg.components_name_out[comp] == "CMB":
+                    self.fwhm_reconstructed[comp] = np.sum(
+                        scalar_acquisition_operators * self.fwhm_tod
+                    ) / (np.sum(scalar_acquisition_operators))
+                if self.preset_fg.components_name_out[comp] == "Dust":
+                    f_dust = c.ModifiedBlackBody(
+                        nu0=self.preset_fg.params_foregrounds["Dust"]["nu0_d"],
+                        beta_d=self.preset_fg.params_foregrounds["Dust"]["beta_d_init"][
+                            0
+                        ],
+                    )
+                    self.fwhm_reconstructed[comp] = np.sum(
+                        scalar_acquisition_operators
+                        * f_dust.eval(self.preset_qubic.joint_out.qubic.allnus)
+                        * self.fwhm_tod
+                    ) / (
+                        np.sum(
+                            scalar_acquisition_operators
+                            * f_dust.eval(self.preset_qubic.joint_out.qubic.allnus)
+                        )
+                    )
+                if self.preset_fg.components_name_out[comp] == "Synchrotron":
+                    f_sync = c.PowerLaw(
+                        nu0=self.preset_fg.params_foregrounds["Synchrotron"]["nu0_s"],
+                        beta_pl=self.preset_fg.params_foregrounds["Synchrotron"][
+                            "beta_s_init"
+                        ][0],
+                    )
+                    self.fwhm_reconstructed[comp] = np.sum(
+                        scalar_acquisition_operators
+                        * f_sync.eval(self.preset_qubic.joint_out.qubic.allnus)
+                        * self.fwhm_tod
+                    ) / (
+                        np.sum(
+                            scalar_acquisition_operators
+                            * f_sync.eval(self.preset_qubic.joint_out.qubic.allnus)
+                        )
+                    )
+
+        elif (
+            not self.preset_qubic.params_qubic["convolution_in"]
+            and not self.preset_qubic.params_qubic["convolution_out"]
+        ):
+            self.fwhm_reconstructed = np.zeros(len(self.preset_fg.components_model_out))
 
         # Print the FWHM values
         self.preset_tools._print_message(f"FWHM for TOD making : {fwhm_tod}")

src/FMM/pipeline.py

@@ -6,6 +6,7 @@
 import numpy as np
 import qubic
 import yaml
+from scipy.optimize import minimize
 from pysimulators.interfaces.healpy import HealpixConvolutionGaussianOperator
 
 ### Local packages
@@ -155,7 +156,7 @@ def __init__(self, comm, file, params):
 
         ### Initial maps
         self.m_nu_in = self.get_input_map()
-        
+
         ### Define reconstructed and TOD operator
         self.get_H()
 
@@ -391,6 +392,25 @@ def get_dict(self, key="in"):
 
         return dict_qubic
 
+    def _get_scalar_acquisition_operator(self):
+        """
+        Function that will compute "scalar acquisition operatord" by applying the acquisition operators to a vector full of ones.
+        These scalar operators will be used to compute the resolutions in the case where we do not add convolutions during reconstruction.
+        """
+        ### Import the acquisition operators
+        acquisition_operators = self.joint.qubic.H
+
+        ### Create the vector full of ones which will be used to compute the scalar operators
+        vector_ones = np.ones(acquisition_operators[0].shapein)
+
+        ### Apply each sub_operator on the vector
+        scalar_acquisition_operators = np.empty(len(self.joint.qubic.allnus))
+        for freq in range(len(self.joint.qubic.allnus)):
+            scalar_acquisition_operators[freq] = np.mean(
+                acquisition_operators[freq](vector_ones)
+            )
+        return scalar_acquisition_operators
+
     def get_convolution(self):
         """QUBIC resolutions.
 
@@ -408,26 +428,25 @@ def get_convolution(self):
         """
 
         ### Define FWHMs
-
         fwhm_in = None
         fwhm_out = None
 
         ### FWHMs during map-making
         if self.params["QUBIC"]["convolution_in"]:
-            fwhm_in = self.joint.qubic.allfwhm.copy()
+            fwhm_in = self.joint_tod.qubic.allfwhm.copy()
         if self.params["QUBIC"]["convolution_out"]:
             fwhm_out = np.array([])
             for irec in range(self.params["QUBIC"]["nrec"]):
                 fwhm_out = np.append(
                     fwhm_out,
                     np.sqrt(
                         self.joint.qubic.allfwhm[
-                            irec * self.fsub : (irec + 1) * self.fsub
+                            irec * self.fsub_out : (irec + 1) * self.fsub_out
                         ]
                         ** 2
                         - np.min(
                             self.joint.qubic.allfwhm[
-                                irec * self.fsub : (irec + 1) * self.fsub
+                                irec * self.fsub_out : (irec + 1) * self.fsub_out
                             ]
                         )
                         ** 2
@@ -445,7 +464,7 @@ def get_convolution(self):
                     fwhm_rec,
                     np.min(
                         self.joint.qubic.allfwhm[
-                            irec * self.fsub : (irec + 1) * self.fsub
+                            irec * self.fsub_out : (irec + 1) * self.fsub_out
                         ]
                     ),
                 )
@@ -455,20 +474,50 @@ def get_convolution(self):
             and self.params["QUBIC"]["convolution_out"] is False
         ):
             fwhm_rec = np.array([])
+            scalar_acquisition_operators = self._get_scalar_acquisition_operator()
+
+            if self.params["Foregrounds"]["Dust"]:
+                f_dust = c.ModifiedBlackBody(nu0=353, beta_d=1.54)
+                weight_factor = f_dust.eval(self.joint.qubic.allnus)
+                fun = lambda nu: np.abs(fraction - f_dust.eval(nu))
+            else:
+                f_cmb = c.CMB()
+                weight_factor = f_cmb.eval(self.joint.qubic.allnus)
+                fun = lambda nu: np.abs(fraction - f_cmb.eval(nu))
+
+            ### Compute expected resolutions and frequencies when not adding convolutions during reconstruction
+            ### See FMM annexe B to understand the computations
+
             for irec in range(self.params["QUBIC"]["nrec"]):
+                numerator_fwhm, denominator_fwhm = 0, 0
+                numerator_nus, denominator_nus = 0, 0
+                for jsub in range(irec * self.fsub_out, (irec + 1) * self.fsub_out):
+                    # Compute the expected reconstructed resolution for sub-acquisition
+                    numerator_fwhm += (
+                        scalar_acquisition_operators[jsub]
+                        * weight_factor[jsub]
+                        * self.fwhm_in[jsub]
+                    )
+                    denominator_fwhm += (
+                        scalar_acquisition_operators[jsub] * weight_factor[jsub]
+                    )
+
+                    # Compute the expected reconstructed frequencies for sub_acquisition
+                    numerator_nus += (
+                        scalar_acquisition_operators[jsub] * weight_factor[jsub]
+                    )
+                    denominator_nus += scalar_acquisition_operators[jsub]
+
+                # Compute the expected resolution
                 fwhm_rec = np.append(
-                    fwhm_rec,
-                    np.mean(
-                        self.joint.qubic.allfwhm[
-                            irec * self.fsub : (irec + 1) * self.fsub
-                        ]
-                    ),
+                    fwhm_rec, np.sum(numerator_fwhm) / np.sum(denominator_fwhm)
                 )
 
-        else:
-            fwhm_rec = np.array([])
-            for irec in range(self.params["QUBIC"]["nrec"]):
-                fwhm_rec = np.append(fwhm_rec, 0)
+                # Compute the expected frequency
+                fraction = np.sum(numerator_nus) / np.sum(denominator_nus)
+                x0 = self.nus_Q[irec]
+                corrected_nu = minimize(fun, x0)
+                self.nus_Q[irec] = corrected_nu["x"]
 
         if self.rank == 0:
             print(f"FWHM for TOD generation : {fwhm_in}")
@@ -492,7 +541,7 @@ def get_input_map(self):
         m_nu_in = np.zeros(
             (self.params["QUBIC"]["nrec"], 12 * self.params["SKY"]["nside"] ** 2, 3)
         )
-        
+
         for i in range(self.params["QUBIC"]["nrec"]):
             m_nu_in[i] = np.mean(
                 self.external_timeline.m_nu[i * self.fsub : (i + 1) * self.fsub], axis=0