Line data Source code
1 : /*
2 : * SingleDishMSFiller.tcc
3 : *
4 : * Created on: Jan 8, 2016
5 : * Author: nakazato
6 : */
7 :
8 : #ifndef SINGLEDISH_FILLER_SINGLEDISHMSFILLER_TCC_
9 : #define SINGLEDISH_FILLER_SINGLEDISHMSFILLER_TCC_
10 :
11 : #include <singledishfiller/Filler/SingleDishMSFiller.h>
12 :
13 : #include <iostream>
14 : #include <map>
15 : #include <memory>
16 :
17 : #include <casacore/casa/Arrays/Cube.h>
18 : #include <casacore/casa/OS/File.h>
19 : #include <casacore/casa/Utilities/Regex.h>
20 : #include <casacore/casa/Utilities/Sort.h>
21 :
22 : #include <casacore/measures/Measures/Stokes.h>
23 :
24 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
25 : #include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
26 : #include <casacore/ms/MeasurementSets/MSPolColumns.h>
27 : #include <casacore/ms/MeasurementSets/MSPointingColumns.h>
28 : #include <casacore/ms/MeasurementSets/MSStateColumns.h>
29 : #include <casacore/ms/MeasurementSets/MSFeedColumns.h>
30 : #include <casacore/ms/MeasurementSets/MSSourceColumns.h>
31 :
32 : #include <casacore/tables/Tables/Table.h>
33 : #include <casacore/tables/Tables/ScalarColumn.h>
34 : #include <casacore/tables/Tables/SetupNewTab.h>
35 : #include <singledishfiller/Filler/AntennaRecord.h>
36 : #include <singledishfiller/Filler/FieldRecord.h>
37 : #include <singledishfiller/Filler/FillerUtil.h>
38 : #include <singledishfiller/Filler/ObservationRecord.h>
39 : #include <singledishfiller/Filler/ProcessorRecord.h>
40 : #include <singledishfiller/Filler/SourceRecord.h>
41 : #include <singledishfiller/Filler/SpectralWindowRecord.h>
42 : #include <singledishfiller/Filler/SysCalRecord.h>
43 : #include <singledishfiller/Filler/WeatherRecord.h>
44 :
45 : using namespace casacore;
46 :
47 : #define ARRAY_BLOCK_SIZE 1024
48 :
49 : namespace {
50 :
51 : template<class _Table, class _Record, class _Reader>
52 0 : inline void fillTable(_Table &table, _Record &record, _Reader const &reader) {
53 : POST_START;
54 :
55 0 : typename _Record::AssociatingColumns columns(table);
56 :
57 0 : size_t irow = 0;
58 0 : record.clear();
59 0 : for (casacore::Bool more_rows = reader(record); more_rows == true;
60 0 : more_rows = reader(record)) {
61 0 : record.add(table, columns);
62 0 : record.fill(irow, columns);
63 0 : ++irow;
64 0 : record.clear();
65 : }
66 :
67 : POST_END;
68 0 : }
69 :
70 : template<class _Table, class _Columns, class _Comparer, class _Updater>
71 0 : inline casacore::Int updateTable(_Table &mytable, _Columns &mycolumns,
72 : _Comparer const &comparer, _Updater const &updater) {
73 : POST_START;
74 :
75 0 : casacore::Int id = -1;
76 0 : if (mycolumns.nrow() >= (casacore::uInt) INT_MAX) {
77 0 : throw casacore::AipsError("Too much row in table");
78 : }
79 0 : for (casacore::uInt i = 0; i < mycolumns.nrow(); ++i) {
80 0 : if (comparer(mycolumns, i)) {
81 0 : id = (casacore::Int) i;
82 : }
83 : }
84 0 : if (id < 0) {
85 0 : id = mycolumns.nrow();
86 0 : mytable.addRow(1, true);
87 0 : updater(mycolumns, id);
88 : }
89 :
90 : POST_END;
91 0 : return id;
92 : }
93 :
94 : template<class _Columns, class _Record>
95 0 : inline void updateSubtable(_Columns &columns, casacore::uInt irow,
96 : _Record const &record) {
97 : // only update timestamp and interval
98 0 : casacore::Double time_org = columns.time()(irow);
99 0 : casacore::Double interval_org = columns.interval()(irow);
100 :
101 0 : casacore::Double time_min_org = time_org - interval_org / 2.0;
102 0 : casacore::Double time_max_org = time_org + interval_org / 2.0;
103 :
104 0 : casacore::Double time_min_in = record.time - record.interval / 2.0;
105 0 : casacore::Double time_max_in = record.time + record.interval / 2.0;
106 :
107 0 : casacore::Double time_min_new = min(time_min_org, time_min_in);
108 0 : casacore::Double time_max_new = max(time_max_org, time_max_in);
109 :
110 0 : if (time_min_new != time_min_org || time_max_new != time_max_org) {
111 0 : casacore::Double time_new = (time_min_new + time_max_new) / 2.0;
112 0 : casacore::Double interval_new = time_max_new - time_min_new;
113 0 : columns.time().put(irow, time_new);
114 0 : columns.interval().put(irow, interval_new);
115 : }
116 0 : }
117 :
118 0 : inline void makeSourceMap(casacore::MSSource const &table, casacore::Record &source_map) {
119 : POST_START;
120 :
121 0 : casacore::ScalarColumn<casacore::String> name_column(table, "NAME");
122 0 : casacore::ScalarColumn<casacore::Int> id_column(table, "SOURCE_ID");
123 0 : casacore::Vector<casacore::Int> id = id_column.getColumn();
124 :
125 0 : casacore::Sort sorter;
126 0 : sorter.sortKey(id.data(), TpInt);
127 0 : casacore::Vector<casacore::uInt> unique_vector;
128 0 : casacore::uInt num_id = sorter.sort(unique_vector, id.nelements(),
129 : casacore::Sort::HeapSort | casacore::Sort::NoDuplicates);
130 0 : for (casacore::uInt i = 0; i < num_id; ++i) {
131 0 : casacore::uInt irow = unique_vector[i];
132 0 : casacore::String const source_name = name_column(irow);
133 0 : casacore::Int const source_id = id[irow];
134 0 : source_map.define(source_name, source_id);
135 0 : }
136 :
137 : POST_END;
138 0 : }
139 :
140 : } // anonymous namespace
141 :
142 : using namespace casacore;
143 :
144 : namespace casa { //# NAMESPACE CASA - BEGIN
145 :
146 : using namespace sdfiller;
147 :
148 : static constexpr ssize_t CONTEXT_BUFFER_SIZE = 10;
149 : static constexpr ssize_t FILLER_STORAGE_SIZE = CONTEXT_BUFFER_SIZE + 2;
150 : typedef sdfiller::ProducerConsumerModelContext<ssize_t, CONTEXT_BUFFER_SIZE> PCMContext;
151 :
152 : extern PCMContext *g_context_p;
153 : extern casa::sdfiller::DataRecord *g_storage_p;
154 :
155 : template<class T>
156 0 : void SingleDishMSFiller<T>::create_context() {
157 : static_assert(0 < FILLER_STORAGE_SIZE && FILLER_STORAGE_SIZE < SSIZE_MAX,
158 : "FILLER_STORAGE_SIZE is wrong value");
159 : static_assert(CONTEXT_BUFFER_SIZE + 2 <= FILLER_STORAGE_SIZE,
160 : "FILLER_STORAGE_SIZE < CONTEXT_BUFFER_SIZE + 2");
161 :
162 0 : constexpr ssize_t END_OF_PRODUCTION = -1;
163 0 : if (!g_context_p) {
164 0 : g_context_p = new PCMContext(END_OF_PRODUCTION);
165 : }
166 0 : if (!g_storage_p) {
167 0 : g_storage_p = new DataRecord[FILLER_STORAGE_SIZE];
168 : }
169 0 : }
170 :
171 : template<class T>
172 0 : void SingleDishMSFiller<T>::destroy_context() {
173 0 : if (g_context_p) {
174 0 : delete g_context_p;
175 0 : g_context_p = nullptr;
176 : }
177 0 : if (g_storage_p) {
178 0 : delete[] g_storage_p;
179 0 : g_storage_p = nullptr;
180 : }
181 0 : }
182 :
183 : template<class T>
184 0 : void *SingleDishMSFiller<T>::consume(void *arg) {
185 : POST_START;
186 :
187 : try {
188 0 : auto filler = reinterpret_cast<SingleDishMSFiller<T> *>(arg);
189 0 : auto reader = filler->reader_.get();
190 :
191 : // std::ostringstream oss;
192 :
193 0 : size_t nrow = reader->getNumberOfRows();
194 0 : assert(nrow < SIZE_MAX);
195 0 : DataAccumulator accumulator;
196 :
197 : // oss << "consume: nrow = " << nrow;
198 : // PCMContext::locked_print(oss.str(), g_context_p);
199 :
200 : ssize_t storage_index;
201 0 : for (size_t irow = 0; irow < nrow + 1; ++irow) {
202 : // oss.str("");
203 : // oss << "consume: start row " << irow;
204 : // PCMContext::locked_print(oss.str(), g_context_p);
205 0 : bool more_products = PCMContext::consume(g_context_p, &storage_index);
206 0 : assert(storage_index < FILLER_STORAGE_SIZE);
207 :
208 : // oss.str("");
209 : // oss << "consume index " << storage_index << " more_products = "
210 : // << more_products;
211 : // PCMContext::locked_print(oss.str(), g_context_p);
212 :
213 0 : if (more_products) {
214 0 : DataRecord *record = &g_storage_p[storage_index];
215 :
216 : // oss.str("");
217 : // oss << "Accumulate record: time = " << record->time << " interval = "
218 : // << record->interval;
219 : // PCMContext::locked_print(oss.str(), g_context_p);
220 :
221 0 : casacore::Bool is_ready = accumulator.queryForGet(record->time);
222 0 : if (is_ready) {
223 0 : filler->flush(accumulator);
224 : }
225 0 : casacore::Bool astatus = accumulator.accumulate(*record);
226 : (void) astatus;
227 :
228 : // oss.str("");
229 : // oss << "astatus = " << astatus;
230 : // PCMContext::locked_print(oss.str(), g_context_p);
231 : } else {
232 0 : break;
233 : }
234 : }
235 :
236 : // PCMContext::locked_print("Final flush", g_context_p);
237 0 : filler->flush(accumulator);
238 0 : } catch (std::runtime_error &e) {
239 0 : std::ostringstream oss;
240 0 : oss << "Exception: " << e.what();
241 0 : PCMContext::locked_print(oss.str(), g_context_p);
242 0 : } catch (...) {
243 0 : PCMContext::locked_print("Unknown exception", g_context_p);
244 : }
245 :
246 : POST_END;
247 0 : pthread_exit(0);
248 : }
249 :
250 : template<class T>
251 0 : void *SingleDishMSFiller<T>::produce(void *arg) {
252 : POST_START;
253 :
254 : try {
255 0 : auto filler = reinterpret_cast<SingleDishMSFiller<T> *>(arg);
256 0 : auto reader = filler->reader_.get();
257 :
258 : // std::ostringstream oss;
259 :
260 0 : size_t nrow = reader->getNumberOfRows();
261 :
262 : // oss << "produce: nrow = " << nrow;
263 : // PCMContext::locked_print(oss.str(), g_context_p);
264 :
265 0 : ssize_t storage_index = 0;
266 :
267 0 : for (size_t irow = 0; irow < nrow; ++irow) {
268 :
269 0 : DataRecord *record = &g_storage_p[storage_index];
270 0 : casacore::Bool status = reader->getData(irow, *record);
271 :
272 : // oss.str("");
273 : // oss << "irow " << irow << " status " << status << std::endl;
274 : // oss << " TIME=" << record->time << " INTERVAL=" << record->interval
275 : // << std::endl;
276 : // oss << "status = " << status;
277 : // PCMContext::locked_print(oss.str(), g_context_p);
278 :
279 0 : if (status) {
280 : // oss.str("");
281 : // oss << "produce index " << storage_index;
282 : // PCMContext::locked_print(oss.str(), g_context_p);
283 :
284 0 : PCMContext::produce(g_context_p, storage_index);
285 : } else {
286 0 : break;
287 : }
288 :
289 0 : storage_index++;
290 0 : storage_index %= FILLER_STORAGE_SIZE;
291 : }
292 :
293 : // PCMContext::locked_print("Done production", g_context_p);
294 :
295 0 : PCMContext::complete_production(g_context_p);
296 :
297 0 : } catch (std::runtime_error &e) {
298 0 : std::ostringstream oss;
299 0 : oss << "Exception: " << e.what();
300 0 : PCMContext::locked_print(oss.str(), g_context_p);
301 :
302 0 : PCMContext::complete_production(g_context_p);
303 0 : } catch (...) {
304 0 : PCMContext::locked_print("Unknown exception", g_context_p);
305 :
306 0 : PCMContext::complete_production(g_context_p);
307 : }
308 :
309 : POST_END;
310 0 : pthread_exit(0);
311 : }
312 :
313 : template<class T>
314 0 : void SingleDishMSFiller<T>::fillMainMT(SingleDishMSFiller<T> *filler) {
315 : POST_START;
316 :
317 : pthread_t consumer_id;
318 : pthread_t producer_id;
319 :
320 : // create context
321 0 : SingleDishMSFiller<T>::create_context();
322 :
323 : try {
324 0 : sdfiller::create_thread(&consumer_id, NULL, SingleDishMSFiller<T>::consume,
325 : filler);
326 0 : sdfiller::create_thread(&producer_id, NULL, SingleDishMSFiller<T>::produce,
327 : filler);
328 :
329 0 : sdfiller::join_thread(&consumer_id, NULL);
330 0 : sdfiller::join_thread(&producer_id, NULL);
331 0 : } catch (...) {
332 0 : SingleDishMSFiller<T>::destroy_context();
333 0 : throw;
334 : }
335 :
336 0 : SingleDishMSFiller<T>::destroy_context();
337 :
338 : POST_END;
339 0 : }
340 :
341 : template<class T>
342 0 : SingleDishMSFiller<T>::SingleDishMSFiller(std::string const &name,
343 : bool parallel) :
344 0 : ms_(), ms_columns_(), data_description_columns_(), feed_columns_(),
345 0 : pointing_columns_(), polarization_columns_(), syscal_columns_(),
346 0 : state_columns_(), weather_columns_(), reader_(new T(name)),
347 0 : is_float_(false), data_key_(), reference_feed_(-1), pointing_time_(),
348 0 : pointing_time_max_(), pointing_time_min_(), num_pointing_time_(),
349 0 : syscal_list_(), subscan_list_(), polarization_type_pool_(), weather_list_(),
350 0 : parallel_(parallel) {
351 0 : }
352 :
353 : template<class T>
354 0 : SingleDishMSFiller<T>::~SingleDishMSFiller() {
355 0 : }
356 :
357 : template<class T>
358 0 : void SingleDishMSFiller<T>::fill() {
359 : POST_START;
360 :
361 : // initialization
362 0 : initialize();
363 :
364 : // Fill tables that can be processed prior to main loop
365 0 : fillPreProcessTables();
366 :
367 : // main loop
368 0 : casacore::LogIO os(casacore::LogOrigin("SingleDishMSFiller", "fill", WHERE));
369 0 : if (parallel_) {
370 0 : os << "Parallel execution of fillMain" << casacore::LogIO::POST;
371 0 : SingleDishMSFiller<T>::fillMainMT(this);
372 : } else {
373 0 : os << "Serial execution of fillMain" << casacore::LogIO::POST;
374 0 : fillMain();
375 : }
376 :
377 : // Fill tables that must be processed after main loop
378 0 : fillPostProcessTables();
379 :
380 : // finalization
381 0 : finalize();
382 :
383 : POST_END;
384 0 : }
385 :
386 : template<class T>
387 0 : void SingleDishMSFiller<T>::initialize() {
388 : POST_START;
389 :
390 : // initialize reader
391 0 : reader_->initialize();
392 :
393 : // query if the data is complex or float
394 0 : is_float_ = reader_->isFloatData();
395 0 : if (is_float_) {
396 : // std::cout << "data column is FLOAT_DATA" << std::endl;
397 0 : data_key_ = "FLOAT_DATA";
398 : } else {
399 : // std::cout << "data column is DATA" << std::endl;
400 0 : data_key_ = "DATA";
401 : }
402 :
403 : // setup MS
404 0 : setupMS();
405 :
406 : // frame information
407 0 : casacore::MDirection::Types direction_frame = reader_->getDirectionFrame();
408 0 : auto mytable = ms_->pointing();
409 0 : casacore::ArrayColumn<casacore::Double> direction_column(mytable, "DIRECTION");
410 0 : casacore::TableRecord &record = direction_column.rwKeywordSet();
411 0 : casacore::Record meas_info = record.asRecord("MEASINFO");
412 0 : casacore::String ref_string = casacore::MDirection::showType(direction_frame);
413 0 : meas_info.define("Ref", ref_string);
414 0 : record.defineRecord("MEASINFO", meas_info);
415 :
416 : POST_END;
417 0 : }
418 :
419 : template<class T>
420 0 : void SingleDishMSFiller<T>::finalize() {
421 : POST_START;
422 :
423 : // finalize reader
424 0 : reader_->finalize();
425 :
426 : POST_END;
427 0 : }
428 :
429 : template<class T>
430 0 : void SingleDishMSFiller<T>::save(std::string const &name) {
431 : POST_START;
432 :
433 : #ifdef SINGLEDISHMSFILLER_DEBUG
434 : std::cout << "Saving MS as \"" << name << "\"" << std::endl;
435 : std::cout << "current working directory is \"" << casacore::Path().absoluteName()
436 : << "\"" << std::endl;
437 : #endif
438 :
439 0 : ms_->deepCopy(name, casacore::Table::New);
440 :
441 : #ifdef SINGLEDISHMSFILLER_DEBUG
442 : casacore::File file(name);
443 : casacore::Bool name_exists = file.exists();
444 : if (name_exists) {
445 : std::cout << "file successfully created" << std::endl;
446 : } else {
447 : std::cout << "failed to create file" << std::endl;
448 : }
449 : #endif
450 :
451 : POST_END;
452 0 : }
453 :
454 : template<class T>
455 0 : void SingleDishMSFiller<T>::setupMS() {
456 : // std::cout << "Start " << __PRETTY_FUNCTION__ << std::endl;
457 :
458 0 : casacore::String const dunit = reader_->getDataUnit();
459 :
460 0 : casacore::TableDesc ms_main_description = casacore::MeasurementSet::requiredTableDesc();
461 0 : if (is_float_) {
462 0 : casacore::MeasurementSet::addColumnToDesc(ms_main_description,
463 : casacore::MSMainEnums::FLOAT_DATA, 2);
464 0 : ms_main_description.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().define("UNIT", dunit);
465 : } else {
466 0 : casacore::MeasurementSet::addColumnToDesc(ms_main_description, casacore::MSMainEnums::DATA, 2);
467 0 : ms_main_description.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().define("UNIT", dunit);
468 : }
469 :
470 0 : casacore::String const scratch_table_name = casacore::File::newUniqueName(".",
471 : "SingleDishMSFillerTemp").originalName();
472 0 : casacore::SetupNewTable newtab(scratch_table_name, ms_main_description, casacore::Table::Scratch);
473 :
474 0 : ms_.reset(new casacore::MeasurementSet(newtab));
475 :
476 : // create subtables
477 0 : casacore::TableDesc ms_antenna_description = casacore::MSAntenna::requiredTableDesc();
478 0 : casacore::SetupNewTable ms_antenna_table(ms_->antennaTableName(),
479 : ms_antenna_description, casacore::Table::Scratch);
480 0 : ms_->rwKeywordSet().defineTable(
481 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::ANTENNA),
482 0 : casacore::Table(ms_antenna_table));
483 :
484 0 : casacore::TableDesc ms_data_desc_description = casacore::MSDataDescription::requiredTableDesc();
485 0 : casacore::SetupNewTable ms_data_desc_table(ms_->dataDescriptionTableName(),
486 : ms_data_desc_description, casacore::Table::Scratch);
487 0 : ms_->rwKeywordSet().defineTable(
488 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::DATA_DESCRIPTION),
489 0 : casacore::Table(ms_data_desc_table));
490 :
491 0 : casacore::TableDesc ms_doppler_description = casacore::MSDoppler::requiredTableDesc();
492 0 : casacore::SetupNewTable ms_doppler_table(ms_->dopplerTableName(),
493 : ms_doppler_description, casacore::Table::Scratch);
494 0 : ms_->rwKeywordSet().defineTable(
495 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::DOPPLER),
496 0 : casacore::Table(ms_doppler_table));
497 :
498 0 : casacore::TableDesc ms_feed_description = casacore::MSFeed::requiredTableDesc();
499 0 : casacore::SetupNewTable ms_feed_table(ms_->feedTableName(), ms_feed_description,
500 : casacore::Table::Scratch);
501 0 : ms_->rwKeywordSet().defineTable(
502 0 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FEED), casacore::Table(ms_feed_table));
503 :
504 0 : casacore::TableDesc ms_field_description = casacore::MSField::requiredTableDesc();
505 0 : casacore::SetupNewTable ms_field_table(ms_->fieldTableName(), ms_field_description,
506 : casacore::Table::Scratch);
507 0 : ms_->rwKeywordSet().defineTable(
508 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FIELD),
509 0 : casacore::Table(ms_field_table));
510 :
511 0 : casacore::TableDesc ms_flag_cmd_description = casacore::MSFlagCmd::requiredTableDesc();
512 0 : casacore::SetupNewTable ms_flag_cmd_table(ms_->flagCmdTableName(),
513 : ms_flag_cmd_description, casacore::Table::Scratch);
514 0 : ms_->rwKeywordSet().defineTable(
515 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FLAG_CMD),
516 0 : casacore::Table(ms_flag_cmd_table));
517 :
518 0 : casacore::TableDesc ms_freq_offset_description = casacore::MSFreqOffset::requiredTableDesc();
519 0 : casacore::SetupNewTable ms_freq_offset_table(ms_->freqOffsetTableName(),
520 : ms_freq_offset_description, casacore::Table::Scratch);
521 0 : ms_->rwKeywordSet().defineTable(
522 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FREQ_OFFSET),
523 0 : casacore::Table(ms_freq_offset_table));
524 :
525 0 : casacore::TableDesc ms_history_description = casacore::MSHistory::requiredTableDesc();
526 0 : casacore::SetupNewTable ms_history_table(ms_->historyTableName(),
527 : ms_history_description, casacore::Table::Scratch);
528 0 : ms_->rwKeywordSet().defineTable(
529 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::HISTORY),
530 0 : casacore::Table(ms_history_table));
531 :
532 0 : casacore::TableDesc ms_observation_description = casacore::MSObservation::requiredTableDesc();
533 0 : casacore::SetupNewTable ms_observation_table(ms_->observationTableName(),
534 : ms_observation_description, casacore::Table::Scratch);
535 0 : ms_->rwKeywordSet().defineTable(
536 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::OBSERVATION),
537 0 : casacore::Table(ms_observation_table));
538 :
539 0 : casacore::TableDesc ms_pointing_description = casacore::MSPointing::requiredTableDesc();
540 0 : casacore::SetupNewTable ms_pointing_table(ms_->pointingTableName(),
541 : ms_pointing_description, casacore::Table::Scratch);
542 0 : ms_->rwKeywordSet().defineTable(
543 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::POINTING),
544 0 : casacore::Table(ms_pointing_table));
545 :
546 0 : casacore::TableDesc ms_polarization_description = casacore::MSPolarization::requiredTableDesc();
547 0 : casacore::SetupNewTable ms_polarization_table(ms_->polarizationTableName(),
548 : ms_polarization_description, casacore::Table::Scratch);
549 0 : ms_->rwKeywordSet().defineTable(
550 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::POLARIZATION),
551 0 : casacore::Table(ms_polarization_table));
552 :
553 0 : casacore::TableDesc ms_processor_description = casacore::MSProcessor::requiredTableDesc();
554 0 : casacore::SetupNewTable ms_processor_table(ms_->processorTableName(),
555 : ms_processor_description, casacore::Table::Scratch);
556 0 : ms_->rwKeywordSet().defineTable(
557 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::PROCESSOR),
558 0 : casacore::Table(ms_processor_table));
559 :
560 0 : casacore::TableDesc ms_source_description = casacore::MSSource::requiredTableDesc();
561 0 : casacore::MSSource::addColumnToDesc(ms_source_description, casacore::MSSourceEnums::TRANSITION,
562 : 1);
563 0 : casacore::MSSource::addColumnToDesc(ms_source_description,
564 : casacore::MSSourceEnums::REST_FREQUENCY, 1);
565 0 : casacore::MSSource::addColumnToDesc(ms_source_description, casacore::MSSourceEnums::SYSVEL, 1);
566 0 : casacore::SetupNewTable ms_source_table(ms_->sourceTableName(), ms_source_description,
567 : casacore::Table::Scratch);
568 0 : ms_->rwKeywordSet().defineTable(
569 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SOURCE),
570 0 : casacore::Table(ms_source_table));
571 :
572 0 : casacore::TableDesc ms_spectral_window_description =
573 : casacore::MSSpectralWindow::requiredTableDesc();
574 0 : casacore::SetupNewTable ms_spectral_window_table(ms_->spectralWindowTableName(),
575 : ms_spectral_window_description, casacore::Table::Scratch);
576 0 : ms_->rwKeywordSet().defineTable(
577 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SPECTRAL_WINDOW),
578 0 : casacore::Table(ms_spectral_window_table));
579 :
580 0 : casacore::TableDesc ms_state_description = casacore::MSState::requiredTableDesc();
581 0 : casacore::SetupNewTable ms_state_table(ms_->stateTableName(), ms_state_description,
582 : casacore::Table::Scratch);
583 0 : ms_->rwKeywordSet().defineTable(
584 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::STATE),
585 0 : casacore::Table(ms_state_table));
586 :
587 0 : casacore::TableDesc ms_syscal_description = casacore::MSSysCal::requiredTableDesc();
588 0 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TCAL_SPECTRUM,
589 : 2);
590 0 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TCAL, 1);
591 0 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TSYS_SPECTRUM,
592 : 2);
593 0 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TSYS, 1);
594 0 : casacore::SetupNewTable ms_syscal_table(ms_->sysCalTableName(), ms_syscal_description,
595 : casacore::Table::Scratch);
596 0 : ms_->rwKeywordSet().defineTable(
597 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SYSCAL),
598 0 : casacore::Table(ms_syscal_table));
599 :
600 0 : casacore::TableDesc ms_weather_description = casacore::MSWeather::requiredTableDesc();
601 0 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
602 : casacore::MSWeatherEnums::TEMPERATURE);
603 0 : casacore::MSWeather::addColumnToDesc(ms_weather_description, casacore::MSWeatherEnums::PRESSURE);
604 0 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
605 : casacore::MSWeatherEnums::REL_HUMIDITY);
606 0 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
607 : casacore::MSWeatherEnums::WIND_SPEED);
608 0 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
609 : casacore::MSWeatherEnums::WIND_DIRECTION);
610 0 : casacore::SetupNewTable ms_weather_table(ms_->weatherTableName(),
611 : ms_weather_description, casacore::Table::Scratch);
612 0 : ms_->rwKeywordSet().defineTable(
613 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::WEATHER),
614 0 : casacore::Table(ms_weather_table));
615 :
616 0 : ms_->initRefs();
617 :
618 : // Set up MSMainColumns
619 0 : ms_columns_.reset(new casacore::MSMainColumns(*ms_));
620 :
621 : // Set up MSDataDescColumns
622 0 : data_description_columns_.reset(
623 0 : new casacore::MSDataDescColumns(ms_->dataDescription()));
624 :
625 : // Set up MSFeedColumns
626 0 : feed_columns_.reset(new casacore::MSFeedColumns(ms_->feed()));
627 :
628 : // Set up MSPointingColumns
629 0 : pointing_columns_.reset(new casacore::MSPointingColumns(ms_->pointing()));
630 :
631 : // Set up MSPolarizationColumns
632 0 : polarization_columns_.reset(new casacore::MSPolarizationColumns(ms_->polarization()));
633 :
634 : // Set up MSSysCalColumns
635 0 : syscal_columns_.reset(new casacore::MSSysCalColumns(ms_->sysCal()));
636 :
637 : // Set up MSStateColumns
638 0 : state_columns_.reset(new casacore::MSStateColumns(ms_->state()));
639 :
640 : // Set up MSWeatherColumns
641 0 : weather_columns_.reset(new casacore::MSWeatherColumns(ms_->weather()));
642 :
643 : // std::cout << "End " << __PRETTY_FUNCTION__ << std::endl;
644 0 : }
645 :
646 : template<class T>
647 0 : void SingleDishMSFiller<T>::fillPreProcessTables() {
648 : POST_START;
649 :
650 : // fill OBSERVATION table
651 0 : fillObservation();
652 :
653 : // fill ANTENNA table
654 0 : fillAntenna();
655 :
656 : // fill PROCESSOR table
657 0 : fillProcessor();
658 :
659 : // fill SOURCE table
660 0 : fillSource();
661 :
662 : // fill FIELD table
663 0 : fillField();
664 :
665 : // fill SPECTRAL_WINDOW table
666 0 : fillSpectralWindow();
667 :
668 : // Generate optional tables if necessary
669 0 : OptionalTables::Generate(*ms_, *reader_);
670 :
671 : POST_END;
672 0 : }
673 :
674 : template<class T>
675 0 : void SingleDishMSFiller<T>::fillPostProcessTables() {
676 : POST_START;
677 :
678 : // fill HISTORY table
679 0 : fillHistory();
680 :
681 : // flush POINTING entry
682 0 : sortPointing();
683 :
684 : POST_END;
685 0 : }
686 :
687 : template<class T>
688 0 : void SingleDishMSFiller<T>::fillMain() {
689 : POST_START;
690 :
691 0 : size_t nrow = reader_->getNumberOfRows();
692 0 : DataAccumulator accumulator;
693 0 : DataRecord record;
694 : // std::cout << "nrow = " << nrow << std::endl;
695 0 : for (size_t irow = 0; irow < nrow; ++irow) {
696 0 : casacore::Bool status = reader_->getData(irow, record);
697 : // std::cout << "irow " << irow << " status " << status << std::endl;
698 : // std::cout << " TIME=" << record.time << " INTERVAL=" << record.interval
699 : // << std::endl;
700 : // std::cout << "status = " << status << std::endl;
701 0 : if (status) {
702 0 : casacore::Bool is_ready = accumulator.queryForGet(record.time);
703 0 : if (is_ready) {
704 0 : flush(accumulator);
705 : }
706 0 : casacore::Bool astatus = accumulator.accumulate(record);
707 : (void) astatus;
708 : // std::cout << "astatus = " << astatus << std::endl;
709 : }
710 : }
711 :
712 0 : flush(accumulator);
713 :
714 : POST_END;
715 0 : }
716 :
717 : template<class T>
718 0 : void SingleDishMSFiller<T>::fillAntenna() {
719 : POST_START;
720 :
721 0 : auto mytable = ms_->antenna();
722 0 : AntennaRecord record;
723 0 : ::fillTable(mytable, record,
724 0 : [&](AntennaRecord &record) {return reader_->getAntennaRow(record);});
725 :
726 : // initialize POINTING table related stuff
727 0 : casacore::uInt nrow = mytable.nrow();
728 0 : num_pointing_time_.resize(nrow);
729 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
730 0 : pointing_time_[i] = casacore::Vector<casacore::Double>(ARRAY_BLOCK_SIZE, -1.0);
731 0 : pointing_time_min_[i] = -1.0;
732 0 : pointing_time_max_[i] = 1.0e30;
733 0 : num_pointing_time_[i] = 0;
734 : }
735 :
736 : POST_END;
737 0 : }
738 :
739 : template<class T>
740 0 : void SingleDishMSFiller<T>::fillObservation() {
741 : POST_START;
742 :
743 0 : ObservationRecord record;
744 0 : ::fillTable(ms_->observation(), record,
745 0 : [&](ObservationRecord &record) {return reader_->getObservationRow(record);});
746 :
747 : POST_END;
748 0 : }
749 :
750 : template<class T>
751 0 : void SingleDishMSFiller<T>::fillProcessor() {
752 : POST_START;
753 :
754 0 : ProcessorRecord record;
755 0 : ::fillTable(ms_->processor(), record,
756 0 : [&](ProcessorRecord &record) {return reader_->getProcessorRow(record);});
757 :
758 : POST_END;
759 0 : }
760 :
761 : template<class T>
762 0 : void SingleDishMSFiller<T>::fillSource() {
763 : POST_START;
764 :
765 0 : SourceRecord record;
766 0 : ::fillTable(ms_->source(), record,
767 0 : [&](SourceRecord &record) {return reader_->getSourceRow(record);});
768 :
769 : POST_END;
770 0 : }
771 :
772 : template<class T>
773 0 : void SingleDishMSFiller<T>::fillField() {
774 : POST_START;
775 :
776 0 : auto mytable = ms_->field();
777 0 : FieldRecord record;
778 0 : record.table = mytable;
779 :
780 : // make (name,id) map for SOURCE table
781 0 : ::makeSourceMap(ms_->source(), record.source_map);
782 :
783 0 : ::fillTable(mytable, record,
784 0 : [&](FieldRecord &record) {return reader_->getFieldRow(record);});
785 :
786 : POST_END;
787 0 : }
788 :
789 : template<class T>
790 0 : void SingleDishMSFiller<T>::fillHistory() {
791 : POST_START;
792 :
793 : // HISTORY table should be filled by upper-level
794 : // application command (e.g. importscantable)
795 : // ms_->history().addRow(1, true);
796 : // casacore::Vector<casacore::String> cols(2);
797 : // cols[0] = "APP_PARAMS";
798 : // cols[1] = "CLI_COMMAND";
799 : // casacore::TableRow row(ms_->history(), cols, true);
800 : // // TODO: fill HISTORY row here
801 : // casacore::TableRecord record = row.record();
802 : // record.print(std::cout);
803 : // row.put(0, record);
804 :
805 : POST_END;
806 0 : }
807 :
808 : template<class T>
809 0 : void SingleDishMSFiller<T>::fillSpectralWindow() {
810 : POST_START;
811 :
812 0 : auto mytable = ms_->spectralWindow();
813 0 : SpectralWindowRecord record;
814 :
815 0 : ::fillTable(mytable, record,
816 0 : [&](SpectralWindowRecord &record) {return reader_->getSpectralWindowRow(record);});
817 :
818 : POST_END;
819 0 : }
820 :
821 : template<class T>
822 0 : casacore::Int SingleDishMSFiller<T>::updatePolarization(casacore::Vector<casacore::Int> const &corr_type,
823 : casacore::Int const &num_pol) {
824 0 : casacore::uInt num_corr = corr_type.size();
825 0 : if (num_pol < 1 || num_corr != (uInt) num_pol) {
826 0 : throw casacore::AipsError("Internal inconsistency in number of correlations");
827 : }
828 0 : casacore::MSPolarization &mytable = ms_->polarization();
829 : //casacore::MSPolarizationColumns mycolumns(mytable);
830 0 : casacore::Matrix<casacore::Int> const corr_product(2, num_pol, 0);
831 0 : auto comparer = [&](casacore::MSPolarizationColumns const &columns, casacore::uInt i) {
832 0 : casacore::Bool match = allEQ(columns.corrType()(i), corr_type);
833 0 : return match;
834 : };
835 0 : auto updater = [&](casacore::MSPolarizationColumns &columns, casacore::uInt i) {
836 0 : columns.numCorr().put(i, num_pol);
837 0 : columns.corrType().put(i, corr_type);
838 0 : columns.corrProduct().put(i, corr_product);
839 : };
840 0 : casacore::Int polarization_id = ::updateTable(mytable, *(polarization_columns_.get()),
841 : comparer, updater);
842 0 : return polarization_id;
843 0 : }
844 :
845 : template<class T>
846 0 : Int SingleDishMSFiller<T>::updateDataDescription(casacore::Int const &polarization_id,
847 : casacore::Int const &spw_id) {
848 0 : if (polarization_id < 0 || spw_id < 0) {
849 0 : throw casacore::AipsError("Invalid ids for DATA_DESCRIPTION");
850 : }
851 0 : casacore::MSDataDescription &mytable = ms_->dataDescription();
852 : //casacore::MSDataDescColumns mycolumns(mytable);
853 0 : auto comparer = [&](casacore::MSDataDescColumns const &columns, casacore::uInt i) {
854 0 : casacore::Bool match = (columns.polarizationId()(i) == polarization_id)
855 0 : && (columns.spectralWindowId()(i) == spw_id);
856 0 : return match;
857 : };
858 0 : auto updater = [&](casacore::MSDataDescColumns &columns, casacore::uInt i) {
859 0 : columns.polarizationId().put(i, polarization_id);
860 0 : columns.spectralWindowId().put(i, spw_id);
861 : };
862 0 : casacore::Int data_desc_id = ::updateTable(mytable, *(data_description_columns_.get()),
863 : comparer, updater);
864 :
865 0 : return data_desc_id;
866 : }
867 :
868 : template<class T>
869 0 : Int SingleDishMSFiller<T>::updateState(casacore::Int const &subscan,
870 : casacore::String const &obs_mode) {
871 0 : casacore::MSState &mytable = ms_->state();
872 0 : static casacore::Regex const regex("^OBSERVE_TARGET#ON_SOURCE");
873 : //static std::vector<casacore::Int> subscan_list;
874 0 : auto comparer =
875 0 : [&](casacore::MSStateColumns &columns, casacore::uInt i) {
876 0 : casacore::Bool match = (subscan == subscan_list_[i]) && (obs_mode == columns.obsMode()(i));
877 0 : return match;
878 : };
879 0 : auto updater = [&](casacore::MSStateColumns &columns, casacore::uInt i) {
880 0 : columns.subScan().put(i, subscan);
881 0 : columns.obsMode().put(i, obs_mode);
882 0 : casacore::Bool is_signal = obs_mode.matches(regex);
883 0 : columns.sig().put(i, is_signal);
884 0 : columns.ref().put(i, !is_signal);
885 :
886 0 : subscan_list_.push_back(subscan);
887 : };
888 0 : casacore::Int state_id = ::updateTable(mytable, *(state_columns_.get()), comparer,
889 : updater);
890 0 : return state_id;
891 : }
892 :
893 : template<class T>
894 0 : Int SingleDishMSFiller<T>::updateFeed(casacore::Int const &feed_id, casacore::Int const &spw_id,
895 : casacore::String const &pol_type) {
896 0 : constexpr casacore::Int num_receptors = 2;
897 0 : static casacore::String const linear_type_arr[2] = { "X", "Y" };
898 0 : static casacore::Vector<casacore::String> linear_type(linear_type_arr, 2, SHARE);
899 0 : static casacore::String const circular_type_arr[2] = { "R", "L" };
900 0 : static casacore::Vector<casacore::String> circular_type(circular_type_arr, 2, SHARE);
901 0 : static casacore::Matrix<casacore::Complex> const pol_response(num_receptors, num_receptors,
902 0 : casacore::Complex(0));
903 0 : casacore::Vector<casacore::String> *polarization_type = nullptr;
904 0 : if (pol_type == "linear") {
905 0 : polarization_type = &linear_type;
906 0 : } else if (pol_type == "circular") {
907 0 : polarization_type = &circular_type;
908 : } else {
909 0 : polarization_type = &linear_type;
910 : }
911 : //static std::vector< casacore::Vector<casacore::String> *> polarization_type_pool;
912 :
913 0 : casacore::String polarization_type_arr[2] = { "X", "Y" };
914 0 : casacore::Vector<casacore::String> polarization_type_storage(polarization_type_arr, 2, SHARE);
915 0 : casacore::Matrix<casacore::Double> const beam_offset(2, num_receptors, 0.0);
916 0 : casacore::Vector<casacore::Double> const receptor_angle(num_receptors, 0.0);
917 0 : casacore::Vector<casacore::Double> const position(3, 0.0);
918 0 : auto comparer = [&](casacore::MSFeedColumns &columns, casacore::uInt i) {
919 0 : casacore::Vector<casacore::String> *current_polarization_type = polarization_type_pool_[i];
920 0 : casacore::Bool match = allEQ(*polarization_type, *current_polarization_type) &&
921 0 : (feed_id == columns.feedId()(i)) &&
922 0 : (spw_id == columns.spectralWindowId()(i));
923 0 : return match;
924 : };
925 0 : auto updater = [&](casacore::MSFeedColumns &columns, casacore::uInt i) {
926 : // TODO: 2016/01/26 TN
927 : // Here I regard "multi-beam receiver" as multi-feed, single-beam receivers
928 : // in casacore::MS v2 notation. So, BEAM_ID is always 0 and FEED_ID is beam identifier
929 : // for "multi-beam" receivers so far. However, if we decide to import beams
930 : // as separated virtual antennas, FEED_ID must always be 0.
931 0 : columns.feedId().put(i, feed_id);
932 0 : columns.beamId().put(i, 0);
933 0 : columns.spectralWindowId().put(i, spw_id);
934 0 : columns.polarizationType().put(i, *polarization_type);
935 0 : columns.beamOffset().put(i, beam_offset);
936 0 : columns.receptorAngle().put(i, receptor_angle);
937 0 : columns.position().put(i, position);
938 0 : columns.polResponse().put(i, pol_response);
939 :
940 0 : polarization_type_pool_.push_back(polarization_type);
941 : };
942 0 : casacore::Int feed_row = ::updateTable(ms_->feed(), *(feed_columns_.get()), comparer,
943 : updater);
944 :
945 : // reference feed
946 0 : if (reference_feed_ < 0) {
947 0 : reference_feed_ = feed_id;
948 : }
949 :
950 0 : return feed_row;
951 0 : }
952 :
953 : template<class T>
954 0 : Int SingleDishMSFiller<T>::updatePointing(casacore::Int const &antenna_id,
955 : casacore::Int const &feed_id, casacore::Double const &time, casacore::Double const &interval,
956 : casacore::Matrix<casacore::Double> const &direction) {
957 : POST_START;
958 :
959 0 : if (reference_feed_ != feed_id) {
960 0 : return -1;
961 : }
962 :
963 0 : auto &mytable = ms_->pointing();
964 0 : casacore::uInt nrow = mytable.nrow();
965 :
966 0 : casacore::uInt *n = &num_pointing_time_[antenna_id];
967 0 : casacore::Double *time_max = &pointing_time_max_[antenna_id];
968 0 : casacore::Double *time_min = &pointing_time_min_[antenna_id];
969 0 : casacore::Vector<casacore::Double> *time_list = &pointing_time_[antenna_id];
970 :
971 0 : auto addPointingRow =
972 0 : [&]() {
973 0 : mytable.addRow(1, true);
974 0 : pointing_columns_->time().put(nrow, time);
975 0 : pointing_columns_->interval().put(nrow, interval);
976 0 : pointing_columns_->antennaId().put(nrow, antenna_id);
977 0 : if (direction.ncolumn() == 1 || allEQ(direction.column(1), 0.0)) {
978 0 : pointing_columns_->numPoly().put(nrow, 0);
979 0 : pointing_columns_->direction().put(nrow, direction(casacore::IPosition(2,0,0), casacore::IPosition(2,1,0)));
980 : } else {
981 0 : pointing_columns_->direction().put(nrow, direction);
982 0 : casacore::Int num_poly = direction.shape()[1] - 1;
983 0 : pointing_columns_->numPoly().put(nrow, num_poly);
984 : }
985 : // add timestamp to the list
986 0 : casacore::uInt nelem = time_list->nelements();
987 0 : if (nelem <= *n) {
988 0 : time_list->resize(nelem + ARRAY_BLOCK_SIZE, true);
989 : }
990 0 : (*time_list)[*n] = time;
991 : // increment number of pointing entry
992 0 : *n += 1;
993 : };
994 :
995 0 : if (*n == 0) {
996 0 : addPointingRow();
997 0 : *time_min = time;
998 0 : *time_max = time;
999 0 : } else if (time < *time_min) {
1000 0 : addPointingRow();
1001 0 : *time_min = time;
1002 0 : } else if (*time_max < time) {
1003 0 : addPointingRow();
1004 0 : *time_max = time;
1005 0 : } else if (allNE(*time_list, time)) {
1006 0 : addPointingRow();
1007 : }
1008 :
1009 : POST_END;
1010 :
1011 0 : return -1;
1012 : }
1013 :
1014 : template<class T>
1015 0 : void SingleDishMSFiller<T>::updateWeather(casacore::Int const &antenna_id,
1016 : casacore::Double const &time, casacore::Double const &interval,
1017 : MSDataRecord const &data_record) {
1018 : WeatherRecord record;
1019 0 : record.clear();
1020 0 : record.antenna_id = antenna_id;
1021 0 : record.time = time;
1022 0 : record.interval = interval;
1023 0 : record.temperature = data_record.temperature;
1024 0 : record.pressure = data_record.pressure;
1025 0 : record.rel_humidity = data_record.rel_humidity;
1026 0 : record.wind_speed = data_record.wind_speed;
1027 0 : record.wind_direction = data_record.wind_direction;
1028 0 : auto &mytable = ms_->weather();
1029 0 : auto pos = std::find(weather_list_.begin(), weather_list_.end(), record);
1030 0 : if (pos == weather_list_.end()) {
1031 0 : weather_list_.push_back(record);
1032 0 : casacore::uInt irow = mytable.nrow();
1033 0 : mytable.addRow(1, true);
1034 0 : record.fill(irow, *(weather_columns_.get()));
1035 : } else {
1036 0 : auto irow = std::distance(weather_list_.begin(), pos);
1037 0 : updateWeather(*(weather_columns_.get()), irow, record);
1038 : }
1039 0 : }
1040 :
1041 : template<class T>
1042 0 : void SingleDishMSFiller<T>::updateWeather(casacore::MSWeatherColumns &columns, casacore::uInt irow,
1043 : WeatherRecord const &record) {
1044 0 : ::updateSubtable(columns, irow, record);
1045 0 : }
1046 :
1047 : template<class T>
1048 0 : void SingleDishMSFiller<T>::updateSysCal(casacore::Int const &antenna_id,
1049 : casacore::Int const &feed_id, casacore::Int const &spw_id, casacore::Double const &time,
1050 : casacore::Double const &interval, MSDataRecord const &data_record) {
1051 : POST_START;
1052 :
1053 0 : SysCalRecord record;
1054 0 : record.clear();
1055 0 : record.antenna_id = antenna_id;
1056 0 : record.feed_id = feed_id;
1057 0 : record.spw_id = spw_id;
1058 0 : record.time = time;
1059 0 : record.interval = interval;
1060 :
1061 : //casacore::Bool tcal_empty = false;
1062 0 : casacore::Bool tsys_empty = false;
1063 :
1064 0 : if (data_record.tcal.empty() || allEQ(data_record.tcal, 1.0f)
1065 0 : || allEQ(data_record.tcal, 0.0f)) {
1066 : //tcal_empty = true;
1067 : } else {
1068 : // std::cout << "tcal seems to be valid " << data_record.tcal << std::endl;
1069 0 : if (data_record.float_data.shape() == data_record.tcal.shape()) {
1070 0 : record.tcal_spectrum.assign(data_record.tcal);
1071 : } else {
1072 0 : casacore::Matrix<casacore::Float> tcal = data_record.tcal;
1073 0 : if (!tcal.empty()) {
1074 0 : record.tcal.assign(tcal.column(0));
1075 : }
1076 0 : }
1077 : }
1078 0 : if (data_record.tsys.empty() || allEQ(data_record.tsys, 1.0f)
1079 0 : || allEQ(data_record.tsys, 0.0f)) {
1080 0 : tsys_empty = true;
1081 : } else {
1082 0 : if (data_record.float_data.shape() == data_record.tsys.shape()) {
1083 0 : record.tsys_spectrum.assign(data_record.tsys);
1084 : } else {
1085 0 : casacore::Matrix<casacore::Float> tsys = data_record.tsys;
1086 0 : if (!tsys.empty()) {
1087 0 : record.tsys.assign(tsys.column(0));
1088 : }
1089 0 : }
1090 : }
1091 :
1092 : // do not add entry if Tsys is empty
1093 : //if (tcal_empty && tsys_empty) {
1094 0 : if (tsys_empty) {
1095 0 : return;
1096 : }
1097 :
1098 0 : auto &mytable = ms_->sysCal();
1099 0 : auto pos = std::find(syscal_list_.begin(), syscal_list_.end(), record);
1100 0 : if (pos == syscal_list_.end()) {
1101 0 : casacore::uInt irow = mytable.nrow();
1102 0 : mytable.addRow(1, true);
1103 0 : record.fill(irow, *(syscal_columns_.get()));
1104 0 : syscal_list_.push_back(SysCalTableRecord(ms_.get(), irow, record));
1105 : } else {
1106 0 : auto irow = std::distance(syscal_list_.begin(), pos);
1107 0 : updateSysCal(*(syscal_columns_.get()), irow, record);
1108 : }
1109 :
1110 : POST_END;
1111 0 : }
1112 :
1113 : template<class T>
1114 0 : void SingleDishMSFiller<T>::updateSysCal(casacore::MSSysCalColumns &columns, casacore::uInt irow,
1115 : SysCalRecord const &record) {
1116 0 : ::updateSubtable(columns, irow, record);
1117 0 : }
1118 :
1119 : template<class T>
1120 0 : void SingleDishMSFiller<T>::updateMain(casacore::Int const &antenna_id, casacore::Int field_id,
1121 : casacore::Int feedId, casacore::Int dataDescriptionId, casacore::Int stateId, casacore::Int const &scan_number,
1122 : casacore::Double const &time, MSDataRecord const &dataRecord) {
1123 : POST_START;
1124 :
1125 : // constant stuff
1126 0 : static casacore::Vector<casacore::Double> const uvw(3, 0.0);
1127 0 : static casacore::Array<casacore::Bool> const flagCategory(casacore::IPosition(3, 0, 0, 0));
1128 :
1129 : // target row id
1130 0 : casacore::uInt irow = ms_->nrow();
1131 :
1132 : // add new row
1133 : //ms_->addRow(1, true);
1134 0 : ms_->addRow(1, false);
1135 :
1136 0 : ms_columns_->uvw().put(irow, uvw);
1137 0 : ms_columns_->flagCategory().put(irow, flagCategory);
1138 0 : ms_columns_->antenna1().put(irow, antenna_id);
1139 0 : ms_columns_->antenna2().put(irow, antenna_id);
1140 0 : ms_columns_->fieldId().put(irow, field_id);
1141 0 : ms_columns_->feed1().put(irow, feedId);
1142 0 : ms_columns_->feed2().put(irow, feedId);
1143 0 : ms_columns_->dataDescId().put(irow, dataDescriptionId);
1144 0 : ms_columns_->stateId().put(irow, stateId);
1145 0 : ms_columns_->scanNumber().put(irow, scan_number);
1146 0 : ms_columns_->time().put(irow, time);
1147 0 : ms_columns_->timeCentroid().put(irow, time);
1148 0 : casacore::Double const &interval = dataRecord.interval;
1149 0 : ms_columns_->interval().put(irow, interval);
1150 0 : ms_columns_->exposure().put(irow, interval);
1151 :
1152 0 : if (is_float_) {
1153 0 : casacore::Matrix<casacore::Float> floatData;
1154 0 : if (dataRecord.isFloat()) {
1155 0 : floatData.reference(dataRecord.float_data);
1156 : } else {
1157 0 : floatData.assign(real(dataRecord.complex_data));
1158 : }
1159 0 : ms_columns_->floatData().put(irow, floatData);
1160 0 : } else {
1161 0 : casacore::Matrix<casacore::Complex> data;
1162 0 : if (dataRecord.isFloat()) {
1163 0 : data.assign(
1164 : makeComplex(dataRecord.float_data,
1165 0 : casacore::Matrix<casacore::Float>(dataRecord.float_data.shape(), 0.0f)));
1166 : } else {
1167 0 : data.reference(dataRecord.complex_data);
1168 : }
1169 0 : ms_columns_->data().put(irow, data);
1170 0 : }
1171 :
1172 0 : ms_columns_->flag().put(irow, dataRecord.flag);
1173 0 : ms_columns_->flagRow().put(irow, dataRecord.flag_row);
1174 0 : ms_columns_->sigma().put(irow, dataRecord.sigma);
1175 0 : ms_columns_->weight().put(irow, dataRecord.weight);
1176 :
1177 : POST_END;
1178 0 : }
1179 :
1180 : template<class T>
1181 0 : void SingleDishMSFiller<T>::flush(DataAccumulator &accumulator) {
1182 : POST_START;
1183 :
1184 0 : size_t nchunk = accumulator.getNumberOfChunks();
1185 : // std::cout << "nchunk = " << nchunk << std::endl;
1186 :
1187 0 : if (nchunk == 0) {
1188 0 : return;
1189 : }
1190 :
1191 0 : for (size_t ichunk = 0; ichunk < nchunk; ++ichunk) {
1192 0 : casacore::Bool status = accumulator.get(ichunk, record_);
1193 : // std::cout << "accumulator status = " << std::endl;
1194 0 : if (status) {
1195 0 : casacore::Double time = record_.time;
1196 0 : casacore::Int antenna_id = record_.antenna_id;
1197 0 : casacore::Int spw_id = record_.spw_id;
1198 0 : casacore::Int feed_id = record_.feed_id;
1199 0 : casacore::Int field_id = record_.field_id;
1200 0 : casacore::Int scan = record_.scan;
1201 0 : casacore::Int subscan = record_.subscan;
1202 0 : casacore::String pol_type = record_.pol_type;
1203 0 : casacore::String obs_mode = record_.intent;
1204 0 : casacore::Int num_pol = record_.num_pol;
1205 0 : casacore::Vector<casacore::Int> &corr_type = record_.corr_type;
1206 0 : casacore::Int polarization_id = updatePolarization(corr_type, num_pol);
1207 0 : updateFeed(feed_id, spw_id, pol_type);
1208 0 : casacore::Int data_desc_id = updateDataDescription(polarization_id, spw_id);
1209 0 : casacore::Int state_id = updateState(subscan, obs_mode);
1210 0 : casacore::Matrix<casacore::Double> &direction = record_.direction;
1211 0 : casacore::Double interval = record_.interval;
1212 :
1213 : // updatePointing must be called after updateFeed
1214 0 : updatePointing(antenna_id, feed_id, time, interval, direction);
1215 :
1216 0 : updateSysCal(antenna_id, feed_id, spw_id, time, interval, record_);
1217 :
1218 0 : updateWeather(antenna_id, time, interval, record_);
1219 :
1220 0 : updateMain(antenna_id, field_id, feed_id, data_desc_id, state_id, scan,
1221 0 : time, record_);
1222 0 : }
1223 : }
1224 : // std::cout << "clear accumulator" << std::endl;
1225 0 : accumulator.clear();
1226 :
1227 : POST_END;
1228 : }
1229 :
1230 : template<class T>
1231 0 : void SingleDishMSFiller<T>::sortPointing() {
1232 : POST_START;
1233 :
1234 : // deallocate POINTING table related stuff
1235 0 : pointing_time_.clear();
1236 0 : pointing_time_min_.clear();
1237 0 : pointing_time_max_.clear();
1238 0 : num_pointing_time_.resize();
1239 :
1240 0 : auto mytable = ms_->pointing();
1241 0 : casacore::MSPointingColumns mycolumns(mytable);
1242 0 : casacore::uInt nrow = mytable.nrow();
1243 0 : casacore::Vector<casacore::Int> antenna_id_list = mycolumns.antennaId().getColumn();
1244 0 : casacore::Vector<casacore::Double> time_list = mycolumns.time().getColumn();
1245 0 : casacore::Sort sorter;
1246 0 : sorter.sortKey(antenna_id_list.data(), TpInt);
1247 0 : sorter.sortKey(time_list.data(), TpDouble);
1248 0 : casacore::Vector<casacore::uInt> index_vector;
1249 0 : sorter.sort(index_vector, nrow);
1250 :
1251 0 : size_t storage_size = nrow * 2 * sizeof(casacore::Double);
1252 0 : std::unique_ptr<void, sdfiller::Deleter> storage(malloc(storage_size));
1253 :
1254 : // sort TIME
1255 : {
1256 0 : casacore::Vector<casacore::Double> sorted(casacore::IPosition(1, nrow),
1257 0 : reinterpret_cast<casacore::Double *>(storage.get()), SHARE);
1258 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1259 0 : sorted[i] = time_list[index_vector[i]];
1260 : }
1261 0 : mycolumns.time().putColumn(sorted);
1262 0 : }
1263 : // sort ANTENNA_ID
1264 : {
1265 0 : casacore::Vector<casacore::Int> sorted(casacore::IPosition(1, nrow),
1266 0 : reinterpret_cast<casacore::Int *>(storage.get()), SHARE);
1267 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1268 0 : sorted[i] = antenna_id_list[index_vector[i]];
1269 : }
1270 0 : mycolumns.antennaId().putColumn(sorted);
1271 0 : }
1272 :
1273 : // sort NUM_POLY
1274 : {
1275 0 : casacore::Vector<casacore::Int> num_poly_list(antenna_id_list);
1276 0 : mycolumns.numPoly().getColumn(num_poly_list);
1277 0 : casacore::Vector<casacore::Int> sorted(casacore::IPosition(1, nrow),
1278 0 : reinterpret_cast<casacore::Int *>(storage.get()), SHARE);
1279 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1280 0 : sorted[i] = antenna_id_list[index_vector[i]];
1281 : }
1282 0 : mycolumns.numPoly().putColumn(sorted);
1283 0 : }
1284 :
1285 : // sort DIRECTION
1286 : {
1287 0 : std::map<casacore::uInt, casacore::Matrix<casacore::Double> > direction;
1288 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1289 0 : direction[i] = mycolumns.direction()(i);
1290 : }
1291 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1292 0 : mycolumns.direction().put(i, direction[index_vector[i]]);
1293 : }
1294 0 : }
1295 :
1296 : // sort INTERVAL
1297 : {
1298 0 : casacore::Vector<casacore::Double> interval_list = mycolumns.interval().getColumn();
1299 0 : for (casacore::uInt i = 0; i < nrow; ++i) {
1300 0 : mycolumns.interval().put(i, interval_list[index_vector[i]]);
1301 : }
1302 0 : }
1303 : POST_END;
1304 0 : }
1305 :
1306 : } //# NAMESPACE CASA - END
1307 :
1308 : #endif /* SINGLEDISH_FILLER_SINGLEDISHMSFILLER_TCC_ */
|
