Normalize beams to unit intensity

python/libconviqt_wrapper/__init__.py

@@ -71,6 +71,9 @@ def encode_comm(comm):
 _conviqt.conviqt_beam_normalize.restype = ct.c_double
 _conviqt.conviqt_beam_normalize.argtypes = [ct.c_void_p]
 
+_conviqt.conviqt_beam_normalized.restype = ct.c_int
+_conviqt.conviqt_beam_normalized.argtypes = [ct.c_void_p]
+
 
 class Beam(object):
     """
@@ -105,6 +108,13 @@ def normalize(self):
             raise RuntimeError("Failed to normalize beam. scale = {}".format(scale))
         return scale
 
+    def normalized(self):
+        result = _conviqt.conviqt_beam_normalized(self._beam)
+        if result < 0:
+            raise RuntimeError("Failed to query beam normalization status")
+        else:
+            return result != 0
+
 
 # Sky functions
 

src/cconviqt.cpp

@@ -63,6 +63,19 @@ extern "C" {
         }
     }
 
+    int conviqt_beam_normalized(void *ptr) {
+        try {
+            conviqt::beam *ref = reinterpret_cast< conviqt::beam * >(ptr);
+            if (ref->normalized()) {
+                return 1;
+            } else {
+                return 0;
+            }
+        } catch (...) {
+            return -1;
+        }
+    }
+
     void *conviqt_sky_new() { return new(std::nothrow) conviqt::sky; }
 
     int conviqt_sky_del(void *ptr) {

src/cconviqt.h

@@ -23,6 +23,7 @@ extern "C" {
 
   double conviqt_beam_normalize(void *ptr);
 
+  int conviqt_beam_normalized(void *ptr);
   // Sky
 
   void *conviqt_sky_new();

src/conviqt.hpp

@@ -37,7 +37,7 @@ namespace conviqt {
 
 class beam {
 public :
-    beam() {}
+    beam() { normalized_ = false; }
     beam(long beamlmax, long beammmax, bool pol,
          std::string infile_beam, MPI_Comm comm);
     int read(long beamlmax, long beammmax, bool pol,
@@ -49,8 +49,10 @@ public :
     int get_lmax(void) { return lmax; }
     int get_mmax(void) { return mmax; }
     double normalize(void);
+    bool normalized() { return normalized_; }
 private :
     Alm< xcomplex<float> > blmT_, blmG_, blmC_;
+    bool normalized_;
     long lmax, mmax;
     bool pol;
     int verbosity = 0;

src/conviqt_math.cpp

@@ -56,20 +56,25 @@ namespace conviqt {
 beam::beam(long beamlmax, long beammmax, bool pol,
            std::string infile_beam,
            MPI_Comm comm=MPI_COMM_WORLD) {
+    normalized_ = false;
     read(beamlmax, beammmax, pol, infile_beam, comm);
 }
 
+
 double beam::normalize() {
     /*
-      Scale the beam expansion to integrate to 0.5 over the sphere
+      Scale the intensity beam expansion to integrate to unity over the sphere
     */
+    if (normalized_) {
+        return 1.0;
+    }
     double scale = 0;
     if (blmT_.Lmax() >= 0) {
         double b00 = blmT_(0, 0).re;
         double current_norm = 2 * b00 / sqrt(1 / pi);
-        scale = 0.5 / current_norm;
+        scale = 1.0 / current_norm;
         if (verbosity > 1) {
-            std::cerr << "Normalizing beam from " << current_norm << " to 0.5 with "
+            std::cerr << "Normalizing beam from " << current_norm << " to 1.0 with "
                       << scale << std::endl;
         }
         blmT_.Scale(scale);
@@ -78,6 +83,7 @@ double beam::normalize() {
             blmC_.Scale(scale);
         }
     }
+    normalized_ = true;
     return scale;
 }
 
@@ -1644,8 +1650,12 @@ int convolver::convolve(pointing &pntarr, bool calibrate) {
         ++n_sort;
     }
 
-    if (totsize != 0 && calibrate) {
+    if (totsize != 0 && calibrate && !b->normalized()) {
         double calibration = 2. / (1. + d->get_epsilon());
+        if (mpiMgr.master() and verbosity > 0) {
+            std::cerr << "Calibrating TOD with " << calibration
+                      << " for unit intensity response" << std::endl;
+        }
         for (long ii = 0; ii < totsize; ++ii) {
             // Insert convolved TOD into the output array
             pntarr[5 * ii + 3] *= calibration;