/** * internal data representation * @author Tobias Weber * @date 1-June-2018 * @license see 'LICENSE' file */ #include "data.h" #include "globals.h" #include "libs/instr.h" #include "libs/mat.h" #include "libs/algos.h" using t_real = t_real_dat; /** * convert an instrument data file to the internal data format */ std::tuple Dataset::convert_instr_file(const char* pcFile) { Dataset dataset; // load instrument data file std::unique_ptr> pInstr(tl2::FileInstrBase::LoadInstr(pcFile)); if(!pInstr) return std::make_tuple(false, dataset); const auto &colnames = pInstr->GetColNames(); const auto &filedata = pInstr->GetData(); if(!pInstr || !colnames.size()) // only valid files with a non-zero column count return std::make_tuple(false, dataset); // process polarisation data pInstr->SetPolNames("p1", "p2", "i11", "i10"); pInstr->ParsePolData(); // get scan axis indices std::vector scan_idx; for(const auto& scanvar : pInstr->GetScannedVars()) { std::size_t idx = 0; pInstr->GetCol(scanvar, &idx); if(idx < colnames.size()) scan_idx.push_back(idx); } // try first axis if none found if(scan_idx.size() == 0) scan_idx.push_back(0); // get counter column index std::vector ctr_idx; { std::size_t idx = 0; pInstr->GetCol(pInstr->GetCountVar(), &idx); if(idx < colnames.size()) ctr_idx.push_back(idx); } // try second axis if none found if(ctr_idx.size() == 0) ctr_idx.push_back(1); // get monitor column index std::vector mon_idx; { std::size_t idx = 0; pInstr->GetCol(pInstr->GetMonVar(), &idx); if(idx < colnames.size()) mon_idx.push_back(idx); } std::size_t numpolstates = pInstr->NumPolChannels(); if(numpolstates == 0) numpolstates = 1; // iterate through all (polarisation) subplots for(std::size_t polstate=0; polstate thedat; tl2::copy_interleave(filedata[idx].begin(), filedata[idx].end(), std::back_inserter(thedat), numpolstates, polstate); data.AddAxis(thedat, colnames[idx]); } // get counters for(std::size_t idx : ctr_idx) { std::vector thedat, theerr; tl2::copy_interleave(filedata[idx].begin(), filedata[idx].end(), std::back_inserter(thedat), numpolstates, polstate); std::transform(thedat.begin(), thedat.end(), std::back_inserter(theerr), [](t_real y) -> t_real { if(tl2::float_equal(y, 0)) return 1; return std::sqrt(y); }); data.AddCounter(thedat, theerr); } // get monitors for(std::size_t idx : mon_idx) { std::vector thedat, theerr; tl2::copy_interleave(filedata[idx].begin(), filedata[idx].end(), std::back_inserter(thedat), numpolstates, polstate); std::transform(thedat.begin(), thedat.end(), std::back_inserter(theerr), [](t_real y) -> t_real { if(tl2::float_equal(y, 0)) return 1; return std::sqrt(y); }); data.AddMonitor(thedat, theerr); } dataset.AddChannel(std::move(data)); } return std::make_tuple(true, dataset); } // ---------------------------------------------------------------------------- // data operators // ---------------------------------------------------------------------------- Data Data::add_pointwise(const Data& dat1, const Data& dat2) { // check if x axes and dimensions are equal constexpr const t_real_dat eps = 0.01; bool compatible = true; compatible = compatible && (dat1.m_x_names == dat2.m_x_names); compatible = compatible && (dat1.m_x.size() == dat2.m_x.size()); if(compatible) { for(std::size_t i=0; i= GetNumMonitors()) { print_err("Invalid monitor selected."); return *this; } Data datret = *this; // normalise all counters for(std::size_t detidx=0; detidx= 1) xlabel = dat.GetAxisName(0); ofstr << "$dat_" << ch << " << ENDDATA\n"; // iterate all x axis names ofstr << "# "; for(std::size_t iX = 0; iX < Nx; ++iX) ofstr << std::setw(iX==0 ? iPrec*1.5 - 2 : iPrec*1.5) << std::left << dat.GetAxisName(iX) << " "; // iterate all counter names for(std::size_t iC = 0; iC < Nc; ++iC) ofstr << std::setw(iPrec*1.5) << std::left << "counter " << " " << std::setw(iPrec*1.5) << std::left << "error" << " "; // iterate all monitors names for(std::size_t iM = 0; iM < Nm; ++iM) ofstr << std::setw(iPrec*1.5) << std::left << "monitor " << " " << std::setw(iPrec*1.5) << std::left << "error" << " "; ofstr << "\n"; // iterate data rows for(std::size_t iRow = 0; iRow < rows; ++iRow) { // iterate all x axes for(std::size_t iX = 0; iX < Nx; ++iX) { const auto& vec = dat.GetAxis(iX); ofstr << std::setw(iPrec*1.5) << std::left << vec[iRow] << " "; } // iterate all counters for(std::size_t iC = 0; iC < Nc; ++iC) { const auto& vec = dat.GetCounter(iC); const auto& vecErr = dat.GetCounterErrors(iC); ofstr << std::setw(iPrec*1.5) << std::left << vec[iRow] << " " << std::setw(iPrec*1.5) << std::left << vecErr[iRow] << " "; } // iterate all monitors for(std::size_t iM = 0; iM < Nm; ++iM) { const auto& vec = dat.GetMonitor(iM); const auto& vecErr = dat.GetMonitorErrors(iM); ofstr << std::setw(iPrec*1.5) << std::left << vec[iRow] << " " << std::setw(iPrec*1.5) << std::left << vecErr[iRow]; } ofstr << "\n"; } ofstr << "ENDDATA\n\n"; } ofstr << "\nset xlabel \"" << xlabel << "\"\n"; ofstr << "set ylabel \"Counts\"\n\n"; ofstr << "\nplot \\\n"; for(std::size_t ch=0; ch