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 24 : inline void fillTable(_Table &table, _Record &record, _Reader const &reader) {
53 : POST_START;
54 :
55 24 : typename _Record::AssociatingColumns columns(table);
56 :
57 24 : size_t irow = 0;
58 24 : record.clear();
59 80 : for (casacore::Bool more_rows = reader(record); more_rows == true;
60 56 : more_rows = reader(record)) {
61 56 : record.add(table, columns);
62 56 : record.fill(irow, columns);
63 56 : ++irow;
64 56 : record.clear();
65 : }
66 :
67 : POST_END;
68 24 : }
69 :
70 : template<class _Table, class _Columns, class _Comparer, class _Updater>
71 1888 : inline casacore::Int updateTable(_Table &mytable, _Columns &mycolumns,
72 : _Comparer const &comparer, _Updater const &updater) {
73 : POST_START;
74 :
75 1888 : casacore::Int id = -1;
76 1888 : if (mycolumns.nrow() >= (casacore::uInt) INT_MAX) {
77 0 : throw casacore::AipsError("Too much row in table");
78 : }
79 7648 : for (casacore::uInt i = 0; i < mycolumns.nrow(); ++i) {
80 5760 : if (comparer(mycolumns, i)) {
81 1824 : id = (casacore::Int) i;
82 : }
83 : }
84 1888 : if (id < 0) {
85 64 : id = mycolumns.nrow();
86 64 : mytable.addRow(1, true);
87 64 : updater(mycolumns, id);
88 : }
89 :
90 : POST_END;
91 1888 : return id;
92 : }
93 :
94 : template<class _Columns, class _Record>
95 528 : inline void updateSubtable(_Columns &columns, casacore::uInt irow,
96 : _Record const &record) {
97 : // only update timestamp and interval
98 528 : casacore::Double time_org = columns.time()(irow);
99 528 : casacore::Double interval_org = columns.interval()(irow);
100 :
101 528 : casacore::Double time_min_org = time_org - interval_org / 2.0;
102 528 : casacore::Double time_max_org = time_org + interval_org / 2.0;
103 :
104 528 : casacore::Double time_min_in = record.time - record.interval / 2.0;
105 528 : casacore::Double time_max_in = record.time + record.interval / 2.0;
106 :
107 528 : casacore::Double time_min_new = min(time_min_org, time_min_in);
108 528 : casacore::Double time_max_new = max(time_max_org, time_max_in);
109 :
110 528 : if (time_min_new != time_min_org || time_max_new != time_max_org) {
111 356 : casacore::Double time_new = (time_min_new + time_max_new) / 2.0;
112 356 : casacore::Double interval_new = time_max_new - time_min_new;
113 356 : columns.time().put(irow, time_new);
114 356 : columns.interval().put(irow, interval_new);
115 : }
116 528 : }
117 :
118 4 : inline void makeSourceMap(casacore::MSSource const &table, casacore::Record &source_map) {
119 : POST_START;
120 :
121 4 : casacore::ScalarColumn<casacore::String> name_column(table, "NAME");
122 4 : casacore::ScalarColumn<casacore::Int> id_column(table, "SOURCE_ID");
123 4 : casacore::Vector<casacore::Int> id = id_column.getColumn();
124 :
125 4 : casacore::Sort sorter;
126 4 : sorter.sortKey(id.data(), TpInt);
127 4 : casacore::Vector<casacore::uInt> unique_vector;
128 4 : casacore::uInt num_id = sorter.sort(unique_vector, id.nelements(),
129 : casacore::Sort::HeapSort | casacore::Sort::NoDuplicates);
130 8 : for (casacore::uInt i = 0; i < num_id; ++i) {
131 4 : casacore::uInt irow = unique_vector[i];
132 4 : casacore::String const source_name = name_column(irow);
133 4 : casacore::Int const source_id = id[irow];
134 4 : source_map.define(source_name, source_id);
135 4 : }
136 :
137 : POST_END;
138 4 : }
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 4 : SingleDishMSFiller<T>::SingleDishMSFiller(std::string const &name,
343 : bool parallel) :
344 4 : ms_(), ms_columns_(), data_description_columns_(), feed_columns_(),
345 4 : pointing_columns_(), polarization_columns_(), syscal_columns_(),
346 4 : state_columns_(), weather_columns_(), reader_(new T(name)),
347 4 : is_float_(false), data_key_(), reference_feed_(-1), pointing_time_(),
348 4 : pointing_time_max_(), pointing_time_min_(), num_pointing_time_(),
349 4 : syscal_list_(), subscan_list_(), polarization_type_pool_(), weather_list_(),
350 4 : parallel_(parallel) {
351 4 : }
352 :
353 : template<class T>
354 4 : SingleDishMSFiller<T>::~SingleDishMSFiller() {
355 4 : }
356 :
357 : template<class T>
358 4 : void SingleDishMSFiller<T>::fill() {
359 : POST_START;
360 :
361 : // initialization
362 4 : initialize();
363 :
364 : // Fill tables that can be processed prior to main loop
365 4 : fillPreProcessTables();
366 :
367 : // main loop
368 8 : casacore::LogIO os(casacore::LogOrigin("SingleDishMSFiller", "fill", WHERE));
369 4 : if (parallel_) {
370 0 : os << "Parallel execution of fillMain" << casacore::LogIO::POST;
371 0 : SingleDishMSFiller<T>::fillMainMT(this);
372 : } else {
373 4 : os << "Serial execution of fillMain" << casacore::LogIO::POST;
374 4 : fillMain();
375 : }
376 :
377 : // Fill tables that must be processed after main loop
378 4 : fillPostProcessTables();
379 :
380 : // finalization
381 4 : finalize();
382 :
383 : POST_END;
384 4 : }
385 :
386 : template<class T>
387 4 : void SingleDishMSFiller<T>::initialize() {
388 : POST_START;
389 :
390 : // initialize reader
391 4 : reader_->initialize();
392 :
393 : // query if the data is complex or float
394 4 : is_float_ = reader_->isFloatData();
395 4 : if (is_float_) {
396 : // std::cout << "data column is FLOAT_DATA" << std::endl;
397 4 : 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 4 : setupMS();
405 :
406 : // frame information
407 4 : casacore::MDirection::Types direction_frame = reader_->getDirectionFrame();
408 4 : auto mytable = ms_->pointing();
409 4 : casacore::ArrayColumn<casacore::Double> direction_column(mytable, "DIRECTION");
410 4 : casacore::TableRecord &record = direction_column.rwKeywordSet();
411 4 : casacore::Record meas_info = record.asRecord("MEASINFO");
412 4 : casacore::String ref_string = casacore::MDirection::showType(direction_frame);
413 4 : meas_info.define("Ref", ref_string);
414 4 : record.defineRecord("MEASINFO", meas_info);
415 :
416 : POST_END;
417 4 : }
418 :
419 : template<class T>
420 4 : void SingleDishMSFiller<T>::finalize() {
421 : POST_START;
422 :
423 : // finalize reader
424 4 : reader_->finalize();
425 :
426 : POST_END;
427 4 : }
428 :
429 : template<class T>
430 4 : 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 4 : 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 4 : }
453 :
454 : template<class T>
455 4 : void SingleDishMSFiller<T>::setupMS() {
456 : // std::cout << "Start " << __PRETTY_FUNCTION__ << std::endl;
457 :
458 4 : casacore::String const dunit = reader_->getDataUnit();
459 :
460 4 : casacore::TableDesc ms_main_description = casacore::MeasurementSet::requiredTableDesc();
461 4 : if (is_float_) {
462 4 : casacore::MeasurementSet::addColumnToDesc(ms_main_description,
463 : casacore::MSMainEnums::FLOAT_DATA, 2);
464 4 : 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 8 : casacore::String const scratch_table_name = casacore::File::newUniqueName(".",
471 : "SingleDishMSFillerTemp").originalName();
472 4 : casacore::SetupNewTable newtab(scratch_table_name, ms_main_description, casacore::Table::Scratch);
473 :
474 4 : ms_.reset(new casacore::MeasurementSet(newtab));
475 :
476 : // create subtables
477 4 : casacore::TableDesc ms_antenna_description = casacore::MSAntenna::requiredTableDesc();
478 4 : casacore::SetupNewTable ms_antenna_table(ms_->antennaTableName(),
479 : ms_antenna_description, casacore::Table::Scratch);
480 8 : ms_->rwKeywordSet().defineTable(
481 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::ANTENNA),
482 8 : casacore::Table(ms_antenna_table));
483 :
484 4 : casacore::TableDesc ms_data_desc_description = casacore::MSDataDescription::requiredTableDesc();
485 4 : casacore::SetupNewTable ms_data_desc_table(ms_->dataDescriptionTableName(),
486 : ms_data_desc_description, casacore::Table::Scratch);
487 8 : ms_->rwKeywordSet().defineTable(
488 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::DATA_DESCRIPTION),
489 8 : casacore::Table(ms_data_desc_table));
490 :
491 4 : casacore::TableDesc ms_doppler_description = casacore::MSDoppler::requiredTableDesc();
492 4 : casacore::SetupNewTable ms_doppler_table(ms_->dopplerTableName(),
493 : ms_doppler_description, casacore::Table::Scratch);
494 8 : ms_->rwKeywordSet().defineTable(
495 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::DOPPLER),
496 8 : casacore::Table(ms_doppler_table));
497 :
498 4 : casacore::TableDesc ms_feed_description = casacore::MSFeed::requiredTableDesc();
499 4 : casacore::SetupNewTable ms_feed_table(ms_->feedTableName(), ms_feed_description,
500 : casacore::Table::Scratch);
501 8 : ms_->rwKeywordSet().defineTable(
502 8 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FEED), casacore::Table(ms_feed_table));
503 :
504 4 : casacore::TableDesc ms_field_description = casacore::MSField::requiredTableDesc();
505 4 : casacore::SetupNewTable ms_field_table(ms_->fieldTableName(), ms_field_description,
506 : casacore::Table::Scratch);
507 8 : ms_->rwKeywordSet().defineTable(
508 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FIELD),
509 8 : casacore::Table(ms_field_table));
510 :
511 4 : casacore::TableDesc ms_flag_cmd_description = casacore::MSFlagCmd::requiredTableDesc();
512 4 : casacore::SetupNewTable ms_flag_cmd_table(ms_->flagCmdTableName(),
513 : ms_flag_cmd_description, casacore::Table::Scratch);
514 8 : ms_->rwKeywordSet().defineTable(
515 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FLAG_CMD),
516 8 : casacore::Table(ms_flag_cmd_table));
517 :
518 4 : casacore::TableDesc ms_freq_offset_description = casacore::MSFreqOffset::requiredTableDesc();
519 4 : casacore::SetupNewTable ms_freq_offset_table(ms_->freqOffsetTableName(),
520 : ms_freq_offset_description, casacore::Table::Scratch);
521 8 : ms_->rwKeywordSet().defineTable(
522 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::FREQ_OFFSET),
523 8 : casacore::Table(ms_freq_offset_table));
524 :
525 4 : casacore::TableDesc ms_history_description = casacore::MSHistory::requiredTableDesc();
526 4 : casacore::SetupNewTable ms_history_table(ms_->historyTableName(),
527 : ms_history_description, casacore::Table::Scratch);
528 8 : ms_->rwKeywordSet().defineTable(
529 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::HISTORY),
530 8 : casacore::Table(ms_history_table));
531 :
532 4 : casacore::TableDesc ms_observation_description = casacore::MSObservation::requiredTableDesc();
533 4 : casacore::SetupNewTable ms_observation_table(ms_->observationTableName(),
534 : ms_observation_description, casacore::Table::Scratch);
535 8 : ms_->rwKeywordSet().defineTable(
536 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::OBSERVATION),
537 8 : casacore::Table(ms_observation_table));
538 :
539 4 : casacore::TableDesc ms_pointing_description = casacore::MSPointing::requiredTableDesc();
540 4 : casacore::SetupNewTable ms_pointing_table(ms_->pointingTableName(),
541 : ms_pointing_description, casacore::Table::Scratch);
542 8 : ms_->rwKeywordSet().defineTable(
543 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::POINTING),
544 8 : casacore::Table(ms_pointing_table));
545 :
546 4 : casacore::TableDesc ms_polarization_description = casacore::MSPolarization::requiredTableDesc();
547 4 : casacore::SetupNewTable ms_polarization_table(ms_->polarizationTableName(),
548 : ms_polarization_description, casacore::Table::Scratch);
549 8 : ms_->rwKeywordSet().defineTable(
550 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::POLARIZATION),
551 8 : casacore::Table(ms_polarization_table));
552 :
553 4 : casacore::TableDesc ms_processor_description = casacore::MSProcessor::requiredTableDesc();
554 4 : casacore::SetupNewTable ms_processor_table(ms_->processorTableName(),
555 : ms_processor_description, casacore::Table::Scratch);
556 8 : ms_->rwKeywordSet().defineTable(
557 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::PROCESSOR),
558 8 : casacore::Table(ms_processor_table));
559 :
560 4 : casacore::TableDesc ms_source_description = casacore::MSSource::requiredTableDesc();
561 4 : casacore::MSSource::addColumnToDesc(ms_source_description, casacore::MSSourceEnums::TRANSITION,
562 : 1);
563 4 : casacore::MSSource::addColumnToDesc(ms_source_description,
564 : casacore::MSSourceEnums::REST_FREQUENCY, 1);
565 4 : casacore::MSSource::addColumnToDesc(ms_source_description, casacore::MSSourceEnums::SYSVEL, 1);
566 4 : casacore::SetupNewTable ms_source_table(ms_->sourceTableName(), ms_source_description,
567 : casacore::Table::Scratch);
568 8 : ms_->rwKeywordSet().defineTable(
569 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SOURCE),
570 8 : casacore::Table(ms_source_table));
571 :
572 4 : casacore::TableDesc ms_spectral_window_description =
573 : casacore::MSSpectralWindow::requiredTableDesc();
574 4 : casacore::SetupNewTable ms_spectral_window_table(ms_->spectralWindowTableName(),
575 : ms_spectral_window_description, casacore::Table::Scratch);
576 8 : ms_->rwKeywordSet().defineTable(
577 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SPECTRAL_WINDOW),
578 8 : casacore::Table(ms_spectral_window_table));
579 :
580 4 : casacore::TableDesc ms_state_description = casacore::MSState::requiredTableDesc();
581 4 : casacore::SetupNewTable ms_state_table(ms_->stateTableName(), ms_state_description,
582 : casacore::Table::Scratch);
583 8 : ms_->rwKeywordSet().defineTable(
584 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::STATE),
585 8 : casacore::Table(ms_state_table));
586 :
587 4 : casacore::TableDesc ms_syscal_description = casacore::MSSysCal::requiredTableDesc();
588 4 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TCAL_SPECTRUM,
589 : 2);
590 4 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TCAL, 1);
591 4 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TSYS_SPECTRUM,
592 : 2);
593 4 : casacore::MSSysCal::addColumnToDesc(ms_syscal_description, casacore::MSSysCalEnums::TSYS, 1);
594 4 : casacore::SetupNewTable ms_syscal_table(ms_->sysCalTableName(), ms_syscal_description,
595 : casacore::Table::Scratch);
596 8 : ms_->rwKeywordSet().defineTable(
597 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::SYSCAL),
598 8 : casacore::Table(ms_syscal_table));
599 :
600 4 : casacore::TableDesc ms_weather_description = casacore::MSWeather::requiredTableDesc();
601 4 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
602 : casacore::MSWeatherEnums::TEMPERATURE);
603 4 : casacore::MSWeather::addColumnToDesc(ms_weather_description, casacore::MSWeatherEnums::PRESSURE);
604 4 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
605 : casacore::MSWeatherEnums::REL_HUMIDITY);
606 4 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
607 : casacore::MSWeatherEnums::WIND_SPEED);
608 4 : casacore::MSWeather::addColumnToDesc(ms_weather_description,
609 : casacore::MSWeatherEnums::WIND_DIRECTION);
610 4 : casacore::SetupNewTable ms_weather_table(ms_->weatherTableName(),
611 : ms_weather_description, casacore::Table::Scratch);
612 8 : ms_->rwKeywordSet().defineTable(
613 : casacore::MeasurementSet::keywordName(casacore::MeasurementSet::WEATHER),
614 8 : casacore::Table(ms_weather_table));
615 :
616 4 : ms_->initRefs();
617 :
618 : // Set up MSMainColumns
619 4 : ms_columns_.reset(new casacore::MSMainColumns(*ms_));
620 :
621 : // Set up MSDataDescColumns
622 4 : data_description_columns_.reset(
623 4 : new casacore::MSDataDescColumns(ms_->dataDescription()));
624 :
625 : // Set up MSFeedColumns
626 4 : feed_columns_.reset(new casacore::MSFeedColumns(ms_->feed()));
627 :
628 : // Set up MSPointingColumns
629 4 : pointing_columns_.reset(new casacore::MSPointingColumns(ms_->pointing()));
630 :
631 : // Set up MSPolarizationColumns
632 4 : polarization_columns_.reset(new casacore::MSPolarizationColumns(ms_->polarization()));
633 :
634 : // Set up MSSysCalColumns
635 4 : syscal_columns_.reset(new casacore::MSSysCalColumns(ms_->sysCal()));
636 :
637 : // Set up MSStateColumns
638 4 : state_columns_.reset(new casacore::MSStateColumns(ms_->state()));
639 :
640 : // Set up MSWeatherColumns
641 4 : weather_columns_.reset(new casacore::MSWeatherColumns(ms_->weather()));
642 :
643 : // std::cout << "End " << __PRETTY_FUNCTION__ << std::endl;
644 4 : }
645 :
646 : template<class T>
647 4 : void SingleDishMSFiller<T>::fillPreProcessTables() {
648 : POST_START;
649 :
650 : // fill OBSERVATION table
651 4 : fillObservation();
652 :
653 : // fill ANTENNA table
654 4 : fillAntenna();
655 :
656 : // fill PROCESSOR table
657 4 : fillProcessor();
658 :
659 : // fill SOURCE table
660 4 : fillSource();
661 :
662 : // fill FIELD table
663 4 : fillField();
664 :
665 : // fill SPECTRAL_WINDOW table
666 4 : fillSpectralWindow();
667 :
668 : // Generate optional tables if necessary
669 4 : OptionalTables::Generate(*ms_, *reader_);
670 :
671 : POST_END;
672 4 : }
673 :
674 : template<class T>
675 4 : void SingleDishMSFiller<T>::fillPostProcessTables() {
676 : POST_START;
677 :
678 : // fill HISTORY table
679 4 : fillHistory();
680 :
681 : // flush POINTING entry
682 4 : sortPointing();
683 :
684 : POST_END;
685 4 : }
686 :
687 : template<class T>
688 4 : void SingleDishMSFiller<T>::fillMain() {
689 : POST_START;
690 :
691 4 : size_t nrow = reader_->getNumberOfRows();
692 4 : DataAccumulator accumulator;
693 4 : DataRecord record;
694 : // std::cout << "nrow = " << nrow << std::endl;
695 796 : for (size_t irow = 0; irow < nrow; ++irow) {
696 792 : 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 792 : if (status) {
702 792 : casacore::Bool is_ready = accumulator.queryForGet(record.time);
703 792 : if (is_ready) {
704 300 : flush(accumulator);
705 : }
706 792 : casacore::Bool astatus = accumulator.accumulate(record);
707 : (void) astatus;
708 : // std::cout << "astatus = " << astatus << std::endl;
709 : }
710 : }
711 :
712 4 : flush(accumulator);
713 :
714 : POST_END;
715 4 : }
716 :
717 : template<class T>
718 4 : void SingleDishMSFiller<T>::fillAntenna() {
719 : POST_START;
720 :
721 4 : auto mytable = ms_->antenna();
722 4 : AntennaRecord record;
723 4 : ::fillTable(mytable, record,
724 8 : [&](AntennaRecord &record) {return reader_->getAntennaRow(record);});
725 :
726 : // initialize POINTING table related stuff
727 4 : casacore::uInt nrow = mytable.nrow();
728 4 : num_pointing_time_.resize(nrow);
729 8 : for (casacore::uInt i = 0; i < nrow; ++i) {
730 4 : pointing_time_[i] = casacore::Vector<casacore::Double>(ARRAY_BLOCK_SIZE, -1.0);
731 4 : pointing_time_min_[i] = -1.0;
732 4 : pointing_time_max_[i] = 1.0e30;
733 4 : num_pointing_time_[i] = 0;
734 : }
735 :
736 : POST_END;
737 4 : }
738 :
739 : template<class T>
740 4 : void SingleDishMSFiller<T>::fillObservation() {
741 : POST_START;
742 :
743 4 : ObservationRecord record;
744 4 : ::fillTable(ms_->observation(), record,
745 8 : [&](ObservationRecord &record) {return reader_->getObservationRow(record);});
746 :
747 : POST_END;
748 4 : }
749 :
750 : template<class T>
751 4 : void SingleDishMSFiller<T>::fillProcessor() {
752 : POST_START;
753 :
754 4 : ProcessorRecord record;
755 4 : ::fillTable(ms_->processor(), record,
756 8 : [&](ProcessorRecord &record) {return reader_->getProcessorRow(record);});
757 :
758 : POST_END;
759 4 : }
760 :
761 : template<class T>
762 4 : void SingleDishMSFiller<T>::fillSource() {
763 : POST_START;
764 :
765 4 : SourceRecord record;
766 4 : ::fillTable(ms_->source(), record,
767 24 : [&](SourceRecord &record) {return reader_->getSourceRow(record);});
768 :
769 : POST_END;
770 4 : }
771 :
772 : template<class T>
773 4 : void SingleDishMSFiller<T>::fillField() {
774 : POST_START;
775 :
776 4 : auto mytable = ms_->field();
777 4 : FieldRecord record;
778 4 : record.table = mytable;
779 :
780 : // make (name,id) map for SOURCE table
781 4 : ::makeSourceMap(ms_->source(), record.source_map);
782 :
783 4 : ::fillTable(mytable, record,
784 8 : [&](FieldRecord &record) {return reader_->getFieldRow(record);});
785 :
786 : POST_END;
787 4 : }
788 :
789 : template<class T>
790 4 : 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 4 : }
807 :
808 : template<class T>
809 4 : void SingleDishMSFiller<T>::fillSpectralWindow() {
810 : POST_START;
811 :
812 4 : auto mytable = ms_->spectralWindow();
813 4 : SpectralWindowRecord record;
814 :
815 4 : ::fillTable(mytable, record,
816 24 : [&](SpectralWindowRecord &record) {return reader_->getSpectralWindowRow(record);});
817 :
818 : POST_END;
819 4 : }
820 :
821 : template<class T>
822 472 : casacore::Int SingleDishMSFiller<T>::updatePolarization(casacore::Vector<casacore::Int> const &corr_type,
823 : casacore::Int const &num_pol) {
824 472 : casacore::uInt num_corr = corr_type.size();
825 472 : if (num_pol < 1 || num_corr != (uInt) num_pol) {
826 0 : throw casacore::AipsError("Internal inconsistency in number of correlations");
827 : }
828 472 : casacore::MSPolarization &mytable = ms_->polarization();
829 : //casacore::MSPolarizationColumns mycolumns(mytable);
830 472 : casacore::Matrix<casacore::Int> const corr_product(2, num_pol, 0);
831 1404 : auto comparer = [&](casacore::MSPolarizationColumns const &columns, casacore::uInt i) {
832 932 : casacore::Bool match = allEQ(columns.corrType()(i), corr_type);
833 932 : return match;
834 : };
835 512 : auto updater = [&](casacore::MSPolarizationColumns &columns, casacore::uInt i) {
836 8 : columns.numCorr().put(i, num_pol);
837 8 : columns.corrType().put(i, corr_type);
838 8 : columns.corrProduct().put(i, corr_product);
839 : };
840 472 : casacore::Int polarization_id = ::updateTable(mytable, *(polarization_columns_.get()),
841 : comparer, updater);
842 472 : return polarization_id;
843 472 : }
844 :
845 : template<class T>
846 472 : Int SingleDishMSFiller<T>::updateDataDescription(casacore::Int const &polarization_id,
847 : casacore::Int const &spw_id) {
848 472 : if (polarization_id < 0 || spw_id < 0) {
849 0 : throw casacore::AipsError("Invalid ids for DATA_DESCRIPTION");
850 : }
851 472 : casacore::MSDataDescription &mytable = ms_->dataDescription();
852 : //casacore::MSDataDescColumns mycolumns(mytable);
853 5260 : auto comparer = [&](casacore::MSDataDescColumns const &columns, casacore::uInt i) {
854 1852 : casacore::Bool match = (columns.polarizationId()(i) == polarization_id)
855 1852 : && (columns.spectralWindowId()(i) == spw_id);
856 1852 : return match;
857 : };
858 552 : auto updater = [&](casacore::MSDataDescColumns &columns, casacore::uInt i) {
859 20 : columns.polarizationId().put(i, polarization_id);
860 20 : columns.spectralWindowId().put(i, spw_id);
861 : };
862 472 : casacore::Int data_desc_id = ::updateTable(mytable, *(data_description_columns_.get()),
863 : comparer, updater);
864 :
865 472 : return data_desc_id;
866 : }
867 :
868 : template<class T>
869 472 : Int SingleDishMSFiller<T>::updateState(casacore::Int const &subscan,
870 : casacore::String const &obs_mode) {
871 472 : casacore::MSState &mytable = ms_->state();
872 472 : static casacore::Regex const regex("^OBSERVE_TARGET#ON_SOURCE");
873 : //static std::vector<casacore::Int> subscan_list;
874 472 : auto comparer =
875 1580 : [&](casacore::MSStateColumns &columns, casacore::uInt i) {
876 1124 : casacore::Bool match = (subscan == subscan_list_[i]) && (obs_mode == columns.obsMode()(i));
877 1124 : return match;
878 : };
879 488 : auto updater = [&](casacore::MSStateColumns &columns, casacore::uInt i) {
880 16 : columns.subScan().put(i, subscan);
881 16 : columns.obsMode().put(i, obs_mode);
882 16 : casacore::Bool is_signal = obs_mode.matches(regex);
883 16 : columns.sig().put(i, is_signal);
884 16 : columns.ref().put(i, !is_signal);
885 :
886 16 : subscan_list_.push_back(subscan);
887 : };
888 472 : casacore::Int state_id = ::updateTable(mytable, *(state_columns_.get()), comparer,
889 : updater);
890 472 : return state_id;
891 : }
892 :
893 : template<class T>
894 472 : Int SingleDishMSFiller<T>::updateFeed(casacore::Int const &feed_id, casacore::Int const &spw_id,
895 : casacore::String const &pol_type) {
896 472 : constexpr casacore::Int num_receptors = 2;
897 472 : static casacore::String const linear_type_arr[2] = { "X", "Y" };
898 472 : static casacore::Vector<casacore::String> linear_type(linear_type_arr, 2, SHARE);
899 472 : static casacore::String const circular_type_arr[2] = { "R", "L" };
900 472 : static casacore::Vector<casacore::String> circular_type(circular_type_arr, 2, SHARE);
901 472 : static casacore::Matrix<casacore::Complex> const pol_response(num_receptors, num_receptors,
902 1 : casacore::Complex(0));
903 472 : casacore::Vector<casacore::String> *polarization_type = nullptr;
904 472 : if (pol_type == "linear") {
905 472 : 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 1888 : casacore::String polarization_type_arr[2] = { "X", "Y" };
914 472 : casacore::Vector<casacore::String> polarization_type_storage(polarization_type_arr, 2, SHARE);
915 472 : casacore::Matrix<casacore::Double> const beam_offset(2, num_receptors, 0.0);
916 472 : casacore::Vector<casacore::Double> const receptor_angle(num_receptors, 0.0);
917 472 : casacore::Vector<casacore::Double> const position(3, 0.0);
918 7880 : auto comparer = [&](casacore::MSFeedColumns &columns, casacore::uInt i) {
919 1852 : casacore::Vector<casacore::String> *current_polarization_type = polarization_type_pool_[i];
920 1852 : casacore::Bool match = allEQ(*polarization_type, *current_polarization_type) &&
921 3704 : (feed_id == columns.feedId()(i)) &&
922 1852 : (spw_id == columns.spectralWindowId()(i));
923 1852 : return match;
924 : };
925 492 : 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 20 : columns.feedId().put(i, feed_id);
932 20 : columns.beamId().put(i, 0);
933 20 : columns.spectralWindowId().put(i, spw_id);
934 20 : columns.polarizationType().put(i, *polarization_type);
935 20 : columns.beamOffset().put(i, beam_offset);
936 20 : columns.receptorAngle().put(i, receptor_angle);
937 20 : columns.position().put(i, position);
938 20 : columns.polResponse().put(i, pol_response);
939 :
940 20 : polarization_type_pool_.push_back(polarization_type);
941 : };
942 472 : casacore::Int feed_row = ::updateTable(ms_->feed(), *(feed_columns_.get()), comparer,
943 : updater);
944 :
945 : // reference feed
946 472 : if (reference_feed_ < 0) {
947 4 : reference_feed_ = feed_id;
948 : }
949 :
950 472 : return feed_row;
951 1888 : }
952 :
953 : template<class T>
954 472 : 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 472 : if (reference_feed_ != feed_id) {
960 0 : return -1;
961 : }
962 :
963 472 : auto &mytable = ms_->pointing();
964 472 : casacore::uInt nrow = mytable.nrow();
965 :
966 472 : casacore::uInt *n = &num_pointing_time_[antenna_id];
967 472 : casacore::Double *time_max = &pointing_time_max_[antenna_id];
968 472 : casacore::Double *time_min = &pointing_time_min_[antenna_id];
969 472 : casacore::Vector<casacore::Double> *time_list = &pointing_time_[antenna_id];
970 :
971 472 : auto addPointingRow =
972 5472 : [&]() {
973 304 : mytable.addRow(1, true);
974 304 : pointing_columns_->time().put(nrow, time);
975 304 : pointing_columns_->interval().put(nrow, interval);
976 304 : pointing_columns_->antennaId().put(nrow, antenna_id);
977 304 : if (direction.ncolumn() == 1 || allEQ(direction.column(1), 0.0)) {
978 304 : pointing_columns_->numPoly().put(nrow, 0);
979 304 : 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 304 : casacore::uInt nelem = time_list->nelements();
987 304 : if (nelem <= *n) {
988 0 : time_list->resize(nelem + ARRAY_BLOCK_SIZE, true);
989 : }
990 304 : (*time_list)[*n] = time;
991 : // increment number of pointing entry
992 304 : *n += 1;
993 : };
994 :
995 472 : if (*n == 0) {
996 4 : addPointingRow();
997 4 : *time_min = time;
998 4 : *time_max = time;
999 468 : } else if (time < *time_min) {
1000 0 : addPointingRow();
1001 0 : *time_min = time;
1002 468 : } else if (*time_max < time) {
1003 300 : addPointingRow();
1004 300 : *time_max = time;
1005 168 : } else if (allNE(*time_list, time)) {
1006 0 : addPointingRow();
1007 : }
1008 :
1009 : POST_END;
1010 :
1011 472 : return -1;
1012 : }
1013 :
1014 : template<class T>
1015 472 : 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 472 : record.clear();
1020 472 : record.antenna_id = antenna_id;
1021 472 : record.time = time;
1022 472 : record.interval = interval;
1023 472 : record.temperature = data_record.temperature;
1024 472 : record.pressure = data_record.pressure;
1025 472 : record.rel_humidity = data_record.rel_humidity;
1026 472 : record.wind_speed = data_record.wind_speed;
1027 472 : record.wind_direction = data_record.wind_direction;
1028 472 : auto &mytable = ms_->weather();
1029 472 : auto pos = std::find(weather_list_.begin(), weather_list_.end(), record);
1030 472 : if (pos == weather_list_.end()) {
1031 20 : weather_list_.push_back(record);
1032 20 : casacore::uInt irow = mytable.nrow();
1033 20 : mytable.addRow(1, true);
1034 20 : record.fill(irow, *(weather_columns_.get()));
1035 : } else {
1036 452 : auto irow = std::distance(weather_list_.begin(), pos);
1037 452 : updateWeather(*(weather_columns_.get()), irow, record);
1038 : }
1039 472 : }
1040 :
1041 : template<class T>
1042 452 : void SingleDishMSFiller<T>::updateWeather(casacore::MSWeatherColumns &columns, casacore::uInt irow,
1043 : WeatherRecord const &record) {
1044 452 : ::updateSubtable(columns, irow, record);
1045 452 : }
1046 :
1047 : template<class T>
1048 472 : 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 472 : SysCalRecord record;
1054 472 : record.clear();
1055 472 : record.antenna_id = antenna_id;
1056 472 : record.feed_id = feed_id;
1057 472 : record.spw_id = spw_id;
1058 472 : record.time = time;
1059 472 : record.interval = interval;
1060 :
1061 : //casacore::Bool tcal_empty = false;
1062 472 : casacore::Bool tsys_empty = false;
1063 :
1064 944 : if (data_record.tcal.empty() || allEQ(data_record.tcal, 1.0f)
1065 944 : || 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 80 : if (data_record.float_data.shape() == data_record.tcal.shape()) {
1070 80 : 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 944 : if (data_record.tsys.empty() || allEQ(data_record.tsys, 1.0f)
1079 944 : || allEQ(data_record.tsys, 0.0f)) {
1080 392 : tsys_empty = true;
1081 : } else {
1082 80 : if (data_record.float_data.shape() == data_record.tsys.shape()) {
1083 80 : 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 472 : if (tsys_empty) {
1095 392 : return;
1096 : }
1097 :
1098 80 : auto &mytable = ms_->sysCal();
1099 80 : auto pos = std::find(syscal_list_.begin(), syscal_list_.end(), record);
1100 80 : if (pos == syscal_list_.end()) {
1101 4 : casacore::uInt irow = mytable.nrow();
1102 4 : mytable.addRow(1, true);
1103 4 : record.fill(irow, *(syscal_columns_.get()));
1104 4 : syscal_list_.push_back(SysCalTableRecord(ms_.get(), irow, record));
1105 : } else {
1106 76 : auto irow = std::distance(syscal_list_.begin(), pos);
1107 76 : updateSysCal(*(syscal_columns_.get()), irow, record);
1108 : }
1109 :
1110 : POST_END;
1111 472 : }
1112 :
1113 : template<class T>
1114 76 : void SingleDishMSFiller<T>::updateSysCal(casacore::MSSysCalColumns &columns, casacore::uInt irow,
1115 : SysCalRecord const &record) {
1116 76 : ::updateSubtable(columns, irow, record);
1117 76 : }
1118 :
1119 : template<class T>
1120 472 : 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 472 : static casacore::Vector<casacore::Double> const uvw(3, 0.0);
1127 472 : static casacore::Array<casacore::Bool> const flagCategory(casacore::IPosition(3, 0, 0, 0));
1128 :
1129 : // target row id
1130 472 : casacore::uInt irow = ms_->nrow();
1131 :
1132 : // add new row
1133 : //ms_->addRow(1, true);
1134 472 : ms_->addRow(1, false);
1135 :
1136 472 : ms_columns_->uvw().put(irow, uvw);
1137 472 : ms_columns_->flagCategory().put(irow, flagCategory);
1138 472 : ms_columns_->antenna1().put(irow, antenna_id);
1139 472 : ms_columns_->antenna2().put(irow, antenna_id);
1140 472 : ms_columns_->fieldId().put(irow, field_id);
1141 472 : ms_columns_->feed1().put(irow, feedId);
1142 472 : ms_columns_->feed2().put(irow, feedId);
1143 472 : ms_columns_->dataDescId().put(irow, dataDescriptionId);
1144 472 : ms_columns_->stateId().put(irow, stateId);
1145 472 : ms_columns_->scanNumber().put(irow, scan_number);
1146 472 : ms_columns_->time().put(irow, time);
1147 472 : ms_columns_->timeCentroid().put(irow, time);
1148 472 : casacore::Double const &interval = dataRecord.interval;
1149 472 : ms_columns_->interval().put(irow, interval);
1150 472 : ms_columns_->exposure().put(irow, interval);
1151 :
1152 472 : if (is_float_) {
1153 472 : casacore::Matrix<casacore::Float> floatData;
1154 472 : if (dataRecord.isFloat()) {
1155 472 : floatData.reference(dataRecord.float_data);
1156 : } else {
1157 0 : floatData.assign(real(dataRecord.complex_data));
1158 : }
1159 472 : ms_columns_->floatData().put(irow, floatData);
1160 472 : } 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 472 : ms_columns_->flag().put(irow, dataRecord.flag);
1173 472 : ms_columns_->flagRow().put(irow, dataRecord.flag_row);
1174 472 : ms_columns_->sigma().put(irow, dataRecord.sigma);
1175 472 : ms_columns_->weight().put(irow, dataRecord.weight);
1176 :
1177 : POST_END;
1178 472 : }
1179 :
1180 : template<class T>
1181 304 : void SingleDishMSFiller<T>::flush(DataAccumulator &accumulator) {
1182 : POST_START;
1183 :
1184 304 : size_t nchunk = accumulator.getNumberOfChunks();
1185 : // std::cout << "nchunk = " << nchunk << std::endl;
1186 :
1187 304 : if (nchunk == 0) {
1188 0 : return;
1189 : }
1190 :
1191 2536 : for (size_t ichunk = 0; ichunk < nchunk; ++ichunk) {
1192 2232 : casacore::Bool status = accumulator.get(ichunk, record_);
1193 : // std::cout << "accumulator status = " << std::endl;
1194 2232 : if (status) {
1195 472 : casacore::Double time = record_.time;
1196 472 : casacore::Int antenna_id = record_.antenna_id;
1197 472 : casacore::Int spw_id = record_.spw_id;
1198 472 : casacore::Int feed_id = record_.feed_id;
1199 472 : casacore::Int field_id = record_.field_id;
1200 472 : casacore::Int scan = record_.scan;
1201 472 : casacore::Int subscan = record_.subscan;
1202 472 : casacore::String pol_type = record_.pol_type;
1203 472 : casacore::String obs_mode = record_.intent;
1204 472 : casacore::Int num_pol = record_.num_pol;
1205 472 : casacore::Vector<casacore::Int> &corr_type = record_.corr_type;
1206 472 : casacore::Int polarization_id = updatePolarization(corr_type, num_pol);
1207 472 : updateFeed(feed_id, spw_id, pol_type);
1208 472 : casacore::Int data_desc_id = updateDataDescription(polarization_id, spw_id);
1209 472 : casacore::Int state_id = updateState(subscan, obs_mode);
1210 472 : casacore::Matrix<casacore::Double> &direction = record_.direction;
1211 472 : casacore::Double interval = record_.interval;
1212 :
1213 : // updatePointing must be called after updateFeed
1214 472 : updatePointing(antenna_id, feed_id, time, interval, direction);
1215 :
1216 472 : updateSysCal(antenna_id, feed_id, spw_id, time, interval, record_);
1217 :
1218 472 : updateWeather(antenna_id, time, interval, record_);
1219 :
1220 472 : updateMain(antenna_id, field_id, feed_id, data_desc_id, state_id, scan,
1221 472 : time, record_);
1222 472 : }
1223 : }
1224 : // std::cout << "clear accumulator" << std::endl;
1225 304 : accumulator.clear();
1226 :
1227 : POST_END;
1228 : }
1229 :
1230 : template<class T>
1231 4 : void SingleDishMSFiller<T>::sortPointing() {
1232 : POST_START;
1233 :
1234 : // deallocate POINTING table related stuff
1235 4 : pointing_time_.clear();
1236 4 : pointing_time_min_.clear();
1237 4 : pointing_time_max_.clear();
1238 4 : num_pointing_time_.resize();
1239 :
1240 4 : auto mytable = ms_->pointing();
1241 4 : casacore::MSPointingColumns mycolumns(mytable);
1242 4 : casacore::uInt nrow = mytable.nrow();
1243 4 : casacore::Vector<casacore::Int> antenna_id_list = mycolumns.antennaId().getColumn();
1244 4 : casacore::Vector<casacore::Double> time_list = mycolumns.time().getColumn();
1245 4 : casacore::Sort sorter;
1246 4 : sorter.sortKey(antenna_id_list.data(), TpInt);
1247 4 : sorter.sortKey(time_list.data(), TpDouble);
1248 4 : casacore::Vector<casacore::uInt> index_vector;
1249 4 : sorter.sort(index_vector, nrow);
1250 :
1251 4 : size_t storage_size = nrow * 2 * sizeof(casacore::Double);
1252 4 : std::unique_ptr<void, sdfiller::Deleter> storage(malloc(storage_size));
1253 :
1254 : // sort TIME
1255 : {
1256 4 : casacore::Vector<casacore::Double> sorted(casacore::IPosition(1, nrow),
1257 4 : reinterpret_cast<casacore::Double *>(storage.get()), SHARE);
1258 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1259 304 : sorted[i] = time_list[index_vector[i]];
1260 : }
1261 4 : mycolumns.time().putColumn(sorted);
1262 4 : }
1263 : // sort ANTENNA_ID
1264 : {
1265 4 : casacore::Vector<casacore::Int> sorted(casacore::IPosition(1, nrow),
1266 4 : reinterpret_cast<casacore::Int *>(storage.get()), SHARE);
1267 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1268 304 : sorted[i] = antenna_id_list[index_vector[i]];
1269 : }
1270 4 : mycolumns.antennaId().putColumn(sorted);
1271 4 : }
1272 :
1273 : // sort NUM_POLY
1274 : {
1275 4 : casacore::Vector<casacore::Int> num_poly_list(antenna_id_list);
1276 4 : mycolumns.numPoly().getColumn(num_poly_list);
1277 4 : casacore::Vector<casacore::Int> sorted(casacore::IPosition(1, nrow),
1278 4 : reinterpret_cast<casacore::Int *>(storage.get()), SHARE);
1279 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1280 304 : sorted[i] = antenna_id_list[index_vector[i]];
1281 : }
1282 4 : mycolumns.numPoly().putColumn(sorted);
1283 4 : }
1284 :
1285 : // sort DIRECTION
1286 : {
1287 4 : std::map<casacore::uInt, casacore::Matrix<casacore::Double> > direction;
1288 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1289 304 : direction[i] = mycolumns.direction()(i);
1290 : }
1291 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1292 304 : mycolumns.direction().put(i, direction[index_vector[i]]);
1293 : }
1294 4 : }
1295 :
1296 : // sort INTERVAL
1297 : {
1298 4 : casacore::Vector<casacore::Double> interval_list = mycolumns.interval().getColumn();
1299 308 : for (casacore::uInt i = 0; i < nrow; ++i) {
1300 304 : mycolumns.interval().put(i, interval_list[index_vector[i]]);
1301 : }
1302 4 : }
1303 : POST_END;
1304 4 : }
1305 :
1306 : } //# NAMESPACE CASA - END
1307 :
1308 : #endif /* SINGLEDISH_FILLER_SINGLEDISHMSFILLER_TCC_ */
|
