Tweaks arising from the BICEP/Keck calibration campaign (#367)

src/libtoast/src/toast_atm_sim.cpp

@@ -591,11 +591,11 @@ cholmod_sparse * toast::atm_sim_build_sparse_covariance(
 
                 // If the covariance exceeds the threshold, add it to the
                 // sparse matrix
-                if (val * val > 1e-6 * diagonal[icol] * diagonal[irow]) {
-                    myrows.push_back(irow);
-                    mycols.push_back(icol);
-                    myvals.push_back(val);
-                }
+                //if (val * val > 1e-6 * diagonal[icol] * diagonal[irow]) {
+                myrows.push_back(irow);
+                mycols.push_back(icol);
+                myvals.push_back(val);
+                //}
             }
         }
         # pragma omp critical
@@ -698,7 +698,7 @@ cholmod_sparse * toast::atm_sim_sqrt_sparse_covariance(
 
     cholmod_factor * factorization;
 
-    const int ntry = 4;
+    const int ntry = 2;
 
     for (int itry = 0; itry < ntry; ++itry) {
         factorization = cholmod_analyze(cov, chol.chcommon);
@@ -729,12 +729,11 @@ cholmod_sparse * toast::atm_sim_sqrt_sparse_covariance(
                     cholmod_print_sparse(cov, "Covariance matrix", chol.chcommon);
 
                     // DEBUG begin
-                    if (itry > 2) {
+                    if (itry == ntry - 2) {
                         FILE * covfile = fopen("failed_covmat.mtx", "w");
                         cholmod_write_sparse(covfile, cov, NULL, NULL,
                                              chol.chcommon);
                         fclose(covfile);
-                        exit(-1);
                     }
 
                     // DEBUG end

src/libtoast/src/toast_fod_psd.cpp

@@ -64,16 +64,19 @@ void toast::fod_crosssums(int64_t n, double const * x, double const * y,
     #pragma \
     omp parallel for default(none) shared(n, gd, lagmax, xgood, ygood, sums, hits) schedule(static, 100)
     for (int64_t lag = 0; lag < lagmax; ++lag) {
-        int64_t j = lag;
+        int64_t i, j;
         double lagsum = 0.0;
         int64_t hitsum = 0;
-        for (int64_t i = 0; i < (n - lag); ++i) {
+        for (i = 0, j = lag; i < (n - lag); ++i, ++j) {
             lagsum += xgood[i] * ygood[j];
             hitsum += gd[i] * gd[j];
-            j++;
+        }
+        // Use symmetry to double the statistics
+        for (i = 0, j = lag; i < (n - lag); ++i, ++j) {
+            lagsum += xgood[j] * ygood[i];
         }
         sums[lag] = lagsum;
-        hits[lag] = hitsum;
+        hits[lag] = 2 * hitsum;
     }
 
     return;

src/toast/pipeline_tools/atm.py

@@ -415,23 +415,33 @@ def scale_atmosphere_by_frequency(
             absorption = np.hstack(todcomm.allgather(my_absorption))
             loading = np.hstack(todcomm.allgather(my_loading))
         for det in tod.local_dets:
-            try:
-                # Use detector bandpass from the focalplane
-                center = focalplane[det]["bandcenter_ghz"]
-                width = focalplane[det]["bandwidth_ghz"]
-            except Exception:
-                # Use default values for the entire focalplane
-                if freq is None:
-                    raise RuntimeError(
-                        "You must supply the nominal frequency if bandpasses "
-                        "are not available"
-                    )
-                center = freq
-                width = 0.2 * freq
-            nstep = 101
-            # Interpolate the absorption coefficient to do a top hat
-            # integral across the bandpass
-            det_freqs = np.linspace(center - width / 2, center + width / 2, nstep)
+            if "bandpass_transmission" in focalplane[det]:
+                # We have full bandpasses for the detector
+                bandpass_freqs = focalplane[det]["bandpass_freq_ghz"]
+                bandpass = focalplane[det]["bandpass_transmission"]
+            else:
+                if "bandcenter_ghz" in focalplane[det]:
+                    # Use detector bandpass from the focalplane
+                    center = focalplane[det]["bandcenter_ghz"]
+                    width = focalplane[det]["bandwidth_ghz"]
+                else:
+                    # Use default values for the entire focalplane
+                    if freq is None:
+                        raise RuntimeError(
+                            "You must supply the nominal frequency if bandpasses "
+                            "are not available"
+                        )
+                    center = freq
+                    width = 0.2 * freq
+                bandpass_freqs = np.array([center - width / 2, center + width / 2])
+                bandpass = np.ones(2)
+            # Normalize and interpolate the bandpass
+            nstep = 1001
+            fmin, fmax = bandpass_freqs[0], bandpass_freqs[-1]
+            det_freqs = np.linspace(fmin, fmax, nstep)
+            det_bandpass = np.interp(det_freqs, bandpass_freqs, bandpass)
+            det_bandpass /= np.sum(det_bandpass)
+            # Interpolate absorption and loading to bandpass frequencies
             absorption_det = np.interp(det_freqs, freqs, absorption)
             loading_det = np.interp(det_freqs, freqs, loading)
             # From brightness to thermodynamic units
@@ -441,8 +451,8 @@ def scale_atmosphere_by_frequency(
             rj2cmb *= 0.5763279042527544
             absorption_det *= rj2cmb
             # Average across the bandpass
-            absorption_det = np.mean(absorption_det)
-            loading_det = np.mean(loading_det)
+            absorption_det = np.sum(absorption_det * det_bandpass)
+            loading_det = np.sum(loading_det * det_bandpass)
             cachename = "{}_{}".format(cache_name, det)
             ref = tod.cache.reference(cachename)
             ref *= absorption_det

src/toast/todmap/atm.py

@@ -625,6 +625,9 @@ def draw(self):
         # Wind is parallel to surface. Rotate to a frame where the scan
         # is across the X-axis.
 
+        #self._w *= 0.5  # DEBUG
+        #self._wdir += np.pi  # DEBUG
+
         eastward_wind = self._w * np.cos(self._wdir)
         northward_wind = self._w * np.sin(self._wdir)
 

src/toast/todmap/groundfilter.py

@@ -163,6 +163,14 @@ def fit_templates(self, tod, det, templates, ref, good):
             rank = comm.rank
         detranks, sampranks = tod.grid_size
 
+        if good is not None:
+            ngood = np.sum(good)
+        else:
+            ngood = 0
+        ngood_tot = comm.allreduce(ngood)
+        if ngood_tot == 0:
+            return None
+
         ntemplate = len(templates)
         invcov = np.zeros([ntemplate, ntemplate])
         proj = np.zeros(ntemplate)
@@ -243,7 +251,7 @@ def exec(self, data):
                     good = None
 
                 coeff = self.fit_templates(tod, det, templates, ref, good)
-                if ref is not None:
+                if coeff is not None:
                     self.subtract_templates(ref, good, coeff, cheby_trend, cheby_filter)
 
                 del ref