diff --git a/mflike/foreground.py b/mflike/foreground.py
index 9e64825..e499cee 100644
--- a/mflike/foreground.py
+++ b/mflike/foreground.py
@@ -493,8 +493,8 @@ def initialize(self):
         super().initialize()
         if self.bands is None:
             self.bands = {
-                f"{exp}_s0": {"nu": [self.bandint_freqs_T[iexp]], "bandpass": [1.0]}
-                for iexp, exp in enumerate(self.experiments)}
+                f"{exp}_{spin}": {"nu": [self.bandint_freqs_T[iexp]], "bandpass": [1.0]}
+                for spin in ["s0","s2"] for iexp, exp in enumerate(self.experiments)}
         self._initialized = False
         self.init_bandpowers()
 
@@ -516,7 +516,10 @@ def init_bandpowers(self):
                     self.log, "One band has width = 0, set a positive width and run again"
                 )
         
-        self.use_beam_profile = bool(self.beam_profile)
+        if not self._initialized:
+            self.use_beam_profile = False
+        else:
+            self.use_beam_profile = bool(self.beam_profile)
         if self.use_beam_profile:
             if not self.beam_profile.get("Gaussian_beam"):
                 self.beam_file = self.beam_profile.get("beam_from_file")
@@ -527,6 +530,7 @@ def init_bandpowers(self):
             # reading the possible dictionary with beam profiles associated to bandpass shifts
             # this has to be present if we want to propagate bandpass shifts to the chromatic beams
             # otherwise they are left unchanged 
+
             self.bandsh_beams_path = self.beam_profile.get("Bandpass_shifted_beams")
             if self.bandsh_beams_path:
                 self.bandpass_shifted_beams = self._read_yaml_file(self.bandsh_beams_path)
@@ -750,48 +754,40 @@ def _init_beam_from_file(self):
         """
 
         if not self.beam_file:
-            # option to read beam from sacc.
-            try:
-                bool(mflike.beams)
-            except:
+            # option to read beam from sacc
+            if not bool(self.beams):
                 raise ValueError("Beams not stored in sacc files!")
-            else:
-                self.beams = mflike.beams
         else:
             self.beams = self._read_yaml_file(self.beam_file)
-
-        if self._initialized:
-            #checking that the freq array is compatible with the bandpass one
-            for exp in self.experiments:
-                # checking nu is the same as the bandpass one
-                if not np.allclose(self.beams[f"{exp}_s0"]['nu'], self.bands[f"{exp}_s0"]['nu'], atol = 1e-5):
-                    raise LoggedError(self.log, f"Frequency array for beam {exp}_s0 is not the same as the bandpass one!")
-                if not np.allclose(self.beams[f"{exp}_s2"]['nu'], self.bands[f"{exp}_s2"]['nu'], atol = 1e-5):
-                    raise LoggedError(self.log, f"Frequency array for beam {exp}_s2 is not the same as the bandpass one!")
+        
+        #checking that the freq array is compatible with the bandpass one
+        from itertools import product
+        for exp, spin in product(self.experiments, ["s0", "s2"]):
+            key =  f"{exp}_{spin}"
+            # checking nu is the same as the bandpass one
+            if not np.allclose(self.beams[key]['nu'], self.bands[key]['nu'], atol = 1e-5):
+                raise LoggedError(self.log, f"Frequency array for beam {key} is not the same as the bandpass one!")
+            # checking beam shape is correct
+            if not self.beams[key]["beams"].shape[0] == self.bands[key]['nu'].size:
+                shape_b = self.beams[key]["beams"].shape[0]
+                shape_n = self.bands[key]['nu'].size
+                raise LoggedError(self.log, f"beam {key} has a wrong shape! It is shape ({shape_b}, ells), but shoule be ({shape_n}, ells)")
 
     def _init_gauss_beams(self):
         """
         Computes the dictionary of beams for each frequency of self.experiments
         """
         self.beams = {}
-        for iexp, exp in enumerate(self.experiments):
-            gauss_pars = self.gaussian_params[f"{exp}_s0"]
+        from itertools import product
+        for exp, spin in product(self.experiments, ["s0", "s2"]):
+            key = f"{exp}_{spin}"
+            gauss_pars = self.gaussian_params[key]
             FWHM0 = np.asarray(gauss_pars["FWHM_0"])
-            #using the same freq array as the bandpass one
-            nu = np.asarray(self.bands[f"{exp}_s0"]['nu'])
-            nu0 = np.asarray(gauss_pars["nu_0"])
-            alpha = np.asarray(gauss_pars["alpha"])
-            # computing temperature beam for exp
-            self.beams[f"{exp}_s0"] = {"nu": nu, "beams": self.gauss_beams(FWHM0, nu, nu0, alpha, False)}
-            # doing the same for polarization
-            gauss_pars = self.gaussian_params[f"{exp}_s2"]
-            FWHM0 = np.asarray(gauss_pars["FWHM_0"])
-            # nu pol should be the same as the T one, I'll comment it for now
-            #nu = np.asarray(self.bands[f"{exp}_s2"]['nu'])
+            nu = np.asarray(self.bands[key]['nu'])
             nu0 = np.asarray(gauss_pars["nu_0"])
             alpha = np.asarray(gauss_pars["alpha"])
             # checking nu is the same as the bandpass one
-            self.beams[f"{exp}_s2"] = {"nu": nu, "beams": self.gauss_beams(FWHM0, nu, nu0, alpha, True)}
+            self.beams[key] = {"nu": nu, "beams": self.gauss_beams(FWHM0, nu, nu0, alpha, pol=spin=="s2")} 
 
     def gauss_beams(self, fwhm0, nu, nu0, alpha, pol):
         r"""
@@ -810,7 +806,6 @@ def gauss_beams(self, fwhm0, nu, nu0, alpha, pol):
         :return: a :math:`b^{Gauss.}_{\ell}(\nu)` = ``array(freqs, lmax +2)`` with Gaussian beam
                  profiles for each frequency in :math:`\nu` (from :math:`\ell = 0`)
         """
-        from astropy import constants, units
         import healpy as hp
 
         fwhm = fwhm0 * (nu / nu0)**(-alpha/2.)
diff --git a/mflike/mflike.py b/mflike/mflike.py
index 0fb34e4..fcfa40c 100644
--- a/mflike/mflike.py
+++ b/mflike/mflike.py
@@ -410,10 +410,10 @@ def get_sacc_names(pol, exp_1, exp_2):
         self.beams = {}
         for name, tracer in s.tracers.items():
             self.bands[name] = {"nu": tracer.nu, "bandpass": tracer.bandpass}
-            # trying to read beams, if present
-            if hasattr(tracer, "beam"):
+            # trying to read beams, if present, and check if it is empty
+            if hasattr(tracer, "beam") and bool(tracer.beam):
                 self.beams[name] = {"nu": tracer.nu, "beams": tracer.beam}
-
+        
         # Put lcuts in a format that is recognisable by CAMB.
         self.lcuts = {k.lower(): c for k, c in self.lcuts.items()}
         if "et" in self.lcuts:
diff --git a/mflike/tests/test_mflike.py b/mflike/tests/test_mflike.py
index e39c16f..ba3f680 100644
--- a/mflike/tests/test_mflike.py
+++ b/mflike/tests/test_mflike.py
@@ -364,7 +364,7 @@ def compute_FWHM(nu):
         # generating the dictionary needed for the bandpass shift case
         test_path = os.path.dirname(__file__)
         import subprocess
-        subprocess.run("python "+os.path.join(test_path, "../../mflike_utils/utils.py"), shell=True, check=True) 
+        subprocess.run("python "+os.path.join(test_path, "../../scripts/generate_beams_w_bandpass_shifts.py"), shell=True, check=True) 
 
         model.close()
         
diff --git a/mflike_utils/utils.py b/scripts/generate_beams_w_bandpass_shifts.py
similarity index 96%
rename from mflike_utils/utils.py
rename to scripts/generate_beams_w_bandpass_shifts.py
index eafc1a4..a71fbc4 100644
--- a/mflike_utils/utils.py
+++ b/scripts/generate_beams_w_bandpass_shifts.py
@@ -1,5 +1,5 @@
 r"""
-Simple code to generate the bandpass shift dictionary needed in the presence of bandpass shifts different from 0.
+Simple code to generate the beam dictionary needed in the presence of bandpass shifts different from 0.
 We compute simple gaussian beams for :math:`\nu_0 + \Delta \nu`, assuming a diffraction limited experiment.
 This is thought to be used in the absence of data coming from the planets beams measurements.
 """