source: trunk/src/MSFiller.cpp@ 2784

Last change on this file since 2784 was 2775, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on indexing of WEATHER_ID when more than one weather station
is included in input MS.


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