fix polarized flux conversion, revert polarization changes

toltec-astro · May 17, 2022 · 2bc3d42 · 2bc3d42
1 parent 7861b4b
commit 2bc3d42
Show file tree

Hide file tree

Showing 5 changed files with 37 additions and 18 deletions.
diff --git a/include/citlali/core/engine/engine.h b/include/citlali/core/engine/engine.h
@@ -394,17 +394,6 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
         get_config(crpix2,std::tuple{"mapmaking","crpix2"});
         get_config(cunit,std::tuple{"mapmaking","cunit"},{"MJy/Sr","mJy/beam"});
 
-        cflux.resize(toltec_io.barea_keys.size());
-
-        for (Eigen::Index i=0; i<toltec_io.barea_keys.size(); i++) {
-            if (cunit == "mJy/beam") {
-                cflux(i) = toltec_io.barea_keys[i]*MJY_SR_TO_mJY_ASEC;
-            }
-            else if (cunit == "MJy/Sr") {
-                cflux(i) = 1;
-            }
-        }
-
         // get beammap config options
         get_config(cutoff,std::tuple{"beammap","iter_tolerance"});
         get_config(max_iterations,std::tuple{"beammap","iter_max"});

diff --git a/include/citlali/core/engine/lali.h b/include/citlali/core/engine/lali.h
@@ -201,15 +201,22 @@ auto Lali::run() {
 
         for (auto const& stokes_params: polarization.stokes_params) {
             TCData<TCDataKind::RTC,Eigen::MatrixXd> in2;
+
+            SPDLOG_INFO("in scans {}",in.scans.data);
+
             auto [map_index_vector, det_index_vector] =  polarization.create_rtc(in, in2, stokes_params.first, this);
 
+            SPDLOG_INFO("in2 scans {}",in2.scans.data);
+
             /*Stage 1: RTCProc*/
             RTCProc rtcproc;
             {
                 tula::logging::scoped_timeit timer("rtcproc.run()");
                 rtcproc.run(in2, out, map_index_vector, det_index_vector, this);
             }
 
+            SPDLOG_INFO("out scans {}",out.scans.data);
+
             if (run_tod_output) {
                 if (ts_chunk_type == "rtc") {
 

diff --git a/include/citlali/core/timestream/rtc/calibrate.h b/include/citlali/core/timestream/rtc/calibrate.h
@@ -7,13 +7,13 @@ namespace timestream {
 template <typename DerivedA, typename DerivedB, typename DerivedC, typename DerivedD>
 void calibrate(Eigen::DenseBase<DerivedA> &in, Eigen::DenseBase<DerivedB> &flxscale,
                Eigen::DenseBase<DerivedC> &map_index_vector, Eigen::DenseBase<DerivedC> &det_index_vector,
-               Eigen::DenseBase<DerivedD> &cscale) {
+               Eigen::DenseBase<DerivedD> &cflux) {
 
     for (Eigen::Index det=0; det<in.cols(); det++) {
         Eigen::Index di = det_index_vector(det);
         Eigen::Index mi = map_index_vector(det);
 
-        in.col(det) = in.col(det)*flxscale(di)*cscale(mi);
+        in.col(det) = in.col(det)*flxscale(di)*cflux(mi);
     }
 }
 

diff --git a/include/citlali/core/timestream/rtc/polarization.h b/include/citlali/core/timestream/rtc/polarization.h
@@ -100,7 +100,7 @@ class Polarization {
         else {
             Eigen::MatrixXd data;
             Eigen::Index ori;
-            auto pa2 = in.tel_meta_data.data["ParAng"].array();// - pi;
+            auto pa2 = in.tel_meta_data.data["ParAng"].array() - pi;
 
             if (sp == "Q") {
                 SPDLOG_INFO("creating Q timestream");
@@ -152,15 +152,15 @@ class Polarization {
                 Eigen::Index di = det_index_vector(i);
                 // rotate by PA
                 if (sp == "Q") {
-                    auto qs0 = cos(-2*pa2.array())*data.col(i).array();
-                    auto us0 = sin(-2*pa2.array())*data.col(i).array();
+                    auto qs0 = cos(2*pa2.array())*data.col(i).array();
+                    auto us0 = sin(2*pa2.array())*data.col(i).array();
 
                     derotate_detector(qs0, us0, i, out, nsamples, ndet, di, engine->calib_data["y_t"](di), engine->run_hwp);
                 }
 
                 else if (sp == "U") {
-                    auto qs0 = -sin(-2*pa2.array())*data.col(i).array();
-                    auto us0 = cos(-2*pa2.array())*data.col(i).array();
+                    auto qs0 = -sin(2*pa2.array())*data.col(i).array();
+                    auto us0 = cos(2*pa2.array())*data.col(i).array();
 
                     derotate_detector(qs0, us0, i, out, nsamples, ndet, di, engine->calib_data["y_t"](di), engine->run_hwp);
                 }

diff --git a/src/citlali/cli/main.cpp b/src/citlali/cli/main.cpp
@@ -1181,6 +1181,7 @@ int run(const rc_t &rc) {
                     // load telescope file
                     todproc.engine().get_telescope(rawobs.teldata().filepath());
 
+                    // coadded exposure time
                     todproc.engine().c_t_exp += todproc.engine().tel_meta_params["t_exp"];
 
                     // calculate physical pointing vectors
@@ -1363,6 +1364,28 @@ int run(const rc_t &rc) {
                     todproc.engine().array_indices = array_indices.at(i);
                     todproc.engine().det_indices = det_indices.at(i);
 
+                    todproc.engine().cflux.resize(map_counts[i]);
+
+                    Eigen::Index k = 0;
+                    Eigen::Index l = 0;
+                    for (Eigen::Index j=0; j<todproc.engine().cflux.size(); j++) {
+                        if (todproc.engine().cunit == "mJy/beam") {
+                             todproc.engine().cflux(j) = toltec_io.barea_keys[l]*MJY_SR_TO_mJY_ASEC;
+                            if (k == toltec_io.barea_keys.size() - 1) {
+                                k = 0;
+                                l++;
+                             }
+                            else {
+                                k++;
+                            }
+                        }
+                        else if (todproc.engine().cunit == "MJy/Sr") {
+                             todproc.engine().cflux(j) = 1;
+                        }
+                    }
+
+                    SPDLOG_INFO("cflux {}",todproc.engine().cflux);
+
                     // do general setup that is only run once per rawobs before grppi pipeline
                     {
                         tula::logging::scoped_timeit timer("engine setup()");