source: trunk/src/MSFiller.cpp@ 2748

Last change on this file since 2748 was 2746, 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...

Fixed a known bug that SYSCAL table is not properly handled when
import/export.


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