diff --git a/src/CMM/preset/preset_acquisition.py b/src/CMM/preset/preset_acquisition.py index 4ebfef4..37ac28f 100644 --- a/src/CMM/preset/preset_acquisition.py +++ b/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}") diff --git a/src/FMM/pipeline.py b/src/FMM/pipeline.py index 9ccd704..dbe88f1 100644 --- a/src/FMM/pipeline.py +++ b/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,13 +428,12 @@ 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"]): @@ -422,12 +441,12 @@ def get_convolution(self): 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