source: trunk/src/MSFiller.cpp@ 2884

Last change on this file since 2884 was 2869, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5841

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave, test_importasdm_sd, all sd regressions

Put in Release Notes: Yes

Module(s): sd

Description: Describe your changes here...

Filler/writer are changed so that SCANNO is consistent with
SCAN_NUMBER in MS and scanNumber in ASDM.

It affects several regressions and unit tests so that tests
are necessary to be updated. This change should be verified
by the updated unit tests/regression tests.


File size: 68.6 KB
Line 
1//
2// C++ Interface: MSFiller
3//
4// Description:
5//
6// This class is specific filler for MS format
7// New version that is implemented using TableVisitor instead of TableIterator
8//
9// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
10//
11// Copyright: See COPYING file that comes with this distribution
12//
13//
14
15#include <assert.h>
16#include <iostream>
17#include <map>
18#include <set>
19
20#include <tables/Tables/ExprNode.h>
21#include <tables/Tables/TableIter.h>
22#include <tables/Tables/TableColumn.h>
23#include <tables/Tables/ScalarColumn.h>
24#include <tables/Tables/ArrayColumn.h>
25#include <tables/Tables/TableParse.h>
26#include <tables/Tables/TableRow.h>
27
28#include <casa/Containers/RecordField.h>
29#include <casa/Logging/LogIO.h>
30#include <casa/Arrays/Slicer.h>
31#include <casa/Quanta/MVTime.h>
32#include <casa/OS/Path.h>
33
34#include <measures/Measures/Stokes.h>
35#include <measures/Measures/MEpoch.h>
36#include <measures/Measures/MCEpoch.h>
37#include <measures/Measures/MFrequency.h>
38#include <measures/Measures/MCFrequency.h>
39#include <measures/Measures/MPosition.h>
40#include <measures/Measures/MCPosition.h>
41#include <measures/Measures/MDirection.h>
42#include <measures/Measures/MCDirection.h>
43#include <measures/Measures/MeasConvert.h>
44#include <measures/TableMeasures/ScalarMeasColumn.h>
45#include <measures/TableMeasures/ArrayMeasColumn.h>
46#include <measures/TableMeasures/ScalarQuantColumn.h>
47#include <measures/TableMeasures/ArrayQuantColumn.h>
48
49#include <ms/MeasurementSets/MSAntennaIndex.h>
50
51#include <atnf/PKSIO/SrcType.h>
52
53#include "MSFiller.h"
54#include "STHeader.h"
55
56#include "MathUtils.h"
57
58using namespace casa ;
59using namespace std ;
60
61namespace asap {
62
63class BaseMSFillerVisitor: public TableVisitor {
64 uInt lastRecordNo ;
65 Int lastObservationId ;
66 Int lastFeedId ;
67 Int lastFieldId ;
68 Int lastDataDescId ;
69 Int lastScanNo ;
70 Int lastStateId ;
71 Double lastTime ;
72protected:
73 const Table &table;
74 uInt count;
75public:
76 BaseMSFillerVisitor(const Table &table)
77 : table(table)
78 {
79 count = 0;
80 }
81
82 virtual void enterObservationId(const uInt /*recordNo*/, Int /*columnValue*/) { }
83 virtual void leaveObservationId(const uInt /*recordNo*/, Int /*columnValue*/) { }
84 virtual void enterFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
85 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
86 virtual void enterFieldId(const uInt /*recordNo*/, Int /*columnValue*/) { }
87 virtual void leaveFieldId(const uInt /*recordNo*/, Int /*columnValue*/) { }
88 virtual void enterDataDescId(const uInt /*recordNo*/, Int /*columnValue*/) { }
89 virtual void leaveDataDescId(const uInt /*recordNo*/, Int /*columnValue*/) { }
90 virtual void enterScanNo(const uInt /*recordNo*/, Int /*columnValue*/) { }
91 virtual void leaveScanNo(const uInt /*recordNo*/, Int /*columnValue*/) { }
92 virtual void enterStateId(const uInt /*recordNo*/, Int /*columnValue*/) { }
93 virtual void leaveStateId(const uInt /*recordNo*/, Int /*columnValue*/) { }
94 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
95 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
96
97 virtual Bool visitRecord(const uInt /*recordNo*/,
98 const Int /*ObservationId*/,
99 const Int /*feedId*/,
100 const Int /*fieldId*/,
101 const Int /*dataDescId*/,
102 const Int /*scanNo*/,
103 const Int /*stateId*/,
104 const Double /*time*/) { return True ; }
105
106 virtual Bool visit(Bool isFirst, const uInt recordNo,
107 const uInt nCols, void const *const colValues[]) {
108 Int observationId, feedId, fieldId, dataDescId, scanNo, stateId;
109 Double time;
110 { // prologue
111 uInt i = 0;
112 {
113 const Int *col = (const Int *)colValues[i++];
114 observationId = col[recordNo];
115 }
116 {
117 const Int *col = (const Int *)colValues[i++];
118 feedId = col[recordNo];
119 }
120 {
121 const Int *col = (const Int *)colValues[i++];
122 fieldId = col[recordNo];
123 }
124 {
125 const Int *col = (const Int *)colValues[i++];
126 dataDescId = col[recordNo];
127 }
128 {
129 const Int *col = (const Int *)colValues[i++];
130 scanNo = col[recordNo];
131 }
132 {
133 const Int *col = (const Int *)colValues[i++];
134 stateId = col[recordNo];
135 }
136 {
137 const Double *col = (const Double *)colValues[i++];
138 time = col[recordNo];
139 }
140 assert(nCols == i);
141 }
142
143 if (isFirst) {
144 enterObservationId(recordNo, observationId);
145 enterFeedId(recordNo, feedId);
146 enterFieldId(recordNo, fieldId);
147 enterDataDescId(recordNo, dataDescId);
148 enterScanNo(recordNo, scanNo);
149 enterStateId(recordNo, stateId);
150 enterTime(recordNo, time);
151 } else {
152 if (lastObservationId != observationId) {
153 leaveTime(lastRecordNo, lastTime);
154 leaveStateId(lastRecordNo, lastStateId);
155 leaveScanNo(lastRecordNo, lastScanNo);
156 leaveDataDescId(lastRecordNo, lastDataDescId);
157 leaveFieldId(lastRecordNo, lastFieldId);
158 leaveFeedId(lastRecordNo, lastFeedId);
159 leaveObservationId(lastRecordNo, lastObservationId);
160
161 enterObservationId(recordNo, observationId);
162 enterFeedId(recordNo, feedId);
163 enterFieldId(recordNo, fieldId);
164 enterDataDescId(recordNo, dataDescId);
165 enterScanNo(recordNo, scanNo);
166 enterStateId(recordNo, stateId);
167 enterTime(recordNo, time);
168 } else if (lastFeedId != feedId) {
169 leaveTime(lastRecordNo, lastTime);
170 leaveStateId(lastRecordNo, lastStateId);
171 leaveScanNo(lastRecordNo, lastScanNo);
172 leaveDataDescId(lastRecordNo, lastDataDescId);
173 leaveFieldId(lastRecordNo, lastFieldId);
174 leaveFeedId(lastRecordNo, lastFeedId);
175
176 enterFeedId(recordNo, feedId);
177 enterFieldId(recordNo, fieldId);
178 enterDataDescId(recordNo, dataDescId);
179 enterScanNo(recordNo, scanNo);
180 enterStateId(recordNo, stateId);
181 enterTime(recordNo, time);
182 } else if (lastFieldId != fieldId) {
183 leaveTime(lastRecordNo, lastTime);
184 leaveStateId(lastRecordNo, lastStateId);
185 leaveScanNo(lastRecordNo, lastScanNo);
186 leaveDataDescId(lastRecordNo, lastDataDescId);
187 leaveFieldId(lastRecordNo, lastFieldId);
188
189 enterFieldId(recordNo, fieldId);
190 enterDataDescId(recordNo, dataDescId);
191 enterScanNo(recordNo, scanNo);
192 enterStateId(recordNo, stateId);
193 enterTime(recordNo, time);
194 } else if (lastDataDescId != dataDescId) {
195 leaveTime(lastRecordNo, lastTime);
196 leaveStateId(lastRecordNo, lastStateId);
197 leaveScanNo(lastRecordNo, lastScanNo);
198 leaveDataDescId(lastRecordNo, lastDataDescId);
199
200 enterDataDescId(recordNo, dataDescId);
201 enterScanNo(recordNo, scanNo);
202 enterStateId(recordNo, stateId);
203 enterTime(recordNo, time);
204 } else if (lastScanNo != scanNo) {
205 leaveTime(lastRecordNo, lastTime);
206 leaveStateId(lastRecordNo, lastStateId);
207 leaveScanNo(lastRecordNo, lastScanNo);
208
209 enterScanNo(recordNo, scanNo);
210 enterStateId(recordNo, stateId);
211 enterTime(recordNo, time);
212 } else if (lastStateId != stateId) {
213 leaveTime(lastRecordNo, lastTime);
214 leaveStateId(lastRecordNo, lastStateId);
215
216 enterStateId(recordNo, stateId);
217 enterTime(recordNo, time);
218 } else if (lastTime != time) {
219 leaveTime(lastRecordNo, lastTime);
220 enterTime(recordNo, time);
221 }
222 }
223 count++;
224 Bool result = visitRecord(recordNo, observationId, feedId, fieldId, dataDescId,
225 scanNo, stateId, time);
226
227 { // epilogue
228 lastRecordNo = recordNo;
229
230 lastObservationId = observationId;
231 lastFeedId = feedId;
232 lastFieldId = fieldId;
233 lastDataDescId = dataDescId;
234 lastScanNo = scanNo;
235 lastStateId = stateId;
236 lastTime = time;
237 }
238 return result ;
239 }
240
241 virtual void finish() {
242 if (count > 0) {
243 leaveTime(lastRecordNo, lastTime);
244 leaveStateId(lastRecordNo, lastStateId);
245 leaveScanNo(lastRecordNo, lastScanNo);
246 leaveDataDescId(lastRecordNo, lastDataDescId);
247 leaveFieldId(lastRecordNo, lastFieldId);
248 leaveFeedId(lastRecordNo, lastFeedId);
249 leaveObservationId(lastRecordNo, lastObservationId);
250 }
251 }
252};
253
254class MSFillerVisitor: public BaseMSFillerVisitor, public MSFillerUtils {
255public:
256 MSFillerVisitor(const Table &from, Scantable &to)
257 : BaseMSFillerVisitor(from),
258 scantable(to)
259 {
260 antennaId = 0 ;
261 rowidx = 0 ;
262 tablerow = TableRow( scantable.table() ) ;
263 feedEntry = Vector<Int>( 64, -1 ) ;
264 nbeam = 0 ;
265 ifmap.clear() ;
266 const TableDesc &desc = table.tableDesc() ;
267 if ( desc.isColumn( "DATA" ) )
268 dataColumnName = "DATA" ;
269 else if ( desc.isColumn( "FLOAT_DATA" ) )
270 dataColumnName = "FLOAT_DATA" ;
271 getpt = False ;
272 isWeather_ = False ;
273 isSysCal = False ;
274 isTcal = False ;
275 cycleNo = 0 ;
276 numSysCalRow = 0 ;
277 header = scantable.getHeader() ;
278 fluxUnit( header.fluxunit ) ;
279
280 // MS subtables
281 const TableRecord &hdr = table.keywordSet();
282 obstab = hdr.asTable( "OBSERVATION" ) ;
283 spwtab = hdr.asTable( "SPECTRAL_WINDOW" ) ;
284 statetab = hdr.asTable( "STATE" ) ;
285 ddtab = hdr.asTable( "DATA_DESCRIPTION" ) ;
286 poltab = hdr.asTable( "POLARIZATION" ) ;
287 fieldtab = hdr.asTable( "FIELD" ) ;
288 anttab = hdr.asTable( "ANTENNA" ) ;
289 if ( hdr.isDefined( "SYSCAL" ) )
290 sctab = hdr.asTable( "SYSCAL" ) ;
291 if ( hdr.isDefined( "SOURCE" ) )
292 srctab = hdr.asTable( "SOURCE" ) ;
293
294 // attach to columns
295 // MS MAIN
296 intervalCol.attach( table, "INTERVAL" ) ;
297 flagRowCol.attach( table, "FLAG_ROW" ) ;
298 flagCol.attach( table, "FLAG" ) ;
299 if ( dataColumnName.compare( "DATA" ) == 0 )
300 dataCol.attach( table, dataColumnName ) ;
301 else
302 floatDataCol.attach( table, dataColumnName ) ;
303
304 // set dummy epoch
305 mf.set( currentTime ) ;
306
307 //
308 // add rows to scantable
309 //
310 // number of polarization is up to 4
311 uInt addrow = table.nrow() * maxNumPol() ;
312 scantable.table().addRow( addrow ) ;
313
314 // attach to columns
315 // Scantable MAIN
316 TableRecord &r = tablerow.record() ;
317 timeRF.attachToRecord( r, "TIME" ) ;
318 intervalRF.attachToRecord( r, "INTERVAL" ) ;
319 directionRF.attachToRecord( r, "DIRECTION" ) ;
320 azimuthRF.attachToRecord( r, "AZIMUTH" ) ;
321 elevationRF.attachToRecord( r, "ELEVATION" ) ;
322 scanRateRF.attachToRecord( r, "SCANRATE" ) ;
323 weatherIdRF.attachToRecord( r, "WEATHER_ID" ) ;
324 cycleNoRF.attachToRecord( r, "CYCLENO" ) ;
325 flagRowRF.attachToRecord( r, "FLAGROW" ) ;
326 polNoRF.attachToRecord( r, "POLNO" ) ;
327 tcalIdRF.attachToRecord( r, "TCAL_ID" ) ;
328 spectraRF.attachToRecord( r, "SPECTRA" ) ;
329 flagtraRF.attachToRecord( r, "FLAGTRA" ) ;
330 tsysRF.attachToRecord( r, "TSYS" ) ;
331 beamNoRF.attachToRecord( r, "BEAMNO" ) ;
332 ifNoRF.attachToRecord( r, "IFNO" ) ;
333 freqIdRF.attachToRecord( r, "FREQ_ID" ) ;
334 moleculeIdRF.attachToRecord( r, "MOLECULE_ID" ) ;
335 sourceNameRF.attachToRecord( r, "SRCNAME" ) ;
336 sourceProperMotionRF.attachToRecord( r, "SRCPROPERMOTION" ) ;
337 sourceDirectionRF.attachToRecord( r, "SRCDIRECTION" ) ;
338 sourceVelocityRF.attachToRecord( r, "SRCVELOCITY" ) ;
339 focusIdRF.attachToRecord( r, "FOCUS_ID" ) ;
340 fieldNameRF.attachToRecord( r, "FIELDNAME" ) ;
341 sourceTypeRF.attachToRecord( r, "SRCTYPE" ) ;
342 scanNoRF.attachToRecord( r, "SCANNO" ) ;
343
344 // put values
345 RecordFieldPtr<Int> refBeamNoRF( r, "REFBEAMNO" ) ;
346 *refBeamNoRF = -1 ;
347 RecordFieldPtr<Int> fitIdRF( r, "FIT_ID" ) ;
348 *fitIdRF = -1 ;
349 RecordFieldPtr<Float> opacityRF( r, "OPACITY" ) ;
350 *opacityRF = 0.0 ;
351 }
352
353 virtual void enterObservationId(const uInt /*recordNo*/, Int columnValue) {
354 //printf("%u: ObservationId: %d\n", recordNo, columnValue);
355 // update header
356 if ( header.observer.empty() )
357 getScalar( String("OBSERVER"), (uInt)columnValue, obstab, header.observer ) ;
358 if ( header.project.empty() )
359 getScalar( "PROJECT", (uInt)columnValue, obstab, header.project ) ;
360 if ( header.utc == 0.0 ) {
361 Vector<MEpoch> amp ;
362 getArrayMeas( "TIME_RANGE", (uInt)columnValue, obstab, amp ) ;
363 obsEpoch = amp[0];
364 header.utc = obsEpoch.get( "d" ).getValue() ;
365 }
366 if ( header.antennaname.empty() )
367 getScalar( "TELESCOPE_NAME", (uInt)columnValue, obstab, header.antennaname ) ;
368 }
369 virtual void leaveObservationId(const uInt /*recordNo*/, Int /*columnValue*/) {
370 // update header
371 header.nbeam = max( header.nbeam, (Int)nbeam ) ;
372
373 nbeam = 0 ;
374 feedEntry = -1 ;
375 }
376 virtual void enterFeedId(const uInt /*recordNo*/, Int columnValue) {
377 //printf("%u: FeedId: %d\n", recordNo, columnValue);
378
379 // update feed entry
380 if ( allNE( feedEntry, columnValue ) ) {
381 feedEntry[nbeam] = columnValue ;
382 nbeam++ ;
383 }
384
385 // put values
386 *beamNoRF = (uInt)columnValue ;
387 *focusIdRF = (uInt)0 ;
388 }
389 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) {
390 uInt nelem = feedEntry.nelements() ;
391 if ( nbeam > nelem ) {
392 feedEntry.resize( nelem+64, True ) ;
393 Slicer slice( IPosition( 1, nelem ), IPosition( 1, feedEntry.nelements()-1 ) ) ;
394 feedEntry( slice ) = -1 ;
395 }
396 }
397 virtual void enterFieldId(const uInt /*recordNo*/, Int columnValue) {
398 //printf("%u: FieldId: %d\n", recordNo, columnValue);
399 // update sourceId and fieldName
400 getScalar( "SOURCE_ID", (uInt)columnValue, fieldtab, sourceId ) ;
401 String fieldName ;
402 getScalar( "NAME", (uInt)columnValue, fieldtab, fieldName ) ;
403 fieldName += "__" + String::toString( columnValue ) ;
404
405 // put values
406 *fieldNameRF = fieldName ;
407 }
408 virtual void leaveFieldId(const uInt /*recordNo*/, Int /*columnValue*/) {
409 sourceId = -1 ;
410 }
411 virtual void enterDataDescId(const uInt /*recordNo*/, Int columnValue) {
412 //printf("%u: DataDescId: %d\n", recordNo, columnValue);
413 // update polarization and spectral window ids
414 getScalar( "POLARIZATION_ID", (uInt)columnValue, ddtab, polId ) ;
415 getScalar( "SPECTRAL_WINDOW_ID", (uInt)columnValue, ddtab, spwId ) ;
416
417 // polarization setup
418 getScalar( "NUM_CORR", (uInt)polId, poltab, npol ) ;
419 Vector<Int> corrtype ;
420 getArray( "CORR_TYPE", (uInt)polId, poltab, corrtype ) ;
421 polnos = getPolNos( corrtype ) ;
422
423 // process SOURCE table
424 String sourceName ;
425 Vector<Double> sourcePM, restFreqs, sysVels ;
426 Vector<String> transition ;
427 processSource( sourceId, spwId, sourceName, sourceDir, sourcePM,
428 restFreqs, transition, sysVels ) ;
429
430 // spectral setup
431 uInt freqId ;
432 Double reffreq, bandwidth ;
433 String freqref ;
434 getScalar( "NUM_CHAN", (uInt)spwId, spwtab, nchan ) ;
435 Bool iswvr = (Bool)(nchan == 4) ;
436 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
437 if ( iter == ifmap.end() ) {
438 //MEpoch me ;
439 //getScalarMeas( "TIME", recordNo, table, me ) ;
440 //spectralSetup( spwId, me, antpos, sourceDir,
441 spectralSetup(spwId, obsEpoch, antpos, sourceDir,
442 freqId, nchan,
443 freqref, reffreq, bandwidth);
444 ifmap.insert( pair<Int,uInt>(spwId,freqId) ) ;
445 }
446 else {
447 freqId = iter->second ;
448 }
449 sp.resize( npol, nchan ) ;
450 fl.resize( npol, nchan ) ;
451
452
453 // molecular setup
454 STMolecules mtab = scantable.molecules() ;
455 uInt molId = mtab.addEntry( restFreqs, transition, transition ) ;
456
457 // process SYSCAL table
458 if ( isSysCal )
459 processSysCal( spwId ) ;
460
461 // update header
462 if ( !iswvr ) {
463 header.nchan = max( header.nchan, nchan ) ;
464 header.bandwidth = max( header.bandwidth, bandwidth ) ;
465 if ( header.reffreq == -1.0 )
466 header.reffreq = reffreq ;
467 header.npol = max( header.npol, npol ) ;
468 if ( header.poltype.empty() )
469 header.poltype = getPolType( corrtype[0] ) ;
470 if ( header.freqref.empty() )
471 header.freqref = freqref ;
472 }
473
474 // put values
475 *ifNoRF = (uInt)spwId ;
476 *freqIdRF = freqId ;
477 *moleculeIdRF = molId ;
478 *sourceNameRF = sourceName ;
479 sourceProperMotionRF.define( sourcePM ) ;
480 Vector<Double> srcD = sourceDir.getAngle().getValue( "rad" ) ;
481 sourceDirectionRF.define( srcD ) ;
482 if ( !sysVels.empty() )
483 *sourceVelocityRF = sysVels[0] ;
484 else {
485 *sourceVelocityRF = (Double)0.0 ;
486 }
487 }
488 virtual void leaveDataDescId(const uInt /*recordNo*/, Int /*columnValue*/) {
489 npol = 0 ;
490 nchan = 0 ;
491 numSysCalRow = 0 ;
492 }
493 virtual void enterScanNo(const uInt /*recordNo*/, Int columnValue) {
494 //printf("%u: ScanNo: %d\n", recordNo, columnValue);
495 // put value
496 // CAS-5841: SCANNO should be consistent with MS SCAN_NUMBER
497 *scanNoRF = (uInt)columnValue ;
498 }
499 virtual void leaveScanNo(const uInt /*recordNo*/, Int /*columnValue*/) {
500 cycleNo = 0 ;
501 }
502 virtual void enterStateId(const uInt /*recordNo*/, Int columnValue) {
503 //printf("%u: StateId: %d\n", recordNo, columnValue);
504 // SRCTYPE
505 Int srcType = getSrcType( columnValue ) ;
506
507 // update header
508 if ( header.obstype.empty() )
509 getScalar( "OBS_MODE", (uInt)columnValue, statetab, header.obstype ) ;
510
511 // put value
512 *sourceTypeRF = srcType ;
513 }
514 virtual void leaveStateId(const uInt /*recordNo*/, Int /*columnValue*/) { }
515 virtual void enterTime(const uInt recordNo, Double columnValue) {
516 //printf("%u: Time: %f\n", recordNo, columnValue);
517 currentTime = MEpoch( Quantity( columnValue, "s" ), MEpoch::UTC ) ;
518
519 // DIRECTION, AZEL, and SCANRATE
520 Vector<Double> direction, azel ;
521 Vector<Double> scanrate( 2, 0.0 ) ;
522 if ( getpt )
523 getDirection( direction, azel, scanrate ) ;
524 else
525 getSourceDirection( direction, azel, scanrate ) ;
526
527 // INTERVAL
528 Double interval = intervalCol.asdouble( recordNo ) ;
529
530 // WEATHER_ID
531 uInt wid = 0 ;
532 if ( isWeather_ )
533 wid = getWeatherId() ;
534
535 // put value
536 Double t = currentTime.get( "d" ).getValue() ;
537 *timeRF = t ;
538 *intervalRF = interval ;
539 directionRF.define( direction ) ;
540 *azimuthRF = (Float)azel[0] ;
541 *elevationRF = (Float)azel[1] ;
542 scanRateRF.define( scanrate ) ;
543 *weatherIdRF = wid ;
544 }
545 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
546 virtual Bool visitRecord(const uInt recordNo,
547 const Int /*observationId*/,
548 const Int /*feedId*/,
549 const Int /*fieldId*/,
550 const Int /*dataDescId*/,
551 const Int /*scanNo*/,
552 const Int /*stateId*/,
553 const Double /*time*/)
554 {
555 //printf("%u: %d, %d, %d, %d, %d, %d, %f\n", recordNo,
556 //observationId, feedId, fieldId, dataDescId, scanNo, stateId, time);
557
558 // SPECTRA and FLAGTRA
559 //Matrix<Float> sp;
560 //Matrix<uChar> fl;
561 spectraAndFlagtra( recordNo, sp, fl ) ;
562
563 // FLAGROW
564 Bool flr = flagRowCol.asBool( recordNo ) ;
565
566 // TSYS
567 Matrix<Float> tsys ;
568 uInt scIdx = getSysCalIndex() ;
569 if ( numSysCalRow > 0 ) {
570 tsys = sysCalTsysCol( syscalRow[scIdx] ) ;
571 }
572 else {
573 tsys.resize( npol, 1 ) ;
574 tsys = 1.0 ;
575 }
576
577 // TCAL_ID
578 Block<uInt> tcalids( npol, 0 ) ;
579 if ( numSysCalRow > 0 ) {
580 tcalids = getTcalId( syscalTime[scIdx] ) ;
581 }
582 else {
583 tcalids = getDummyTcalId( spwId ) ;
584 }
585
586 // put value
587 *cycleNoRF = cycleNo ;
588 *flagRowRF = (uInt)flr ;
589
590 // for each polarization component
591 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
592 // put value depending on polarization component
593 *polNoRF = polnos[ipol] ;
594 *tcalIdRF = tcalids[ipol] ;
595 spectraRF.define( sp.row( ipol ) ) ;
596 flagtraRF.define( fl.row( ipol ) ) ;
597 tsysRF.define( tsys.row( ipol ) ) ;
598
599 // commit row
600 tablerow.put( rowidx ) ;
601 rowidx++ ;
602 }
603
604 // increment CYCLENO
605 cycleNo++ ;
606
607 return True ;
608 }
609 virtual void finish()
610 {
611 BaseMSFillerVisitor::finish();
612 //printf("Total: %u\n", count);
613 // remove redundant rows
614 //cout << "filled " << rowidx << " rows out of " << scantable.nrow() << " rows" << endl ;
615 if ( scantable.nrow() > (Int)rowidx ) {
616 uInt numRemove = scantable.nrow() - rowidx ;
617 //cout << "numRemove = " << numRemove << endl ;
618 Vector<uInt> rows( numRemove ) ;
619 indgen( rows, rowidx ) ;
620 scantable.table().removeRow( rows ) ;
621 }
622
623 // antenna name and station name
624 String antennaName ;
625 getScalar( "NAME", (uInt)antennaId, anttab, antennaName ) ;
626 String stationName ;
627 getScalar( "STATION", (uInt)antennaId, anttab, stationName ) ;
628
629 // update header
630 header.nif = ifmap.size() ;
631 header.antennaposition = antpos.get( "m" ).getValue() ;
632 if ( header.antennaname.empty() || header.antennaname == antennaName )
633 header.antennaname = antennaName ;
634 else
635 header.antennaname += "//" + antennaName ;
636 if ( !stationName.empty() && stationName != antennaName )
637 header.antennaname += "@" + stationName ;
638 if ( header.fluxunit.empty() || header.fluxunit == "CNTS" )
639 header.fluxunit = "K" ;
640 header.epoch = "UTC" ;
641 header.equinox = 2000.0 ;
642 if (header.freqref == "TOPO") {
643 header.freqref = "TOPOCENT";
644 } else if (header.freqref == "GEO") {
645 header.freqref = "GEOCENTR";
646 } else if (header.freqref == "BARY") {
647 header.freqref = "BARYCENT";
648 } else if (header.freqref == "GALACTO") {
649 header.freqref = "GALACTOC";
650 } else if (header.freqref == "LGROUP") {
651 header.freqref = "LOCALGRP";
652 } else if (header.freqref == "CMB") {
653 header.freqref = "CMBDIPOL";
654 } else if (header.freqref == "REST") {
655 header.freqref = "SOURCE";
656 }
657 scantable.setHeader( header ) ;
658 }
659 void setAntenna( Int id )
660 {
661 antennaId = id ;
662
663 Vector< Quantum<Double> > pos ;
664 getArrayQuant( "POSITION", (uInt)antennaId, anttab, pos ) ;
665 antpos = MPosition( MVPosition( pos ), MPosition::ITRF ) ;
666 mf.set( antpos ) ;
667 }
668 void setPointingTable( const Table &tab, String columnToUse="DIRECTION" )
669 {
670 // input POINTING table must be
671 // 1) selected by antenna
672 // 2) sorted by TIME
673 ROScalarColumn<Double> tcol( tab, "TIME" ) ;
674 ROArrayColumn<Double> dcol( tab, columnToUse ) ;
675 tcol.getColumn( pointingTime ) ;
676 dcol.getColumn( pointingDirection ) ;
677 const TableRecord &rec = dcol.keywordSet() ;
678 String pointingRef = rec.asRecord( "MEASINFO" ).asString( "Ref" ) ;
679 MDirection::getType( dirType, pointingRef ) ;
680 getpt = True ;
681
682 // initialize toj2000 and toazel
683 initConvert() ;
684 }
685 void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it,
686 const Vector<uInt> &idx )
687 {
688 isWeather_ = True ;
689 weatherTime_ = t ;
690 weatherInterval_ = it ;
691 weatherIndex_ = idx;
692 }
693 void setSysCalRecord( const Record &r )
694 //void setSysCalRecord( const map< String,Vector<uInt> > &r )
695 {
696 isSysCal = True ;
697 isTcal = True ;
698 syscalRecord = r ;
699 if ( syscalRecord.nfields() == 0 )
700 isTcal = False ;
701
702 const TableDesc &desc = sctab.tableDesc() ;
703 uInt nrow = sctab.nrow() ;
704 syscalRow.resize( nrow ) ;
705 syscalTime.resize( nrow ) ;
706 syscalInterval.resize( nrow ) ;
707 String tsysCol = "NONE" ;
708 Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ;
709 for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) {
710 if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) )
711 tsysCol = tsysCols[i] ;
712 }
713 sysCalTsysCol.attach( sctab, tsysCol ) ;
714 }
715 STHeader getHeader() { return header ; }
716 uInt getNumBeam() { return nbeam ; }
717 uInt getFilledRowNum() { return rowidx ; }
718private:
719 void initConvert()
720 {
721 toj2000 = MDirection::Convert( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
722 toazel = MDirection::Convert( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
723 }
724
725 void fluxUnit( String &u )
726 {
727 ROTableColumn col( table, dataColumnName ) ;
728 const TableRecord &rec = col.keywordSet() ;
729 if ( rec.isDefined( "UNIT" ) )
730 u = rec.asString( "UNIT" ) ;
731 else if ( rec.isDefined( "QuantumUnits" ) )
732 u = rec.asString( "QuantumUnits" ) ;
733 if ( u.empty() )
734 u = "K" ;
735 }
736 void processSource( Int sourceId, Int spwId,
737 String &name, MDirection &dir, Vector<Double> &pm,
738 Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel )
739 {
740 // find row
741 uInt nrow = srctab.nrow() ;
742 Int idx = -1 ;
743 ROTableRow row( srctab ) ;
744 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
745 const TableRecord &r = row.get( irow ) ;
746 if ( r.asInt( "SOURCE_ID" ) == sourceId ) {
747 Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ;
748 if ( tmpSpwId == spwId || tmpSpwId == -1 ) {
749 idx = (Int)irow ;
750 break ;
751 }
752 }
753 }
754
755 // fill
756 Int numLines = 0 ;
757 if ( idx != -1 ) {
758 const TableRecord &r = row.get( idx ) ;
759 name = r.asString( "NAME" ) ;
760 getScalarMeas( "DIRECTION", idx, srctab, dir ) ;
761 pm = r.toArrayDouble( "PROPER_MOTION" ) ;
762 numLines = r.asInt( "NUM_LINES" ) ;
763 }
764 else {
765 name = "" ;
766 pm = Vector<Double>( 2, 0.0 ) ;
767 dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
768 }
769 if ( !getpt ) {
770 String ref = dir.getRefString() ;
771 MDirection::getType( dirType, ref ) ;
772
773 // initialize toj2000 and toazel
774 initConvert() ;
775 }
776
777 rf.resize( numLines ) ;
778 trans.resize( numLines ) ;
779 vel.resize( numLines ) ;
780 if ( numLines > 0 ) {
781 Block<Bool> isDefined = row.getDefined() ;
782 Vector<String> colNames = row.columnNames() ;
783 Vector<Int> indexes( 3, -1 ) ;
784 Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ;
785 for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) {
786 if ( anyEQ( indexes, -1 ) ) {
787 for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) {
788 if ( colNames[icol] == cols[jcol] )
789 indexes[jcol] = icol ;
790 }
791 }
792 }
793 if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) {
794 Vector< Quantum<Double> > qrf ;
795 getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ;
796 for ( int i = 0 ; i < numLines ; i++ )
797 rf[i] = qrf[i].getValue( "Hz" ) ;
798 }
799 if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) {
800 getArray( "TRANSITION", idx, srctab, trans ) ;
801 }
802 if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) {
803 Vector< Quantum<Double> > qsv ;
804 getArrayQuant( "SYSVEL", idx, srctab, qsv ) ;
805 for ( int i = 0 ; i < numLines ; i++ )
806 vel[i] = qsv[i].getValue( "m/s" ) ;
807 }
808 }
809 }
810 void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md,
811 uInt &freqId, Int &nchan,
812 String &freqref, Double &reffreq, Double &bandwidth )
813 {
814 // fill
815 Int measFreqRef ;
816 getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ;
817 MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ;
818 //freqref = MFrequency::showType( freqRef ) ;
819 //freqref = "LSRK" ;
820 freqref = "TOPO";
821 Quantum<Double> q ;
822 getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ;
823 bandwidth = q.getValue( "Hz" ) ;
824 getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ;
825 reffreq = q.getValue( "Hz" ) ;
826 Double refpix = 0.5 * ( (Double)nchan-1.0 ) ;
827 Int refchan = ( nchan - 1 ) / 2 ;
828 Bool even = (Bool)( nchan % 2 == 0 ) ;
829 Vector< Quantum<Double> > qa ;
830 getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ;
831// Double increment = qa[refchan].getValue( "Hz" ) ;
832 Double increment = abs(qa[refchan].getValue( "Hz" )) ;
833 getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ;
834 if ( nchan == 1 ) {
835 Int netSideband ;
836 getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ;
837 if ( netSideband == 1 ) increment *= -1.0 ;
838 }
839 else {
840 if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) )
841 increment *= -1.0 ;
842 }
843 Double refval = qa[refchan].getValue( "Hz" ) ;
844 if ( even )
845 refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ;
846
847 // add new row to FREQUENCIES
848 Table ftab = scantable.frequencies().table() ;
849 freqId = ftab.nrow() ;
850 ftab.addRow() ;
851 TableRow row( ftab ) ;
852 TableRecord &r = row.record() ;
853 RecordFieldPtr<uInt> idRF( r, "ID" ) ;
854 *idRF = freqId ;
855 RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ;
856 RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ;
857 RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ;
858 *refpixRF = refpix ;
859 *refvalRF = refval ;
860 *incrRF = increment ;
861 row.put( freqId ) ;
862 }
863 void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl )
864 {
865 Matrix<Bool> b = flagCol( recordNo ) ;
866 if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) {
867 sp = floatDataCol( recordNo ) ;
868 convertArray( fl, b ) ;
869 }
870 else {
871 Bool notyet = True ;
872 Matrix<Complex> c = dataCol( recordNo ) ;
873 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
874 if ( ( header.poltype == "linear" || header.poltype == "circular" )
875 && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) {
876 if ( notyet ) {
877 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
878 IPosition start( 1, 0 ) ;
879 IPosition end( 1, 2*nchan-1 ) ;
880 IPosition inc( 1, 2 ) ;
881 if ( polnos[ipol] == 2 ) {
882 sp.row( ipol ) = tmp( start, end, inc ) ;
883 Vector<Bool> br = b.row( ipol ) ;
884 Vector<uChar> flr = fl.row( ipol ) ;
885 convertArray( flr, br ) ;
886 start = IPosition( 1, 1 ) ;
887 Int jpol = ipol+1 ;
888 while( polnos[jpol] != 3 && jpol < npol )
889 jpol++ ;
890 sp.row( jpol ) = tmp( start, end, inc ) ;
891 flr.reference( fl.row( jpol ) ) ;
892 convertArray( flr, br ) ;
893 }
894 else if ( polnos[ipol] == 3 ) {
895 sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ;
896 Int jpol = ipol+1 ;
897 while( polnos[jpol] != 2 && jpol < npol )
898 jpol++ ;
899 Vector<Bool> br = b.row( ipol ) ;
900 Vector<uChar> flr = fl.row( jpol ) ;
901 sp.row( jpol ) = tmp( start, end, inc ) ;
902 convertArray( flr, br ) ;
903 start = IPosition( 1, 1 ) ;
904 sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ;
905 flr.reference( fl.row( ipol ) ) ;
906 convertArray( flr, br ) ;
907 }
908 notyet = False ;
909 }
910 }
911 else {
912 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
913 IPosition start( 1, 0 ) ;
914 IPosition end( 1, 2*nchan-1 ) ;
915 IPosition inc( 1, 2 ) ;
916 sp.row( ipol ) = tmp( start, end, inc ) ;
917 Vector<Bool> br = b.row( ipol ) ;
918 Vector<uChar> flr = fl.row( ipol ) ;
919 convertArray( flr, br ) ;
920 }
921 }
922 }
923 }
924 uInt binarySearch( Vector<Double> &timeList, Double target )
925 {
926 Int low = 0 ;
927 Int high = timeList.nelements() ;
928 uInt idx = 0 ;
929
930 while ( low <= high ) {
931 idx = (Int)( 0.5 * ( low + high ) ) ;
932 Double t = timeList[idx] ;
933 if ( t < target )
934 low = idx + 1 ;
935 else if ( t > target )
936 high = idx - 1 ;
937 else {
938 return idx ;
939 }
940 }
941
942 idx = max( 0, min( low, high ) ) ;
943 return idx ;
944 }
945 void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
946 {
947 // @todo At the moment, do binary search every time
948 // if this is bottleneck, frequency of binary search must be reduced
949 Double t = currentTime.get( "s" ).getValue() ;
950 uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ;
951 Matrix<Double> d ;
952 if ( pointingTime[idx] == t )
953 d = pointingDirection.xyPlane( idx ) ;
954 else if ( pointingTime[idx] < t ) {
955 if ( idx == pointingTime.nelements()-1 )
956 d = pointingDirection.xyPlane( idx ) ;
957 else
958 d = interp( pointingTime[idx], pointingTime[idx+1], t,
959 pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ;
960 }
961 else {
962 if ( idx == 0 )
963 d = pointingDirection.xyPlane( idx ) ;
964 else
965 d = interp( pointingTime[idx-1], pointingTime[idx], t,
966 pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
967 }
968 mf.set( currentTime ) ;
969 Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
970 if ( dirType != MDirection::J2000 ) {
971 dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
972 }
973 else {
974 dir = d.column( 0 ) ;
975 }
976 if ( dirType != MDirection::AZELGEO ) {
977 azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
978 }
979 else {
980 azel = d.column( 0 ) ;
981 }
982 if ( d.ncolumn() > 1 )
983 srate = d.column( 1 ) ;
984 }
985 void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &/*srate*/ )
986 {
987 dir = sourceDir.getAngle( "rad" ).getValue() ;
988 mf.set( currentTime ) ;
989 azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
990 if ( dirType != MDirection::J2000 ) {
991 dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
992 }
993 }
994 String detectSeparator( String &s )
995 {
996 String tmp = s.substr( 0, s.find_first_of( "," ) ) ;
997 Char *separators[] = { ":", "#", ".", "_" } ;
998 uInt nsep = 4 ;
999 for ( uInt i = 0 ; i < nsep ; i++ ) {
1000 if ( tmp.find( separators[i] ) != String::npos )
1001 return separators[i] ;
1002 }
1003 return "" ;
1004 }
1005 Int getSrcType( Int stateId )
1006 {
1007 // get values
1008 Bool sig ;
1009 getScalar( "SIG", stateId, statetab, sig ) ;
1010 Bool ref ;
1011 getScalar( "REF", stateId, statetab, ref ) ;
1012 Double cal ;
1013 getScalar( "CAL", stateId, statetab, cal ) ;
1014 String obsmode ;
1015 getScalar( "OBS_MODE", stateId, statetab, obsmode ) ;
1016 String sep = detectSeparator( obsmode ) ;
1017
1018 Int srcType = SrcType::NOTYPE ;
1019 if ( sep == ":" )
1020 srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ;
1021 else if ( sep == "." || sep == "#" )
1022 srcTypeALMA( srcType, sep, obsmode ) ;
1023 else if ( sep == "_" )
1024 srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ;
1025 else
1026 srcTypeDefault( srcType, sig, ref ) ;
1027
1028 return srcType ;
1029 }
1030 void srcTypeDefault( Int &st, Bool &sig, Bool &ref )
1031 {
1032 if ( sig ) st = SrcType::SIG ;
1033 else if ( ref ) st = SrcType::REF ;
1034 }
1035 void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal )
1036 {
1037 Int epos = mode.find_first_of( sep ) ;
1038 Int nextpos = mode.find_first_of( sep, epos+1 ) ;
1039 String m1 = mode.substr( 0, epos ) ;
1040 String m2 = mode.substr( epos+1, nextpos-epos-1 ) ;
1041 if ( m1 == "Nod" ) {
1042 st = SrcType::NOD ;
1043 }
1044 else if ( m1 == "OffOn" ) {
1045 if ( m2 == "PSWITCHON" ) st = SrcType::PSON ;
1046 if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ;
1047 }
1048 else {
1049 if ( m2 == "FSWITCH" ) {
1050 if ( sig ) st = SrcType::FSON ;
1051 else if ( ref ) st = SrcType::FSOFF ;
1052 }
1053 }
1054 if ( cal > 0.0 ) {
1055 if ( st == SrcType::NOD )
1056 st = SrcType::NODCAL ;
1057 else if ( st == SrcType::PSON )
1058 st = SrcType::PONCAL ;
1059 else if ( st == SrcType::PSOFF )
1060 st = SrcType::POFFCAL ;
1061 else if ( st == SrcType::FSON )
1062 st = SrcType::FONCAL ;
1063 else if ( st == SrcType::FSOFF )
1064 st = SrcType::FOFFCAL ;
1065 else
1066 st = SrcType::CAL ;
1067 }
1068 }
1069 void srcTypeALMA( Int &st, String &sep, String &mode )
1070 {
1071 Int epos = mode.find_first_of( "," ) ;
1072 String first = mode.substr( 0, epos ) ;
1073 epos = first.find_first_of( sep ) ;
1074 Int nextpos = first.find_first_of( sep, epos+1 ) ;
1075 String m1 = first.substr( 0, epos ) ;
1076 String m2 = first.substr( epos+1, nextpos-epos-1 ) ;
1077 if ( m1.find( "CALIBRATE_" ) == 0 ) {
1078 if ( m2.find( "ON_SOURCE" ) == 0 )
1079 st = SrcType::PONCAL ;
1080 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1081 st = SrcType::POFFCAL ;
1082 }
1083 else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) {
1084 if ( m2.find( "ON_SOURCE" ) == 0 )
1085 st = SrcType::PSON ;
1086 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1087 st = SrcType::PSOFF ;
1088 }
1089 }
1090 void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref )
1091 {
1092 Int epos = mode.find_first_of( "," ) ;
1093 String first = mode.substr( 0, epos ) ;
1094 string substr[4] ;
1095 int numSubstr = split( first, substr, 4, sep ) ;
1096 String m1( substr[0] ) ;
1097 String m2( substr[2] ) ;
1098 if ( numSubstr == 4 ) {
1099 if ( m1.find( "CALIBRATE" ) == 0 ) {
1100 if ( m2.find( "ON" ) == 0 )
1101 st = SrcType::PONCAL ;
1102 else if ( m2.find( "OFF" ) == 0 )
1103 st = SrcType::POFFCAL ;
1104 }
1105 else if ( m1.find( "OBSERVE" ) == 0 ) {
1106 if ( m2.find( "ON" ) == 0 )
1107 st = SrcType::PSON ;
1108 else if ( m2.find( "OFF" ) == 0 )
1109 st = SrcType::PSOFF ;
1110 }
1111 }
1112 else {
1113 if ( sig ) st = SrcType::SIG ;
1114 else if ( ref ) st = SrcType::REF ;
1115 }
1116 }
1117 Block<uInt> getPolNos( Vector<Int> &corr )
1118 {
1119 Block<uInt> polnos( npol ) ;
1120 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
1121 if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX )
1122 polnos[ipol] = 0 ;
1123 else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY )
1124 polnos[ipol] = 1 ;
1125 else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY )
1126 polnos[ipol] = 2 ;
1127 else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX )
1128 polnos[ipol] = 3 ;
1129 }
1130 return polnos ;
1131 }
1132 String getPolType( Int &corr )
1133 {
1134 String poltype = "" ;
1135 if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V )
1136 poltype = "stokes" ;
1137 else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX )
1138 poltype = "linear" ;
1139 else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR )
1140 poltype = "circular" ;
1141 else if ( corr == Stokes::Plinear || corr == Stokes::Pangle )
1142 poltype = "linpol" ;
1143 return poltype ;
1144 }
1145 uInt getWeatherId()
1146 {
1147 // if only one row, return 0
1148 if ( weatherTime_.nelements() == 1 )
1149 return 0 ;
1150
1151 // @todo At the moment, do binary search every time
1152 // if this is bottleneck, frequency of binary search must be reduced
1153 Double t = currentTime.get( "s" ).getValue() ;
1154 uInt idx = min( binarySearch( weatherTime_, t ), weatherTime_.nelements()-1 ) ;
1155 if ( weatherTime_[idx] < t ) {
1156 if ( idx != weatherTime_.nelements()-1 ) {
1157 if ( weatherTime_[idx+1] - t < 0.5 * weatherInterval_[idx+1] )
1158 idx++ ;
1159 }
1160 }
1161 else if ( weatherTime_[idx] > t ) {
1162 if ( idx != 0 ) {
1163 if ( weatherTime_[idx] - t > 0.5 * weatherInterval_[idx] )
1164 idx-- ;
1165 }
1166 }
1167 return weatherIndex_[idx] ;
1168 }
1169 void processSysCal( Int &spwId )
1170 {
1171 // get feedId from row
1172 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1173
1174 uInt nrow = sctab.nrow() ;
1175 ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ;
1176 Vector<Int> aids = col.getColumn() ;
1177 col.attach( sctab, "FEED_ID" ) ;
1178 Vector<Int> fids = col.getColumn() ;
1179 col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ;
1180 Vector<Int> sids = col.getColumn() ;
1181 ROScalarColumn<Double> timeCol( sctab, "TIME" ) ;
1182 ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ;
1183 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1184 if ( aids[irow] == antennaId
1185 && fids[irow] == feedId
1186 && sids[irow] == spwId ) {
1187 syscalRow[numSysCalRow] = irow ;
1188 syscalTime[numSysCalRow] = timeCol( irow ) ;
1189 syscalInterval[numSysCalRow] = intCol( irow ) ;
1190 numSysCalRow++ ;
1191 }
1192 }
1193 }
1194 uInt getSysCalIndex()
1195 {
1196 // if only one row, return 0
1197 if ( numSysCalRow == 1 || !isSysCal )
1198 return 0 ;
1199
1200 // @todo At the moment, do binary search every time
1201 // if this is bottleneck, frequency of binary search must be reduced
1202 Double t = currentTime.get( "s" ).getValue() ;
1203 Vector<Double> tslice = syscalTime( Slice(0, numSysCalRow) ) ;
1204 uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ;
1205 if ( syscalTime[idx] < t ) {
1206 if ( idx != numSysCalRow-1 ) {
1207 if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] )
1208 idx++ ;
1209 }
1210 }
1211 else if ( syscalTime[idx] > t ) {
1212 if ( idx != 0 ) {
1213 if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] )
1214 idx-- ;
1215 }
1216 }
1217 return idx ;
1218 }
1219 Block<uInt> getTcalId( Double &t )
1220 {
1221 // return 0 if no SysCal table
1222 if ( !isSysCal or !isTcal ) {
1223 return Block<uInt>( 4, 0 ) ;
1224 }
1225
1226 // get feedId from row
1227 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1228
1229 // key
1230 String key = keyTcal( feedId, spwId, t ) ;
1231
1232 // retrieve ids
1233 Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ;
1234 //Vector<uInt> ids = syscalRecord[key] ;
1235 uInt np = ids[1] - ids[0] + 1 ;
1236 Block<uInt> tcalids( np ) ;
1237 if ( np > 0 ) {
1238 tcalids[0] = ids[0] ;
1239 if ( np > 1 ) {
1240 tcalids[1] = ids[1] ;
1241 for ( uInt ip = 2 ; ip < np ; ip++ )
1242 tcalids[ip] = ids[0] + ip - 1 ;
1243 }
1244 }
1245 return tcalids ;
1246 }
1247 Block<uInt> getDummyTcalId( Int spwId )
1248 {
1249 Block<uInt> idList(4, 0);
1250 uInt nfields = syscalRecord.nfields();
1251 Int idx = -1;
1252 for (uInt i = 0; i< nfields ; i++ ) {
1253 String spw = "SPW" + String::toString(spwId);
1254 if (syscalRecord.name(i).find(spw) != String::npos) {
1255 idx = i;
1256 break;
1257 }
1258 }
1259 if ( idx > -1) {
1260 Vector<uInt> tmp = syscalRecord.asArrayuInt(idx);
1261 for (uInt j = 0 ; j < 4 ; j++) {
1262 idList[j] = tmp[0];
1263 }
1264 }
1265 return idList;
1266 }
1267 uInt maxNumPol()
1268 {
1269 ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ;
1270 return max( numCorrCol.getColumn() ) ;
1271 }
1272
1273 Scantable &scantable;
1274 Int antennaId;
1275 uInt rowidx;
1276 String dataColumnName;
1277 TableRow tablerow;
1278 STHeader header;
1279 Vector<Int> feedEntry;
1280 uInt nbeam;
1281 Int npol;
1282 Int nchan;
1283 Int sourceId;
1284 Int polId;
1285 Int spwId;
1286 uInt cycleNo;
1287 MDirection sourceDir;
1288 MPosition antpos;
1289 MEpoch currentTime;
1290 MEpoch obsEpoch;
1291 MeasFrame mf;
1292 MDirection::Convert toj2000;
1293 MDirection::Convert toazel;
1294 map<Int,uInt> ifmap;
1295 Block<uInt> polnos;
1296 Bool getpt;
1297 Vector<Double> pointingTime;
1298 Cube<Double> pointingDirection;
1299 MDirection::Types dirType;
1300 Bool isWeather_;
1301 Vector<Double> weatherTime_;
1302 Vector<Double> weatherInterval_;
1303 Vector<uInt> weatherIndex_;
1304 Bool isSysCal;
1305 Bool isTcal;
1306 Record syscalRecord;
1307 //map< String,Vector<uInt> > syscalRecord;
1308 uInt numSysCalRow ;
1309 Vector<uInt> syscalRow;
1310 Vector<Double> syscalTime;
1311 Vector<Double> syscalInterval;
1312 //String tsysCol;
1313 //String tcalCol;
1314
1315 // MS subtables
1316 Table obstab;
1317 Table sctab;
1318 Table spwtab;
1319 Table statetab;
1320 Table ddtab;
1321 Table poltab;
1322 Table fieldtab;
1323 Table anttab;
1324 Table srctab;
1325 Matrix<Float> sp;
1326 Matrix<uChar> fl;
1327
1328 // MS MAIN columns
1329 ROTableColumn intervalCol;
1330 ROTableColumn flagRowCol;
1331 ROArrayColumn<Float> floatDataCol;
1332 ROArrayColumn<Complex> dataCol;
1333 ROArrayColumn<Bool> flagCol;
1334
1335 // MS SYSCAL columns
1336 ROArrayColumn<Float> sysCalTsysCol;
1337
1338 // Scantable MAIN columns
1339 RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF;
1340 RecordFieldPtr< Vector<Double> > directionRF,scanRateRF,
1341 sourceProperMotionRF,sourceDirectionRF;
1342 RecordFieldPtr<Float> azimuthRF,elevationRF;
1343 RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF,
1344 ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF;
1345 RecordFieldPtr< Vector<Float> > spectraRF,tsysRF;
1346 RecordFieldPtr< Vector<uChar> > flagtraRF;
1347 RecordFieldPtr<String> sourceNameRF,fieldNameRF;
1348 RecordFieldPtr<Int> sourceTypeRF;
1349};
1350
1351class BaseTcalVisitor: public TableVisitor {
1352 uInt lastRecordNo ;
1353 Int lastAntennaId ;
1354 Int lastFeedId ;
1355 Int lastSpwId ;
1356 Double lastTime ;
1357protected:
1358 const Table &table;
1359 uInt count;
1360public:
1361 BaseTcalVisitor(const Table &table)
1362 : table(table)
1363 {
1364 count = 0;
1365 }
1366
1367 virtual void enterAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1368 virtual void leaveAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1369 virtual void enterFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1370 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1371 virtual void enterSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1372 virtual void leaveSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1373 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1374 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1375
1376 virtual Bool visitRecord(const uInt /*recordNo*/,
1377 const Int /*antennaId*/,
1378 const Int /*feedId*/,
1379 const Int /*spwId*/,
1380 const Double /*time*/) { return True ; }
1381
1382 virtual Bool visit(Bool isFirst, const uInt recordNo,
1383 const uInt nCols, void const *const colValues[]) {
1384 Int antennaId, feedId, spwId;
1385 Double time;
1386 { // prologue
1387 uInt i = 0;
1388 {
1389 const Int *col = (const Int *)colValues[i++];
1390 antennaId = col[recordNo];
1391 }
1392 {
1393 const Int *col = (const Int *)colValues[i++];
1394 feedId = col[recordNo];
1395 }
1396 {
1397 const Int *col = (const Int *)colValues[i++];
1398 spwId = col[recordNo];
1399 }
1400 {
1401 const Double *col = (const Double *)colValues[i++];
1402 time = col[recordNo];
1403 }
1404 assert(nCols == i);
1405 }
1406
1407 if (isFirst) {
1408 enterAntennaId(recordNo, antennaId);
1409 enterFeedId(recordNo, feedId);
1410 enterSpwId(recordNo, spwId);
1411 enterTime(recordNo, time);
1412 } else {
1413 if ( lastAntennaId != antennaId ) {
1414 leaveTime(lastRecordNo, lastTime);
1415 leaveSpwId(lastRecordNo, lastSpwId);
1416 leaveFeedId(lastRecordNo, lastFeedId);
1417 leaveAntennaId(lastRecordNo, lastAntennaId);
1418
1419 enterAntennaId(recordNo, antennaId);
1420 enterFeedId(recordNo, feedId);
1421 enterSpwId(recordNo, spwId);
1422 enterTime(recordNo, time);
1423 }
1424 else if (lastFeedId != feedId) {
1425 leaveTime(lastRecordNo, lastTime);
1426 leaveSpwId(lastRecordNo, lastSpwId);
1427 leaveFeedId(lastRecordNo, lastFeedId);
1428
1429 enterFeedId(recordNo, feedId);
1430 enterSpwId(recordNo, spwId);
1431 enterTime(recordNo, time);
1432 } else if (lastSpwId != spwId) {
1433 leaveTime(lastRecordNo, lastTime);
1434 leaveSpwId(lastRecordNo, lastSpwId);
1435
1436 enterSpwId(recordNo, spwId);
1437 enterTime(recordNo, time);
1438 } else if (lastTime != time) {
1439 leaveTime(lastRecordNo, lastTime);
1440 enterTime(recordNo, time);
1441 }
1442 }
1443 count++;
1444 Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time);
1445
1446 { // epilogue
1447 lastRecordNo = recordNo;
1448
1449 lastAntennaId = antennaId;
1450 lastFeedId = feedId;
1451 lastSpwId = spwId;
1452 lastTime = time;
1453 }
1454 return result ;
1455 }
1456
1457 virtual void finish() {
1458 if (count > 0) {
1459 leaveTime(lastRecordNo, lastTime);
1460 leaveSpwId(lastRecordNo, lastSpwId);
1461 leaveFeedId(lastRecordNo, lastFeedId);
1462 leaveAntennaId(lastRecordNo, lastAntennaId);
1463 }
1464 }
1465};
1466
1467class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
1468public:
1469 TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid )
1470 //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid )
1471 : BaseTcalVisitor( table ),
1472 tcal(tcaltab),
1473 rec(r),
1474 antenna(aid)
1475 {
1476 process = False ;
1477 rowidx = 0 ;
1478
1479 // attach to SYSCAL columns
1480 timeCol.attach( table, "TIME" ) ;
1481
1482 // add rows
1483 uInt addrow = table.nrow() * 4 ;
1484 tcal.addRow( addrow ) ;
1485
1486 // attach to TCAL columns
1487 row = TableRow( tcal ) ;
1488 TableRecord &trec = row.record() ;
1489 idRF.attachToRecord( trec, "ID" ) ;
1490 timeRF.attachToRecord( trec, "TIME" ) ;
1491 tcalRF.attachToRecord( trec, "TCAL" ) ;
1492 }
1493
1494 virtual void enterAntennaId(const uInt /*recordNo*/, Int columnValue) {
1495 if ( columnValue == antenna )
1496 process = True ;
1497 }
1498 virtual void leaveAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) {
1499 process = False ;
1500 }
1501 virtual void enterFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1502 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1503 virtual void enterSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1504 virtual void leaveSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1505 virtual void enterTime(const uInt recordNo, Double /*columnValue*/) {
1506 qtime = timeCol( recordNo ) ;
1507 }
1508 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1509 virtual Bool visitRecord(const uInt recordNo,
1510 const Int /*antennaId*/,
1511 const Int feedId,
1512 const Int spwId,
1513 const Double /*time*/)
1514 {
1515 //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ;
1516 if ( process ) {
1517 String sTime = MVTime( qtime ).string( MVTime::YMD ) ;
1518 *timeRF = sTime ;
1519 uInt oldidx = rowidx ;
1520 Matrix<Float> subtcal = tcalCol( recordNo ) ;
1521 Vector<uInt> idminmax( 2 ) ;
1522 for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) {
1523 *idRF = rowidx ;
1524 tcalRF.define( subtcal.row( ipol ) ) ;
1525
1526 // commit row
1527 row.put( rowidx ) ;
1528 rowidx++ ;
1529 }
1530
1531 idminmax[0] = oldidx ;
1532 idminmax[1] = rowidx - 1 ;
1533
1534 String key = keyTcal( feedId, spwId, sTime ) ;
1535 rec.define( key, idminmax ) ;
1536 //rec[key] = idminmax ;
1537 }
1538 return True ;
1539 }
1540 virtual void finish()
1541 {
1542 BaseTcalVisitor::finish() ;
1543
1544 if ( tcal.nrow() > rowidx ) {
1545 uInt numRemove = tcal.nrow() - rowidx ;
1546 //cout << "numRemove = " << numRemove << endl ;
1547 Vector<uInt> rows( numRemove ) ;
1548 indgen( rows, rowidx ) ;
1549 tcal.removeRow( rows ) ;
1550 }
1551
1552 }
1553 void setTcalColumn( String &col )
1554 {
1555 //colName = col ;
1556 tcalCol.attach( table, col ) ;
1557 }
1558private:
1559 Table &tcal;
1560 Record &rec;
1561 //map< String,Vector<uInt> > &rec;
1562 Int antenna;
1563 uInt rowidx;
1564 Bool process;
1565 Quantum<Double> qtime;
1566 TableRow row;
1567 String colName;
1568
1569 // MS SYSCAL columns
1570 ROScalarQuantColumn<Double> timeCol;
1571 ROArrayColumn<Float> tcalCol;
1572
1573 // TCAL columns
1574 RecordFieldPtr<uInt> idRF;
1575 RecordFieldPtr<String> timeRF;
1576 RecordFieldPtr< Vector<Float> > tcalRF;
1577};
1578
1579MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
1580 : table_( stable ),
1581 tablename_( "" ),
1582 antenna_( -1 ),
1583 antennaStr_(""),
1584 getPt_( True ),
1585 isFloatData_( False ),
1586 isData_( False ),
1587 isDoppler_( False ),
1588 isFlagCmd_( False ),
1589 isFreqOffset_( False ),
1590 isHistory_( False ),
1591 isProcessor_( False ),
1592 isSysCal_( False ),
1593 isWeather_( False ),
1594 colTsys_( "TSYS_SPECTRUM" ),
1595 colTcal_( "TCAL_SPECTRUM" )
1596{
1597 os_ = LogIO() ;
1598 os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
1599}
1600
1601MSFiller::~MSFiller()
1602{
1603 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
1604}
1605
1606bool MSFiller::open( const std::string &filename, const casa::Record &rec )
1607{
1608 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
1609 //double startSec = mathutil::gettimeofday_sec() ;
1610 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
1611 //os_ << " filename = " << filename << endl ;
1612
1613 // parsing MS options
1614 if ( rec.isDefined( "ms" ) ) {
1615 Record msrec = rec.asRecord( "ms" ) ;
1616 if ( msrec.isDefined( "getpt" ) ) {
1617 getPt_ = msrec.asBool( "getpt" ) ;
1618 }
1619 if ( msrec.isDefined( "antenna" ) ) {
1620 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
1621 antenna_ = msrec.asInt( "antenna" ) ;
1622 }
1623 else {
1624 //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
1625 antennaStr_ = msrec.asString( "antenna" ) ;
1626 }
1627 }
1628 else {
1629 antenna_ = 0 ;
1630 }
1631 }
1632
1633 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
1634 tablename_ = tmpMS->tableName() ;
1635 if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
1636 MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
1637 Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
1638 if ( id.size() > 0 )
1639 antenna_ = id[0] ;
1640 else {
1641 delete tmpMS ;
1642 //throw( AipsError( "Antenna " + antennaStr_ + " doesn't exist." ) ) ;
1643 os_ << LogIO::SEVERE << "Antenna " << antennaStr_ << " doesn't exist." << LogIO::POST ;
1644 return False ;
1645 }
1646 }
1647
1648 os_ << "Parsing MS options" << endl ;
1649 os_ << " getPt = " << (getPt_ ? "True" : "False") << endl ;
1650 os_ << " antenna = " << antenna_ << endl ;
1651 os_ << " antennaStr = " << antennaStr_ << LogIO::POST;
1652
1653 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
1654 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
1655
1656 delete tmpMS ;
1657
1658 // check which data column exists
1659 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
1660 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
1661
1662 //double endSec = mathutil::gettimeofday_sec() ;
1663 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1664 return true ;
1665}
1666
1667void MSFiller::fill()
1668{
1669 //double startSec = mathutil::gettimeofday_sec() ;
1670 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
1671
1672 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
1673
1674 // Initialize header
1675 STHeader sdh ;
1676 initHeader( sdh ) ;
1677 table_->setHeader( sdh ) ;
1678
1679 // check if optional table exists
1680 const TableRecord &msrec = mstable_.keywordSet() ;
1681 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
1682 if ( isDoppler_ )
1683 if ( mstable_.doppler().nrow() == 0 )
1684 isDoppler_ = False ;
1685 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
1686 if ( isFlagCmd_ )
1687 if ( mstable_.flagCmd().nrow() == 0 )
1688 isFlagCmd_ = False ;
1689 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
1690 if ( isFreqOffset_ )
1691 if ( mstable_.freqOffset().nrow() == 0 )
1692 isFreqOffset_ = False ;
1693 isHistory_ = msrec.isDefined( "HISTORY" ) ;
1694 if ( isHistory_ )
1695 if ( mstable_.history().nrow() == 0 )
1696 isHistory_ = False ;
1697 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
1698 if ( isProcessor_ )
1699 if ( mstable_.processor().nrow() == 0 )
1700 isProcessor_ = False ;
1701 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
1702 if ( isSysCal_ )
1703 if ( mstable_.sysCal().nrow() == 0 )
1704 isSysCal_ = False ;
1705 isWeather_ = msrec.isDefined( "WEATHER" ) ;
1706 if ( isWeather_ )
1707 if ( mstable_.weather().nrow() == 0 )
1708 isWeather_ = False ;
1709
1710 // column name for Tsys and Tcal
1711 if ( isSysCal_ ) {
1712 const MSSysCal &caltab = mstable_.sysCal() ;
1713 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
1714 colTcal_ = "TCAL" ;
1715 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1716 colTcal_ = "NONE" ;
1717 }
1718 if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
1719 colTsys_ = "TSYS" ;
1720 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1721 colTsys_ = "NONE" ;
1722 }
1723 }
1724 else {
1725 colTcal_ = "NONE" ;
1726 colTsys_ = "NONE" ;
1727 }
1728
1729 // Access to MS subtables
1730 //MSField &fieldtab = mstable_.field() ;
1731 //MSPolarization &poltab = mstable_.polarization() ;
1732 //MSDataDescription &ddtab = mstable_.dataDescription() ;
1733 //MSObservation &obstab = mstable_.observation() ;
1734 //MSSource &srctab = mstable_.source() ;
1735 //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
1736 //MSSysCal &caltab = mstable_.sysCal() ;
1737 MSPointing &pointtab = mstable_.pointing() ;
1738 //MSState &stattab = mstable_.state() ;
1739 //MSAntenna &anttab = mstable_.antenna() ;
1740
1741 // SUBTABLES: FREQUENCIES
1742 //string freqFrame = getFrame() ;
1743 string baseFrame = frameFromSpwTable() ;
1744 table_->frequencies().setFrame( baseFrame ) ;
1745 table_->frequencies().setFrame( baseFrame, True ) ;
1746
1747 // SUBTABLES: WEATHER
1748 fillWeather() ;
1749
1750 // SUBTABLES: FOCUS
1751 fillFocus() ;
1752
1753 // SUBTABLES: TCAL
1754 fillTcal() ;
1755
1756 // SUBTABLES: FIT
1757 //fillFit() ;
1758
1759 // SUBTABLES: HISTORY
1760 //fillHistory() ;
1761
1762 /***
1763 * Start iteration using TableVisitor
1764 ***/
1765 Table stab = table_->table() ;
1766 {
1767 static const char *cols[] = {
1768 "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
1769 "STATE_ID", "TIME",
1770 NULL
1771 };
1772 static const TypeManagerImpl<Int> tmInt;
1773 static const TypeManagerImpl<Double> tmDouble;
1774 static const TypeManager *const tms[] = {
1775 &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
1776 };
1777 //double t0 = mathutil::gettimeofday_sec() ;
1778 MSFillerVisitor myVisitor(mstable_, *table_ );
1779 //double t1 = mathutil::gettimeofday_sec() ;
1780 //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1781 myVisitor.setAntenna( antenna_ ) ;
1782 //myVisitor.setHeader( sdh ) ;
1783 if ( getPt_ ) {
1784 Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
1785 myVisitor.setPointingTable( ptsel ) ;
1786 }
1787 if ( isWeather_ )
1788 myVisitor.setWeatherTime( mwTime_, mwInterval_, mwIndex_ ) ;
1789 if ( isSysCal_ )
1790 myVisitor.setSysCalRecord( tcalrec_ ) ;
1791
1792 //double t2 = mathutil::gettimeofday_sec() ;
1793 traverseTable(mstable_, cols, tms, &myVisitor);
1794 //double t3 = mathutil::gettimeofday_sec() ;
1795 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1796
1797 sdh = myVisitor.getHeader() ;
1798 }
1799 /***
1800 * End iteration using TableVisitor
1801 ***/
1802
1803 // set header
1804 //sdh = myVisitor.getHeader() ;
1805 //table_->setHeader( sdh ) ;
1806
1807 // save path to POINTING table
1808 // 2011/07/06 TN
1809 // Path to POINTING table in original MS will not be written
1810 // if getPt_ is True
1811 Path datapath( tablename_ ) ;
1812 if ( !getPt_ ) {
1813 String pTabName = datapath.absoluteName() + "/POINTING" ;
1814 stab.rwKeywordSet().define( "POINTING", pTabName ) ;
1815 }
1816
1817 // for GBT
1818 if ( sdh.antennaname.contains( "GBT" ) ) {
1819 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
1820 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
1821 }
1822
1823 // for MS created from ASDM
1824 const TableRecord &msKeys = mstable_.keywordSet() ;
1825 uInt nfields = msKeys.nfields() ;
1826 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1827 String name = msKeys.name( ifield ) ;
1828 //os_ << "name = " << name << LogIO::POST ;
1829 if ( name.find( "ASDM" ) != String::npos ) {
1830 String asdmpath = msKeys.asTable( ifield ).tableName() ;
1831 os_ << "ASDM table: " << asdmpath << LogIO::POST ;
1832 stab.rwKeywordSet().define( name, asdmpath ) ;
1833 }
1834 }
1835
1836 //double endSec = mathutil::gettimeofday_sec() ;
1837 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1838}
1839
1840void MSFiller::close()
1841{
1842 //tablesel_.closeSubTables() ;
1843 mstable_.closeSubTables() ;
1844 //tablesel_.unlock() ;
1845 mstable_.unlock() ;
1846}
1847
1848void MSFiller::fillWeather()
1849{
1850 //double startSec = mathutil::gettimeofday_sec() ;
1851 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1852
1853 if ( !isWeather_ ) {
1854 // add dummy row
1855 table_->weather().table().addRow(1,True) ;
1856 return ;
1857 }
1858
1859 Table mWeather = mstable_.weather() ;
1860 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1861 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1862 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1863 if ( mWeatherSel.nrow() == 0 ) {
1864 os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1865 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1866 if ( mWeatherSel.nrow() == 0 ) {
1867 os_ << "No rows in WEATHER table" << LogIO::POST ;
1868 }
1869 }
1870 uInt wnrow = mWeatherSel.nrow() ;
1871 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1872
1873 if ( wnrow == 0 )
1874 return ;
1875
1876 Table wtab = table_->weather().table() ;
1877 wtab.addRow( wnrow ) ;
1878
1879 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
1880 Int stationId = -1 ;
1881 if ( stationInfoExists ) {
1882 // determine which station is closer
1883 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1884 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
1885 Vector<Int> stationIds = stationCol.getColumn() ;
1886 Vector<Int> stationIdList( 0 ) ;
1887 Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
1888 uInt numStation = 0 ;
1889 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
1890 if ( !anyEQ( stationIdList, stationIds[i] ) ) {
1891 numStation++ ;
1892 stationIdList.resize( numStation, True ) ;
1893 stationIdList[numStation-1] = stationIds[i] ;
1894 stationPosList.resize( numStation, 3, True ) ;
1895 stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
1896 }
1897 }
1898 //os_ << "staionIdList = " << stationIdList << endl ;
1899 Table mAntenna = mstable_.antenna() ;
1900 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
1901 Vector<Double> antpos = antposCol( antenna_ ) ;
1902 Double minDiff = -1.0 ;
1903 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
1904 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
1905 if ( minDiff < 0.0 || minDiff > diff ) {
1906 minDiff = diff ;
1907 stationId = stationIdList[i] ;
1908 }
1909 }
1910 }
1911 //os_ << "stationId = " << stationId << endl ;
1912
1913 ScalarColumn<Float> *fCol ;
1914 ROScalarColumn<Float> *sharedFloatCol ;
1915 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1916 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1917 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1918 fCol->putColumn( *sharedFloatCol ) ;
1919 delete sharedFloatCol ;
1920 delete fCol ;
1921 }
1922 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1923 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1924 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1925 fCol->putColumn( *sharedFloatCol ) ;
1926 delete sharedFloatCol ;
1927 delete fCol ;
1928 }
1929 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1930 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1931 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1932 fCol->putColumn( *sharedFloatCol ) ;
1933 delete sharedFloatCol ;
1934 delete fCol ;
1935 }
1936 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1937 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1938 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1939 fCol->putColumn( *sharedFloatCol ) ;
1940 delete sharedFloatCol ;
1941 delete fCol ;
1942 }
1943 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1944 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1945 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1946 fCol->putColumn( *sharedFloatCol ) ;
1947 delete sharedFloatCol ;
1948 delete fCol ;
1949 }
1950 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1951 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1952 idCol.put( irow, irow ) ;
1953
1954 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1955 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1956 String tUnit = tqCol.getUnits() ;
1957 Vector<Double> mwTime = tCol.getColumn() ;
1958 if ( tUnit == "d" )
1959 mwTime *= 86400.0 ;
1960 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1961 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1962 String iUnit = tqCol.getUnits() ;
1963 Vector<Double> mwInterval = tCol.getColumn() ;
1964 if ( iUnit == "d" )
1965 mwInterval *= 86400.0 ;
1966
1967 if ( stationId > 0 ) {
1968 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1969 Vector<Int> stationVec = stationCol.getColumn() ;
1970 uInt wsnrow = ntrue( stationVec == stationId ) ;
1971 mwTime_.resize( wsnrow ) ;
1972 mwInterval_.resize( wsnrow ) ;
1973 mwIndex_.resize( wsnrow ) ;
1974 uInt wsidx = 0 ;
1975 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
1976 if ( stationId == stationVec[irow] ) {
1977 mwTime_[wsidx] = mwTime[irow] ;
1978 mwInterval_[wsidx] = mwInterval[irow] ;
1979 mwIndex_[wsidx] = irow ;
1980 wsidx++ ;
1981 }
1982 }
1983 }
1984 else {
1985 mwTime_ = mwTime ;
1986 mwInterval_ = mwInterval ;
1987 mwIndex_.resize( mwTime_.size() ) ;
1988 indgen( mwIndex_ ) ;
1989 }
1990 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1991 //os_ << "mwIndex_=" << mwIndex_ << LogIO::POST;
1992 //double endSec = mathutil::gettimeofday_sec() ;
1993 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1994}
1995
1996void MSFiller::fillFocus()
1997{
1998 //double startSec = mathutil::gettimeofday_sec() ;
1999 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
2000 // tentative
2001 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
2002 //double endSec = mathutil::gettimeofday_sec() ;
2003 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2004}
2005
2006void MSFiller::fillTcal()
2007{
2008 //double startSec = mathutil::gettimeofday_sec() ;
2009 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
2010
2011 if ( !isSysCal_ ) {
2012 // add dummy row
2013 os_ << "No SYSCAL rows" << LogIO::POST ;
2014 table_->tcal().table().addRow(1,True) ;
2015 Vector<Float> defaultTcal( 1, 1.0 ) ;
2016 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2017 tcalCol.put( 0, defaultTcal ) ;
2018 return ;
2019 }
2020
2021 if ( colTcal_ == "NONE" ) {
2022 // add dummy row
2023 os_ << "No TCAL column" << LogIO::POST ;
2024 table_->tcal().table().addRow(1,True) ;
2025 Vector<Float> defaultTcal( 1, 1.0 ) ;
2026 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2027 tcalCol.put( 0, defaultTcal ) ;
2028 return ;
2029 }
2030
2031 Table &sctab = mstable_.sysCal() ;
2032 if ( sctab.nrow() == 0 ) {
2033 os_ << "No SYSCAL rows" << LogIO::POST ;
2034 return ;
2035 }
2036 ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
2037 Vector<Int> ant = antCol.getColumn() ;
2038 if ( allNE( ant, antenna_ ) ) {
2039 os_ << "No SYSCAL rows" << LogIO::POST ;
2040 return ;
2041 }
2042 ROTableColumn tcalCol( sctab, colTcal_ ) ;
2043 Bool notDefined = False ;
2044 for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
2045 if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
2046 notDefined = True ;
2047 break ;
2048 }
2049 }
2050 if ( notDefined ) {
2051 os_ << "No TCAL value" << LogIO::POST ;
2052 table_->tcal().table().addRow(1,True) ;
2053 Vector<Float> defaultTcal( 1, 1.0 ) ;
2054 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2055 tcalCol.put( 0, defaultTcal ) ;
2056 return ;
2057 }
2058
2059 static const char *cols[] = {
2060 "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
2061 NULL
2062 };
2063 static const TypeManagerImpl<Int> tmInt;
2064 static const TypeManagerImpl<Double> tmDouble;
2065 static const TypeManager *const tms[] = {
2066 &tmInt, &tmInt, &tmInt, &tmDouble, NULL
2067 };
2068 Table tab = table_->tcal().table() ;
2069 TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
2070 visitor.setTcalColumn( colTcal_ ) ;
2071
2072 traverseTable(sctab, cols, tms, &visitor);
2073
2074 infillTcal();
2075
2076 //tcalrec_.print( std::cout ) ;
2077 //double endSec = mathutil::gettimeofday_sec() ;
2078 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2079}
2080
2081void MSFiller::infillTcal()
2082{
2083 uInt nfields = tcalrec_.nfields() ;
2084 set<Int> spwAvailable;
2085 for (uInt i = 0; i < nfields; i++) {
2086 String name = tcalrec_.name(i);
2087 size_t pos1 = name.find(':') + 4;
2088 size_t pos2 = name.find(':',pos1);
2089 Int spwid = String::toInt(name.substr(pos1,pos2-pos1));
2090 //cout << "spwid=" << spwid << endl;
2091 spwAvailable.insert(spwid);
2092 }
2093 Table spwtab = mstable_.spectralWindow();
2094 Table tcaltab = table_->tcal().table();
2095 ScalarColumn<uInt> idCol(tcaltab, "ID");
2096 ScalarColumn<String> timeCol(tcaltab, "TIME");
2097 ArrayColumn<Float> tcalCol(tcaltab, "TCAL");
2098 ROScalarColumn<Int> numChanCol(spwtab, "NUM_CHAN");
2099 Int numSpw = spwtab.nrow();
2100 Int dummyFeed = 0;
2101 Double dummyTime = 0.0;
2102 Vector<uInt> idminmax(2);
2103 for (Int i = 0; i < numSpw; i++) {
2104 if (spwAvailable.find(i) == spwAvailable.end()) {
2105 String key = keyTcal(dummyFeed, i, dummyTime);
2106 Vector<Float> tcal(numChanCol(i), 1.0);
2107 uInt nrow = tcaltab.nrow();
2108 tcaltab.addRow(1);
2109 idCol.put(nrow, nrow);
2110 timeCol.put(nrow, "");
2111 tcalCol.put(nrow, tcal);
2112 idminmax = nrow;
2113 tcalrec_.define(key, idminmax);
2114 }
2115 }
2116 //tcalrec_.print(cout);
2117}
2118
2119string MSFiller::getFrame()
2120{
2121 MFrequency::Types frame = MFrequency::DEFAULT ;
2122 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
2123 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
2124 uInt nrow = numChanCol.nrow() ;
2125 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
2126 uInt nref = 0 ;
2127 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2128 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
2129 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
2130 nref++ ;
2131 }
2132 }
2133 if ( nref > 0 )
2134 frame = (MFrequency::Types)measFreqRef[0] ;
2135
2136 return MFrequency::showType( frame ) ;
2137}
2138
2139void MSFiller::initHeader( STHeader &header )
2140{
2141 header.nchan = 0 ;
2142 header.npol = 0 ;
2143 header.nif = 0 ;
2144 header.nbeam = 0 ;
2145 header.observer = "" ;
2146 header.project = "" ;
2147 header.obstype = "" ;
2148 header.antennaname = "" ;
2149 header.antennaposition.resize( 3 ) ;
2150 header.equinox = 0.0 ;
2151 header.freqref = "" ;
2152 header.reffreq = -1.0 ;
2153 header.bandwidth = 0.0 ;
2154 header.utc = 0.0 ;
2155 header.fluxunit = "" ;
2156 header.epoch = "" ;
2157 header.poltype = "" ;
2158}
2159
2160string MSFiller::frameFromSpwTable()
2161{
2162 string frameString;
2163 Table tab = mstable_.spectralWindow();
2164 ROScalarColumn<Int> mfrCol(tab, "MEAS_FREQ_REF");
2165 Vector<Int> mfr = mfrCol.getColumn();
2166 if (allEQ(mfr,mfr[0])) {
2167 frameString = MFrequency::showType(mfr[0]);
2168 //cout << "all rows have same frame: " << frameString << endl;
2169 }
2170 else {
2171 mfrCol.attach(tab, "NUM_CHAN");
2172 for (uInt i = 0; i < tab.nrow(); i++) {
2173 if (mfrCol(i) != 4) {
2174 frameString = MFrequency::showType(mfr[i]);
2175 break;
2176 }
2177 }
2178 if (frameString.size() == 0) {
2179 frameString = "TOPO";
2180 }
2181 }
2182
2183 //cout << "frameString = " << frameString << endl;
2184
2185 return frameString;
2186}
2187
2188};
Note: See TracBrowser for help on using the repository browser.