PSD noise ratio cuts

sotodlib/mapmaking/demod_mapmaker.py

@@ -624,6 +624,16 @@ def make_demod_map(context, obslist, noise_model, info,
     n_dets = comm.allreduce(n_dets)
     for subinfo in info:
         subinfo['number_dets'] = n_dets
+
+        # Blindly add *all* scalars in this aman to sub_info.
+        # The ones in the AtomicInfo class will be automatically added to the AtomicInfo, the others ignored.
+        info_aman_name = 'preprocess.split_flags'
+        if obs is not None and info_aman_name in obs:
+            info_aman = obs[info_aman_name]
+            for info_entry in info_aman.keys():
+                if np.isscalar(info_aman[info_entry]):
+                    if info_entry not in subinfo:  # Don't overwrite existing entries
+                        subinfo[info_entry] = info_aman[info_entry]
     # if we skip all the obs then we return error and output
     if nobs_kept == 0:
         return errors, outputs, None

sotodlib/mapmaking/utils.py

@@ -787,6 +787,7 @@ class AtomicInfo(Base):
     moon_distance: Mapped[Optional[float]]
     wind_speed: Mapped[Optional[float]]
     wind_direction: Mapped[Optional[float]]
+    rqu_avg: Mapped[Optional[float]]
 
     def __init__(self, obs_id, telescope, freq_channel, wafer, ctime, split_label):
         self.obs_id = obs_id

sotodlib/obs_ops/splits.py

@@ -89,7 +89,7 @@ def get_split_flags(aman, proc_aman=None, split_cfg=None):
     # Set default set of splits
     default_cfg = {'high_gain': 0.115, 'high_tau': 1.5e-3,
                    'det_A': 'A', 'pol_angle': 35, 'crossover': 'BL', 'right_focal_plane': 0,
-                   'top_focal_plane': 0, 'central_pixels': 0.071 }
+                   'top_focal_plane': 0, 'central_pixels': 0.071, 'high_rq': 1.1, 'high_ru': 1.12}
     if split_cfg is None:
         split_cfg = default_cfg
 
@@ -169,6 +169,29 @@ def get_split_flags(aman, proc_aman=None, split_cfg=None):
         fm.wrap('scan_left', proc_aman.turnaround_flags.left_scan)
         fm.wrap('scan_right', proc_aman.turnaround_flags.right_scan)
 
+    if 'noise_ratio_Q' in proc_aman:
+        # 1/f noise split
+        rq = proc_aman.noise_ratio_Q.rdets
+        if rq.ndim > 1:  # Mean over subscans
+            rq = np.nanmean(rq, axis=-1)
+        fm.wrap_dets('high_rq', rq > split_cfg['high_rq'])
+        fm.wrap_dets('low_rq',  rq <= split_cfg['high_rq'])
+        split_aman.wrap('rq_avg', np.nanmean(proc_aman.noise_ratio_Q.rmean))
+    if 'noise_ratio_U' in proc_aman:
+        ru = proc_aman.noise_ratio_U.rdets
+        if ru.ndim > 1:  # Mean over subscans
+            ru = np.nanmean(ru, axis=-1)
+        fm.wrap_dets('high_ru', ru > split_cfg['high_ru'])
+        fm.wrap_dets('low_ru', ru <= split_cfg['high_ru'])
+        split_aman.wrap('ru_avg', np.nanmean(proc_aman.noise_ratio_U.rmean))
+    if 'noise_ratio_Q' in proc_aman and 'noise_ratio_U' in proc_aman:
+        rqu = np.mean([rq, ru], axis=0)
+        high_rqu = np.mean([split_cfg['high_rq'], split_cfg['high_ru']])
+        fm.wrap_dets('high_rqu', rqu > high_rqu)
+        fm.wrap_dets('low_rqu', rqu <= high_rqu)
+        split_aman.wrap('rqu_avg', np.nanmean([proc_aman.noise_ratio_Q.rmean, proc_aman.noise_ratio_U.rmean]))
+
+
     for k in split_cfg.keys():
         split_aman.wrap(f'{k}_threshold', split_cfg[k])
 

sotodlib/preprocess/processes.py

@@ -455,6 +455,76 @@ def plot(self, aman, proc_aman, filename):
             plot_psd(aman, signal=attrgetter(f"{self.wrap}.Pxx")(proc_aman),
                      xx=attrgetter(f"{self.wrap}.freqs")(proc_aman), filename=filename, **self.plot_cfgs)
 
+class NoiseRatio(_Preprocess):
+    """ Compute ratios of "signal band" to white noise in PSDs.
+
+    Example config block::
+
+    - name: "noise_ratio"
+      psd: "psdQ"
+      wrap: "noise_ratio_Q"
+      subscan: False
+      calc:
+        f_sel: [0.04, 0.14]
+        f_wn: [0.6, 1.0]
+      save: True
+      select:
+        r_max: 1.19
+        select_per_detector: True
+
+    .. autofunction:: sotodlib.tod_ops.fft_ops.noise_ratio
+"""
+    name = "noise_ratio"
+
+    def __init__(self, step_cfgs):
+        self.psd = step_cfgs.get('psd', 'psd')
+        self.wrap = step_cfgs.get('wrap', 'noise_ratio')
+        self.subscan = step_cfgs.get('subscan', False)
+
+        super().__init__(step_cfgs)
+
+    def calc_and_save(self, aman, proc_aman):
+        if self.psd not in proc_aman:
+            raise ValueError("PSD is not saved in Preprocessing AxisManager")
+        psd = proc_aman[self.psd]
+        pxx = psd.Pxx_ss if self.subscan else psd.Pxx
+
+        if self.calc_cfgs is None:
+            self.calc_cfgs = {}
+
+        calc_aman = tod_ops.fft_ops.noise_ratio(proc_aman, pxx, psd.freqs, subscan=self.subscan, **self.calc_cfgs)
+        self.save(proc_aman, calc_aman)
+        return aman, proc_aman
+
+    def save(self, proc_aman, calc_aman):
+        if self.save_cfgs is None:
+            return
+        else:
+            proc_aman.wrap(self.wrap, calc_aman)
+
+    def select(self, meta, proc_aman=None, in_place=True):
+        if self.select_cfgs is None:
+            return meta
+
+        if proc_aman is None:
+            proc_aman = meta.preprocess
+
+        per_detector = self.select_cfgs.get("select_per_detector", True)
+        if per_detector:
+            noise_ratio = proc_aman[self.wrap]["rdets"]
+        else:
+            noise_ratio = proc_aman[self.wrap]["rmean"]
+
+        if self.subscan:
+            noise_ratio = np.nanmean(noise_ratio, axis=-1) # Mean over subscans
+        keep = np.ones_like(meta.dets.vals, dtype=bool)
+        keep &= (noise_ratio <= np.float64(self.select_cfgs["r_max"]))
+        if in_place:
+            meta.restrict("dets", meta.dets.vals[keep])
+            return meta
+        else:
+            return keep
+
 
 class GetStats(_Preprocess):
     """
@@ -2734,6 +2804,7 @@ def process(self, aman, proc_aman, sim=False):
 _Preprocess.register(Jumps)
 _Preprocess.register(FixJumps)
 _Preprocess.register(PSDCalc)
+_Preprocess.register(NoiseRatio)
 _Preprocess.register(Noise)
 _Preprocess.register(Calibrate)
 _Preprocess.register(EstimateHWPSS)
@@ -2767,3 +2838,4 @@ def process(self, aman, proc_aman, sim=False):
 _Preprocess.register(SmurfGapsFlags)
 _Preprocess.register(GetTauHWP)
 _Preprocess.register(Move)
+

sotodlib/tod_ops/fft_ops.py

@@ -503,6 +503,53 @@ def calc_wn(aman, pxx=None, freqs=None, nseg=None, low_f=5, high_f=10):
     wn = np.sqrt(wn2)
     return wn
 
+def noise_ratio(aman, pxx, freqs, f_sig=(0.04, 0.14), f_wn=(0.6, 1.0), subscan=False):
+    """Compute the ratio of the mean PSD in two frequency regions to evaluate the noise.
+
+    Arguments
+    ---------
+        aman: AxisManager
+            Only used for matching the dets and subscans dimensions in the output aman.
+        pxx: np.ndarray[float]
+            Input PSD. Can be [dets, nufreq] or [dets, nufreq, subscans] (NOT just [nufreq]).
+        freqs: np.ndarray[float]
+            frequency information related to the psd.
+        f_sig: tuple
+            2-tuple of frequencies giving the range of the numerator ("signal band")
+        f_wn: tuple
+            2-tuple of frequencies giving the range of the denominator ("white noise")
+        subscan: bool
+            True if the PSD is split by subscans
+
+    Returns
+    -------
+        calc_aman: AxisManager
+            Axis manager with fields "rdets" giving the per-detector ratio and
+            "rmean" giving the ratio for the mean (over detectors) PSD.
+    """
+
+    fselect = np.logical_and(freqs >= f_sig[0], freqs <= f_sig[1])
+    fwn = np.logical_and(freqs >= f_wn[0], freqs <= f_wn[1])
+
+    pxxmean = np.mean(pxx, axis=0)
+    rmean = np.mean(pxxmean[fselect], axis=0) / np.mean(pxxmean[fwn], axis=0)
+    rdets = np.mean(pxx[:, fselect], axis=1) / np.mean(pxx[:, fwn], axis=1)
+    rmean = np.array([rmean]*aman.dets.count)  # Make a dets axis since scalar wrapping is broken
+
+    if not subscan:
+        calc_aman = core.AxisManager(aman.dets)
+        calc_aman.wrap("rdets", rdets, [(0,"dets")])
+        calc_aman.wrap("rmean", rmean, [(0,"dets")])
+    else:
+        try:
+            calc_aman = core.AxisManager(aman.dets, aman.turnaround_flags.subscan_info.subscans)
+        except AttributeError:
+            calc_aman = core.AxisManager(aman.dets, aman.subscan_info.subscans)
+        calc_aman.wrap("rdets", rdets, [(0,"dets"), (1,"subscans")])
+        calc_aman.wrap("rmean", rmean, [(0,"dets"), (1,"subscans")])
+
+    return calc_aman
+
 def noise_model(f, params, **fixed_param):
     """
     Noise model for power spectrum with white noise, and 1/f noise.