Skip to content

Commit

Permalink
add exp time, adjust noise maps, add bool coverage, adjust wiener fil…
Browse files Browse the repository at this point in the history 
…ter, jinc map, beammap adjustments, wcs changes
  • Loading branch information
@mmccrackan
mmccrackan committed May 17, 2022
1 parent fda1a39 commit 7861b4b
Show file tree
Hide file tree
Showing 20 changed files with 581 additions and 149 deletions.
19 changes: 11 additions & 8 deletions include/citlali/core/engine/beammap.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,12 +156,14 @@ auto Beammap::run_timestream() {

in.tel_meta_data.data["SourceEl"] = tel_meta_data["SourceEl"].segment(start_index, scan_length);

Eigen::Matrix<Eigen::Index, Eigen::Dynamic, 1> map_index_vector, det_index_vector;
//Eigen::Matrix<Eigen::Index, Eigen::Dynamic, 1> map_index_vector, det_index_vector;

auto [map_index_vector, det_index_vector] = polarization.create_rtc(in, in, "I", this);

/*Stage 1: RTCProc*/
RTCProc rtcproc;
TCData<TCDataKind::PTC,Eigen::MatrixXd> out;
rtcproc.run(in, out, det_index_vector, this);
rtcproc.run(in, out, map_index_vector, det_index_vector, this);

// move out into the PTCData vector
ptcs0.at(out.index.data - 1) = std::move(out);
Expand Down Expand Up @@ -453,7 +455,6 @@ void Beammap::output(MC &mout, fits_out_vec_t &f_ios, fits_out_vec_t & nf_ios, b
table.row(ti+4) = mout.pfit.row(ci).template cast <float> ();
table.row(ti+4 + 1) = mout.perror.row(ci).template cast <float> ();
ci++;
SPDLOG_INFO("ci {} ti {}", ci, ti);
}

table.row(toltec_io.beammap_apt_header.size()-1) = converge_iter.cast <float> ();
Expand All @@ -471,13 +472,13 @@ void Beammap::output(MC &mout, fits_out_vec_t &f_ios, fits_out_vec_t & nf_ios, b
meta["nw"].push_back("units: N/A");
meta["nw"].push_back("network index");

meta["flxscale"].push_back("units: Mjy/sr");
meta["flxscale"].push_back("units: MJy/Sr");
meta["flxscale"].push_back("flux conversion scale");

meta["amp"].push_back("units: Mjy/sr");
meta["amp"].push_back("units: MJy/Sr");
meta["amp"].push_back("fitted amplitude");

meta["amp_err"].push_back("units: Mjy/sr");
meta["amp_err"].push_back("units: MJy/Sr");
meta["amp_err"].push_back("fitted amplitude error");

meta["x_t"].push_back("units: arcsec");
Expand Down Expand Up @@ -553,8 +554,8 @@ void Beammap::output(MC &mout, fits_out_vec_t &f_ios, fits_out_vec_t & nf_ios, b
pixel_size,source_center,toltec_io.array_freqs[i],
polarization.stokes_params,hdu_name);
// add fit parameters to hdus
hdu->addKey("amp", (float)mout.pfit(0,i),"amplitude (Mjy/sr)");
hdu->addKey("amp_err", (float)mout.perror(0,i),"amplitude error (Mjy/sr)");
hdu->addKey("amp", (float)mout.pfit(0,i),"amplitude (MJy/Sr)");
hdu->addKey("amp_err", (float)mout.perror(0,i),"amplitude error (MJy/Sr)");
hdu->addKey("x_t", (float)mout.pfit(1,i),"az offset (arcsec)");
hdu->addKey("x_t_err", (float)mout.perror(1,i),"az offset error (arcsec)");
hdu->addKey("y_t", (float)mout.pfit(2,i),"alt offset (arcsec)");
Expand Down Expand Up @@ -590,6 +591,8 @@ void Beammap::output(MC &mout, fits_out_vec_t &f_ios, fits_out_vec_t & nf_ios, b
f_ios.at(i).pfits->pHDU().addKey("WAV", toltec_io.name_keys[i], "Array Name");
// add obsnum
f_ios.at(i).pfits->pHDU().addKey("OBSNUM", obsnum, "Observation Number");
// exp time
f_ios.at(i).pfits->pHDU().addKey("t_exptime", tel_meta_params["t_exp"], "Exposure Time (sec)");
// add units
f_ios.at(i).pfits->pHDU().addKey("UNIT", obsnum, "MJy/Sr");
// add conversion
Expand Down
10 changes: 10 additions & 0 deletions include/citlali/core/engine/calib.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
struct ToltecCalib {
// name and scale
std::map<std::string, double> header_keys = {
{"uid", 1},
{"nw", 1},
{"array", 1},
{"flxscale", 1},
Expand Down Expand Up @@ -99,7 +100,16 @@ class Calib: public ToltecCalib {
throw DataIOError{fmt::format(
"failed to load data from netCDF file {}", filepath)};
}
}

void check_missing() {
for (Eigen::Index i=0; i<calib_data["uid"].size(); i++) {
std::ostringstream os;

os << calib_data["uid"](i);
std::string digits = os.str();

}
}
};

48 changes: 41 additions & 7 deletions include/citlali/core/engine/engine.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
#include <citlali/core/timestream/ptc/ptcproc.h>

#include <citlali/core/map/naive_mm_2.h>
#include <citlali/core/map/jinc_mm.h>
#include <citlali/core/map/jinc_mm_2.h>
#include <citlali/core/map/wiener_filter.h>

struct reduControls {
Expand Down Expand Up @@ -132,6 +132,9 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
int max_iterations;
double cutoff;

// jinc map-maker
Eigen::VectorXd radius, jinc_weights;

// control for fitting
bool run_fit;

Expand Down Expand Up @@ -263,6 +266,8 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
engine_config = _c;
SPDLOG_INFO("getting config options");

ToltecIO toltec_io;

get_config(run_tod_output,std::tuple{"timestream","output","enabled"});
if (run_tod_output) {
get_config(ts_format,std::tuple{"timestream","output","format"});
Expand All @@ -276,7 +281,7 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
if (ex_name != "seq") {
SPDLOG_INFO("timestream output requires sequential execution policy. setting policy to seq");
ex_name = "seq";
}
}
}

get_config(nthreads,std::tuple{"runtime","n_threads"});
Expand Down Expand Up @@ -334,6 +339,11 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
if (run_kernel) {
get_config(kernel.filepath,std::tuple{"timestream","kernel","filepath"});
get_config(kernel.kernel_type,std::tuple{"timestream","kernel","type"},{"internal_gaussian","internal_airy","image"});

auto kernel_node = engine_config.get_node(std::tuple{"timestream","kernel","image_ext_name"});
auto kernel_node_size = kernel_node.size();

std::string hdu_name;
get_config(kernel.hdu_name,std::tuple{"timestream","kernel","image_ext_name"});
kernel.setup();
}
Expand All @@ -349,11 +359,26 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
// get mapmaking config options
get_config(map_grouping,std::tuple{"mapmaking","grouping"});
get_config(mapping_method,std::tuple{"mapmaking","method"});

if (mapping_method == "jinc") {

double r_max = 1.5;
double a = 1.1;
double b = 4.75;
double c = 2.;

radius = Eigen::VectorXd::LinSpaced(1000, 1e-10, r_max);
jinc_weights.resize(radius.size());

for (Eigen::Index i=0; i<radius.size(); i++) {
jinc_weights(i) = jinc_func(radius(i),a,b,c,r_max,1);
}
}

get_config(map_type,std::tuple{"mapmaking","pixel_axes"});
get_config(pixel_size,std::tuple{"mapmaking","pixel_size_arcsec"});
get_config(cmb.cov_cut,std::tuple{"coadd","cov_cut"});


// convert pixel size to radians at start
pixel_size *= RAD_ASEC;
// assign pixel size to map buffer and coadded map buffer for convenience
Expand All @@ -367,7 +392,18 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
get_config(y_size_pix,std::tuple{"mapmaking","y_size_pix"});
get_config(crpix1,std::tuple{"mapmaking","crpix1"});
get_config(crpix2,std::tuple{"mapmaking","crpix2"});
get_config(cunit,std::tuple{"mapmaking","cunit"});
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"});
Expand Down Expand Up @@ -404,7 +440,7 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
get_config(coadd_filter_type,std::tuple{"coadd","filtering","type"});

if (run_noise == false) {
SPDLOG_INFO("noise maps are needed for map filtering.");
SPDLOG_ERROR("noise maps are needed for map filtering.");
std::exit(EXIT_FAILURE);
}

Expand All @@ -418,9 +454,7 @@ class EngineBase: public reduControls, public reduClasses, public Telescope, pub
get_config(gaussian_template_fwhm_rad["a2000"],std::tuple{"wiener_filter","gaussian_template_fwhm_arcsec","a2000"});

for (auto const& pair : gaussian_template_fwhm_rad) {
SPDLOG_INFO("gaussian_template_fwhm_rad[pair.first] {}",gaussian_template_fwhm_rad[pair.first]);
gaussian_template_fwhm_rad[pair.first] = gaussian_template_fwhm_rad[pair.first]*RAD_ASEC;
SPDLOG_INFO("gaussian_template_fwhm_rad[pair.first] {}",gaussian_template_fwhm_rad[pair.first]);
}

wiener_filter.run_kernel = run_kernel;
Expand Down
Loading

0 comments on commit 7861b4b

Please sign in to comment.