Added script that solves satv1 pointing model parameters.

sotodlib/coords/fp_containers.py

@@ -141,6 +141,7 @@ class FocalPlane:
     id_strs: NDArray[np.str_]  # (ndet,)
     avg_fp: NDArray[np.floating]  # (ndim, ndet)
     weights: NDArray[np.floating]  # (ndet,)
+    det_boresight: NDArray[np.floating] # (ndim, ndet)
     transformed: NDArray[np.floating]  # (ndet, ndim)
     center: NDArray[np.floating]  # (1, ndim)
     center_transformed: NDArray[np.floating]  # (1, ndim)
@@ -219,6 +220,7 @@ def empty(cls, template, stream_id, wafer_slot, n_aman, config=""):
         tot_weight = np.zeros((len(template.det_ids), 2))
         avg_fp = np.full_like(template.fp, np.nan)
         weight = np.zeros((len(template.det_ids), 2))
+        det_boresight = np.zeros((len(template.det_ids), 3)) + np.nan # az, el, roll
         transformed = template.fp.copy()
         center = template.center.copy()
         center_transformed = template.center.copy()
@@ -232,6 +234,7 @@ def empty(cls, template, stream_id, wafer_slot, n_aman, config=""):
             template.id_strs,
             avg_fp,
             weight,
+            det_boresight,
             transformed,
             center,
             center_transformed,
@@ -260,8 +263,17 @@ def map_by_det_id(self, aman):
         xi = aman.pointing.xi[msk][srt][mapping]
         eta = aman.pointing.eta[msk][srt][mapping]
         r2 = np.nan + np.zeros_like(eta)
+        az = np.deg2rad(aman.obs_info.az_center) * np.ones_like(eta)
+        el = np.deg2rad(aman.obs_info.el_center) * np.ones_like(eta)
+        roll = np.deg2rad(aman.obs_info.roll_center) * np.ones_like(eta)
         if "R2" in aman.pointing:
             r2 = aman.pointing.R2[msk][srt][mapping]
+        if "az" in aman.pointing:
+            az = aman.pointing.az[msk][srt][mapping]
+        if "el" in aman.pointing:
+            el = aman.pointing.el[msk][srt][mapping]
+        if "roll" in aman.pointing:
+            roll = aman.pointing.roll[msk][srt][mapping]
         if "polarization" in aman:
             # name of field just a placeholder for now
             gamma = aman.polarization.polang[msk][srt][mapping]
@@ -270,14 +282,16 @@ def map_by_det_id(self, aman):
         else:
             gamma = np.full(len(xi), np.nan)
         fp = np.column_stack((xi, eta, gamma))
-        return fp, r2, template_msk
+        det_boresight = np.column_stack((az, el, roll))
+        return fp, r2, det_boresight, template_msk
 
-    def add_fp(self, i, fp, weights, template_msk):
-        if self.full_fp is None or self.tot_weight is None:
+    def add_fp(self, i, fp, weights, det_boresight, template_msk):
+        if self.full_fp is None or self.tot_weight is None or self.det_boresight is None:
             raise ValueError("full_fp or tot_weight not initialized")
         self.full_fp[template_msk, :, i] = fp * weights[:, 0][..., None]
         weights = np.nan_to_num(weights)
         self.tot_weight[template_msk] += weights
+        self.det_boresight[template_msk, :] = det_boresight
 
     def save(self, f, db_info, group):
         logger.info("Saving %s", self.stream_id)
@@ -323,6 +337,9 @@ def save(self, f, db_info, group):
             ("gamma_m", np.float32),
             ("weights", np.float32),
             ("r2", np.float32),
+            ("az", np.float32),
+            ("el", np.float32),
+            ("roll", np.float32),
             ("n_point", np.int8),
             ("n_gamma", np.int8),
         ]
@@ -333,6 +350,7 @@ def save(self, f, db_info, group):
                 *(self.transformed.T),
                 *(self.avg_fp.T),
                 *(self.weights.T),
+                *(self.det_boresight.T),
                 self.n_point,
                 self.n_gamma,
             ),
@@ -378,6 +396,14 @@ def load(cls, group, include_cm=None):
                 np.array(fp_full["gamma_m"]),
             )
         )
+
+        det_boresight = np.column_stack(
+            (
+                np.array(fp_full["az"]),
+                np.array(fp_full["el"]),
+                np.array(fp_full["roll"]),
+            )
+        )
         # For backwards compatibility
         weights = np.array(fp_full["weights"])
         if "r2" in fp_full.keys:
@@ -424,6 +450,7 @@ def load(cls, group, include_cm=None):
             np.array(id_strs),
             avg_fp,
             np.array(weights),
+            det_boresight,
             transformed,
             center,
             center_transformed,

sotodlib/core/flagman.py

@@ -172,7 +172,7 @@ def reduce(self, flags=None, method='union', wrap=False, new_flag=None,
             flags: List of flags to collapse together. Uses their names.
                    If flags is None then all flags are reduced
             method: How to collapse the data. Accepts 'union','intersect',
-                        or function.
+                    'except', or function.
             wrap: if True, add reduced flag to self
             new_flag: name of new flag, required if wrap is True
             remove_reduced: if True, remove all reduced flags from self
@@ -198,6 +198,8 @@ def reduce(self, flags=None, method='union', wrap=False, new_flag=None,
             op = lambda x, y: x+y
         elif method == 'intersect':
             op = lambda x, y: x*y
+        elif method == 'except':
+            op = lambda x, y: x*~y
         else:
             op = method
         out = reduce(op, to_reduce)

sotodlib/preprocess/processes.py

@@ -1,3 +1,4 @@
+import pdb
 import numpy as np
 from operator import attrgetter
 import copy
@@ -1832,6 +1833,7 @@ class SourceFlags(_Preprocess):
     .. autofunction:: sotodlib.tod_ops.flags.get_source_flags
     """
     name = "source_flags"
+
     def __init__(self, step_cfgs):
         self.source_flags_name = step_cfgs.get('source_flags_name', 'source_flags')
         self.save_name = self.source_flags_name
@@ -1881,7 +1883,7 @@ def calc_and_save(self, aman, proc_aman):
                     source_aman.wrap(source + '_inv',
                                     RangesMatrix.ones([aman.dets.count, aman.samps.count]),
                                     [(0, 'dets'), (1, 'samps')])
-
+                  
         self.save(proc_aman, source_aman)
 
         return aman, proc_aman
@@ -2102,8 +2104,7 @@ def process(self, aman, proc_aman, sim=False, data_aman=None):
                                     aman.samps.offset + aman.samps.count - trim))
             proc_aman.restrict('samps', (proc_aman.samps.offset + trim,
                                          proc_aman.samps.offset + proc_aman.samps.count - trim))
-        return aman, proc_aman
-
+        return aman, proc_aman  
 
 class DetcalNanCuts(_Preprocess):
     """
@@ -2696,9 +2697,9 @@ def process(self, aman, proc_aman, sim=False, data_aman=None):
         aman['flags'].wrap(self.process_cfgs['total_flags_label'], total_flags)
 
         return aman, proc_aman
-
+            
 class CombineFlags(_Preprocess):
-    """Do the conbine of relevant flags for mapping
+    """Do the combination of relevant flags for mapping
 
 
     Saves results for aman under the "flags.[total_flags_label]" field.
@@ -2709,8 +2710,11 @@ class CombineFlags(_Preprocess):
           process:
             flag_labels: ['glitches.glitch_flags', 'source_flags.jupiter_inv']
             total_flags_label: 'glitch_flags'
-            method: 'union' # You can select a method from ['union', '+', 'intersect', '*'].
-            #method: ['+', '*'] # Or you can pass individual method for each flags as a list. Lentgh must match the length of flag_labels.
+            method: 'union' # You can select a method from ['union', '+', 'intersect', '*', 'except', '-'].
+            #method: ['+', '*'] # Or you can pass individual method for each flags as a list. 
+               # Length of list must match the length of flag_labels.
+               # If a list, the first method must be '+', as if adding the first flag set to an empty flag set. 
+               # Operations are performed strictly from Left to Right, '*' are not performed first.
 
     """
     name = "combine_flags"
@@ -2727,13 +2731,12 @@ def process(self, aman, proc_aman, sim=False, data_aman=None):
         if isinstance(self.process_cfgs['method'], list):
             if len(self.process_cfgs['flag_labels']) != len(self.process_cfgs['method']):
                 raise ValueError("The length of method does not match to the length of flag_labels")
-            elif any(method not in ['+', 'union', '*', 'intersect'] for method in self.process_cfgs['method']):
-                raise ValueError("The method provided does not match one of '+', '*', 'union', or 'intersect'")
-        elif self.process_cfgs['method'] in ['+', 'union', '*', 'intersect']:
+            elif any(method not in ['+', 'union', '*', 'intersect', '-', 'except'] for method in self.process_cfgs['method']):
+                raise ValueError("One or more methods in list are not valid")
+        elif self.process_cfgs['method'] in ['+', 'union', '*', 'intersect', '-', 'except']:
             self.process_cfgs['method'] =  ['+'] + (len(self.process_cfgs['flag_labels']) - 1)*[self.process_cfgs['method']]
         else:
-            raise ValueError("The method matches neither list nor the one of the ['+', 'union', '*', 'intersect']")
-
+            raise ValueError("The method matches neither list nor the one of the valid operations")
         total_flags = RangesMatrix.zeros([proc_aman.dets.count, proc_aman.samps.count]) # get an empty flags with shape (Ndets,Nsamps)
         for i, (method, label) in enumerate(zip(self.process_cfgs['method'], self.process_cfgs['flag_labels'])):
             _label = attrgetter(label)
@@ -2746,6 +2749,8 @@ def process(self, aman, proc_aman, sim=False, data_aman=None):
                     total_flags += _label(proc_aman) # The + operator is the union operator in this case
                 elif method in ['*', 'intersect']:
                     total_flags *= _label(proc_aman) # The * operator is the intersect operator in this case
+                elif method in ['-', 'except']:
+                    total_flags *= ~ _label(proc_aman) # The - operator is the except operator in this case
 
         if 'flags' not in aman._fields:
             from sotodlib.core import FlagManager

sotodlib/site_pipeline/cli.py

@@ -59,6 +59,7 @@ def main(obs_id=None, config_file=None, logger=None):
     update_obsdb_ancil,
     make_coadd_atomic_map,
     make_cosamp_hk,
+    solve_pointing_model,
 )
 
 # Dictionary matching element name to a submodule (which must have
@@ -81,7 +82,8 @@ def main(obs_id=None, config_file=None, logger=None):
     'update-obsdb': update_obsdb,
     'update-obsdb-ancil': update_obsdb_ancil,
     'make-cosamp-hk': make_cosamp_hk,
-    'make-coadd-atomic-map': make_coadd_atomic_map
+    'make-coadd-atomic-map': make_coadd_atomic_map,
+    'solve-pointing-model': solve_pointing_model,
 }
 
 CLI_NAME = 'so-site-pipeline'

sotodlib/site_pipeline/finalize_focal_plane.py

@@ -6,6 +6,10 @@
 from importlib import import_module
 from typing import List, Optional
 
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+
 import git
 import h5py
 import megham.transform as mt
@@ -403,7 +407,7 @@ def _mk_pointing_config(telescope_flavor, tube_slot, wafer_slot, config):
 def _restrict_inliers(aman, focal_plane):
     # TODO: Use gamma as well
     # Map to template
-    fp, _, template_msk = focal_plane.map_by_det_id(aman)
+    fp, _, _, template_msk = focal_plane.map_by_det_id(aman)
     fp = fp[:, :2]
     inliers = np.ones(len(fp), dtype=bool)
 
@@ -779,7 +783,7 @@ def main():
                 _restrict_inliers(aman, focal_plane)
 
                 # Mapping to template
-                fp, r2, template_msk = focal_plane.map_by_det_id(aman)
+                fp, r2, det_boresight, template_msk = focal_plane.map_by_det_id(aman)
                 focal_plane.template.add_wafer_info(aman, template_msk)
 
                 # Try an initial alignment and get weights
@@ -816,7 +820,7 @@ def main():
 
                 # Store weighted values
                 weights = np.column_stack((weights, r2))
-                focal_plane.add_fp(i, fp, weights, template_msk)
+                focal_plane.add_fp(i, fp, weights, det_boresight, template_msk)
 
                 n_obs += 1