diff --git a/src/libtoast/src/toast_tod_filter.cpp b/src/libtoast/src/toast_tod_filter.cpp
index d59847c73..e265661e4 100644
--- a/src/libtoast/src/toast_tod_filter.cpp
+++ b/src/libtoast/src/toast_tod_filter.cpp
@@ -81,6 +81,7 @@ void toast::filter_polynomial(int64_t order, size_t n, uint8_t * flags,
                 last = &full_templates[(iorder - 1) * scanlen];
                 lastlast = &full_templates[(iorder - 2) * scanlen];
                 double orderinv = 1. / iorder;
+
                 #pragma omp simd
                 for (size_t i = 0; i < scanlen; ++i) {
                     const double x = xstart + i * dx;
@@ -130,6 +131,16 @@ void toast::filter_polynomial(int64_t order, size_t n, uint8_t * flags,
                                     singular_values.data(), rcond_limit,
                                     &rank, WORK.data(), LWORK, &info);
 
+        if (info != 0) {
+            auto log = toast::Logger::get();
+            std::ostringstream o;
+            o << "DGELLS:  " << ngood << "/" << scanlen << " good samples, order " <<
+            norder;
+            o << " failed with info " << info;
+            log.error(o.str().c_str(), TOAST_HERE());
+            throw std::runtime_error(o.str().c_str());
+        }
+
         for (int iorder = 0; iorder < norder; ++iorder) {
             double * temp = &full_templates[iorder * scanlen];
             for (int isignal = 0; isignal < nsignal; ++isignal) {
@@ -143,16 +154,6 @@ void toast::filter_polynomial(int64_t order, size_t n, uint8_t * flags,
                 }
             }
         }
-
-        if (info != 0) {
-            auto log = toast::Logger::get();
-            std::ostringstream o;
-            o << "DGELLS:  " << ngood << "/" << scanlen << " good samples, order " <<
-            norder;
-            o << " failed with info " << info;
-            log.error(o.str().c_str(), TOAST_HERE());
-            throw std::runtime_error(o.str().c_str());
-        }
     }
 }
 
diff --git a/src/toast/_libtoast/ops_filterbin.cpp b/src/toast/_libtoast/ops_filterbin.cpp
index db3d4de40..11f6b9977 100644
--- a/src/toast/_libtoast/ops_filterbin.cpp
+++ b/src/toast/_libtoast/ops_filterbin.cpp
@@ -273,7 +273,6 @@ void expand_matrix(py::array_t <double> compressed_matrix,
     }
 }
 
-
 void build_template_covariance(std::vector <int64_t> & starts,
                                std::vector <int64_t> & stops,
                                std::vector <py::array_t <double,
@@ -312,7 +311,6 @@ void build_template_covariance(std::vector <int64_t> & starts,
     }
 }
 
-
 // void build_template_covariance(
 //     py::buffer starts,
 //     py::buffer stops,
@@ -340,12 +338,12 @@ void build_template_covariance(std::vector <int64_t> & starts,
 //         stops, "stops", 1, temp_shape, {n_template}
 //     );
 
-//     // Template covariance    
+//     // Template covariance
 //     double * raw_template_cov = extract_buffer <double> (
-//         template_covariance, 
-//         "template_covariance", 
-//         2, 
-//         temp_shape, 
+//         template_covariance,
+//         "template_covariance",
+//         2,
+//         temp_shape,
 //         {n_template, n_template}
 //     );
 
@@ -377,7 +375,10 @@ void build_template_covariance(std::vector <int64_t> & starts,
 //                 val = 1;
 //             }
 //             for (int64_t isample = istart; isample < istop; ++isample) {
-//                 //std::cerr << "DBG: (" << row << "," << col << "," << isample << ") trow[" << isample << "-" << row_start << "]=" << raw_row[isample - row_start] << " tcol[" << isample << "-" << col_start << "]=" << raw_col[isample - col_start] << " g=" << raw_good[isample] << std::endl;
+//                 //std::cerr << "DBG: (" << row << "," << col << "," << isample << ")
+// trow[" << isample << "-" << row_start << "]=" << raw_row[isample - row_start] << "
+// tcol[" << isample << "-" << col_start << "]=" << raw_col[isample - col_start] << "
+// g=" << raw_good[isample] << std::endl;
 //                 val += raw_row[isample - row_start]
 //                        * raw_col[isample - col_start]
 //                        * raw_good[isample];
diff --git a/src/toast/_libtoast/tod_filter.cpp b/src/toast/_libtoast/tod_filter.cpp
index 17477d8f2..380561d5c 100644
--- a/src/toast/_libtoast/tod_filter.cpp
+++ b/src/toast/_libtoast/tod_filter.cpp
@@ -7,7 +7,9 @@
 
 
 void sum_detectors(py::array_t <int64_t,
-                                py::array::c_style | py::array::forcecast> detectors,
+                                py::array::c_style | py::array::forcecast> det_indx,
+                   py::array_t <int64_t,
+                                py::array::c_style | py::array::forcecast> flag_indx,
                    py::array_t <unsigned char,
                                 py::array::c_style | py::array::forcecast> shared_flags,
                    unsigned char shared_flag_mask,
@@ -20,14 +22,15 @@ void sum_detectors(py::array_t <int64_t,
                                 py::array::c_style | py::array::forcecast> sum_data,
                    py::array_t <int64_t,
                                 py::array::c_style | py::array::forcecast> hits) {
-    auto fast_detectors = detectors.unchecked <1>();
+    auto fast_detindx = det_indx.unchecked <1>();
+    auto fast_flagindx = flag_indx.unchecked <1>();
     auto fast_shared_flags = shared_flags.unchecked <1>();
     auto fast_det_data = det_data.unchecked <2>();
     auto fast_det_flags = det_flags.unchecked <2>();
     auto fast_sum_data = sum_data.mutable_unchecked <1>();
     auto fast_hits = hits.mutable_unchecked <1>();
 
-    size_t ndet = fast_detectors.shape(0);
+    size_t ndet = fast_detindx.shape(0);
     size_t nsample = fast_det_data.shape(1);
 
     size_t buflen = 10000;
@@ -38,11 +41,12 @@ void sum_detectors(py::array_t <int64_t,
         size_t sample_start = ibuf * buflen;
         size_t sample_stop = std::min(sample_start + buflen, nsample);
         for (size_t idet = 0; idet < ndet; ++idet) {
-            int64_t row = fast_detectors(idet);
+            int64_t datarow = fast_detindx(idet);
+            int64_t flagrow = fast_flagindx(idet);
             for (size_t sample = sample_start; sample < sample_stop; ++sample) {
                 if ((fast_shared_flags(sample) & shared_flag_mask) != 0) continue;
-                if ((fast_det_flags(row, sample) & det_flag_mask) != 0) continue;
-                fast_sum_data(sample) += fast_det_data(row, sample);
+                if ((fast_det_flags(flagrow, sample) & det_flag_mask) != 0) continue;
+                fast_sum_data(sample) += fast_det_data(datarow, sample);
                 fast_hits(sample)++;
             }
         }
diff --git a/src/toast/noise_sim.py b/src/toast/noise_sim.py
index 5ad476bb3..befe2cdfe 100644
--- a/src/toast/noise_sim.py
+++ b/src/toast/noise_sim.py
@@ -137,7 +137,7 @@ def NET(self, det):
     def _detector_weight(self, det):
         if self._NET[det] == 0:
             nunit = self._NET[det].unit
-            return 0.0 * (1.0 / (nunit ** 2) / u.Hz)
+            return 0.0 * (1.0 / (nunit**2) / u.Hz)
         else:
             wt = 1.0 / (self._NET[det] ** 2) / self._rate[det]
             return wt.decompose()
diff --git a/src/toast/observation.py b/src/toast/observation.py
index fb27222c5..d67bdf0cd 100644
--- a/src/toast/observation.py
+++ b/src/toast/observation.py
@@ -102,14 +102,14 @@ def set_default_values(values=None):
 
     # composite masks for convenience
     defaults["shared_mask_nonscience"] = (
-        defaults["shared_mask_invalid"] | 
-        defaults["shared_mask_unstable_scanrate"] | 
-        defaults["shared_mask_irregular"]
+        defaults["shared_mask_invalid"]
+        | defaults["shared_mask_unstable_scanrate"]
+        | defaults["shared_mask_irregular"]
     )
     defaults["det_mask_nonscience"] = (
-        defaults["det_mask_invalid"] | 
-        defaults["det_mask_processing"] | 
-        defaults["det_mask_sso"]
+        defaults["det_mask_invalid"]
+        | defaults["det_mask_processing"]
+        | defaults["det_mask_sso"]
     )
 
     if values is not None:
@@ -443,6 +443,7 @@ def select_local_detectors(
             for det in self.local_detectors:
                 if (det in sel_set) and (det in good):
                     dets.append(det)
+        # print(f"SELECT mask {int(flagmask)} {selection}: {dets}", flush=True)
         return dets
 
     # Detector set distribution
diff --git a/src/toast/ops/arithmetic.py b/src/toast/ops/arithmetic.py
index 8704ec362..de7c35fa4 100644
--- a/src/toast/ops/arithmetic.py
+++ b/src/toast/ops/arithmetic.py
@@ -70,8 +70,8 @@ def _exec(self, data, detectors=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
-            if len(dets) == 0:
+            local_dets = ob.select_local_detectors(detectors)
+            if len(local_dets) == 0:
                 # Nothing to do for this observation
                 continue
             if self.first not in ob.detdata:
@@ -90,6 +90,16 @@ def _exec(self, data, detectors=None, **kwargs):
             first_units = ob.detdata[self.first].units
             second_units = ob.detdata[self.second].units
 
+            # Operate on the intersection of detectors
+            dets = list(
+                sorted(
+                    set.intersection(
+                        set(ob.detdata[self.first].detectors),
+                        set(ob.detdata[self.second].detectors),
+                    )
+                )
+            )
+
             if self.result == self.first:
                 result_units = first_units
                 scale_first = 1.0
diff --git a/src/toast/ops/demodulation.py b/src/toast/ops/demodulation.py
index 4c0522367..c0fa71525 100644
--- a/src/toast/ops/demodulation.py
+++ b/src/toast/ops/demodulation.py
@@ -192,21 +192,44 @@ def _exec(self, data, detectors=None, **kwargs):
         timer = Timer()
         timer.start()
         for obs in demodulate_obs:
-            dets = obs.select_local_detectors(detectors, flagmask=self.det_flag_mask)
+            # Get the detectors which are not cut with per-detector flags
+            local_dets = obs.select_local_detectors(
+                detectors, flagmask=self.det_flag_mask
+            )
+            # print(f"DEMOD Proc {obs.comm.group_rank}: local {local_dets}")
+            if obs.comm.comm_group is None:
+                all_dets = local_dets
+                # print(f"DEMOD Proc {obs.comm.group_rank}: comm_group is None")
+            else:
+                proc_dets = obs.comm.comm_group.gather(local_dets, root=0)
+                all_dets = None
+                if obs.comm.comm_group.rank == 0:
+                    # print(f"DEMOD gathered {proc_dets}")
+                    all_dets = set()
+                    for pdets in proc_dets:
+                        for d in pdets:
+                            all_dets.add(d)
+                    all_dets = list(sorted(all_dets))
+                all_dets = obs.comm.comm_group.bcast(all_dets, root=0)
+                # print(f"DEMOD Proc {obs.comm.group_rank}: all {all_dets}")
 
             offset = obs.local_index_offset
             nsample = obs.n_local_samples
 
             fsample = obs.telescope.focalplane.sample_rate
             fmax, hwp_rate = self._get_fmax(obs)
+            # print(
+            #     f"DEMOD:  fsample={fsample}, fmax={fmax}, hwprate={hwp_rate}",
+            #     flush=True,
+            # )
             wkernel = self._get_wkernel(fmax, fsample)
             lowpass = Lowpass(wkernel, fmax, fsample, offset, self.nskip, self.window)
 
             # Create a new observation to hold the demodulated and downsampled data
 
-            demod_telescope = self._demodulate_telescope(obs)
+            demod_telescope = self._demodulate_telescope(obs, all_dets)
             demod_times = self._demodulate_times(obs)
-            demod_detsets = self._demodulate_detsets(obs)
+            demod_detsets = self._demodulate_detsets(obs, all_dets)
             demod_sample_sets = self._demodulate_sample_sets(obs)
             demod_process_rows = obs.dist.process_rows
 
@@ -230,7 +253,7 @@ def _exec(self, data, detectors=None, **kwargs):
             # Allocate storage
 
             demod_dets = []
-            for det in dets:
+            for det in local_dets:
                 for prefix in self.prefixes:
                     demod_dets.append(f"{prefix}_{det}")
             n_local = demod_obs.n_local_samples
@@ -244,6 +267,7 @@ def _exec(self, data, detectors=None, **kwargs):
 
             self._demodulate_shared_data(obs, demod_obs)
 
+            # (f"DEMOD Proc {obs.comm.group_rank} local_dets: {demod_obs.local_detectors}", flush=True)
             exists_data = demod_obs.detdata.ensure(
                 self.det_data,
                 detectors=demod_dets,
@@ -254,10 +278,12 @@ def _exec(self, data, detectors=None, **kwargs):
                 self.det_flags, detectors=demod_dets, dtype=np.uint8
             )
 
-            self._demodulate_flags(obs, demod_obs, dets, wkernel, offset)
-            self._demodulate_signal(data, obs, demod_obs, dets, lowpass)
-            self._demodulate_pointing(data, obs, demod_obs, dets, lowpass, offset)
-            self._demodulate_noise(obs, demod_obs, dets, fsample, hwp_rate, lowpass)
+            self._demodulate_flags(obs, demod_obs, local_dets, wkernel, offset)
+            self._demodulate_signal(data, obs, demod_obs, local_dets, lowpass)
+            self._demodulate_pointing(data, obs, demod_obs, local_dets, lowpass, offset)
+            self._demodulate_noise(
+                obs, demod_obs, local_dets, fsample, hwp_rate, lowpass
+            )
 
             self._demodulate_intervals(obs, demod_obs)
 
@@ -282,7 +308,9 @@ def _exec(self, data, detectors=None, **kwargs):
     def _get_fmax(self, obs):
         times = obs.shared[self.times].data
         hwp_angle = np.unwrap(obs.shared[self.hwp_angle].data)
-        hwp_rate = np.mean(np.diff(hwp_angle) / np.diff(times)) / (2 * np.pi) * u.Hz
+        hwp_rate = np.absolute(
+            np.mean(np.diff(hwp_angle) / np.diff(times)) / (2 * np.pi) * u.Hz
+        )
         if self.fmax is not None:
             fmax = self.fmax
         else:
@@ -300,22 +328,30 @@ def _get_wkernel(self, fmax, fsample):
         return wkernel
 
     @function_timer
-    def _demodulate_telescope(self, obs):
+    def _demodulate_telescope(self, obs, all_dets):
         focalplane = obs.telescope.focalplane
         det_data = focalplane.detector_data
         field_names = det_data.colnames
         # Initialize fields to empty lists
-        fields = dict([(name, []) for name in field_names])
-        for idet, det in enumerate(det_data["name"]):
+        fields = {name: list() for name in field_names}
+        all_set = set(all_dets)
+        for row, det in enumerate(det_data["name"]):
+            if det not in all_set:
+                continue
             for field_name in field_names:
                 # Each detector translates into 3 or 5 new entries
                 for prefix in self.prefixes:
                     if field_name == "name":
                         fields[field_name].append(f"{prefix}_{det}")
                     else:
-                        fields[field_name].append(det_data[field_name][idet])
-        fields = [fields[field_name] for field_name in field_names]
-        demod_det_data = QTable(fields, names=field_names)
+                        fields[field_name].append(det_data[field_name][row])
+        demod_det_data = QTable(
+            [fields[field_name] for field_name in field_names], names=field_names
+        )
+        my_all = list()
+        for name in demod_det_data["name"]:
+            my_all.append(name)
+
         demod_focalplane = Focalplane(
             detector_data=demod_det_data,
             field_of_view=focalplane.field_of_view,
@@ -359,23 +395,40 @@ def _demodulate_shared_data(self, obs, demod_obs):
         return
 
     @function_timer
-    def _demodulate_detsets(self, obs):
-        """Lump all derived detectors into detector sets"""
-        detsets = obs.all_detector_sets
-        demod_detsets = []
-        if detsets is None:
-            for det in obs.all_detectors:
-                demod_detset = []
+    def _demodulate_detsets(self, obs, all_dets):
+        """In order to force local detectors to remain on their original
+        process, we create a detector set for each row of the process
+        grid.
+        """
+        log = Logger.get()
+        if obs.comm_col_size == 1:
+            # One process row
+            detsets = [all_dets]
+        else:
+            local_proc_dets = obs.comm_col.gather(obs.local_detectors, root=0)
+            detsets = None
+            if obs.comm_col_rank == 0:
+                all_set = set(all_dets)
+                detsets = list()
+                for iprow, pdets in enumerate(local_proc_dets):
+                    plocal = list()
+                    for d in pdets:
+                        if d in all_set:
+                            plocal.append(d)
+                    if len(plocal) == 0:
+                        msg = f"obs {obs.name}, process row {iprow} has no"
+                        msg += " unflagged detectors.  This may cause an error."
+                        log.warning(msg)
+                    detsets.append(plocal)
+            detsets = obs.comm_col.bcast(detsets, root=0)
+
+        demod_detsets = list()
+        for dset in detsets:
+            demod_detset = list()
+            for det in dset:
                 for prefix in self.prefixes:
                     demod_detset.append(f"{prefix}_{det}")
-                demod_detsets.append(demod_detset)
-        else:
-            for detset in detsets:
-                demod_detset = []
-                for det in detset:
-                    for prefix in self.prefixes:
-                        demod_detset.append(f"{prefix}_{det}")
-                demod_detsets.append(demod_detset)
+            demod_detsets.append(demod_detset)
         return demod_detsets
 
     @function_timer
@@ -634,7 +687,7 @@ class StokesWeightsDemod(Operator):
     )
 
     weights = Unicode(
-        defaults.weights, help="Observation detdata key for output weights"
+        f"{defaults.weights}_demod", help="Observation detdata key for output weights"
     )
 
     single_precision = Bool(False, help="If True, use 32bit float in output")
diff --git a/src/toast/ops/elevation_noise.py b/src/toast/ops/elevation_noise.py
index a3a158c18..7b30d4590 100644
--- a/src/toast/ops/elevation_noise.py
+++ b/src/toast/ops/elevation_noise.py
@@ -213,11 +213,13 @@ def _exec(self, data, detectors=None, **kwargs):
             noise = obs[self.noise_model]
 
             # We will be collectively building the scale factor for all detectors.
-            # Allocate arrays for communication
+            # Allocate arrays for communication.
 
-            local_net_factors = np.zeros(len(obs.local_detectors), dtype=np.float64)
-            local_tot_factors = np.zeros(len(obs.local_detectors), dtype=np.float64)
-            local_rates = np.zeros(len(obs.local_detectors), dtype=np.float64)
+            local_dets = obs.select_local_detectors(flagmask=defaults.det_mask_invalid)
+
+            local_net_factors = np.zeros(len(local_dets), dtype=np.float64)
+            local_tot_factors = np.zeros(len(local_dets), dtype=np.float64)
+            local_rates = np.zeros(len(local_dets), dtype=np.float64)
             local_weights_in = list()
 
             # We are building up a data product (a noise model) which has values for
@@ -244,7 +246,7 @@ def _exec(self, data, detectors=None, **kwargs):
                         flags = None
                 view_flags.append(flags)
 
-            for idet, det in enumerate(obs.local_detectors):
+            for idet, det in enumerate(local_dets):
                 local_rates[idet] = focalplane.sample_rate.to_value(u.Hz)
                 local_weights_in.append(noise.detector_weight(det))
 
@@ -350,11 +352,11 @@ def _exec(self, data, detectors=None, **kwargs):
 
             # Modify all psds first, since the first call to detector_weight()
             # will trigger the calculation for all detectors.
-            for idet, det in enumerate(obs.local_detectors):
+            for idet, det in enumerate(local_dets):
                 out_noise.psd(det)[:] *= local_tot_factors[idet]
 
             local_weights_out = list()
-            for idet, det in enumerate(obs.local_detectors):
+            for idet, det in enumerate(local_dets):
                 local_weights_out.append(out_noise.detector_weight(det))
 
             if obs.comm_row_rank == 0:
diff --git a/src/toast/ops/gainscrambler.py b/src/toast/ops/gainscrambler.py
index d85f3ed02..093023144 100644
--- a/src/toast/ops/gainscrambler.py
+++ b/src/toast/ops/gainscrambler.py
@@ -77,6 +77,8 @@ def _exec(self, data, detectors=None, **kwargs):
             counter1 = 0
             counter2 = 0
 
+            dets_present = set(obs.detdata[self.det_data].detectors)
+
             for det in dets:
                 # Test the detector pattern
                 if pat is not None and pat.match(det) is None:
@@ -93,8 +95,8 @@ def _exec(self, data, detectors=None, **kwargs):
                 )
 
                 gain = self.center + rngdata[0] * self.sigma
-
-                obs.detdata[self.det_data][det] *= gain
+                if det in dets_present:
+                    obs.detdata[self.det_data][det] *= gain
 
         return
 
diff --git a/src/toast/ops/groundfilter.py b/src/toast/ops/groundfilter.py
index d9a9b6d8c..218ce2084 100644
--- a/src/toast/ops/groundfilter.py
+++ b/src/toast/ops/groundfilter.py
@@ -232,7 +232,6 @@ def build_templates(self, obs):
     def fit_templates(
         self,
         obs,
-        det,
         templates,
         ref,
         good,
@@ -350,27 +349,24 @@ def _exec(self, data, detectors=None, **kwargs):
             last_cov = None
             last_rcond = None
             ndet = len(obs.local_detectors)
-            for idet, det in enumerate(
-                obs.select_local_detectors(detectors, flagmask=self.det_flag_mask)
+
+            for det in obs.select_local_detectors(
+                detectors, flagmask=self.det_flag_mask
             ):
                 if data.comm.group_rank == 0:
-                    msg = (
-                        f"{log_prefix} OpGroundFilter: "
-                        f"Processing detector # {idet + 1} / {ndet}"
-                    )
+                    msg = f"{log_prefix} OpGroundFilter: " f"Processing detector {det}"
                     log.verbose(msg)
 
-                ref = obs.detdata[self.det_data][idet]
+                ref = obs.detdata[self.det_data][det]
                 if self.det_flags is not None:
-                    def_flags = obs.detdata[self.det_flags][idet] & self.det_flag_mask
-                    good = np.logical_and(common_flags == 0, def_flags == 0)
+                    test_flags = obs.detdata[self.det_flags][det] & self.det_flag_mask
+                    good = np.logical_and(common_flags == 0, test_flags == 0)
                 else:
                     good = common_flags == 0
 
                 t1 = time()
                 coeff, last_invcov, last_cov, last_rcond = self.fit_templates(
                     obs,
-                    det,
                     templates,
                     ref,
                     good,
@@ -388,6 +384,11 @@ def _exec(self, data, detectors=None, **kwargs):
                     log.verbose(msg)
 
                 if coeff is None:
+                    # All samples flagged or template fit failed.
+                    curflag = obs.local_detector_flags[det]
+                    obs.update_local_detector_flags(
+                        {det: curflag | self.ground_flag_mask}
+                    )
                     continue
 
                 t1 = time()
diff --git a/src/toast/ops/hwpfilter.py b/src/toast/ops/hwpfilter.py
index 1cbdebaab..cc042a6b6 100644
--- a/src/toast/ops/hwpfilter.py
+++ b/src/toast/ops/hwpfilter.py
@@ -171,7 +171,6 @@ def build_templates(self, obs):
     def fit_templates(
         self,
         obs,
-        det,
         templates,
         ref,
         good,
@@ -225,7 +224,6 @@ def fit_templates(
             if cov is None:
                 cov = get_pseudoinverse(invcov)
         coeff = np.dot(cov, proj)
-
         return coeff, invcov, cov, rcond
 
     @function_timer
@@ -234,11 +232,11 @@ def subtract_templates(self, ref, good, coeff, legendre_trend, fourier_filter):
         if self.detrend:
             trend = np.zeros_like(ref)
             add_templates(trend, legendre_trend, coeff[: self.trend_order + 1])
-            ref -= trend
+            ref[:] -= trend
         # HWP template
         hwptemplate = np.zeros_like(ref)
         add_templates(hwptemplate, fourier_filter, coeff[self.trend_order + 1 :])
-        ref -= hwptemplate
+        ref[:] -= hwptemplate
         return
 
     @function_timer
@@ -289,27 +287,23 @@ def _exec(self, data, detectors=None, **kwargs):
             last_cov = None
             last_rcond = None
             ndet = len(obs.local_detectors)
-            for idet, det in enumerate(
-                obs.select_local_detectors(detectors, flagmask=self.det_flag_mask)
+            for det in obs.select_local_detectors(
+                detectors, flagmask=self.det_flag_mask
             ):
                 if data.comm.group_rank == 0:
-                    msg = (
-                        f"{log_prefix} OpHWPFilter: "
-                        f"Processing detector # {idet + 1} / {ndet}"
-                    )
+                    msg = f"{log_prefix} OpHWPFilter: " f"Processing detector {det}"
                     log.verbose(msg)
 
-                ref = obs.detdata[self.det_data][idet]
+                ref = obs.detdata[self.det_data][det]
                 if self.det_flags is not None:
-                    def_flags = obs.detdata[self.det_flags][idet] & self.det_flag_mask
-                    good = np.logical_and(common_flags == 0, def_flags == 0)
+                    test_flags = obs.detdata[self.det_flags][det] & self.det_flag_mask
+                    good = np.logical_and(common_flags == 0, test_flags == 0)
                 else:
                     good = common_flags == 0
 
                 t1 = time()
                 coeff, last_invcov, last_cov, last_rcond = self.fit_templates(
                     obs,
-                    det,
                     templates,
                     ref,
                     good,
@@ -327,6 +321,9 @@ def _exec(self, data, detectors=None, **kwargs):
                     log.verbose(msg)
 
                 if coeff is None:
+                    # All samples flagged or template fit failed.
+                    curflag = obs.local_detector_flags[det]
+                    obs.update_local_detector_flags({det: curflag | self.hwp_flag_mask})
                     continue
 
                 t1 = time()
diff --git a/src/toast/ops/madam.py b/src/toast/ops/madam.py
index 279803b7d..6e2a0592b 100644
--- a/src/toast/ops/madam.py
+++ b/src/toast/ops/madam.py
@@ -521,7 +521,8 @@ def _prepare(self, params, data, detectors):
         for ob in data.obs:
             # Get the detectors we are using for this observation
             dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
-            all_dets.update(dets)
+            for d in dets:
+                all_dets.add(d)
 
             # Check that the timestamps exist.
             if self.times not in ob.shared:
@@ -588,7 +589,7 @@ def _prepare(self, params, data, detectors):
         nsamp = nsamp_valid
 
         interval_starts = np.array(interval_starts, dtype=np.int64)
-        all_dets = sorted(all_dets)
+        all_dets = list(sorted(all_dets))
         ndet = len(all_dets)
 
         nnz_full = len(self.stokes_weights.mode)
@@ -715,7 +716,8 @@ def _stage_data(
                         nse_scale = float(ob[self.noise_scale])
 
                 for det in all_dets:
-                    if det not in ob.local_detectors:
+                    local_dets = ob.select_local_detectors(flagmask=self.det_flag_mask)
+                    if det not in local_dets:
                         continue
                     psd = nse.psd(det).to_value(u.K**2 * u.second) * nse_scale**2
                     detw = nse.detector_weight(det).to_value(1.0 / u.K**2)
@@ -755,7 +757,7 @@ def _stage_data(
                 None,
                 None,
                 None,
-                None,
+                self.det_flag_mask,
                 do_purge=False,
             )
         else:
@@ -777,7 +779,7 @@ def _stage_data(
                     None,
                     None,
                     None,
-                    None,
+                    self.det_flag_mask,
                 )
             else:
                 # Allocate and copy all at once.
@@ -798,7 +800,7 @@ def _stage_data(
                     None,
                     None,
                     None,
-                    None,
+                    self.det_flag_mask,
                     do_purge=False,
                 )
 
@@ -855,7 +857,7 @@ def _stage_data(
                 None,
                 None,
                 None,
-                None,
+                self.det_flag_mask,
                 operator=self.stokes_weights,
             )
 
@@ -989,6 +991,7 @@ def _unstage_data(
                     self._madam_signal_raw,
                     interval_starts,
                     1,
+                    self.det_flag_mask,
                 )
                 del self._madam_signal
                 del self._madam_signal_raw
@@ -1004,6 +1007,7 @@ def _unstage_data(
                     self._madam_signal,
                     interval_starts,
                     1,
+                    self.det_flag_mask,
                 )
 
             log_time_memory(
diff --git a/src/toast/ops/madam_utils.py b/src/toast/ops/madam_utils.py
index 58b82df7a..64640c60f 100644
--- a/src/toast/ops/madam_utils.py
+++ b/src/toast/ops/madam_utils.py
@@ -79,8 +79,9 @@ def stage_local(
             )
     for ob in data.obs:
         views = ob.view[view]
+        local_dets = set(ob.select_local_detectors(flagmask=det_mask))
         for idet, det in enumerate(dets):
-            if det not in ob.local_detectors:
+            if det not in local_dets:
                 continue
             if operator is not None:
                 # Synthesize data for staging
@@ -185,6 +186,7 @@ def restore_local(
     madam_buffer,
     interval_starts,
     nnz,
+    det_mask,
 ):
     """Helper function to create a detdata buffer from madam data."""
     n_det = len(dets)
@@ -195,6 +197,7 @@ def restore_local(
             ob.detdata.create(detdata_name, dtype=detdata_dtype)
         else:
             ob.detdata.create(detdata_name, dtype=detdata_dtype, sample_shape=(nnz,))
+        local_dets = ob.select_local_detectors(flagmask=det_mask)
         # Loop over views
         views = ob.view[view]
         for ivw, vw in enumerate(views):
@@ -205,23 +208,21 @@ def restore_local(
             else:
                 view_samples = vw.stop - vw.start
             offset = interval_starts[interval]
-            ldet = 0
-            for det in dets:
-                if det not in ob.local_detectors:
+            for idet, det in enumerate(dets):
+                if det not in local_dets:
                     continue
-                idet = ob.local_detectors.index(det)
+                idet = local_dets.index(det)
                 slc = slice(
                     (idet * nsamp + offset) * nnz,
                     (idet * nsamp + offset + view_samples) * nnz,
                     1,
                 )
                 if nnz > 1:
-                    views.detdata[detdata_name][ivw][ldet] = madam_buffer[slc].reshape(
+                    views.detdata[detdata_name][ivw][det] = madam_buffer[slc].reshape(
                         (-1, nnz)
                     )
                 else:
-                    views.detdata[detdata_name][ivw][ldet] = madam_buffer[slc]
-                ldet += 1
+                    views.detdata[detdata_name][ivw][det] = madam_buffer[slc]
             interval += 1
     return
 
@@ -239,6 +240,7 @@ def restore_in_turns(
     madam_buffer_raw,
     interval_starts,
     nnz,
+    det_mask,
 ):
     """When restoring data, take turns copying it."""
     for copying in range(n_copy_groups):
@@ -254,6 +256,7 @@ def restore_in_turns(
                 madam_buffer,
                 interval_starts,
                 nnz,
+                det_mask,
             )
             madam_buffer_raw.clear()
         nodecomm.barrier()
diff --git a/src/toast/ops/mapmaker_binning.py b/src/toast/ops/mapmaker_binning.py
index b5269dd18..adf6cc1c8 100644
--- a/src/toast/ops/mapmaker_binning.py
+++ b/src/toast/ops/mapmaker_binning.py
@@ -213,6 +213,11 @@ def _exec(self, data, detectors=None, **kwargs):
             data[self.binned].reset()
             data[self.binned].update_units(1.0 / self.det_data_units)
 
+        # Use the same detector mask in the pointing
+        self.pixel_pointing.detector_pointing.det_flag_mask = self.det_flag_mask
+        if hasattr(self.stokes_weights, "detector_pointing"):
+            self.stokes_weights.detector_pointing.det_flag_mask = self.det_flag_mask
+
         # Noise weighted map.  We output this to the final binned map location,
         # since we will multiply by the covariance in-place.
 
diff --git a/src/toast/ops/mapmaker_solve.py b/src/toast/ops/mapmaker_solve.py
index 5d5da3f5f..0a3817b28 100644
--- a/src/toast/ops/mapmaker_solve.py
+++ b/src/toast/ops/mapmaker_solve.py
@@ -160,6 +160,7 @@ def _exec(self, data, detectors=None, **kwargs):
             map_key=self.binning.binned,
             det_data=det_temp,
             det_data_units=self.det_data_units,
+            det_flag_mask=self.binning.det_flag_mask,
             subtract=True,
         )
 
@@ -167,6 +168,7 @@ def _exec(self, data, detectors=None, **kwargs):
         noise_weight = NoiseWeight(
             noise_model=self.binning.noise_model,
             det_data=det_temp,
+            det_flag_mask=self.binning.det_flag_mask,
             view=pixels.view,
         )
 
@@ -174,8 +176,12 @@ def _exec(self, data, detectors=None, **kwargs):
         self.template_matrix.transpose = True
         self.template_matrix.det_data = det_temp
         self.template_matrix.det_data_units = self.det_data_units
+        self.template_matrix.det_flag_mask = self.binning.det_flag_mask
         self.template_matrix.view = pixels.view
 
+        # Initialize template matrix for all detectors
+        self.template_matrix.initialize(data)
+
         # Create a pipeline that projects the binned map and applies noise
         # weights and templates.
 
@@ -211,8 +217,7 @@ def _exec(self, data, detectors=None, **kwargs):
 
         log.debug_rank("MapMaker   RHS begin run projection pipeline", comm=comm)
 
-        # NOTE: we set `use_accel=False` as templates cannot be initialized on device
-        proj_pipe.apply(data, detectors=detectors, use_accel=False)
+        proj_pipe.apply(data, detectors=detectors)
 
         log.debug_rank(
             "MapMaker   RHS projection pipeline finished in", comm=comm, timer=timer
@@ -374,8 +379,12 @@ def _exec(self, data, detectors=None, **kwargs):
         self.template_matrix.transpose = False
         self.template_matrix.det_data = self.det_temp
         self.template_matrix.det_data_units = self.det_data_units
+        self.template_matrix.det_flag_mask = self.binning.det_flag_mask
         self.template_matrix.view = pixels.view
 
+        # Initialize template matrix for all detectors
+        self.template_matrix.initialize(data)
+
         self.binning.det_data = self.det_temp
         self.binning.det_data_units = self.det_data_units
 
@@ -432,6 +441,7 @@ def _exec(self, data, detectors=None, **kwargs):
             map_key=self.binning.binned,
             det_data=self.det_temp,
             det_data_units=self.det_data_units,
+            det_flag_mask=self.binning.det_flag_mask,
             subtract=True,
         )
 
@@ -439,6 +449,7 @@ def _exec(self, data, detectors=None, **kwargs):
         noise_weight = NoiseWeight(
             noise_model=self.binning.noise_model,
             det_data=self.det_temp,
+            det_flag_mask=self.binning.det_flag_mask,
             view=pixels.view,
         )
 
diff --git a/src/toast/ops/mapmaker_templates.py b/src/toast/ops/mapmaker_templates.py
index b1b8d6685..8528e2c69 100644
--- a/src/toast/ops/mapmaker_templates.py
+++ b/src/toast/ops/mapmaker_templates.py
@@ -173,6 +173,27 @@ def reset_templates(self):
         """Mark templates to be re-initialized on next call to exec()."""
         self._initialized = False
 
+    @function_timer
+    def initialize(self, data, use_accel=False):
+        if not self._initialized:
+            if use_accel:
+                # fail when a user tries to run the initialization pipeline on GPU
+                raise RuntimeError(
+                    "You cannot currently initialize templates on device (please disable accel for this operator/pipeline)."
+                )
+            for tmpl in self.templates:
+                if not tmpl.enabled:
+                    continue
+                if tmpl.view is None:
+                    tmpl.view = self.view
+                tmpl.det_data_units = self.det_data_units
+                tmpl.det_flags = self.det_flags
+                tmpl.det_flag_mask = self.det_flag_mask
+                # This next line will trigger calculation of the number
+                # of amplitudes within each template.
+                tmpl.data = data
+            self._initialized = True
+
     @function_timer
     def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         log = Logger.get()
@@ -202,26 +223,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             msg += "but does not support accelerators"
             raise RuntimeError(msg)
 
-        # On the first call, we initialize all templates using the Data instance and
-        # the fixed options for view, flagging, etc.
+        # Ensure we have initialized templates with the full set of detectors.
         if not self._initialized:
-            if use_accel:
-                # fail when a user tries to run the initialization pipeline on GPU
-                raise RuntimeError(
-                    "You cannot currently initialize templates on device (please disable accel for this operator/pipeline)."
-                )
-            for tmpl in self.templates:
-                if not tmpl.enabled:
-                    continue
-                if tmpl.view is None:
-                    tmpl.view = self.view
-                tmpl.det_data_units = self.det_data_units
-                tmpl.det_flags = self.det_flags
-                tmpl.det_flag_mask = self.det_flag_mask
-                # This next line will trigger calculation of the number
-                # of amplitudes within each template.
-                tmpl.data = data
-            self._initialized = True
+            raise RuntimeError("You must call initialize() before calling exec()")
 
         # Set the data we are using for this execution
         for tmpl in self.templates:
@@ -238,7 +242,10 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             input_units = 1.0 / self.det_data_units
             for ob in data.obs:
                 if self.det_data not in ob.detdata:
-                    print(f"detector data {self.det_data} not in obs {ob.name}", flush=True)
+                    print(
+                        f"detector data {self.det_data} not in obs {ob.name}",
+                        flush=True,
+                    )
                 if ob.detdata[self.det_data].units != input_units:
                     msg = f"obs {ob.name} detdata {self.det_data}"
                     msg += f" does not have units of {input_units}"
@@ -573,6 +580,16 @@ def _exec(self, data, detectors=None, **kwargs):
         save_tmpl_flags = self.template_matrix.det_flags
         save_tmpl_mask = self.template_matrix.det_flag_mask
 
+        # The pointing matrix used for the solve.  The per-detector flags
+        # are normally reset when the binner is run, but here we set them
+        # explicitly since we will use these pointing matrix operators for
+        # setting up the solver flags below.
+        solve_pixels = self.binning.pixel_pointing
+        solve_weights = self.binning.stokes_weights
+        solve_pixels.detector_pointing.det_flag_mask = save_det_flag_mask
+        if hasattr(solve_weights, "detector_pointing"):
+            solve_weights.detector_pointing.det_flag_mask = save_det_flag_mask
+
         # Output data products, prefixed with the name of the operator and optionally
         # the MC index.
 
@@ -637,7 +654,7 @@ def _exec(self, data, detectors=None, **kwargs):
                     continue
                 # Create the new solver flags
                 exists = ob.detdata.ensure(
-                    self.solver_flags, dtype=np.uint8, detectors=detectors
+                    self.solver_flags, dtype=np.uint8, detectors=dets
                 )
                 # The data views
                 views = ob.view[solve_view]
@@ -681,21 +698,18 @@ def _exec(self, data, detectors=None, **kwargs):
             # pointing, we want to do that later when building the covariance and
             # the pixel distribution.
 
-            # Use the same pointing operator as the binning
-            scan_pointing = self.binning.pixel_pointing
-
             scanner = ScanMask(
                 det_flags=self.solver_flags,
-                pixels=scan_pointing.pixels,
+                det_flag_mask=save_det_flag_mask,
+                pixels=solve_pixels.pixels,
                 view=solve_view,
-                # mask_bits=1,
             )
 
             scanner.det_flags_value = 2
             scanner.mask_key = self.mask
 
             scan_pipe = Pipeline(
-                detector_sets=["SINGLE"], operators=[scan_pointing, scanner]
+                detector_sets=["SINGLE"], operators=[solve_pixels, scanner]
             )
 
             if self.mask is not None:
@@ -743,8 +757,8 @@ def _exec(self, data, detectors=None, **kwargs):
                 det_data_units=det_data_units,
                 det_flags=self.solver_flags,
                 det_flag_mask=255,
-                pixel_pointing=self.binning.pixel_pointing,
-                stokes_weights=self.binning.stokes_weights,
+                pixel_pointing=solve_pixels,
+                stokes_weights=solve_weights,
                 noise_model=self.binning.noise_model,
                 rcond_threshold=self.solve_rcond_threshold,
                 sync_type=self.binning.sync_type,
diff --git a/src/toast/ops/mapmaker_utils/mapmaker_utils.py b/src/toast/ops/mapmaker_utils/mapmaker_utils.py
index 40d522299..df2601417 100644
--- a/src/toast/ops/mapmaker_utils/mapmaker_utils.py
+++ b/src/toast/ops/mapmaker_utils/mapmaker_utils.py
@@ -542,7 +542,7 @@ class BuildNoiseWeighted(Operator):
     If any samples have compromised telescope pointing, those pixel indices should
     have already been set to a negative value by the operator that generated the
     pointing matrix.  Individual detector flags can optionally be applied to
-    timesamples when accumulating data.  The detector mask defaults to cutting 
+    timesamples when accumulating data.  The detector mask defaults to cutting
     "non-science" samples.
 
     """
@@ -1083,6 +1083,11 @@ def _exec(self, data, detectors=None, **kwargs):
                 msg = f"You must set the '{trait}' trait before calling exec()"
                 raise RuntimeError(msg)
 
+        # Set pointing flags
+        self.pixel_pointing.detector_pointing.det_flag_mask = self.det_flag_mask
+        if hasattr(self.stokes_weights, "detector_pointing"):
+            self.stokes_weights.detector_pointing.det_flag_mask = self.det_flag_mask
+
         # Construct the pointing distribution if it does not already exist
 
         if self.pixel_dist not in data:
diff --git a/src/toast/ops/noise_weight/noise_weight.py b/src/toast/ops/noise_weight/noise_weight.py
index 7fb03e418..374237f8f 100644
--- a/src/toast/ops/noise_weight/noise_weight.py
+++ b/src/toast/ops/noise_weight/noise_weight.py
@@ -7,6 +7,7 @@
 
 from ...accelerator import ImplementationType
 from ...noise_sim import AnalyticNoise
+from ...observation import default_values as defaults
 from ...timing import Timer, function_timer
 from ...traits import Int, Unicode, UseEnum, trait_docs
 from ...utils import Environment, Logger
@@ -40,6 +41,11 @@ class NoiseWeight(Operator):
         None, allow_none=True, help="Observation detdata key for the timestream data"
     )
 
+    det_flag_mask = Int(
+        defaults.det_mask_invalid,
+        help="Bit mask value for optional detector flagging",
+    )
+
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
@@ -51,17 +57,17 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
         implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
         for ob in data.obs:
-            # Get the detectors we are using for this observation.  Since
-            # we are modifying the timestream units, we consider all detectors.
-            dets = ob.select_local_detectors(detectors, flagmask=0)
-            if len(dets) == 0:
-                # Nothing to do for this observation
-                continue
-
             data_input_units = ob.detdata[self.det_data].units
             data_invcov_units = 1.0 / data_input_units**2
             data_output_units = 1.0 / data_input_units
 
+            dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
+            if len(dets) == 0:
+                # Nothing to do for this observation, but
+                # update the units of the output
+                ob.detdata[self.det_data].update_units(data_output_units)
+                continue
+
             # Check that the noise model exists
             if self.noise_model not in ob:
                 msg = "Noise model {} does not exist in observation {}".format(
diff --git a/src/toast/ops/operator.py b/src/toast/ops/operator.py
index 2510b5302..486a76421 100644
--- a/src/toast/ops/operator.py
+++ b/src/toast/ops/operator.py
@@ -218,10 +218,3 @@ def translate(cls, props):
         if "API" in props:
             del props["API"]
         return props
-
-    def __str__(self):
-        """
-        Converts the operator into a short human readable string
-        returns only its class name (for example: 'Operator')
-        """
-        return str(self.__class__).split(".")[-1][:-2]
diff --git a/src/toast/ops/pipeline.py b/src/toast/ops/pipeline.py
index e0cf8aa5e..769f1683b 100644
--- a/src/toast/ops/pipeline.py
+++ b/src/toast/ops/pipeline.py
@@ -170,7 +170,11 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         # notify user of device->host data movements introduced by CPU operators
         if (self._unstaged_data is not None) and (not self._unstaged_data.is_empty()):
-            cpu_ops = {str(op) for op in self.operators if not op.supports_accel()}
+            cpu_ops = {
+                op.__class__.__qualname__
+                for op in self.operators
+                if not op.supports_accel()
+            }
             log.debug(
                 f"{pstr} {self}, had to move {self._unstaged_data} back to host as {cpu_ops} do not support accel."
             )
@@ -357,4 +361,4 @@ def __str__(self):
         """
         Converts the pipeline into a human-readable string.
         """
-        return f"Pipeline{[str(op) for op in self.operators]}"
+        return f"Pipeline{[op.__class__.__qualname__ for op in self.operators]}"
diff --git a/src/toast/ops/pixels_healpix/pixels_healpix.py b/src/toast/ops/pixels_healpix/pixels_healpix.py
index 19c68460d..633258f3d 100644
--- a/src/toast/ops/pixels_healpix/pixels_healpix.py
+++ b/src/toast/ops/pixels_healpix/pixels_healpix.py
@@ -162,7 +162,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(
+                detectors, flagmask=self.detector_pointing.det_flag_mask
+            )
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
@@ -222,7 +224,9 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                         # this calculation.  This could eventually be a kernel.
                         ob.detdata[self.pixels].accel_update_host()
                         restore_dev = True
-                    for det in ob.select_local_detectors(detectors):
+                    for det in ob.select_local_detectors(
+                        detectors, flagmask=self.detector_pointing.det_flag_mask
+                    ):
                         for vslice in view_slices:
                             good = ob.detdata[self.pixels][det, vslice] >= 0
                             self._local_submaps[
diff --git a/src/toast/ops/pixels_wcs.py b/src/toast/ops/pixels_wcs.py
index ee66f9225..38287a60c 100644
--- a/src/toast/ops/pixels_wcs.py
+++ b/src/toast/ops/pixels_wcs.py
@@ -395,7 +395,9 @@ def _exec(self, data, detectors=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(
+                detectors, flagmask=self.detector_pointing.det_flag_mask
+            )
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
@@ -442,7 +444,9 @@ def _exec(self, data, detectors=None, **kwargs):
                         # this calculation.  This could eventually be a kernel.
                         ob.detdata[self.pixels].accel_update_host()
                         restore_dev = True
-                    for det in ob.select_local_detectors(detectors):
+                    for det in ob.select_local_detectors(
+                        detectors, flagmask=self.detector_pointing.det_flag_mask
+                    ):
                         for vslice in view_slices:
                             good = ob.detdata[self.pixels][det, vslice] >= 0
                             self._local_submaps[
@@ -475,7 +479,9 @@ def _exec(self, data, detectors=None, **kwargs):
                 center_lonlat = ob.shared[self.center_offset].data
 
             # Process all detectors
-            for det in ob.select_local_detectors(detectors):
+            for det in ob.select_local_detectors(
+                detectors, flagmask=self.detector_pointing.det_flag_mask
+            ):
                 for vslice in view_slices:
                     # Timestream of detector quaternions
                     quats = ob.detdata[quats_name][det][vslice]
diff --git a/src/toast/ops/pointing_detector/pointing_detector.py b/src/toast/ops/pointing_detector/pointing_detector.py
index 74df6b60a..045552cb8 100644
--- a/src/toast/ops/pointing_detector/pointing_detector.py
+++ b/src/toast/ops/pointing_detector/pointing_detector.py
@@ -37,6 +37,11 @@ class PointingDetectorSimple(Operator):
         defaults.shared_mask_invalid, help="Bit mask value for optional flagging"
     )
 
+    det_flag_mask = Int(
+        defaults.det_mask_invalid,
+        help="Bit mask value for optional detector flagging",
+    )
+
     boresight = Unicode(
         defaults.boresight_radec, help="Observation shared key for boresight"
     )
@@ -119,7 +124,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
diff --git a/src/toast/ops/polyfilter/polyfilter.py b/src/toast/ops/polyfilter/polyfilter.py
index 3617f5779..16278d14b 100644
--- a/src/toast/ops/polyfilter/polyfilter.py
+++ b/src/toast/ops/polyfilter/polyfilter.py
@@ -546,27 +546,29 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 shared_flags = np.zeros(obs.n_local_samples, dtype=np.uint8)
 
             signals = []
-            idets = []
+            filter_dets = []
             last_flags = None
             in_place = True
-            for idet, det in enumerate(dets):
+            for det in dets:
                 # Test the detector pattern
                 if pat.match(det) is None:
                     continue
 
-                ref = obs.detdata[self.det_data][idet]
-                signal = ref.astype(np.float64, copy=False)
-                if not signal is ref:
+                ref = obs.detdata[self.det_data][det]
+                if isinstance(ref[0], np.float64):
+                    signal = ref
+                else:
                     in_place = False
+                    signal = np.array(ref, dtype=np.float64)
                 if self.det_flags is not None:
-                    det_flags = obs.detdata[self.det_flags][idet] & self.det_flag_mask
+                    det_flags = obs.detdata[self.det_flags][det] & self.det_flag_mask
                     flags = shared_flags | det_flags
                 else:
                     flags = shared_flags
 
                 if last_flags is None or np.all(last_flags == flags):
+                    filter_dets.append(det)
                     signals.append(signal)
-                    idets.append(idet)
                 else:
                     filter_polynomial(
                         self.order,
@@ -578,10 +580,10 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                         use_accel=use_accel,
                     )
                     if not in_place:
-                        for i, x in zip(idets, signals):
-                            obs.detdata[self.det_data][i] = x
+                        for fdet, x in zip(filter_dets, signals):
+                            obs.detdata[self.det_data][fdet] = x
                     signals = [signal]
-                    idets = [idet]
+                    filter_dets = [det]
                 last_flags = flags.copy()
 
             if len(signals) > 0:
@@ -595,8 +597,8 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                     use_accel=use_accel,
                 )
                 if not in_place:
-                    for i, x in zip(idets, signals):
-                        obs.detdata[self.det_data][i] = x
+                    for fdet, x in zip(filter_dets, signals):
+                        obs.detdata[self.det_data][fdet] = x
 
             # Optionally flag unfiltered data
             if self.shared_flags is not None and self.poly_flag_mask is not None:
@@ -807,7 +809,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             # Loop over all values of the focalplane key
             for value in values:
                 local_dets = []
-                for idet, det in enumerate(temp_ob.local_detectors):
+                for det in temp_ob.local_detectors:
                     if temp_ob.local_detector_flags[det] & self.det_flag_mask:
                         continue
                     if pat.match(det) is None:
@@ -817,8 +819,12 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                         and focalplane[det][self.focalplane_key] != value
                     ):
                         continue
-                    local_dets.append(idet)
-                local_dets = np.array(local_dets)
+                    local_dets.append(det)
+
+                # The indices into the detector data, which may be different than
+                # the index into the full set of local detectors.
+                data_indices = temp_ob.detdata[self.det_data].indices(local_dets)
+                flag_indices = temp_ob.detdata[self.det_flags].indices(local_dets)
 
                 # Average all detectors that match the key
                 template = np.zeros(nsample)
@@ -834,7 +840,8 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                     det_flags = np.zeros([ndet, nsample], dtype=np.uint8)
 
                 sum_detectors(
-                    local_dets,
+                    data_indices,
+                    flag_indices,
                     shared_flags,
                     self.shared_flag_mask,
                     temp_ob.detdata[self.det_data].data,
@@ -850,7 +857,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                     comm.Allreduce(MPI.IN_PLACE, hits, op=MPI.SUM)
 
                 subtract_mean(
-                    local_dets,
+                    data_indices,
                     temp_ob.detdata[self.det_data].data,
                     template,
                     hits,
diff --git a/src/toast/ops/scan_map/scan_map.py b/src/toast/ops/scan_map/scan_map.py
index 9ece11ab8..70e35a12b 100644
--- a/src/toast/ops/scan_map/scan_map.py
+++ b/src/toast/ops/scan_map/scan_map.py
@@ -38,6 +38,11 @@ class ScanMap(Operator):
         defaults.det_data_units, help="Output units if creating detector data"
     )
 
+    det_flag_mask = Int(
+        defaults.det_mask_invalid,
+        help="Bit mask value for optional detector flagging",
+    )
+
     view = Unicode(
         None, allow_none=True, help="Use this view of the data in all observations"
     )
@@ -90,7 +95,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
@@ -209,8 +214,13 @@ class ScanMask(Operator):
     )
 
     det_flags_value = Int(
-        defaults.det_mask_processing, 
-        help="The detector flag value to set where the mask result is non-zero"
+        defaults.det_mask_processing,
+        help="The detector flag value to set where the mask result is non-zero",
+    )
+
+    det_flag_mask = Int(
+        defaults.det_mask_invalid,
+        help="Bit mask value for optional detector flagging",
     )
 
     view = Unicode(
@@ -257,7 +267,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(detectors, self.det_flag_mask)
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
@@ -324,6 +334,11 @@ class ScanScale(Operator):
         None, allow_none=True, help="Observation detdata key for the timestream data"
     )
 
+    det_flag_mask = Int(
+        defaults.det_mask_invalid,
+        help="Bit mask value for optional detector flagging",
+    )
+
     view = Unicode(
         None, allow_none=True, help="Use this view of the data in all observations"
     )
@@ -373,7 +388,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
diff --git a/src/toast/ops/sim_crosstalk.py b/src/toast/ops/sim_crosstalk.py
index 877d09e7c..228d75e15 100644
--- a/src/toast/ops/sim_crosstalk.py
+++ b/src/toast/ops/sim_crosstalk.py
@@ -125,7 +125,8 @@ def invert_xtalk_mat(matdic):
 
 @trait_docs
 class CrossTalk(Operator):
-    """
+    """Simulate readout crosstalk between channels.
+
     1.  The cross talk matrix can just be a dictionary of
     dictionaries of values (i.e. a sparse matrix) on every process.
     It does not need to be a dense matrix loaded from an HDF5 file.
@@ -138,6 +139,7 @@ class CrossTalk(Operator):
     accumulates to its local detectors, and passes it along.
     This continues until every process has accumulated the data
     from the other processes in the column.
+
     """
 
     # Class traits
diff --git a/src/toast/ops/sim_ground_utils.py b/src/toast/ops/sim_ground_utils.py
index e581f85cb..a86e61d97 100644
--- a/src/toast/ops/sim_ground_utils.py
+++ b/src/toast/ops/sim_ground_utils.py
@@ -8,9 +8,9 @@
 import numpy as np
 from astropy import units as u
 
+from ..coordinates import to_DJD
 from ..intervals import IntervalList
 from ..timing import Timer, function_timer
-from ..coordinates import to_DJD
 
 
 @function_timer
diff --git a/src/toast/ops/sim_hwp.py b/src/toast/ops/sim_hwp.py
index 3ccc01867..8cd9a3fb3 100644
--- a/src/toast/ops/sim_hwp.py
+++ b/src/toast/ops/sim_hwp.py
@@ -9,8 +9,8 @@
 from astropy import units as u
 
 from .. import rng
-from ..observation import default_values as defaults
 from ..mpi import MPI
+from ..observation import default_values as defaults
 from ..timing import Timer, function_timer
 from ..traits import Int, Quantity, Unicode, trait_docs
 from ..utils import GlobalTimers, Logger, Timer, dtype_to_aligned, name_UID
diff --git a/src/toast/ops/sim_tod_atm.py b/src/toast/ops/sim_tod_atm.py
index 8adb34e23..2a45fd7b1 100644
--- a/src/toast/ops/sim_tod_atm.py
+++ b/src/toast/ops/sim_tod_atm.py
@@ -243,7 +243,6 @@ def _check_detector_weights(self, proposal):
             # Check that this operator has the traits we expect
             for trt in [
                 "view",
-                "quats",
                 "weights",
                 "mode",
             ]:
diff --git a/src/toast/ops/sss.py b/src/toast/ops/sss.py
index 07be9ef1f..faa7944ce 100644
--- a/src/toast/ops/sss.py
+++ b/src/toast/ops/sss.py
@@ -97,7 +97,9 @@ def _exec(self, data, detectors=None, **kwargs):
         comm = data.comm.comm_group
 
         for obs in data.obs:
-            dets = obs.select_local_detectors(detectors)
+            dets = obs.select_local_detectors(
+                detectors, flagmask=defaults.det_mask_invalid
+            )
             log_prefix = f"{group} : {obs.name} : "
 
             exists = obs.detdata.ensure(
diff --git a/src/toast/ops/statistics.py b/src/toast/ops/statistics.py
index 5502c6de6..a258174b7 100644
--- a/src/toast/ops/statistics.py
+++ b/src/toast/ops/statistics.py
@@ -106,14 +106,29 @@ def _exec(self, data, detectors=None, **kwargs):
                 fname_out = f"{self.name}_{obs.name}.h5"
             if self.output_dir is not None:
                 fname_out = os.path.join(self.output_dir, fname_out)
-            all_dets = list(obs.all_detectors)
-            ndet = len(all_dets)
 
+            # Get the list of all detectors that are not cut
+            obs_dets = obs.select_local_detectors(
+                detectors, flagmask=self.det_flag_mask
+            )
+            if obs.comm.group_size == 1:
+                all_dets = obs_dets
+            else:
+                proc_dets = obs.comm.comm_group.gather(obs_dets, root=0)
+                all_dets = None
+                if obs.comm.group_rank == 0:
+                    all_set = set()
+                    for pdets in proc_dets:
+                        for d in pdets:
+                            all_set.add(d)
+                    all_dets = list(sorted(all_set))
+                all_dets = obs.comm.comm_group.bcast(all_dets, root=0)
+
+            ndet = len(all_dets)
             hits = np.zeros([ndet], dtype=int)
             means = np.zeros([ndet], dtype=float)
             stats = np.zeros([nstat, ndet], dtype=float)
 
-            obs_dets = obs.select_local_detectors(detectors)
             views = obs.view[self.view]
 
             # Measure the mean separately to simplify the math
diff --git a/src/toast/ops/stokes_weights/stokes_weights.py b/src/toast/ops/stokes_weights/stokes_weights.py
index ffdd9ec6d..1a09d3729 100644
--- a/src/toast/ops/stokes_weights/stokes_weights.py
+++ b/src/toast/ops/stokes_weights/stokes_weights.py
@@ -89,12 +89,6 @@ class StokesWeights(Operator):
         defaults.weights, help="Observation detdata key for output weights"
     )
 
-    quats = Unicode(
-        None,
-        allow_none=True,
-        help="Observation detdata key for output quaternions",
-    )
-
     single_precision = Bool(False, help="If True, use 32bit float in output")
 
     cal = Unicode(
@@ -157,13 +151,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                 raise RuntimeError("If using HWP, you must specify the fp_gamma key")
 
         # Expand detector pointing
-        if self.quats is not None:
-            quats_name = self.quats
-        else:
-            if self.detector_pointing.quats is not None:
-                quats_name = self.detector_pointing.quats
-            else:
-                quats_name = "quats"
+        quats_name = self.detector_pointing.quats
 
         view = self.view
         if view is None:
@@ -171,12 +159,13 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
             view = self.detector_pointing.view
 
         # Expand detector pointing
-        self.detector_pointing.quats = quats_name
         self.detector_pointing.apply(data, detectors=detectors, use_accel=use_accel)
 
         for ob in data.obs:
             # Get the detectors we are using for this observation
-            dets = ob.select_local_detectors(detectors)
+            dets = ob.select_local_detectors(
+                detectors, flagmask=self.detector_pointing.det_flag_mask
+            )
             if len(dets) == 0:
                 # Nothing to do for this observation
                 continue
diff --git a/src/toast/templates/fourier2d.py b/src/toast/templates/fourier2d.py
index bf1783cd4..19bce2bdc 100644
--- a/src/toast/templates/fourier2d.py
+++ b/src/toast/templates/fourier2d.py
@@ -124,6 +124,8 @@ def _initialize(self, new_data):
 
             # Build up detector list
             for d in ob.local_detectors:
+                if ob.local_detector_flags[d] & self.det_flag_mask:
+                    continue
                 if d not in all_dets:
                     all_dets[d] = None
 
@@ -265,6 +267,8 @@ def evaluate_template(theta, phi, radius):
                 norms_view = self._norms[norm_slice].reshape((-1, self._nmode))
 
                 for det in ob.local_detectors:
+                    if det not in self._all_dets:
+                        continue
                     detweight = 1.0
                     if noise is not None:
                         detweight = noise.detector_weight(det).to_value(invvar_units)
@@ -327,6 +331,9 @@ def _zeros(self):
         return z
 
     def _add_to_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         for iob, ob in enumerate(self.data.obs):
             if detector not in ob.local_detectors:
                 continue
@@ -347,6 +354,9 @@ def _add_to_signal(self, detector, amplitudes, **kwargs):
                 )
 
     def _project_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         for iob, ob in enumerate(self.data.obs):
             if detector not in ob.local_detectors:
                 continue
diff --git a/src/toast/templates/gaintemplate.py b/src/toast/templates/gaintemplate.py
index e198294ba..e49193ddd 100644
--- a/src/toast/templates/gaintemplate.py
+++ b/src/toast/templates/gaintemplate.py
@@ -65,9 +65,11 @@ def _initialize(self, new_data):
         # but sorted in order of occurrence.
         all_dets = OrderedDict()
 
+        # Build up detector list
         for iob, ob in enumerate(new_data.obs):
-            # Build up detector list
             for d in ob.local_detectors:
+                if ob.local_detector_flags[d] & self.det_flag_mask:
+                    continue
                 if d not in all_dets:
                     all_dets[d] = None
 
@@ -134,6 +136,8 @@ def _initialize(self, new_data):
 
                 self._precond[iob][ivw] = dict()
                 for det in ob.local_detectors:
+                    if det not in self._all_dets:
+                        continue
                     detweight = 1.0
                     if noise is not None:
                         detweight = noise.detector_weight(det).to_value(invvar_units)
@@ -162,11 +166,16 @@ def _zeros(self):
         return z
 
     def _add_to_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         norder = self.order + 1
         offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
             if detector not in ob.local_detectors:
                 continue
+            if detector not in ob.detdata[self.det_data].detectors:
+                continue
             for ivw, vw in enumerate(ob.view[self.view].detdata[self.det_data]):
                 legendre_poly = self._templates[iob][ivw]
                 poly_amps = amplitudes.local[offset : offset + norder]
@@ -176,11 +185,16 @@ def _add_to_signal(self, detector, amplitudes, **kwargs):
                 vw[detector] += gain_fluctuation
 
     def _project_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         norder = self.order + 1
         offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
             if detector not in ob.local_detectors:
                 continue
+            if detector not in ob.detdata[self.det_data].detectors:
+                continue
             for ivw, vw in enumerate(ob.view[self.view].detdata[self.det_data]):
                 legendre_poly = self._templates[iob][ivw]
                 signal_estimate = ob.detdata[self.template_name][detector]
diff --git a/src/toast/templates/offset/kernels_jax.py b/src/toast/templates/offset/kernels_jax.py
index e0c01ddf9..c803b66da 100644
--- a/src/toast/templates/offset/kernels_jax.py
+++ b/src/toast/templates/offset/kernels_jax.py
@@ -474,7 +474,9 @@ def offset_apply_diag_precond_inner(offset_var, amplitudes_in, amplitudes_out):
 
 
 @kernel(impl=ImplementationType.JAX, name="offset_apply_diag_precond")
-def offset_apply_diag_precond_jax(offset_var, amplitudes_in, amplitude_flags, amplitudes_out, use_accel):
+def offset_apply_diag_precond_jax(
+    offset_var, amplitudes_in, amplitude_flags, amplitudes_out, use_accel
+):
     """
     Simple multiplication.
 
diff --git a/src/toast/templates/offset/offset.py b/src/toast/templates/offset/offset.py
index 842322ff1..65c7a7bdc 100644
--- a/src/toast/templates/offset/offset.py
+++ b/src/toast/templates/offset/offset.py
@@ -3,6 +3,7 @@
 # a BSD-style license that can be found in the LICENSE file.
 
 import json
+import os
 from collections import OrderedDict
 
 import h5py
@@ -71,6 +72,12 @@ class Offset(Template):
 
     precond_width = Int(20, help="Preconditioner width in terms of offsets / baselines")
 
+    debug_plots = Unicode(
+        None,
+        allow_none=True,
+        help="If not None, make debugging plots in this directory",
+    )
+
     @traitlets.validate("precond_width")
     def _check_precond_width(self, proposal):
         w = proposal["value"]
@@ -167,13 +174,16 @@ def _initialize(self, new_data):
                 if self.use_noise_prior:
                     obstime = ob.shared[self.times][-1] - ob.shared[self.times][0]
                     tbase = self.step_time.to_value(u.second)
-                    fbase = 1.0 / tbase
                     powmin = np.floor(np.log10(1 / obstime)) - 1
-                    powmax = min(np.ceil(np.log10(1 / tbase)) + 2, self._obs_rate[iob])
+                    powmax = min(
+                        np.ceil(np.log10(1 / tbase)) + 2, np.log10(self._obs_rate[iob])
+                    )
                     self._freq[iob] = np.logspace(powmin, powmax, 1000)
 
             # Build up detector list
             for d in ob.local_detectors:
+                if ob.local_detector_flags[d] & self.det_flag_mask:
+                    continue
                 if d not in all_dets:
                     all_dets[d] = None
 
@@ -264,7 +274,9 @@ def _initialize(self, new_data):
                                 amplen = step_length
                                 if amp == n_amp_view - 1:
                                     amplen = view_samples - voff
-                                self._offsetvar[offset + amp] = 1.0 / (detnoise * amplen)
+                                self._offsetvar[offset + amp] = 1.0 / (
+                                    detnoise * amplen
+                                )
                                 voff += step_length
                         else:
                             flags = views.detdata[self.det_flags][ivw]
@@ -274,9 +286,10 @@ def _initialize(self, new_data):
                                 if amp == n_amp_view - 1:
                                     amplen = view_samples - voff
                                 n_good = amplen - np.count_nonzero(
-                                    flags[det][voff : voff + amplen] & self.det_flag_mask
+                                    flags[det][voff : voff + amplen]
+                                    & self.det_flag_mask
                                 )
-                                if ((n_good / amplen) <= self.good_fraction):
+                                if (n_good / amplen) <= self.good_fraction:
                                     # This detector is cut or too many samples flagged
                                     self._offsetvar[offset + amp] = 0.0
                                     self._amp_flags[offset + amp] = True
@@ -317,6 +330,51 @@ def _initialize(self, new_data):
                         det,
                     )
 
+                    if self.debug_plots is not None:
+                        set_matplotlib_backend()
+                        import matplotlib.pyplot as plt
+
+                        fname = os.path.join(
+                            self.debug_plots, f"{self.name}_{det}_{ob.name}_psd.pdf"
+                        )
+                        psdfreq = ob[self.noise_model].freq(det).to_value(u.Hz)
+                        psd = (
+                            ob[self.noise_model]
+                            .psd(det)
+                            .to_value(self.det_data_units**2 * u.second)
+                        )
+                        corrpsd = self._remove_white_noise(psdfreq, psd)
+
+                        fig = plt.figure(figsize=[12, 12])
+                        ax = fig.add_subplot(2, 1, 1)
+                        ax.loglog(
+                            psdfreq,
+                            psd,
+                            color="black",
+                            label="Original PSD",
+                        )
+                        ax.loglog(
+                            psdfreq,
+                            corrpsd,
+                            color="red",
+                            label="Correlated PSD",
+                        )
+                        ax.set_xlabel("Frequency [Hz]")
+                        ax.set_ylabel("PSD [K$^2$ / Hz]")
+                        ax.legend(loc="best")
+
+                        ax = fig.add_subplot(2, 1, 2)
+                        ax.loglog(
+                            self._freq[iob],
+                            offset_psd,
+                            label=f"Offset PSD",
+                        )
+                        ax.set_xlabel("Frequency [Hz]")
+                        ax.set_ylabel("PSD [K$^2$ / Hz]")
+                        ax.legend(loc="best")
+                        fig.savefig(fname)
+                        plt.close(fig)
+
                     # "Noise weight" (time-domain inverse variance)
                     detnoise = (
                         ob[self.noise_model]
@@ -335,6 +393,18 @@ def _initialize(self, new_data):
                     self._filters[iob][det] = list()
                     self._precond[iob][det] = list()
 
+                    if self.debug_plots is not None:
+                        ffilter = os.path.join(
+                            self.debug_plots, f"{self.name}_{det}_{ob.name}_filters.pdf"
+                        )
+                        fprec = os.path.join(
+                            self.debug_plots, f"{self.name}_{det}_{ob.name}_prec.pdf"
+                        )
+                        figfilter = plt.figure(figsize=[12, 8])
+                        axfilter = figfilter.add_subplot(1, 1, 1)
+                        figprec = plt.figure(figsize=[12, 8])
+                        axprec = figprec.add_subplot(1, 1, 1)
+
                     # Loop over views
                     views = ob.view[self.view]
                     for ivw, vw in enumerate(views):
@@ -370,6 +440,13 @@ def _initialize(self, new_data):
 
                         self._filters[iob][det].append(noisefilter)
 
+                        if self.debug_plots is not None:
+                            axfilter.plot(
+                                np.arange(len(noisefilter)),
+                                noisefilter,
+                                label=f"Noise filter {ivw}",
+                            )
+
                         # Build the preconditioner
                         lower = None
                         preconditioner = None
@@ -388,6 +465,12 @@ def _initialize(self, new_data):
                             icenter = preconditioner.size // 2
                             if detnoise != 0:
                                 preconditioner[icenter] += 1.0 / detnoise
+                            if self.debug_plots is not None:
+                                axprec.plot(
+                                    np.arange(len(preconditioner)),
+                                    preconditioner,
+                                    label=f"Toeplitz preconditioner {ivw}",
+                                )
                         else:
                             # We are using a banded matrix for the preconditioner.
                             # This contains a Toeplitz component from the inverse
@@ -421,17 +504,35 @@ def _initialize(self, new_data):
                             )
                         self._precond[iob][det].append((preconditioner, lower))
                         offset += n_amp_view
+
+                    if self.debug_plots is not None:
+                        axfilter.set_xlabel("Sample Lag")
+                        axfilter.set_ylabel("Amplitude")
+                        axfilter.legend(loc="best")
+                        figfilter.savefig(ffilter)
+                        axprec.set_xlabel("Sample Lag")
+                        axprec.set_ylabel("Amplitude")
+                        axprec.legend(loc="best")
+                        figprec.savefig(fprec)
+                        plt.close(figfilter)
+                        plt.close(figprec)
+
         log.verbose(f"Offset variance = {self._offsetvar}")
         return
 
     # Helper functions for noise / preconditioner calculations
 
     def _interpolate_psd(self, x, lfreq, lpsd):
-        result = np.zeros(x.size)
-        good = np.abs(x) > 1e-10
-        logx = np.log(np.abs(x[good]))
+        # Threshold for zero frequency
+        thresh = 1.0e-6
+        lowf = x < thresh
+        good = np.logical_not(lowf)
+
+        logx = np.empty_like(x)
+        logx[lowf] = np.log(thresh)
+        logx[good] = np.log(x[good])
         logresult = np.interp(logx, lfreq, lpsd)
-        result[good] = np.exp(logresult)
+        result = np.exp(logresult)
         return result
 
     def _truncate(self, noisefilter, lim=1e-4):
@@ -444,6 +545,42 @@ def _truncate(self, noisefilter, lim=1e-4):
         noisefilter = noisefilter[icenter - icut : icenter + icut + 1]
         return noisefilter
 
+    def _remove_white_noise(self, freq, psd):
+        """Remove the white noise component of the PSD."""
+        corrpsd = psd.copy()
+        n_corrpsd = len(corrpsd)
+        plat_off = int(0.8 * n_corrpsd)
+        if n_corrpsd - plat_off < 10:
+            if n_corrpsd < 10:
+                # Crazy spectrum...
+                plat_off = 0
+            else:
+                plat_off = n_corrpsd - 10
+
+        cfreq = np.log(freq[plat_off:])
+        cdata = np.log(corrpsd[plat_off:])
+
+        def lin_func(x, a, b, c):
+            # Line
+            return a * (x - b) + c
+
+        params, params_cov = scipy.optimize.curve_fit(
+            lin_func, cfreq, cdata, p0=[0.0, cfreq[-1], cdata[-1]]
+        )
+
+        cdata = lin_func(cfreq, params[0], params[1], params[2])
+        cdata = np.exp(cdata)
+        plat = cdata[-1]
+
+        # Given the range between the white noise plateau and the maximum
+        # values of the PSD, we set a minimum value for any spectral bins
+        # that are small or negative.
+        corrmax = np.amax(corrpsd)
+        corrthresh = 1.0e-10 * corrmax - plat
+        corrpsd -= plat
+        corrpsd[corrpsd < corrthresh] = corrthresh
+        return corrpsd
+
     def _get_offset_psd(self, noise, freq, step_time, det):
         """Compute the PSD of the baseline offsets."""
         psdfreq = noise.freq(det).to_value(u.Hz)
@@ -451,37 +588,47 @@ def _get_offset_psd(self, noise, freq, step_time, det):
         rate = noise.rate(det).to_value(u.Hz)
 
         # Remove the white noise component from the PSD
-        psd = psd.copy()
-        psd -= np.amin(psd[psdfreq > 1.0])
-        psd[psd < 1e-30] = 1e-30
+        psd = self._remove_white_noise(psdfreq, psd)
 
-        # The calculation of `offset_psd` is from Keihänen, E. et al:
-        # "Making CMB temperature and polarization maps with Madam",
-        # A&A 510:A57, 2010
+        # Log PSD for interpolation
         logfreq = np.log(psdfreq)
         logpsd = np.log(psd)
 
-        def g(x):
-            bad = np.abs(x) < 1e-10
-            good = np.logical_not(bad)
-            arg = np.pi * x[good]
-            result = bad.astype(np.float64)
-            result[good] = (np.sin(arg) / arg) ** 2
-            return result
+        # The calculation of `offset_psd` is based on Keihänen, E. et al:
+        # "Making CMB temperature and polarization maps with Madam",
+        # A&A 510:A57, 2010, with a small algebra correction.
 
+        m_max = 5
         tbase = step_time
         fbase = 1.0 / tbase
 
+        def g(f, m):
+            # The frequencies are constructed without the zero frequency,
+            # so we do not need to handle it here.
+            # result = np.sin(np.pi * f * tbase) ** 2 / (np.pi * (f * tbase + m)) ** 2
+            x = np.pi * (f * tbase + m)
+            bad = np.abs(x) < 1.0e-30
+            good = np.logical_not(bad)
+            result = np.empty_like(x)
+            result[bad] = 1.0
+            result[good] = np.sin(x[good]) ** 2 / x[good] ** 2
+            return result
+
         offset_psd = np.zeros_like(freq)
-        offset_psd[:] += self._interpolate_psd(freq, logfreq, logpsd) * g(freq * tbase)
 
-        for m in range(1, 2):
+        # The m = 0 term
+        offset_psd = self._interpolate_psd(freq, logfreq, logpsd) * g(freq, 0)
+
+        # The remaining terms
+        for m in range(1, m_max):
+            # Positive m
             offset_psd[:] += self._interpolate_psd(
                 freq + m * fbase, logfreq, logpsd
-            ) * g(freq * tbase + m)
+            ) * g(freq, m)
+            # Negative m
             offset_psd[:] += self._interpolate_psd(
                 freq - m * fbase, logfreq, logpsd
-            ) * g(freq * tbase - m)
+            ) * g(freq, -m)
 
         offset_psd *= fbase
         return offset_psd
@@ -501,6 +648,10 @@ def _step_length(self, stime, rate):
     def _add_to_signal(self, detector, amplitudes, use_accel=None, **kwargs):
         log = Logger.get()
 
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
+
         # Kernel selection
         implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
@@ -584,6 +735,10 @@ def _add_to_signal(self, detector, amplitudes, use_accel=None, **kwargs):
     def _project_signal(self, detector, amplitudes, use_accel=None, **kwargs):
         log = Logger.get()
 
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
+
         # Kernel selection
         implementation, use_accel = self.select_kernels(use_accel=use_accel)
 
@@ -656,13 +811,15 @@ def _add_prior(self, amplitudes_in, amplitudes_out, use_accel=None, **kwargs):
             raise NotImplementedError(
                 "offset template add_prior on accelerator not implemented"
             )
+        if self.debug_plots is not None:
+            set_matplotlib_backend()
+            import matplotlib.pyplot as plt
+
         for det in self._all_dets:
             offset = self._det_start[det]
             for iob, ob in enumerate(self.data.obs):
                 if det not in ob.local_detectors:
                     continue
-                if det not in ob.detdata[self.det_data].detectors:
-                    continue
                 for ivw, vw in enumerate(ob.view[self.view].detdata[self.det_data]):
                     n_amp_view = self._obs_views[iob][ivw]
                     amp_slice = slice(offset, offset + n_amp_view, 1)
@@ -672,12 +829,51 @@ def _add_prior(self, amplitudes_in, amplitudes_out, use_accel=None, **kwargs):
                     if det in self._filters[iob]:
                         # There is some contribution from this detector
                         amps_out[:] += scipy.signal.convolve(
-                            amps_in, self._filters[iob][det][ivw], mode="same", method="direct"
+                            amps_in,
+                            self._filters[iob][det][ivw],
+                            mode="same",
+                            method="direct",
                         )
+
+                        if self.debug_plots is not None:
+                            # Find the first unused file name in the sequence
+                            iter = -1
+                            while iter < 0 or os.path.isfile(fname):
+                                iter += 1
+                                fname = os.path.join(
+                                    self.debug_plots,
+                                    f"{self.name}_{det}_{ob.name}_prior_{ivw}_{iter}.pdf",
+                                )
+                            fig = plt.figure(figsize=[12, 8])
+                            ax = fig.add_subplot(1, 1, 1)
+                            ax.plot(
+                                np.arange(len(amps_in)),
+                                amps_in,
+                                color="black",
+                                label="Input Amplitudes",
+                            )
+                            ax.plot(
+                                np.arange(len(amps_in)),
+                                amps_out,
+                                color="red",
+                                label="Output Amplitudes",
+                            )
+
                         amps_out[amp_flags_in != 0] = 0.0
-                        print(f"DBG filter: {self._filters[iob][det][ivw]}")
-                        print(f"DBG amps_in: {amps_in}")
-                        print(f"DBG amps_out: {amps_out}", flush=True)
+
+                        if self.debug_plots is not None:
+                            ax.plot(
+                                np.arange(len(amps_in)),
+                                amps_out,
+                                color="green",
+                                label="Output Amplitudes (flagged)",
+                            )
+                            ax.set_xlabel("Amplitude Index")
+                            ax.set_ylabel("Value")
+                            ax.legend(loc="best")
+                            fig.savefig(fname)
+                            plt.close(fig)
+
                     else:
                         amps_out[:] = 0.0
                     offset += n_amp_view
@@ -721,7 +917,9 @@ def _apply_precond(self, amplitudes_in, amplitudes_out, use_accel=None, **kwargs
                                 # scipy.signal.convolve will use either `convolve` or
                                 # `fftconvolve` depending on the size of the inputs
                                 amps_out = scipy.signal.convolve(
-                                    amps_in, self._precond[iob][det][ivw][0], mode="same"
+                                    amps_in,
+                                    self._precond[iob][det][ivw][0],
+                                    mode="same",
                                 )
                             else:
                                 # Use pre-computed Cholesky decomposition.  Note that this
diff --git a/src/toast/templates/periodic.py b/src/toast/templates/periodic.py
index 56dcc95c4..53478f829 100644
--- a/src/toast/templates/periodic.py
+++ b/src/toast/templates/periodic.py
@@ -154,8 +154,11 @@ def _initialize(self, new_data):
             total_bins += obins
             self._obs_nbins.append(obins)
             self._obs_incr.append(oincr)
+
             # Build up detector list
             for d in ob.local_detectors:
+                if ob.local_detector_flags[d] & self.det_flag_mask:
+                    continue
                 if d not in all_dets:
                     all_dets[d] = None
 
@@ -228,6 +231,10 @@ def _initialize(self, new_data):
                 det_indx = ob.detdata[self.det_data].indices([det])[0]
                 amp_hits = self._amp_hits[amp_offset : amp_offset + nbins]
                 amp_flags = self._amp_flags[amp_offset : amp_offset + nbins]
+                if self.det_flags is not None:
+                    flag_indx = ob.detdata[self.det_flags].indices([det])[0]
+                else:
+                    flag_indx = None
                 for vw in ob.intervals[self.view].data:
                     vw_slc = slice(vw.first, vw.last + 1, 1)
                     if self.is_detdata_key:
@@ -239,6 +246,7 @@ def _initialize(self, new_data):
                         iob,
                         ob,
                         vw,
+                        flag_indx=flag_indx,
                         det_flags=True,
                     )
                     np.add.at(
@@ -259,7 +267,9 @@ def _zeros(self):
         z.local_flags[:] = np.where(self._amp_flags, 1, 0)
         return z
 
-    def _view_flags_and_index(self, det_indx, ob_indx, ob, view, det_flags=False):
+    def _view_flags_and_index(
+        self, det_indx, ob_indx, ob, view, flag_indx=None, det_flags=False
+    ):
         """Get the flags and amplitude indices for one detector and view."""
         vw_slc = slice(view.first, view.last + 1, 1)
         vw_len = view.last - view.first + 1
@@ -281,7 +291,7 @@ def _view_flags_and_index(self, det_indx, ob_indx, ob, view, det_flags=False):
                 bad = np.zeros(vw_len, dtype=bool)
         if det_flags and self.det_flags is not None:
             # We have some det flags
-            bad |= ob.detdata[self.det_flags][det_indx, vw_slc] & self.det_flag_mask
+            bad |= ob.detdata[self.det_flags][flag_indx, vw_slc] & self.det_flag_mask
         good = np.logical_not(bad)
 
         # Find the amplitude index for every good sample
@@ -295,6 +305,10 @@ def _view_flags_and_index(self, det_indx, ob_indx, ob, view, det_flags=False):
         return good, amp_indx
 
     def _add_to_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
+
         amp_offset = self._det_offset[detector]
         for iob, ob in enumerate(self.data.obs):
             if self.is_detdata_key:
@@ -306,6 +320,8 @@ def _add_to_signal(self, detector, amplitudes, **kwargs):
 
             if detector not in ob.local_detectors:
                 continue
+            if detector not in ob.detdata[self.det_data].detectors:
+                continue
             nbins = self._obs_nbins[iob]
             det_indx = ob.detdata[self.det_data].indices([detector])[0]
             amps = amplitudes.local[amp_offset : amp_offset + nbins]
@@ -323,6 +339,10 @@ def _add_to_signal(self, detector, amplitudes, **kwargs):
             amp_offset += nbins
 
     def _project_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
+
         amp_offset = self._det_offset[detector]
         for iob, ob in enumerate(self.data.obs):
             if self.is_detdata_key:
@@ -334,9 +354,16 @@ def _project_signal(self, detector, amplitudes, **kwargs):
 
             if detector not in ob.local_detectors:
                 continue
+            if detector not in ob.detdata[self.det_data].detectors:
+                continue
+
             nbins = self._obs_nbins[iob]
             det_indx = ob.detdata[self.det_data].indices([detector])[0]
             amps = amplitudes.local[amp_offset : amp_offset + nbins]
+            if self.det_flags is not None:
+                flag_indx = ob.detdata[self.det_flags].indices([detector])[0]
+            else:
+                flag_indx = None
             for vw in ob.intervals[self.view].data:
                 vw_slc = slice(vw.first, vw.last + 1, 1)
                 good, amp_indx = self._view_flags_and_index(
@@ -344,6 +371,7 @@ def _project_signal(self, detector, amplitudes, **kwargs):
                     iob,
                     ob,
                     vw,
+                    flag_indx=flag_indx,
                     det_flags=True,
                 )
                 # Accumulate to amplitudes
diff --git a/src/toast/templates/subharmonic.py b/src/toast/templates/subharmonic.py
index 946117a31..ade60ccee 100644
--- a/src/toast/templates/subharmonic.py
+++ b/src/toast/templates/subharmonic.py
@@ -58,9 +58,11 @@ def _initialize(self, new_data):
         # but sorted in order of occurrence.
         all_dets = OrderedDict()
 
+        # Build up detector list
         for iob, ob in enumerate(new_data.obs):
-            # Build up detector list
             for d in ob.local_detectors:
+                if ob.local_detector_flags[d] & self.det_flag_mask:
+                    continue
                 if d not in all_dets:
                     all_dets[d] = None
 
@@ -140,6 +142,8 @@ def _initialize(self, new_data):
 
                 self._precond[iob][ivw] = dict()
                 for det in ob.local_detectors:
+                    if det not in self._all_dets:
+                        continue
                     detweight = 1.0
                     if noise is not None:
                         detweight = noise.detector_weight(det).to_value(invvar_units)
@@ -170,6 +174,9 @@ def _zeros(self):
         return z
 
     def _add_to_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         norder = self.order + 1
         offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
@@ -182,11 +189,16 @@ def _add_to_signal(self, detector, amplitudes, **kwargs):
                 offset += norder
 
     def _project_signal(self, detector, amplitudes, **kwargs):
+        if detector not in self._all_dets:
+            # This must have been cut by per-detector flags during initialization
+            return
         norder = self.order + 1
         offset = self._det_start[detector]
         for iob, ob in enumerate(self.data.obs):
             if detector not in ob.local_detectors:
                 continue
+            if detector not in ob.detdata[self.det_data].detectors:
+                continue
             for ivw, vw in enumerate(ob.view[self.view].detdata[self.det_data]):
                 amp_view = amplitudes.local[offset : offset + norder]
                 for order, template in enumerate(self._templates[iob][ivw]):
diff --git a/src/toast/tests/_helpers.py b/src/toast/tests/_helpers.py
index 1da05d288..39d54eaf2 100644
--- a/src/toast/tests/_helpers.py
+++ b/src/toast/tests/_helpers.py
@@ -3,6 +3,7 @@
 # a BSD-style license that can be found in the LICENSE file.
 
 import os
+import re
 import shutil
 import tempfile
 from datetime import datetime
@@ -24,6 +25,7 @@
 from ..observation import DetectorData, Observation
 from ..observation import default_values as defaults
 from ..pixels import PixelData
+from ..pixels_io_healpix import write_healpix_fits
 from ..schedule import GroundSchedule
 from ..schedule_sim_ground import run_scheduler
 from ..schedule_sim_satellite import create_satellite_schedule
@@ -224,6 +226,7 @@ def create_satellite_data(
     pixel_per_process=1,
     hwp_rpm=9.0,
     single_group=False,
+    flagged_pixels=True,
 ):
     """Create a data object with a simple satellite sim.
 
@@ -245,6 +248,10 @@ def create_satellite_data(
     toastcomm = create_comm(mpicomm, single_group=single_group)
     data = Data(toastcomm)
 
+    if flagged_pixels:
+        # We are going to flag half the pixels
+        pixel_per_process *= 2
+
     tele = create_space_telescope(
         toastcomm.group_size,
         sample_rate=sample_rate,
@@ -281,6 +288,17 @@ def create_satellite_data(
     )
     sim_sat.apply(data)
 
+    if flagged_pixels:
+        det_pat = re.compile(r"D(.*)[AB]-.*")
+        for ob in data.obs:
+            det_flags = dict()
+            for det in ob.local_detectors:
+                det_mat = det_pat.match(det)
+                idet = int(det_mat.group(1))
+                if idet % 2 != 0:
+                    det_flags[det] = defaults.det_mask_invalid
+            ob.update_local_detector_flags(det_flags)
+
     return data
 
 
@@ -496,7 +514,7 @@ def create_healpix_ring_satellite(
     return data
 
 
-def create_fake_sky(data, dist_key, map_key):
+def create_fake_sky(data, dist_key, map_key, hpix_out=None):
     np.random.seed(987654321)
     dist = data[dist_key]
     pix_data = PixelData(dist, np.float64, n_value=3, units=defaults.det_data_units)
@@ -512,6 +530,15 @@ def create_fake_sky(data, dist_key, map_key):
             pix_data.data[off, :, 2] = U_data
             off += 1
     data[map_key] = pix_data
+    if hpix_out is not None:
+        write_healpix_fits(
+            pix_data,
+            hpix_out,
+            nest=True,
+            comm_bytes=10000000,
+            report_memory=False,
+            single_precision=False,
+        )
 
 
 def create_fake_mask(data, dist_key, mask_key):
@@ -713,6 +740,7 @@ def create_ground_data(
     split=False,
     turnarounds_invalid=False,
     single_group=False,
+    flagged_pixels=True,
 ):
     """Create a data object with a simple ground sim.
 
@@ -732,6 +760,10 @@ def create_ground_data(
     toastcomm = create_comm(mpicomm, single_group=single_group)
     data = Data(toastcomm)
 
+    if flagged_pixels:
+        # We are going to flag half the pixels
+        pixel_per_process *= 2
+
     tele = create_ground_telescope(
         toastcomm.group_size,
         sample_rate=sample_rate,
@@ -808,6 +840,17 @@ def create_ground_data(
         sim_ground.turnaround_mask = 2
     sim_ground.apply(data)
 
+    if flagged_pixels:
+        det_pat = re.compile(r"D(.*)[AB]-.*")
+        for ob in data.obs:
+            det_flags = dict()
+            for det in ob.local_detectors:
+                det_mat = det_pat.match(det)
+                idet = int(det_mat.group(1))
+                if idet % 2 != 0:
+                    det_flags[det] = defaults.det_mask_invalid
+            ob.update_local_detector_flags(det_flags)
+
     return data
 
 
diff --git a/src/toast/tests/io_compression.py b/src/toast/tests/io_compression.py
index d10bd27bf..997a83776 100644
--- a/src/toast/tests/io_compression.py
+++ b/src/toast/tests/io_compression.py
@@ -198,7 +198,10 @@ def test_roundtrip_detdata(self):
 
     def create_data(self, pixel_per_process=1, single_group=False):
         data = create_ground_data(
-            self.comm, pixel_per_process=pixel_per_process, single_group=single_group
+            self.comm,
+            pixel_per_process=pixel_per_process,
+            single_group=single_group,
+            flagged_pixels=False,
         )
 
         # Simple detector pointing
diff --git a/src/toast/tests/noise.py b/src/toast/tests/noise.py
index ae021994c..184164c7c 100644
--- a/src/toast/tests/noise.py
+++ b/src/toast/tests/noise.py
@@ -15,6 +15,7 @@
 from ..io import H5File
 from ..noise import Noise
 from ..noise_sim import AnalyticNoise
+from ..observation import default_values as defaults
 from ._helpers import close_data, create_outdir, create_satellite_data
 from .mpi import MPITestCase
 from .ops_noise_estim import plot_noise_estim_compare
@@ -117,7 +118,7 @@ def test_noise_fit(self):
         for ob in data.obs:
             in_model = ob[noise_model.noise_model]
             out_model = ob[noise_fitter.out_model]
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 in_psd = in_model.psd(det)
                 fit_psd = out_model.psd(det)
 
diff --git a/src/toast/tests/ops_crosslinking.py b/src/toast/tests/ops_crosslinking.py
index 8d1fe34ab..efe455e63 100644
--- a/src/toast/tests/ops_crosslinking.py
+++ b/src/toast/tests/ops_crosslinking.py
@@ -50,7 +50,8 @@ def test_crosslinking(self):
         # Check that the total number of samples makes sense
         nsample = 0
         for obs in data.obs:
-            nsample += obs.n_local_samples * len(obs.local_detectors)
+            dets = obs.select_local_detectors(flagmask=crosslinking.det_flag_mask)
+            nsample += obs.n_local_samples * len(dets)
         if self.comm is not None:
             nsample = self.comm.reduce(nsample)
 
diff --git a/src/toast/tests/ops_groundfilter.py b/src/toast/tests/ops_groundfilter.py
index b25665acb..17cb314cb 100644
--- a/src/toast/tests/ops_groundfilter.py
+++ b/src/toast/tests/ops_groundfilter.py
@@ -67,7 +67,7 @@ def test_groundfilter(self):
         groundfilter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = ob.shared[defaults.shared_flags].data & self.shared_flag_mask
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
@@ -107,7 +107,7 @@ def test_groundfilter_split(self):
             leftgoing = (shared_flags & defaults.shared_mask_scan_rightleft) != 0
             az = ob.shared[defaults.azimuth].data * 100
             rms[ob.name] = dict()
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = ob.shared[defaults.shared_flags].data & self.shared_flag_mask
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
@@ -125,7 +125,7 @@ def test_groundfilter_split(self):
             split_template=True,
             det_data=defaults.det_data,
             det_flags=defaults.det_flags,
-            det_flag_mask=255,
+            det_flag_mask=defaults.det_mask_invalid,
             shared_flags=defaults.shared_flags,
             shared_flag_mask=self.shared_flag_mask,
             view=None,
@@ -133,7 +133,7 @@ def test_groundfilter_split(self):
         groundfilter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = ob.shared[defaults.shared_flags].data & self.shared_flag_mask
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
diff --git a/src/toast/tests/ops_hwpfilter.py b/src/toast/tests/ops_hwpfilter.py
index 7f8ca50b3..a56fd2c0f 100644
--- a/src/toast/tests/ops_hwpfilter.py
+++ b/src/toast/tests/ops_hwpfilter.py
@@ -61,7 +61,7 @@ def test_hwpfilter(self):
         hwpfilter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=hwpfilter.det_flag_mask):
                 flags = ob.shared[defaults.shared_flags].data & self.shared_flag_mask
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
diff --git a/src/toast/tests/ops_madam.py b/src/toast/tests/ops_madam.py
index 743e459dc..f0405a73c 100644
--- a/src/toast/tests/ops_madam.py
+++ b/src/toast/tests/ops_madam.py
@@ -136,7 +136,7 @@ def test_madam_det_out(self):
         rms = dict()
         for ob in data.obs:
             rms[ob.name] = dict()
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = np.array(ob.shared[defaults.shared_flags])
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
@@ -203,13 +203,6 @@ def test_madam_det_out(self):
             stokes_weights=weights,
             det_out="destriped",
             noise_model="noise_model",
-            copy_groups=2,
-            purge_det_data=False,
-            restore_det_data=False,
-            shared_flags=defaults.shared_flags,
-            shared_flag_mask=1,
-            det_flags=defaults.det_flags,
-            det_flag_mask=1,
         )
         madam.apply(data)
 
@@ -242,7 +235,7 @@ def test_madam_det_out(self):
         #     plt.close()
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = np.array(ob.shared[defaults.shared_flags])
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
diff --git a/src/toast/tests/ops_mapmaker.py b/src/toast/tests/ops_mapmaker.py
index 357387fb3..91bf7758f 100644
--- a/src/toast/tests/ops_mapmaker.py
+++ b/src/toast/tests/ops_mapmaker.py
@@ -7,6 +7,7 @@
 
 import healpy as hp
 import numpy as np
+import scipy
 from astropy import units as u
 
 from .. import ops as ops
@@ -270,6 +271,9 @@ def test_compare_madam_noprior(self):
             output_dir=testdir,
             save_cleaned=True,
             overwrite_cleaned=False,
+            copy_groups=2,
+            purge_det_data=True,
+            restore_det_data=True,
         )
 
         # Make the map
@@ -309,7 +313,7 @@ def test_compare_madam_noprior(self):
         pars["precond_width_max"] = 1
         pars["use_cgprecond"] = "F"
         pars["use_fprecond"] = "F"
-        pars["info"] = 2
+        pars["info"] = 3
         pars["path_output"] = testdir
 
         madam = ops.Madam(
@@ -337,7 +341,7 @@ def test_compare_madam_noprior(self):
         # Compare local destriped TOD on every process
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 input_signal = ob.detdata["input_signal"][det]
                 madam_signal = ob.detdata["madam_cleaned"][det]
                 toast_signal = ob.detdata["toastmap_cleaned"][det]
@@ -504,7 +508,12 @@ def test_compare_madam_diagpre(self):
         weights.apply(data)
 
         # Create fake polarized sky pixel values locally
-        create_fake_sky(data, "pixel_dist", "fake_map")
+        create_fake_sky(
+            data,
+            "pixel_dist",
+            "fake_map",
+            hpix_out=os.path.join(testdir, "input_map.fits"),
+        )
 
         # Scan map into timestreams
         scanner = ops.ScanMap(
@@ -554,6 +563,7 @@ def test_compare_madam_diagpre(self):
             step_time=step_seconds * u.second,
             use_noise_prior=True,
             precond_width=1,
+            debug_plots=testdir,
         )
 
         tmatrix = ops.TemplateMatrix(templates=[tmpl])
@@ -616,7 +626,8 @@ def test_compare_madam_diagpre(self):
         pars["precond_width_min"] = 1
         pars["precond_width_max"] = 1
         pars["use_cgprecond"] = "F"
-        pars["use_fprecond"] = "T"
+        pars["use_fprecond"] = "F"
+        pars["info"] = 3
         pars["path_output"] = testdir
 
         madam = ops.Madam(
@@ -664,14 +675,33 @@ def test_compare_madam_diagpre(self):
 
             # Compare maps
 
+            input_map = hp.read_map(
+                os.path.join(testdir, "input_map.fits"), field=None, nest=True
+            )
             toast_map = hp.read_map(toast_map_path, field=None, nest=True)
             madam_map = hp.read_map(madam_map_path, field=None, nest=True)
+
             # Set madam unhit pixels to zero
             for stokes, ststr in zip(range(3), ["I", "Q", "U"]):
                 mask = hp.mask_bad(madam_map[stokes])
                 madam_map[stokes][mask] = 0.0
+                input_map[stokes][mask] = 0.0
                 diff_map = toast_map[stokes] - madam_map[stokes]
+                madam_diff_truth = madam_map[stokes] - input_map[stokes]
+                toast_diff_truth = toast_map[stokes] - input_map[stokes]
+
                 print("diff map {} has rms {}".format(ststr, np.std(diff_map)))
+                print(
+                    "toast-truth map {} has rms {}".format(
+                        ststr, np.std(toast_diff_truth)
+                    )
+                )
+                print(
+                    "madam-truth map {} has rms {}".format(
+                        ststr, np.std(madam_diff_truth)
+                    )
+                )
+
                 outfile = os.path.join(testdir, "madam_map_{}.png".format(ststr))
                 hp.mollview(madam_map[stokes], xsize=1600, nest=True)
                 plt.savefig(outfile)
@@ -684,6 +714,14 @@ def test_compare_madam_diagpre(self):
                 hp.mollview(diff_map, xsize=1600, nest=True)
                 plt.savefig(outfile)
                 plt.close()
+                outfile = os.path.join(testdir, "madam_diff_truth_{}.png".format(ststr))
+                hp.mollview(madam_diff_truth, xsize=1600, nest=True)
+                plt.savefig(outfile)
+                plt.close()
+                outfile = os.path.join(testdir, "toast_diff_truth_{}.png".format(ststr))
+                hp.mollview(toast_diff_truth, xsize=1600, nest=True)
+                plt.savefig(outfile)
+                plt.close()
 
                 if not np.allclose(
                     toast_map[stokes], madam_map[stokes], atol=0.1, rtol=0.05
@@ -730,7 +768,12 @@ def test_compare_madam_bandpre(self):
         weights.apply(data)
 
         # Create fake polarized sky pixel values locally
-        create_fake_sky(data, "pixel_dist", "fake_map")
+        create_fake_sky(
+            data,
+            "pixel_dist",
+            "fake_map",
+            hpix_out=os.path.join(testdir, "input_map.fits"),
+        )
 
         # Scan map into timestreams
         scanner = ops.ScanMap(
@@ -782,6 +825,7 @@ def test_compare_madam_bandpre(self):
             step_time=step_seconds * u.second,
             use_noise_prior=True,
             precond_width=10,
+            debug_plots=testdir,
         )
 
         tmatrix = ops.TemplateMatrix(templates=[tmpl])
@@ -843,7 +887,7 @@ def test_compare_madam_bandpre(self):
         pars["kfilter"] = "T"
         pars["precond_width_min"] = 10
         pars["precond_width_max"] = 10
-        pars["use_cgprecond"] = "T"
+        pars["use_cgprecond"] = "F"
         pars["use_fprecond"] = "F"
         pars["info"] = 3
         pars["path_output"] = testdir
@@ -893,13 +937,19 @@ def test_compare_madam_bandpre(self):
 
             # Compare maps
 
+            input_map = hp.read_map(
+                os.path.join(testdir, "input_map.fits"), field=None, nest=True
+            )
             toast_map = hp.read_map(toast_map_path, field=None, nest=True)
             madam_map = hp.read_map(madam_map_path, field=None, nest=True)
             # Set madam unhit pixels to zero
             for stokes, ststr in zip(range(3), ["I", "Q", "U"]):
                 mask = hp.mask_bad(madam_map[stokes])
                 madam_map[stokes][mask] = 0.0
+                input_map[stokes][mask] = 0.0
                 diff_map = toast_map[stokes] - madam_map[stokes]
+                madam_diff_truth = madam_map[stokes] - input_map[stokes]
+                toast_diff_truth = toast_map[stokes] - input_map[stokes]
                 # print("diff map {} has rms {}".format(ststr, np.std(diff_map)))
                 outfile = os.path.join(testdir, "madam_map_{}.png".format(ststr))
                 hp.mollview(madam_map[stokes], xsize=1600, nest=True)
@@ -913,6 +963,14 @@ def test_compare_madam_bandpre(self):
                 hp.mollview(diff_map, xsize=1600, nest=True)
                 plt.savefig(outfile)
                 plt.close()
+                outfile = os.path.join(testdir, "madam_diff_truth_{}.png".format(ststr))
+                hp.mollview(madam_diff_truth, xsize=1600, nest=True)
+                plt.savefig(outfile)
+                plt.close()
+                outfile = os.path.join(testdir, "toast_diff_truth_{}.png".format(ststr))
+                hp.mollview(toast_diff_truth, xsize=1600, nest=True)
+                plt.savefig(outfile)
+                plt.close()
 
                 if not np.allclose(
                     toast_map[stokes], madam_map[stokes], atol=0.1, rtol=0.05
diff --git a/src/toast/tests/ops_mapmaker_solve.py b/src/toast/tests/ops_mapmaker_solve.py
index 21c642c5c..11ffcd592 100644
--- a/src/toast/tests/ops_mapmaker_solve.py
+++ b/src/toast/tests/ops_mapmaker_solve.py
@@ -187,8 +187,8 @@ def test_lhs(self):
 
         tmatrix.amplitudes = "amplitudes"
         tmatrix.det_data = defaults.det_data
-        tmatrix.data = data
         tmatrix.transpose = False
+        tmatrix.initialize(data)
         tmatrix.apply(data)
 
         # Pointing operator
diff --git a/src/toast/tests/ops_mapmaker_utils.py b/src/toast/tests/ops_mapmaker_utils.py
index 3199d9548..56eca6ee3 100644
--- a/src/toast/tests/ops_mapmaker_utils.py
+++ b/src/toast/tests/ops_mapmaker_utils.py
@@ -65,7 +65,7 @@ def test_hits(self):
         # Manual check
         check_hits = PixelData(data["pixel_dist"], np.int64, n_value=1)
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 local_sm, local_pix = data["pixel_dist"].global_pixel_to_submap(
                     ob.detdata["pixels"][det]
                 )
@@ -159,8 +159,7 @@ def test_inv_cov(self):
         for ob in data.obs:
             noise = ob["noise_model"]
             noise_corr = ob["noise_model_corr"]
-
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 detweight = noise.detector_weight(det).to_value(invnpp_units)
                 detweight_corr = noise_corr.detector_weight(det).to_value(invnpp_units)
 
@@ -322,7 +321,7 @@ def test_zmap(self):
         for ob in data.obs:
             noise = ob["noise_model"]
             noise_corr = ob["noise_model_corr"]
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 wt = ob.detdata[defaults.weights][det]
                 local_sm, local_pix = data["pixel_dist"].global_pixel_to_submap(
                     ob.detdata[defaults.pixels][det]
diff --git a/src/toast/tests/ops_noise_estim.py b/src/toast/tests/ops_noise_estim.py
index 5b1810d8f..cb4954fdc 100644
--- a/src/toast/tests/ops_noise_estim.py
+++ b/src/toast/tests/ops_noise_estim.py
@@ -508,7 +508,7 @@ def test_noise_estim_model(self):
                 input_model = ob["noise_model"]
                 estim_model = ob["noise_estimate"]
                 fit_model = ob[noise_fitter.out_model]
-                for det in ob.local_detectors:
+                for det in ob.select_local_detectors(flagmask=estim.det_flag_mask):
                     fname = os.path.join(
                         self.outdir, f"estimate_model_{ob.name}_{det}.pdf"
                     )
@@ -533,7 +533,7 @@ def test_noise_estim_model(self):
             input_model = ob["noise_model"]
             estim_model = ob["noise_estimate"]
             fit_model = ob[noise_fitter.out_model]
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=estim.det_flag_mask):
                 np.testing.assert_almost_equal(
                     np.mean(input_model.psd(det)[-5:]).value,
                     np.mean(fit_model.psd(det)[-5:]).value,
diff --git a/src/toast/tests/ops_pointing_wcs.py b/src/toast/tests/ops_pointing_wcs.py
index 879b09c14..33e7e53bd 100644
--- a/src/toast/tests/ops_pointing_wcs.py
+++ b/src/toast/tests/ops_pointing_wcs.py
@@ -395,7 +395,7 @@ def create_source_data(
             stheta = (90.0 - ob.shared["source"].data[:, 1]) * np.pi / 180.0
             spos = qa.from_iso_angles(stheta, sphi, np.zeros_like(stheta))
             sdir = qa.rotate(spos, zaxis)
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 fwhm = 2 * ob.telescope.focalplane[det]["fwhm"].to_value(u.arcmin)
                 coeff = 1.0 / (fwhm * np.sqrt(2 * np.pi))
                 pre = -0.5 * (1.0 / fwhm) ** 2
diff --git a/src/toast/tests/ops_polyfilter.py b/src/toast/tests/ops_polyfilter.py
index a3691d767..ee4fcfabb 100644
--- a/src/toast/tests/ops_polyfilter.py
+++ b/src/toast/tests/ops_polyfilter.py
@@ -79,7 +79,8 @@ def test_polyfilter(self):
 
         for ob in data.obs:
             times = np.array(ob.shared[defaults.times])
-            for det in ob.local_detectors:
+
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = (
                     np.array(ob.shared[defaults.shared_flags])
                     & defaults.shared_mask_invalid
@@ -109,7 +110,7 @@ def test_polyfilter(self):
         polyfilter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = (
                     np.array(ob.shared[defaults.shared_flags])
                     & defaults.shared_mask_invalid
@@ -171,7 +172,7 @@ def test_polyfilter_trend(self):
 
         for ob in data.obs:
             times = np.array(ob.shared[defaults.times])
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = (
                     np.array(ob.shared[defaults.shared_flags])
                     & defaults.shared_mask_invalid
@@ -201,7 +202,7 @@ def test_polyfilter_trend(self):
         polyfilter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = (
                     np.array(ob.shared[defaults.shared_flags])
                     & defaults.shared_mask_invalid
@@ -297,7 +298,7 @@ def test_polyfilter2D(self):
         # Add 2D polynomial modes.
         coeff = np.arange(norder**2)
         for obs in data.obs:
-            for det in obs.local_detectors:
+            for det in obs.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 det_quat = obs.telescope.focalplane[det]["quat"]
                 x, y, z = qa.rotate(det_quat, ZAXIS)
                 theta, phi = np.arcsin([x, y])
@@ -332,7 +333,9 @@ def test_polyfilter2D(self):
             fig = plt.figure(figsize=[18, 12])
             ax = fig.add_subplot(1, 3, 1)
             ob = data.obs[0]
-            for idet, det in enumerate(ob.local_detectors):
+            for idet, det in enumerate(
+                ob.select_local_detectors(flagmask=defaults.det_mask_invalid)
+            ):
                 good = np.logical_and(
                     ob.detdata[defaults.det_flags][det, :] == 0,
                     ob.shared[defaults.shared_flags].data == 0,
@@ -366,7 +369,9 @@ def test_polyfilter2D(self):
 
         if data.comm.world_rank == 0:
             ax = fig.add_subplot(1, 3, 2)
-            for idet, det in enumerate(ob.local_detectors):
+            for idet, det in enumerate(
+                ob.select_local_detectors(flagmask=defaults.det_mask_invalid)
+            ):
                 good = np.logical_and(
                     ob.detdata["pyflags"][det, :] == 0,
                     ob.shared[defaults.shared_flags].data == 0,
@@ -387,7 +392,9 @@ def test_polyfilter2D(self):
 
         if data.comm.world_rank == 0:
             ax = fig.add_subplot(1, 3, 3)
-            for idet, det in enumerate(ob.local_detectors):
+            for idet, det in enumerate(
+                ob.select_local_detectors(flagmask=defaults.det_mask_invalid)
+            ):
                 good = np.logical_and(
                     ob.detdata[defaults.det_flags][det, :] == 0,
                     ob.shared[defaults.shared_flags].data == 0,
@@ -402,7 +409,7 @@ def test_polyfilter2D(self):
 
         # Check for consistency
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 check = np.allclose(
                     ob.detdata[defaults.det_data][det], ob.detdata["pyfilter"][det]
                 )
@@ -474,7 +481,7 @@ def test_common_mode_filter(self):
         for ob in data.obs:
             rms[ob.name] = dict()
             times = ob.shared[defaults.times]
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = np.array(ob.shared[defaults.shared_flags])
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
@@ -495,7 +502,7 @@ def test_common_mode_filter(self):
         common_filter.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 flags = np.array(ob.shared[defaults.shared_flags])
                 flags |= ob.detdata[defaults.det_flags][det]
                 good = flags == 0
diff --git a/src/toast/tests/ops_scan_map.py b/src/toast/tests/ops_scan_map.py
index 6c3f3b143..3bdd17a25 100644
--- a/src/toast/tests/ops_scan_map.py
+++ b/src/toast/tests/ops_scan_map.py
@@ -60,7 +60,9 @@ def test_scan(self):
 
             # Manual check of the projection of map values to timestream
             for ob in data.obs:
-                for det in ob.local_detectors:
+                for det in ob.select_local_detectors(
+                    flagmask=defaults.det_mask_invalid
+                ):
                     wt = ob.detdata[weights.weights][det]
                     local_sm, local_pix = data["pixel_dist"].global_pixel_to_submap(
                         ob.detdata[pixels.pixels][det]
@@ -179,7 +181,7 @@ def test_mask(self):
 
         mask_data = data["fake_mask"]
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 local_sm, local_pix = data["pixel_dist"].global_pixel_to_submap(
                     ob.detdata[pixels.pixels][det]
                 )
diff --git a/src/toast/tests/ops_sim_tod_atm.py b/src/toast/tests/ops_sim_tod_atm.py
index 5d5c2249c..b874c2359 100644
--- a/src/toast/tests/ops_sim_tod_atm.py
+++ b/src/toast/tests/ops_sim_tod_atm.py
@@ -255,7 +255,9 @@ def test_sim_interp(self):
 
         if rank == 0:
             for obs in data.obs:
-                for det in obs.local_detectors:
+                for det in obs.select_local_detectors(
+                    flagmask=defaults.det_mask_invalid
+                ):
                     sig0 = obs.detdata["full_signal"][det]
                     sig1 = obs.detdata["interpolated_signal"][det]
                     assert np.std(sig1 - sig0) / np.std(sig0) < 1e-2
@@ -481,7 +483,7 @@ def test_loading(self):
         # a non-zero signal
 
         for obs in data.obs:
-            for det in obs.local_detectors:
+            for det in obs.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 sig = obs.detdata[defaults.det_data][det]
                 assert np.std(sig) != 0
 
@@ -547,7 +549,7 @@ def test_bandpass(self):
 
         # Check that the atmospheric fluctuations are stronger at higher frequency
         for obs in data.obs:
-            for det in obs.local_detectors:
+            for det in obs.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 new_rms = np.std(obs.detdata[defaults.det_data][det])
                 assert new_rms > 1.1 * old_rms[obs.name][det]
 
diff --git a/src/toast/tests/ops_statistics.py b/src/toast/tests/ops_statistics.py
index c22e0d7ff..06f591a27 100644
--- a/src/toast/tests/ops_statistics.py
+++ b/src/toast/tests/ops_statistics.py
@@ -60,7 +60,7 @@ def test_statistics(self):
                 skew = f["skewness"][:].copy()
                 kurt = f["kurtosis"][:].copy()
             # Test the statistics
-            for det in obs.local_detectors:
+            for det in obs.select_local_detectors(flagmask=statistics.det_flag_mask):
                 idet = detectors.index(det)
                 sig = obs.detdata[defaults.det_data][det]
                 # Test variance
diff --git a/src/toast/tests/template_gain.py b/src/toast/tests/template_gain.py
index 0457d015b..6bdfcef19 100644
--- a/src/toast/tests/template_gain.py
+++ b/src/toast/tests/template_gain.py
@@ -9,6 +9,7 @@
 from astropy import units as u
 
 from .. import ops
+from ..observation import default_values as defaults
 from ..templates import GainTemplate
 from ..utils import rate_from_times
 from ._helpers import close_data, create_outdir, create_satellite_data
@@ -79,7 +80,7 @@ def test_gainfit(self):
         calibration.apply(data)
 
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 np.testing.assert_allclose(
                     ob.detdata["calibrated"][det], np.ones(ob.n_local_samples)
                 )
diff --git a/src/toast/tests/template_offset.py b/src/toast/tests/template_offset.py
index acc16f5a8..6b5049829 100644
--- a/src/toast/tests/template_offset.py
+++ b/src/toast/tests/template_offset.py
@@ -60,7 +60,7 @@ def test_projection(self):
 
         # Verify
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 np.testing.assert_equal(ob.detdata[defaults.det_data][det], 1.0)
 
         # Accumulate amplitudes
@@ -84,7 +84,7 @@ def test_projection(self):
             slices.append(slice((n_step - 1) * step_samples, ob.n_local_samples, 1))
             sizes.append(ob.n_local_samples - (n_step - 1) * step_samples)
 
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 for slc, sz in zip(slices, sizes):
                     np.testing.assert_equal(
                         np.sum(ob.detdata[defaults.det_data][det, slc]), 1.0 * sz
@@ -168,7 +168,7 @@ def test_accel(self):
             slices.append(slice((n_step - 1) * step_samples, ob.n_local_samples, 1))
             sizes.append(ob.n_local_samples - (n_step - 1) * step_samples)
 
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 for slc, sz in zip(slices, sizes):
                     np.testing.assert_equal(
                         np.sum(ob.detdata[defaults.det_data][det, slc]), 1.0 * sz
@@ -263,7 +263,9 @@ def test_accel_pipeline(self):
         # input amplitudes (with starting value 1.0) without clearing those,
         # so we add 1.0 to the expected values.
         offset = 0
-        all_dets = list(data.obs[0].local_detectors)
+        all_dets = list(
+            data.obs[0].select_local_detectors(flagmask=defaults.det_data_mask)
+        )
         expected = list()
         for det in all_dets:
             for iob, ob in enumerate(data.obs):
diff --git a/src/toast/tests/template_periodic.py b/src/toast/tests/template_periodic.py
index 7234a674a..ee3ac6c31 100644
--- a/src/toast/tests/template_periodic.py
+++ b/src/toast/tests/template_periodic.py
@@ -37,7 +37,7 @@ def setUp(self):
         self.outdir = create_outdir(self.comm, fixture_name)
         self.nside = 64
         np.random.seed(123456)
-        if not (
+        if (
             ("CONDA_BUILD" in os.environ)
             or ("CIBUILDWHEEL" in os.environ)
             or ("CI" in os.environ)
@@ -560,14 +560,15 @@ def test_satellite_hwp(self):
 
         # Compare the cleaned timestreams to the original
         for ob in data.obs:
-            for det in ob.local_detectors:
+            for det in ob.select_local_detectors(flagmask=defaults.det_mask_invalid):
                 dof = np.sqrt(ob.n_local_samples)
                 in_rms = np.std(ob.detdata["input"][det])
                 thresh = hwpss_scale * (in_rms / dof)
                 out_rms = np.std(ob.detdata[defaults.det_data][det])
-                print(f"Using threshold {thresh}, input rms = {in_rms}")
                 if out_rms > thresh:
-                    print(f"{ob.name}:{det} output rms = {out_rms} exceeds threshold")
+                    msg = f"{ob.name}:{det} output rms = {out_rms} exceeds threshold"
+                    msg += f" ({thresh}) for input rms = {in_rms}"
+                    print(msg)
                     raise RuntimeError("Failed satellite hwpss template regression")
 
         close_data(data)