Latest processes for maps

sotodlib/coords/demod.py

@@ -182,20 +182,24 @@ def from_map(tod, signal_map, cuts=None, flip_gamma=True, wrap=False, modulated=
         return signal_sim
 
 
-def rotate_demodQU(tod, update_focal_plane=True):
+def rotate_demodQU(tod, sign=1, offset=0, update_focal_plane=True):
     """
     Apply detectors' polarization angle calibration to the HWP demodulated Q and U timestreams to
     place all detectors' Q and U timestreams in a common telescope frame. This updates tod.demodQ
     and tod.demodU, in place.
 
     Args:
         tod : an axisManager object
-            update_focal_plane (bool, optional): Whether to set focal_plane.gamma angles to zero,
+        update_focal_plane (bool, optional): Whether to set focal_plane.gamma angles to zero,
             consistent with new coordinate reference. Make this true for polarization mapmaking
             using make_map.
+        offset : float, optional
+            The rotation angle in degrees to apply (default is 0).
+        sign : int, optional
+            A sign factor to control the direction of the rotation (default is +1).
 
     """
-    demodC = ((tod.demodQ + 1j*tod.demodU).T * np.exp(-2j*tod.focal_plane.gamma)).T
+    demodC = ((tod.demodQ + 1j*tod.demodU).T * np.exp( sign*(-2j*tod.focal_plane.gamma + 1j*np.deg2rad(offset)) )).T
     tod.demodQ = demodC.real
     tod.demodU = demodC.imag
     del demodC

sotodlib/preprocess/processes.py

@@ -9,10 +9,11 @@
 
 from sotodlib.core.flagman import (has_any_cuts, has_all_cut,
                                    count_cuts,
-                                    sparse_to_ranges_matrix)
+                                   sparse_to_ranges_matrix)
 
 from .pcore import _Preprocess, _FracFlaggedMixIn
 from .. import flag_utils
+from ..core import AxisManager
 
 
 class FFTTrim(_Preprocess):
@@ -21,7 +22,7 @@ class FFTTrim(_Preprocess):
 
     .. autofunction:: sotodlib.tod_ops.fft_trim
     """
-    name = "fft_trim"    
+    name = "fft_trim"
     def process(self, aman, proc_aman):
         tod_ops.fft_trim(aman, **self.process_cfgs)
 
@@ -113,7 +114,6 @@ def __init__(self, step_cfgs):
         self.signal = step_cfgs.get('signal', 'signal')
 
         super().__init__(step_cfgs)
-
     def calc_and_save(self, aman, proc_aman):
         _, trend_aman = tod_ops.flags.get_trending_flags(
             aman, merge=False, full_output=True,
@@ -551,15 +551,28 @@ def select(self, meta, proc_aman=None):
         if proc_aman is None:
             proc_aman = meta.preprocess
 
-        self.select_cfgs['name'] = self.select_cfgs.get('name','noise')
+        if 'wrap_name' in self.save_cfgs:
+            self.select_cfgs['name'] = self.select_cfgs.get('name', self.save_cfgs['wrap_name'])
+        else:
+            self.select_cfgs['name'] = self.select_cfgs.get('name', 'noise')
 
         if self.fit:
             wn = proc_aman[self.select_cfgs['name']].fit[:,1]
+            fk = proc_aman[self.select_cfgs['name']].fit[:,0]
         else:
             wn = proc_aman[self.select_cfgs['name']].white_noise
+            fk = None 
         if self.subscan:
             wn = np.nanmean(wn, axis=-1) # Mean over subscans
-        keep = wn <= np.float64(self.select_cfgs["max_noise"])
+            if fk is not None:
+                fk = np.nanmean(fk, axis=-1) # Mean over subscans
+        keep = np.ones_like(wn, dtype=bool)
+        if "max_noise" in self.select_cfgs.keys():
+            keep = (wn <= np.float64(self.select_cfgs["max_noise"]))
+        if "min_noise" in self.select_cfgs.keys():
+            keep &= (wn >= np.float64(self.select_cfgs["min_noise"]))
+        if fk is not None and "max_fknee" in self.select_cfgs.keys():
+            keep &= (fk <= np.float64(self.select_cfgs["max_fknee"]))
         meta.restrict("dets", meta.dets.vals[keep])
         return meta
 
@@ -865,7 +878,7 @@ class FlagTurnarounds(_Preprocess):
     .. autofunction:: sotodlib.tod_ops.flags.get_turnaround_flags
     """
     name = 'flag_turnarounds'
-    
+
     def calc_and_save(self, aman, proc_aman):
         if self.calc_cfgs is None:
             self.calc_cfgs = {}
@@ -879,7 +892,7 @@ def calc_and_save(self, aman, proc_aman):
             calc_aman.wrap('turnarounds', ta, [(0, 'dets'), (1, 'samps')])
             calc_aman.wrap('left_scan', left, [(0, 'dets'), (1, 'samps')])
             calc_aman.wrap('right_scan', right, [(0, 'dets'), (1, 'samps')])
-            
+
         if self.calc_cfgs['method'] == 'az':
             ta = tod_ops.flags.get_turnaround_flags(aman, **self.calc_cfgs)
             calc_aman = core.AxisManager(aman.dets, aman.samps)
@@ -898,7 +911,7 @@ def save(self, proc_aman, turn_aman):
 
     def process(self, aman, proc_aman):
         tod_ops.flags.get_turnaround_flags(aman, **self.process_cfgs)
-        
+
 class SubPolyf(_Preprocess):
     """Fit TOD in each subscan with polynominal of given order and subtract it.
         All process configs go to `sotodlib.tod_ops.sub_polyf`.
@@ -1095,7 +1108,7 @@ def save(self, proc_aman, hwp_angle_aman):
             return
         if self.save_cfgs:
             proc_aman.wrap("hwp_angle", hwp_angle_aman)
-        
+
 
 class FourierFilter(_Preprocess):
     """
@@ -1126,6 +1139,7 @@ class FourierFilter(_Preprocess):
     See :ref:`fourier-filters` documentation for more details.
     """
     name = 'fourier_filter'
+
     def __init__(self, step_cfgs):
         self.signal_name = step_cfgs.get('signal_name', 'signal')
         # By default signal is overwritted by the filtered signal
@@ -1196,6 +1210,7 @@ class PCARelCal(_Preprocess):
     See :ref:`pca-background` for more details on the method.
     """
     name = 'pca_relcal'
+
     def __init__(self, step_cfgs):
         self.signal = step_cfgs.get('signal', 'signal')
         self.run = step_cfgs.get('pca_run', 'run1')
@@ -1272,7 +1287,7 @@ def select(self, meta, proc_aman=None):
         keep = ~proc_aman[self.run_name]['pca_det_mask']
         meta.restrict("dets", meta.dets.vals[keep])
         return meta
-    
+
     def plot(self, aman, proc_aman, filename):
         if self.plot_cfgs is None:
             return
@@ -1315,13 +1330,13 @@ def calc_and_save(self, aman, proc_aman):
         ptp_aman = core.AxisManager(aman.dets, aman.samps)
         ptp_aman.wrap('ptp_flags', mskptps, [(0, 'dets'), (1, 'samps')])
         self.save(proc_aman, ptp_aman)
-    
-    def save(self, proc_aman, dark_aman):
+
+    def save(self, proc_aman, calc_aman):
         if self.save_cfgs is None:
             return
         if self.save_cfgs:
-            proc_aman.wrap("ptp_flags", dark_aman)
-    
+            proc_aman.wrap("ptp_flags", calc_aman)
+
     def select(self, meta, proc_aman=None):
         if self.select_cfgs is None:
             return meta
@@ -1355,13 +1370,13 @@ def calc_and_save(self, aman, proc_aman):
         ptp_aman = core.AxisManager(aman.dets, aman.samps)
         ptp_aman.wrap('inv_var_flags', mskptps, [(0, 'dets'), (1, 'samps')])
         self.save(proc_aman, ptp_aman)
-    
+
     def save(self, proc_aman, dark_aman):
         if self.save_cfgs is None:
             return
         if self.save_cfgs:
             proc_aman.wrap("inv_var_flags", dark_aman)
-    
+
     def select(self, meta, proc_aman=None):
         if self.select_cfgs is None:
             return meta
@@ -1370,7 +1385,7 @@ def select(self, meta, proc_aman=None):
         keep = ~has_all_cut(proc_aman.inv_var_flags.inv_var_flags)
         meta.restrict("dets", meta.dets.vals[keep])
         return meta
-    
+
 class EstimateT2P(_Preprocess):
     """Estimate T to P leakage coefficients.
 
@@ -1398,7 +1413,7 @@ class EstimateT2P(_Preprocess):
     def calc_and_save(self, aman, proc_aman):
         t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs(aman, **self.calc_cfgs)
         self.save(proc_aman, t2p_aman)
-    
+
     def save(self, proc_aman, t2p_aman):
         if self.save_cfgs is None:
             return
@@ -1422,6 +1437,7 @@ class SubtractT2P(_Preprocess):
     def process(self, aman, proc_aman):
         tod_ops.t2pleakage.subtract_t2p(aman, proc_aman['t2p'],
                                         **self.process_cfgs)
+
 class SplitFlags(_Preprocess):
     """Get flags used for map splitting/bundling.
 
@@ -1465,7 +1481,7 @@ class UnionFlags(_Preprocess):
 
     Saves results for aman under the "flags.[total_flags_label]" field.
 
-     Example config block::
+     Example config block:
 
         - name : "union_flags"
           process:
@@ -1489,6 +1505,150 @@ def process(self, aman, proc_aman):
             aman['flags'].move(self.process_cfgs['total_flags_label'], None)
         aman['flags'].wrap(self.process_cfgs['total_flags_label'], total_flags)
 
+class RotateQU(_Preprocess):
+    """Rotate Q and U components to/from telescope coordinates.
+
+    Example config block::
+
+        - name : "rotate_qu"
+            process:
+            sign: 1 
+            offset: 0 
+            update_focal_plane: True
+
+    .. autofunction:: sotodlib.coords.demod.rotate_demodQU
+    """
+    name = "rotate_qu"
+
+    def process(self, aman, proc_aman):
+        from sotodlib.coords import demod
+        demod.rotate_demodQU(aman, **self.process_cfgs)
+
+class SubtractQUCommonMode(_Preprocess):
+    """Subtract Q and U common mode.
+
+    Example config block::
+
+        - name : "subtract_qu_common_mode"
+          wrap_name: "commonmode"
+          calc: True
+          save: True
+
+    .. autofunction:: sotodlib.tod_ops.deproject.medQU_correct
+    """
+    name = "subtract_qu_common_mode"
+
+    def __init__(self, step_cfgs):
+        self.signal_name_Q = step_cfgs.get('signal_Q', 'demodQ')
+        self.signal_name_U = step_cfgs.get('signal_U', 'demodU')
+        self.wrap_name = step_cfgs.get('wrap_name', 'commonmode')
+        self.run_name_Q = f'{self.signal_name_Q}_{self.wrap_name}'
+        self.run_name_U = f'{self.signal_name_U}_{self.wrap_name}'
+        super().__init__(step_cfgs)
+
+    def calc_and_save(self, aman, proc_aman):
+        for signal_name, run_name in [(self.signal_name_Q, self.run_name_Q), \
+            (self.signal_name_U, self.run_name_U)]:
+            coeffs_aman, med_aman = self._calc_and_save_signal(aman, proc_aman, signal_name, run_name)
+            self.save(proc_aman, med_aman, f'{run_name}_med')
+            self.save(proc_aman, coeffs_aman, f'{run_name}_coeffs')
+
+    def _calc_and_save_signal(self, aman, proc_aman, signal_name, run_name):
+        coeffs, med = tod_ops.deproject.deprojection(aman, signal_name, **self.calc_cfgs)
+        coeffs_aman = core.AxisManager(aman.dets)
+        coeffs_aman.wrap('coeffs', coeffs)
+        med_aman = core.AxisManager(aman.samps)
+        med_aman.wrap('med', med, [(0, 'samps')])
+        return coeffs_aman, med_aman
+
+    def save(self, proc_aman, calc_aman, run_name):
+        if self.save_cfgs is None:
+            return
+        if self.save_cfgs:
+            proc_aman.wrap(run_name, calc_aman)
+
+class FocalplaneNanFlags(_Preprocess):
+    """Find additional detectors which have nans 
+       in their focal plane coordinates.
+
+    Saves results in proc_aman under the "fp_flags" field. 
+
+     Example config block::
+
+        - name : "fp_flags"
+          signal: "signal" # optional
+          calc: True
+          save: True
+          select: True
+
+    .. autofunction:: sotodlib.tod_ops.flags.get_focalplane_flags
+    """
+    name = "fp_flags"
+
+    def calc_and_save(self, aman, proc_aman):
+        mskfp = tod_ops.flags.get_focalplane_flags(aman, **self.calc_cfgs)
+        fp_aman = core.AxisManager(aman.dets, aman.samps)
+        fp_aman.wrap('fp_nans', mskfp, [(0, 'dets'), (1, 'samps')])
+        self.save(proc_aman, fp_aman)
+
+    def save(self, proc_aman, fp_aman):
+        if self.save_cfgs is None:
+            return
+        if self.save_cfgs:
+            proc_aman.wrap("fp_flags", fp_aman)
+
+    def select(self, meta, proc_aman=None):
+        if self.select_cfgs is None:
+            return meta
+        if proc_aman is None:
+            proc_aman = meta.preprocess
+        keep = ~has_all_cut(proc_aman.fp_flags.fp_nans)
+        meta.restrict("dets", meta.dets.vals[keep])
+        return meta
+
+class NoiseFlags(_Preprocess):
+    """Find detectors with anomalous white noise / fknee.
+
+    Saves results in proc_aman under the "wn_flags" field. 
+
+     Example config block::
+
+        - name: "wn_fk_flags"
+          calc:
+            low_wn: 5
+            high_wn: 60
+            high_fk: 6
+          save: True
+          select: True
+
+    .. autofunction:: sotodlib.tod_ops.flags.noise_fit_flags
+    """
+    name = "wn_fk_flags"
+
+    def calc_and_save(self, aman, proc_aman):
+        mskwn, mskfk = tod_ops.flags.noise_fit_flags(aman, **self.calc_cfgs)
+        calc_aman = core.AxisManager(aman.dets, aman.samps)
+        calc_aman.wrap('wn_flags', mskwn)
+        calc_aman.wrap('fknee_flags', mskfk)
+        self.save(proc_aman, calc_aman)
+
+    def save(self, proc_aman, calc_aman):
+        if self.save_cfgs is None:
+            return
+        if self.save_cfgs:
+            proc_aman.wrap("noise_flags", calc_aman)
+
+    def select(self, meta, proc_aman=None):
+        if self.select_cfgs is None:
+            return meta
+        if proc_aman is None:
+            proc_aman = meta.preprocess
+        if "wn_flags" in proc_aman:
+            meta.restrict('dets', proc_aman.dets.vals[proc_aman.noise_flags.wn_flags])
+        if "fknee_flags" in proc_aman:
+            meta.restrict('dets', proc_aman.dets.vals[proc_aman.noise_flags.fknee_flags])
+        return meta
+
 class PointingModel(_Preprocess):
     """Apply pointing model to the TOD.
 
@@ -1507,7 +1667,7 @@ def process(self, aman, proc_aman):
         from sotodlib.coords import pointing_model
         if self.process_cfgs:
             pointing_model.apply_pointing_model(aman)
-            
+
 _Preprocess.register(SplitFlags)
 _Preprocess.register(SubtractT2P)
 _Preprocess.register(EstimateT2P)
@@ -1539,4 +1699,8 @@ def process(self, aman, proc_aman):
 _Preprocess.register(HWPAngleModel)
 _Preprocess.register(GetStats)
 _Preprocess.register(UnionFlags)
-_Preprocess.register(PointingModel)
+_Preprocess.register(RotateQU) 
+_Preprocess.register(SubtractQUCommonMode) 
+_Preprocess.register(FocalplaneNanFlags) 
+_Preprocess.register(NoiseFlags) 
+_Preprocess.register(PointingModel)