source: trunk/src/MSFiller.cpp@ 2765

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

New Development: No

JIRA Issue: Yes CSV-2532 (may be related to CSV-1908 and CSV-2161)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: In tool level, added parameter 'freq_tolsr' to

scantable constructor and function sd.splitant.

Test Programs: test_sdsave, test_importasdm_sd

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

In importing MS to Scantable, frequency frame information is
imported as is by default, i.e., base frame in Scantable is
TOPO for ALMA data, which is forcibly converted to LSRK with
wrong time and direction reference.

Some functions have a boolean parameter 'freq_tolsr' that controls
the above behavior. If freq_tolsr is False (default), frequency
is imported as is, while frequency is converted to LSRK (wrongly)
when it is True.


File size: 69.2 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 {
689 isWeather = True ;
690 weatherTime = t ;
691 weatherInterval = it ;
692 }
693 void setSysCalRecord( const Record &r )
694 //void setSysCalRecord( const map< String,Vector<uInt> > &r )
695 {
696 isSysCal = True ;
697 isTcal = True ;
698 syscalRecord = r ;
699 if ( syscalRecord.nfields() == 0 )
700 isTcal = False ;
701
702 const TableDesc &desc = sctab.tableDesc() ;
703 uInt nrow = sctab.nrow() ;
704 syscalRow.resize( nrow ) ;
705 syscalTime.resize( nrow ) ;
706 syscalInterval.resize( nrow ) ;
707 String tsysCol = "NONE" ;
708 Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ;
709 for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) {
710 if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) )
711 tsysCol = tsysCols[i] ;
712 }
713 sysCalTsysCol.attach( sctab, tsysCol ) ;
714 }
715 void setFreqToLsr(Bool b) { freqToLsr_ = b; }
716 STHeader getHeader() { return header ; }
717 uInt getNumBeam() { return nbeam ; }
718 uInt getFilledRowNum() { return rowidx ; }
719private:
720 void initConvert()
721 {
722 toj2000 = MDirection::Convert( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
723 toazel = MDirection::Convert( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
724 }
725
726 void fluxUnit( String &u )
727 {
728 ROTableColumn col( table, dataColumnName ) ;
729 const TableRecord &rec = col.keywordSet() ;
730 if ( rec.isDefined( "UNIT" ) )
731 u = rec.asString( "UNIT" ) ;
732 else if ( rec.isDefined( "QuantumUnits" ) )
733 u = rec.asString( "QuantumUnits" ) ;
734 if ( u.empty() )
735 u = "K" ;
736 }
737 void processSource( Int sourceId, Int spwId,
738 String &name, MDirection &dir, Vector<Double> &pm,
739 Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel )
740 {
741 // find row
742 uInt nrow = srctab.nrow() ;
743 Int idx = -1 ;
744 ROTableRow row( srctab ) ;
745 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
746 const TableRecord &r = row.get( irow ) ;
747 if ( r.asInt( "SOURCE_ID" ) == sourceId ) {
748 Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ;
749 if ( tmpSpwId == spwId || tmpSpwId == -1 ) {
750 idx = (Int)irow ;
751 break ;
752 }
753 }
754 }
755
756 // fill
757 Int numLines = 0 ;
758 if ( idx != -1 ) {
759 const TableRecord &r = row.get( idx ) ;
760 name = r.asString( "NAME" ) ;
761 getScalarMeas( "DIRECTION", idx, srctab, dir ) ;
762 pm = r.toArrayDouble( "PROPER_MOTION" ) ;
763 numLines = r.asInt( "NUM_LINES" ) ;
764 }
765 else {
766 name = "" ;
767 pm = Vector<Double>( 2, 0.0 ) ;
768 dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
769 }
770 if ( !getpt ) {
771 String ref = dir.getRefString() ;
772 MDirection::getType( dirType, ref ) ;
773
774 // initialize toj2000 and toazel
775 initConvert() ;
776 }
777
778 rf.resize( numLines ) ;
779 trans.resize( numLines ) ;
780 vel.resize( numLines ) ;
781 if ( numLines > 0 ) {
782 Block<Bool> isDefined = row.getDefined() ;
783 Vector<String> colNames = row.columnNames() ;
784 Vector<Int> indexes( 3, -1 ) ;
785 Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ;
786 for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) {
787 if ( anyEQ( indexes, -1 ) ) {
788 for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) {
789 if ( colNames[icol] == cols[jcol] )
790 indexes[jcol] = icol ;
791 }
792 }
793 }
794 if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) {
795 Vector< Quantum<Double> > qrf ;
796 getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ;
797 for ( int i = 0 ; i < numLines ; i++ )
798 rf[i] = qrf[i].getValue( "Hz" ) ;
799 }
800 if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) {
801 getArray( "TRANSITION", idx, srctab, trans ) ;
802 }
803 if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) {
804 Vector< Quantum<Double> > qsv ;
805 getArrayQuant( "SYSVEL", idx, srctab, qsv ) ;
806 for ( int i = 0 ; i < numLines ; i++ )
807 vel[i] = qsv[i].getValue( "m/s" ) ;
808 }
809 }
810 }
811 void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md,
812 uInt &freqId, Int &nchan,
813 String &freqref, Double &reffreq, Double &bandwidth )
814 {
815 // fill
816 Int measFreqRef ;
817 getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ;
818 MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ;
819 //freqref = MFrequency::showType( freqRef ) ;
820 //freqref = "LSRK" ;
821 freqref = "TOPO";
822 Quantum<Double> q ;
823 getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ;
824 bandwidth = q.getValue( "Hz" ) ;
825 getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ;
826 reffreq = q.getValue( "Hz" ) ;
827 Double refpix = 0.5 * ( (Double)nchan-1.0 ) ;
828 Int refchan = ( nchan - 1 ) / 2 ;
829 Bool even = (Bool)( nchan % 2 == 0 ) ;
830 Vector< Quantum<Double> > qa ;
831 getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ;
832// Double increment = qa[refchan].getValue( "Hz" ) ;
833 Double increment = abs(qa[refchan].getValue( "Hz" )) ;
834 getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ;
835 if ( nchan == 1 ) {
836 Int netSideband ;
837 getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ;
838 if ( netSideband == 1 ) increment *= -1.0 ;
839 }
840 else {
841 if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) )
842 increment *= -1.0 ;
843 }
844 Double refval = qa[refchan].getValue( "Hz" ) ;
845 if ( even )
846 refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ;
847 if ( freqToLsr_ && freqRef != MFrequency::LSRK ) {
848 //cout << "do conversion to LSRK" << endl;
849 MeasFrame mframe( me, mp, md ) ;
850 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mframe ) ) ;
851 refval = tolsr( Quantum<Double>( refval, "Hz" ) ).get( "Hz" ).getValue() ;
852 }
853
854 // add new row to FREQUENCIES
855 Table ftab = scantable.frequencies().table() ;
856 freqId = ftab.nrow() ;
857 ftab.addRow() ;
858 TableRow row( ftab ) ;
859 TableRecord &r = row.record() ;
860 RecordFieldPtr<uInt> idRF( r, "ID" ) ;
861 *idRF = freqId ;
862 RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ;
863 RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ;
864 RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ;
865 *refpixRF = refpix ;
866 *refvalRF = refval ;
867 *incrRF = increment ;
868 row.put( freqId ) ;
869 }
870 void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl )
871 {
872 Matrix<Bool> b = flagCol( recordNo ) ;
873 if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) {
874 sp = floatDataCol( recordNo ) ;
875 convertArray( fl, b ) ;
876 }
877 else {
878 Bool notyet = True ;
879 Matrix<Complex> c = dataCol( recordNo ) ;
880 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
881 if ( ( header.poltype == "linear" || header.poltype == "circular" )
882 && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) {
883 if ( notyet ) {
884 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
885 IPosition start( 1, 0 ) ;
886 IPosition end( 1, 2*nchan-1 ) ;
887 IPosition inc( 1, 2 ) ;
888 if ( polnos[ipol] == 2 ) {
889 sp.row( ipol ) = tmp( start, end, inc ) ;
890 Vector<Bool> br = b.row( ipol ) ;
891 Vector<uChar> flr = fl.row( ipol ) ;
892 convertArray( flr, br ) ;
893 start = IPosition( 1, 1 ) ;
894 Int jpol = ipol+1 ;
895 while( polnos[jpol] != 3 && jpol < npol )
896 jpol++ ;
897 sp.row( jpol ) = tmp( start, end, inc ) ;
898 flr.reference( fl.row( jpol ) ) ;
899 convertArray( flr, br ) ;
900 }
901 else if ( polnos[ipol] == 3 ) {
902 sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ;
903 Int jpol = ipol+1 ;
904 while( polnos[jpol] != 2 && jpol < npol )
905 jpol++ ;
906 Vector<Bool> br = b.row( ipol ) ;
907 Vector<uChar> flr = fl.row( jpol ) ;
908 sp.row( jpol ) = tmp( start, end, inc ) ;
909 convertArray( flr, br ) ;
910 start = IPosition( 1, 1 ) ;
911 sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ;
912 flr.reference( fl.row( ipol ) ) ;
913 convertArray( flr, br ) ;
914 }
915 notyet = False ;
916 }
917 }
918 else {
919 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
920 IPosition start( 1, 0 ) ;
921 IPosition end( 1, 2*nchan-1 ) ;
922 IPosition inc( 1, 2 ) ;
923 sp.row( ipol ) = tmp( start, end, inc ) ;
924 Vector<Bool> br = b.row( ipol ) ;
925 Vector<uChar> flr = fl.row( ipol ) ;
926 convertArray( flr, br ) ;
927 }
928 }
929 }
930 }
931 uInt binarySearch( Vector<Double> &timeList, Double target )
932 {
933 Int low = 0 ;
934 Int high = timeList.nelements() ;
935 uInt idx = 0 ;
936
937 while ( low <= high ) {
938 idx = (Int)( 0.5 * ( low + high ) ) ;
939 Double t = timeList[idx] ;
940 if ( t < target )
941 low = idx + 1 ;
942 else if ( t > target )
943 high = idx - 1 ;
944 else {
945 return idx ;
946 }
947 }
948
949 idx = max( 0, min( low, high ) ) ;
950 return idx ;
951 }
952 void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
953 {
954 // @todo At the moment, do binary search every time
955 // if this is bottleneck, frequency of binary search must be reduced
956 Double t = currentTime.get( "s" ).getValue() ;
957 uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ;
958 Matrix<Double> d ;
959 if ( pointingTime[idx] == t )
960 d = pointingDirection.xyPlane( idx ) ;
961 else if ( pointingTime[idx] < t ) {
962 if ( idx == pointingTime.nelements()-1 )
963 d = pointingDirection.xyPlane( idx ) ;
964 else
965 d = interp( pointingTime[idx], pointingTime[idx+1], t,
966 pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ;
967 }
968 else {
969 if ( idx == 0 )
970 d = pointingDirection.xyPlane( idx ) ;
971 else
972 d = interp( pointingTime[idx-1], pointingTime[idx], t,
973 pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
974 }
975 mf.set( currentTime ) ;
976 Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
977 if ( dirType != MDirection::J2000 ) {
978 dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
979 }
980 else {
981 dir = d.column( 0 ) ;
982 }
983 if ( dirType != MDirection::AZELGEO ) {
984 azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
985 }
986 else {
987 azel = d.column( 0 ) ;
988 }
989 if ( d.ncolumn() > 1 )
990 srate = d.column( 1 ) ;
991 }
992 void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &/*srate*/ )
993 {
994 dir = sourceDir.getAngle( "rad" ).getValue() ;
995 mf.set( currentTime ) ;
996 azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
997 if ( dirType != MDirection::J2000 ) {
998 dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
999 }
1000 }
1001 String detectSeparator( String &s )
1002 {
1003 String tmp = s.substr( 0, s.find_first_of( "," ) ) ;
1004 Char *separators[] = { ":", "#", ".", "_" } ;
1005 uInt nsep = 4 ;
1006 for ( uInt i = 0 ; i < nsep ; i++ ) {
1007 if ( tmp.find( separators[i] ) != String::npos )
1008 return separators[i] ;
1009 }
1010 return "" ;
1011 }
1012 Int getSrcType( Int stateId )
1013 {
1014 // get values
1015 Bool sig ;
1016 getScalar( "SIG", stateId, statetab, sig ) ;
1017 Bool ref ;
1018 getScalar( "REF", stateId, statetab, ref ) ;
1019 Double cal ;
1020 getScalar( "CAL", stateId, statetab, cal ) ;
1021 String obsmode ;
1022 getScalar( "OBS_MODE", stateId, statetab, obsmode ) ;
1023 String sep = detectSeparator( obsmode ) ;
1024
1025 Int srcType = SrcType::NOTYPE ;
1026 if ( sep == ":" )
1027 srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ;
1028 else if ( sep == "." || sep == "#" )
1029 srcTypeALMA( srcType, sep, obsmode ) ;
1030 else if ( sep == "_" )
1031 srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ;
1032 else
1033 srcTypeDefault( srcType, sig, ref ) ;
1034
1035 return srcType ;
1036 }
1037 void srcTypeDefault( Int &st, Bool &sig, Bool &ref )
1038 {
1039 if ( sig ) st = SrcType::SIG ;
1040 else if ( ref ) st = SrcType::REF ;
1041 }
1042 void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal )
1043 {
1044 Int epos = mode.find_first_of( sep ) ;
1045 Int nextpos = mode.find_first_of( sep, epos+1 ) ;
1046 String m1 = mode.substr( 0, epos ) ;
1047 String m2 = mode.substr( epos+1, nextpos-epos-1 ) ;
1048 if ( m1 == "Nod" ) {
1049 st = SrcType::NOD ;
1050 }
1051 else if ( m1 == "OffOn" ) {
1052 if ( m2 == "PSWITCHON" ) st = SrcType::PSON ;
1053 if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ;
1054 }
1055 else {
1056 if ( m2 == "FSWITCH" ) {
1057 if ( sig ) st = SrcType::FSON ;
1058 else if ( ref ) st = SrcType::FSOFF ;
1059 }
1060 }
1061 if ( cal > 0.0 ) {
1062 if ( st == SrcType::NOD )
1063 st = SrcType::NODCAL ;
1064 else if ( st == SrcType::PSON )
1065 st = SrcType::PONCAL ;
1066 else if ( st == SrcType::PSOFF )
1067 st = SrcType::POFFCAL ;
1068 else if ( st == SrcType::FSON )
1069 st = SrcType::FONCAL ;
1070 else if ( st == SrcType::FSOFF )
1071 st = SrcType::FOFFCAL ;
1072 else
1073 st = SrcType::CAL ;
1074 }
1075 }
1076 void srcTypeALMA( Int &st, String &sep, String &mode )
1077 {
1078 Int epos = mode.find_first_of( "," ) ;
1079 String first = mode.substr( 0, epos ) ;
1080 epos = first.find_first_of( sep ) ;
1081 Int nextpos = first.find_first_of( sep, epos+1 ) ;
1082 String m1 = first.substr( 0, epos ) ;
1083 String m2 = first.substr( epos+1, nextpos-epos-1 ) ;
1084 if ( m1.find( "CALIBRATE_" ) == 0 ) {
1085 if ( m2.find( "ON_SOURCE" ) == 0 )
1086 st = SrcType::PONCAL ;
1087 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1088 st = SrcType::POFFCAL ;
1089 }
1090 else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) {
1091 if ( m2.find( "ON_SOURCE" ) == 0 )
1092 st = SrcType::PSON ;
1093 else if ( m2.find( "OFF_SOURCE" ) == 0 )
1094 st = SrcType::PSOFF ;
1095 }
1096 }
1097 void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref )
1098 {
1099 Int epos = mode.find_first_of( "," ) ;
1100 String first = mode.substr( 0, epos ) ;
1101 string substr[4] ;
1102 int numSubstr = split( first, substr, 4, sep ) ;
1103 String m1( substr[0] ) ;
1104 String m2( substr[2] ) ;
1105 if ( numSubstr == 4 ) {
1106 if ( m1.find( "CALIBRATE" ) == 0 ) {
1107 if ( m2.find( "ON" ) == 0 )
1108 st = SrcType::PONCAL ;
1109 else if ( m2.find( "OFF" ) == 0 )
1110 st = SrcType::POFFCAL ;
1111 }
1112 else if ( m1.find( "OBSERVE" ) == 0 ) {
1113 if ( m2.find( "ON" ) == 0 )
1114 st = SrcType::PSON ;
1115 else if ( m2.find( "OFF" ) == 0 )
1116 st = SrcType::PSOFF ;
1117 }
1118 }
1119 else {
1120 if ( sig ) st = SrcType::SIG ;
1121 else if ( ref ) st = SrcType::REF ;
1122 }
1123 }
1124 Block<uInt> getPolNos( Vector<Int> &corr )
1125 {
1126 Block<uInt> polnos( npol ) ;
1127 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
1128 if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX )
1129 polnos[ipol] = 0 ;
1130 else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY )
1131 polnos[ipol] = 1 ;
1132 else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY )
1133 polnos[ipol] = 2 ;
1134 else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX )
1135 polnos[ipol] = 3 ;
1136 }
1137 return polnos ;
1138 }
1139 String getPolType( Int &corr )
1140 {
1141 String poltype = "" ;
1142 if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V )
1143 poltype = "stokes" ;
1144 else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX )
1145 poltype = "linear" ;
1146 else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR )
1147 poltype = "circular" ;
1148 else if ( corr == Stokes::Plinear || corr == Stokes::Pangle )
1149 poltype = "linpol" ;
1150 return poltype ;
1151 }
1152 uInt getWeatherId()
1153 {
1154 // if only one row, return 0
1155 if ( weatherTime.nelements() == 1 )
1156 return 0 ;
1157
1158 // @todo At the moment, do binary search every time
1159 // if this is bottleneck, frequency of binary search must be reduced
1160 Double t = currentTime.get( "s" ).getValue() ;
1161 uInt idx = min( binarySearch( weatherTime, t ), weatherTime.nelements()-1 ) ;
1162 if ( weatherTime[idx] < t ) {
1163 if ( idx != weatherTime.nelements()-1 ) {
1164 if ( weatherTime[idx+1] - t < 0.5 * weatherInterval[idx+1] )
1165 idx++ ;
1166 }
1167 }
1168 else if ( weatherTime[idx] > t ) {
1169 if ( idx != 0 ) {
1170 if ( weatherTime[idx] - t > 0.5 * weatherInterval[idx] )
1171 idx-- ;
1172 }
1173 }
1174 return idx ;
1175 }
1176 void processSysCal( Int &spwId )
1177 {
1178 // get feedId from row
1179 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1180
1181 uInt nrow = sctab.nrow() ;
1182 ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ;
1183 Vector<Int> aids = col.getColumn() ;
1184 col.attach( sctab, "FEED_ID" ) ;
1185 Vector<Int> fids = col.getColumn() ;
1186 col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ;
1187 Vector<Int> sids = col.getColumn() ;
1188 ROScalarColumn<Double> timeCol( sctab, "TIME" ) ;
1189 ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ;
1190 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1191 if ( aids[irow] == antennaId
1192 && fids[irow] == feedId
1193 && sids[irow] == spwId ) {
1194 syscalRow[numSysCalRow] = irow ;
1195 syscalTime[numSysCalRow] = timeCol( irow ) ;
1196 syscalInterval[numSysCalRow] = intCol( irow ) ;
1197 numSysCalRow++ ;
1198 }
1199 }
1200 }
1201 uInt getSysCalIndex()
1202 {
1203 // if only one row, return 0
1204 if ( numSysCalRow == 1 || !isSysCal )
1205 return 0 ;
1206
1207 // @todo At the moment, do binary search every time
1208 // if this is bottleneck, frequency of binary search must be reduced
1209 Double t = currentTime.get( "s" ).getValue() ;
1210 Vector<Double> tslice = syscalTime( Slice(0, numSysCalRow) ) ;
1211 uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ;
1212 if ( syscalTime[idx] < t ) {
1213 if ( idx != numSysCalRow-1 ) {
1214 if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] )
1215 idx++ ;
1216 }
1217 }
1218 else if ( syscalTime[idx] > t ) {
1219 if ( idx != 0 ) {
1220 if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] )
1221 idx-- ;
1222 }
1223 }
1224 return idx ;
1225 }
1226 Block<uInt> getTcalId( Double &t )
1227 {
1228 // return 0 if no SysCal table
1229 if ( !isSysCal or !isTcal ) {
1230 return Block<uInt>( 4, 0 ) ;
1231 }
1232
1233 // get feedId from row
1234 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
1235
1236 // key
1237 String key = keyTcal( feedId, spwId, t ) ;
1238
1239 // retrieve ids
1240 Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ;
1241 //Vector<uInt> ids = syscalRecord[key] ;
1242 uInt np = ids[1] - ids[0] + 1 ;
1243 Block<uInt> tcalids( np ) ;
1244 if ( np > 0 ) {
1245 tcalids[0] = ids[0] ;
1246 if ( np > 1 ) {
1247 tcalids[1] = ids[1] ;
1248 for ( uInt ip = 2 ; ip < np ; ip++ )
1249 tcalids[ip] = ids[0] + ip - 1 ;
1250 }
1251 }
1252 return tcalids ;
1253 }
1254 Block<uInt> getDummyTcalId( Int spwId )
1255 {
1256 Block<uInt> idList(4, 0);
1257 uInt nfields = syscalRecord.nfields();
1258 Int idx = -1;
1259 for (uInt i = 0; i< nfields ; i++ ) {
1260 String spw = "SPW" + String::toString(spwId);
1261 if (syscalRecord.name(i).find(spw) != String::npos) {
1262 idx = i;
1263 break;
1264 }
1265 }
1266 if ( idx > -1) {
1267 Vector<uInt> tmp = syscalRecord.asArrayuInt(idx);
1268 for (uInt j = 0 ; j < 4 ; j++) {
1269 idList[j] = tmp[0];
1270 }
1271 }
1272 return idList;
1273 }
1274 uInt maxNumPol()
1275 {
1276 ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ;
1277 return max( numCorrCol.getColumn() ) ;
1278 }
1279
1280 Scantable &scantable;
1281 Int antennaId;
1282 uInt rowidx;
1283 String dataColumnName;
1284 TableRow tablerow;
1285 STHeader header;
1286 Vector<Int> feedEntry;
1287 uInt nbeam;
1288 Int npol;
1289 Int nchan;
1290 Int sourceId;
1291 Int polId;
1292 Int spwId;
1293 uInt cycleNo;
1294 MDirection sourceDir;
1295 MPosition antpos;
1296 MEpoch currentTime;
1297 MEpoch obsEpoch;
1298 MeasFrame mf;
1299 MDirection::Convert toj2000;
1300 MDirection::Convert toazel;
1301 map<Int,uInt> ifmap;
1302 Block<uInt> polnos;
1303 Bool getpt;
1304 Vector<Double> pointingTime;
1305 Cube<Double> pointingDirection;
1306 MDirection::Types dirType;
1307 Bool isWeather;
1308 Vector<Double> weatherTime;
1309 Vector<Double> weatherInterval;
1310 Bool isSysCal;
1311 Bool isTcal;
1312 Record syscalRecord;
1313 //map< String,Vector<uInt> > syscalRecord;
1314 uInt numSysCalRow ;
1315 Vector<uInt> syscalRow;
1316 Vector<Double> syscalTime;
1317 Vector<Double> syscalInterval;
1318 //String tsysCol;
1319 //String tcalCol;
1320 Bool freqToLsr_;
1321
1322 // MS subtables
1323 Table obstab;
1324 Table sctab;
1325 Table spwtab;
1326 Table statetab;
1327 Table ddtab;
1328 Table poltab;
1329 Table fieldtab;
1330 Table anttab;
1331 Table srctab;
1332 Matrix<Float> sp;
1333 Matrix<uChar> fl;
1334
1335 // MS MAIN columns
1336 ROTableColumn intervalCol;
1337 ROTableColumn flagRowCol;
1338 ROArrayColumn<Float> floatDataCol;
1339 ROArrayColumn<Complex> dataCol;
1340 ROArrayColumn<Bool> flagCol;
1341
1342 // MS SYSCAL columns
1343 ROArrayColumn<Float> sysCalTsysCol;
1344
1345 // Scantable MAIN columns
1346 RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF;
1347 RecordFieldPtr< Vector<Double> > directionRF,scanRateRF,
1348 sourceProperMotionRF,sourceDirectionRF;
1349 RecordFieldPtr<Float> azimuthRF,elevationRF;
1350 RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF,
1351 ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF;
1352 RecordFieldPtr< Vector<Float> > spectraRF,tsysRF;
1353 RecordFieldPtr< Vector<uChar> > flagtraRF;
1354 RecordFieldPtr<String> sourceNameRF,fieldNameRF;
1355 RecordFieldPtr<Int> sourceTypeRF;
1356};
1357
1358class BaseTcalVisitor: public TableVisitor {
1359 uInt lastRecordNo ;
1360 Int lastAntennaId ;
1361 Int lastFeedId ;
1362 Int lastSpwId ;
1363 Double lastTime ;
1364protected:
1365 const Table &table;
1366 uInt count;
1367public:
1368 BaseTcalVisitor(const Table &table)
1369 : table(table)
1370 {
1371 count = 0;
1372 }
1373
1374 virtual void enterAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1375 virtual void leaveAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1376 virtual void enterFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1377 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1378 virtual void enterSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1379 virtual void leaveSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1380 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1381 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1382
1383 virtual Bool visitRecord(const uInt /*recordNo*/,
1384 const Int /*antennaId*/,
1385 const Int /*feedId*/,
1386 const Int /*spwId*/,
1387 const Double /*time*/) { return True ; }
1388
1389 virtual Bool visit(Bool isFirst, const uInt recordNo,
1390 const uInt nCols, void const *const colValues[]) {
1391 Int antennaId, feedId, spwId;
1392 Double time;
1393 { // prologue
1394 uInt i = 0;
1395 {
1396 const Int *col = (const Int *)colValues[i++];
1397 antennaId = col[recordNo];
1398 }
1399 {
1400 const Int *col = (const Int *)colValues[i++];
1401 feedId = col[recordNo];
1402 }
1403 {
1404 const Int *col = (const Int *)colValues[i++];
1405 spwId = col[recordNo];
1406 }
1407 {
1408 const Double *col = (const Double *)colValues[i++];
1409 time = col[recordNo];
1410 }
1411 assert(nCols == i);
1412 }
1413
1414 if (isFirst) {
1415 enterAntennaId(recordNo, antennaId);
1416 enterFeedId(recordNo, feedId);
1417 enterSpwId(recordNo, spwId);
1418 enterTime(recordNo, time);
1419 } else {
1420 if ( lastAntennaId != antennaId ) {
1421 leaveTime(lastRecordNo, lastTime);
1422 leaveSpwId(lastRecordNo, lastSpwId);
1423 leaveFeedId(lastRecordNo, lastFeedId);
1424 leaveAntennaId(lastRecordNo, lastAntennaId);
1425
1426 enterAntennaId(recordNo, antennaId);
1427 enterFeedId(recordNo, feedId);
1428 enterSpwId(recordNo, spwId);
1429 enterTime(recordNo, time);
1430 }
1431 else if (lastFeedId != feedId) {
1432 leaveTime(lastRecordNo, lastTime);
1433 leaveSpwId(lastRecordNo, lastSpwId);
1434 leaveFeedId(lastRecordNo, lastFeedId);
1435
1436 enterFeedId(recordNo, feedId);
1437 enterSpwId(recordNo, spwId);
1438 enterTime(recordNo, time);
1439 } else if (lastSpwId != spwId) {
1440 leaveTime(lastRecordNo, lastTime);
1441 leaveSpwId(lastRecordNo, lastSpwId);
1442
1443 enterSpwId(recordNo, spwId);
1444 enterTime(recordNo, time);
1445 } else if (lastTime != time) {
1446 leaveTime(lastRecordNo, lastTime);
1447 enterTime(recordNo, time);
1448 }
1449 }
1450 count++;
1451 Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time);
1452
1453 { // epilogue
1454 lastRecordNo = recordNo;
1455
1456 lastAntennaId = antennaId;
1457 lastFeedId = feedId;
1458 lastSpwId = spwId;
1459 lastTime = time;
1460 }
1461 return result ;
1462 }
1463
1464 virtual void finish() {
1465 if (count > 0) {
1466 leaveTime(lastRecordNo, lastTime);
1467 leaveSpwId(lastRecordNo, lastSpwId);
1468 leaveFeedId(lastRecordNo, lastFeedId);
1469 leaveAntennaId(lastRecordNo, lastAntennaId);
1470 }
1471 }
1472};
1473
1474class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
1475public:
1476 TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid )
1477 //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid )
1478 : BaseTcalVisitor( table ),
1479 tcal(tcaltab),
1480 rec(r),
1481 antenna(aid)
1482 {
1483 process = False ;
1484 rowidx = 0 ;
1485
1486 // attach to SYSCAL columns
1487 timeCol.attach( table, "TIME" ) ;
1488
1489 // add rows
1490 uInt addrow = table.nrow() * 4 ;
1491 tcal.addRow( addrow ) ;
1492
1493 // attach to TCAL columns
1494 row = TableRow( tcal ) ;
1495 TableRecord &trec = row.record() ;
1496 idRF.attachToRecord( trec, "ID" ) ;
1497 timeRF.attachToRecord( trec, "TIME" ) ;
1498 tcalRF.attachToRecord( trec, "TCAL" ) ;
1499 }
1500
1501 virtual void enterAntennaId(const uInt /*recordNo*/, Int columnValue) {
1502 if ( columnValue == antenna )
1503 process = True ;
1504 }
1505 virtual void leaveAntennaId(const uInt /*recordNo*/, Int /*columnValue*/) {
1506 process = False ;
1507 }
1508 virtual void enterFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1509 virtual void leaveFeedId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1510 virtual void enterSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1511 virtual void leaveSpwId(const uInt /*recordNo*/, Int /*columnValue*/) { }
1512 virtual void enterTime(const uInt recordNo, Double /*columnValue*/) {
1513 qtime = timeCol( recordNo ) ;
1514 }
1515 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1516 virtual Bool visitRecord(const uInt recordNo,
1517 const Int /*antennaId*/,
1518 const Int feedId,
1519 const Int spwId,
1520 const Double /*time*/)
1521 {
1522 //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ;
1523 if ( process ) {
1524 String sTime = MVTime( qtime ).string( MVTime::YMD ) ;
1525 *timeRF = sTime ;
1526 uInt oldidx = rowidx ;
1527 Matrix<Float> subtcal = tcalCol( recordNo ) ;
1528 Vector<uInt> idminmax( 2 ) ;
1529 for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) {
1530 *idRF = rowidx ;
1531 tcalRF.define( subtcal.row( ipol ) ) ;
1532
1533 // commit row
1534 row.put( rowidx ) ;
1535 rowidx++ ;
1536 }
1537
1538 idminmax[0] = oldidx ;
1539 idminmax[1] = rowidx - 1 ;
1540
1541 String key = keyTcal( feedId, spwId, sTime ) ;
1542 rec.define( key, idminmax ) ;
1543 //rec[key] = idminmax ;
1544 }
1545 return True ;
1546 }
1547 virtual void finish()
1548 {
1549 BaseTcalVisitor::finish() ;
1550
1551 if ( tcal.nrow() > rowidx ) {
1552 uInt numRemove = tcal.nrow() - rowidx ;
1553 //cout << "numRemove = " << numRemove << endl ;
1554 Vector<uInt> rows( numRemove ) ;
1555 indgen( rows, rowidx ) ;
1556 tcal.removeRow( rows ) ;
1557 }
1558
1559 }
1560 void setTcalColumn( String &col )
1561 {
1562 //colName = col ;
1563 tcalCol.attach( table, col ) ;
1564 }
1565private:
1566 Table &tcal;
1567 Record &rec;
1568 //map< String,Vector<uInt> > &rec;
1569 Int antenna;
1570 uInt rowidx;
1571 Bool process;
1572 Quantum<Double> qtime;
1573 TableRow row;
1574 String colName;
1575
1576 // MS SYSCAL columns
1577 ROScalarQuantColumn<Double> timeCol;
1578 ROArrayColumn<Float> tcalCol;
1579
1580 // TCAL columns
1581 RecordFieldPtr<uInt> idRF;
1582 RecordFieldPtr<String> timeRF;
1583 RecordFieldPtr< Vector<Float> > tcalRF;
1584};
1585
1586MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
1587 : table_( stable ),
1588 tablename_( "" ),
1589 antenna_( -1 ),
1590 antennaStr_(""),
1591 getPt_( True ),
1592 freqToLsr_( False ),
1593 isFloatData_( False ),
1594 isData_( False ),
1595 isDoppler_( False ),
1596 isFlagCmd_( False ),
1597 isFreqOffset_( False ),
1598 isHistory_( False ),
1599 isProcessor_( False ),
1600 isSysCal_( False ),
1601 isWeather_( False ),
1602 colTsys_( "TSYS_SPECTRUM" ),
1603 colTcal_( "TCAL_SPECTRUM" )
1604{
1605 os_ = LogIO() ;
1606 os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
1607}
1608
1609MSFiller::~MSFiller()
1610{
1611 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
1612}
1613
1614bool MSFiller::open( const std::string &filename, const casa::Record &rec )
1615{
1616 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
1617 //double startSec = mathutil::gettimeofday_sec() ;
1618 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
1619 //os_ << " filename = " << filename << endl ;
1620
1621 // parsing MS options
1622 if ( rec.isDefined( "ms" ) ) {
1623 Record msrec = rec.asRecord( "ms" ) ;
1624 if ( msrec.isDefined( "getpt" ) ) {
1625 getPt_ = msrec.asBool( "getpt" ) ;
1626 }
1627 if ( msrec.isDefined( "antenna" ) ) {
1628 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
1629 antenna_ = msrec.asInt( "antenna" ) ;
1630 }
1631 else {
1632 //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
1633 antennaStr_ = msrec.asString( "antenna" ) ;
1634 }
1635 }
1636 else {
1637 antenna_ = 0 ;
1638 }
1639 if ( msrec.isDefined( "freq_tolsr" ) ) {
1640 freqToLsr_ = msrec.asBool( "freq_tolsr" ) ;
1641 }
1642 }
1643
1644 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
1645 tablename_ = tmpMS->tableName() ;
1646 if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
1647 MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
1648 Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
1649 if ( id.size() > 0 )
1650 antenna_ = id[0] ;
1651 else {
1652 delete tmpMS ;
1653 //throw( AipsError( "Antenna " + antennaStr_ + " doesn't exist." ) ) ;
1654 os_ << LogIO::SEVERE << "Antenna " << antennaStr_ << " doesn't exist." << LogIO::POST ;
1655 return False ;
1656 }
1657 }
1658
1659 os_ << "Parsing MS options" << endl ;
1660 os_ << " getPt = " << (getPt_ ? "True" : "False") << endl ;
1661 os_ << " antenna = " << antenna_ << endl ;
1662 os_ << " antennaStr = " << antennaStr_ << endl ;
1663 os_ << " freqToLsr = " << (freqToLsr_ ? "True" : "False") << LogIO::POST;
1664
1665 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
1666 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
1667
1668 delete tmpMS ;
1669
1670 // check which data column exists
1671 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
1672 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
1673
1674 //double endSec = mathutil::gettimeofday_sec() ;
1675 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1676 return true ;
1677}
1678
1679void MSFiller::fill()
1680{
1681 //double startSec = mathutil::gettimeofday_sec() ;
1682 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
1683
1684 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
1685
1686 // Initialize header
1687 STHeader sdh ;
1688 initHeader( sdh ) ;
1689 table_->setHeader( sdh ) ;
1690
1691 // check if optional table exists
1692 const TableRecord &msrec = mstable_.keywordSet() ;
1693 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
1694 if ( isDoppler_ )
1695 if ( mstable_.doppler().nrow() == 0 )
1696 isDoppler_ = False ;
1697 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
1698 if ( isFlagCmd_ )
1699 if ( mstable_.flagCmd().nrow() == 0 )
1700 isFlagCmd_ = False ;
1701 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
1702 if ( isFreqOffset_ )
1703 if ( mstable_.freqOffset().nrow() == 0 )
1704 isFreqOffset_ = False ;
1705 isHistory_ = msrec.isDefined( "HISTORY" ) ;
1706 if ( isHistory_ )
1707 if ( mstable_.history().nrow() == 0 )
1708 isHistory_ = False ;
1709 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
1710 if ( isProcessor_ )
1711 if ( mstable_.processor().nrow() == 0 )
1712 isProcessor_ = False ;
1713 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
1714 if ( isSysCal_ )
1715 if ( mstable_.sysCal().nrow() == 0 )
1716 isSysCal_ = False ;
1717 isWeather_ = msrec.isDefined( "WEATHER" ) ;
1718 if ( isWeather_ )
1719 if ( mstable_.weather().nrow() == 0 )
1720 isWeather_ = False ;
1721
1722 // column name for Tsys and Tcal
1723 if ( isSysCal_ ) {
1724 const MSSysCal &caltab = mstable_.sysCal() ;
1725 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
1726 colTcal_ = "TCAL" ;
1727 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1728 colTcal_ = "NONE" ;
1729 }
1730 if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
1731 colTsys_ = "TSYS" ;
1732 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1733 colTsys_ = "NONE" ;
1734 }
1735 }
1736 else {
1737 colTcal_ = "NONE" ;
1738 colTsys_ = "NONE" ;
1739 }
1740
1741 // Access to MS subtables
1742 //MSField &fieldtab = mstable_.field() ;
1743 //MSPolarization &poltab = mstable_.polarization() ;
1744 //MSDataDescription &ddtab = mstable_.dataDescription() ;
1745 //MSObservation &obstab = mstable_.observation() ;
1746 //MSSource &srctab = mstable_.source() ;
1747 //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
1748 //MSSysCal &caltab = mstable_.sysCal() ;
1749 MSPointing &pointtab = mstable_.pointing() ;
1750 //MSState &stattab = mstable_.state() ;
1751 //MSAntenna &anttab = mstable_.antenna() ;
1752
1753 // SUBTABLES: FREQUENCIES
1754 //string freqFrame = getFrame() ;
1755 string freqFrame = "LSRK" ;
1756 table_->frequencies().setFrame( freqFrame ) ;
1757 if ( freqToLsr_ ) {
1758 table_->frequencies().setFrame( freqFrame, True ) ;
1759 }
1760 else {
1761 string baseFrame = frameFromSpwTable() ;
1762 table_->frequencies().setFrame( baseFrame, True ) ;
1763 }
1764
1765 // SUBTABLES: WEATHER
1766 fillWeather() ;
1767
1768 // SUBTABLES: FOCUS
1769 fillFocus() ;
1770
1771 // SUBTABLES: TCAL
1772 fillTcal() ;
1773
1774 // SUBTABLES: FIT
1775 //fillFit() ;
1776
1777 // SUBTABLES: HISTORY
1778 //fillHistory() ;
1779
1780 /***
1781 * Start iteration using TableVisitor
1782 ***/
1783 Table stab = table_->table() ;
1784 {
1785 static const char *cols[] = {
1786 "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
1787 "STATE_ID", "TIME",
1788 NULL
1789 };
1790 static const TypeManagerImpl<Int> tmInt;
1791 static const TypeManagerImpl<Double> tmDouble;
1792 static const TypeManager *const tms[] = {
1793 &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
1794 };
1795 //double t0 = mathutil::gettimeofday_sec() ;
1796 MSFillerVisitor myVisitor(mstable_, *table_ );
1797 //double t1 = mathutil::gettimeofday_sec() ;
1798 //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1799 myVisitor.setAntenna( antenna_ ) ;
1800 //myVisitor.setHeader( sdh ) ;
1801 if ( getPt_ ) {
1802 Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
1803 myVisitor.setPointingTable( ptsel ) ;
1804 }
1805 if ( isWeather_ )
1806 myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ;
1807 if ( isSysCal_ )
1808 myVisitor.setSysCalRecord( tcalrec_ ) ;
1809 myVisitor.setFreqToLsr( freqToLsr_ ) ;
1810
1811 //double t2 = mathutil::gettimeofday_sec() ;
1812 traverseTable(mstable_, cols, tms, &myVisitor);
1813 //double t3 = mathutil::gettimeofday_sec() ;
1814 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1815
1816 sdh = myVisitor.getHeader() ;
1817 }
1818 /***
1819 * End iteration using TableVisitor
1820 ***/
1821
1822 // set header
1823 //sdh = myVisitor.getHeader() ;
1824 //table_->setHeader( sdh ) ;
1825
1826 // save path to POINTING table
1827 // 2011/07/06 TN
1828 // Path to POINTING table in original MS will not be written
1829 // if getPt_ is True
1830 Path datapath( tablename_ ) ;
1831 if ( !getPt_ ) {
1832 String pTabName = datapath.absoluteName() + "/POINTING" ;
1833 stab.rwKeywordSet().define( "POINTING", pTabName ) ;
1834 }
1835
1836 // for GBT
1837 if ( sdh.antennaname.contains( "GBT" ) ) {
1838 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
1839 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
1840 }
1841
1842 // for MS created from ASDM
1843 const TableRecord &msKeys = mstable_.keywordSet() ;
1844 uInt nfields = msKeys.nfields() ;
1845 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1846 String name = msKeys.name( ifield ) ;
1847 //os_ << "name = " << name << LogIO::POST ;
1848 if ( name.find( "ASDM" ) != String::npos ) {
1849 String asdmpath = msKeys.asTable( ifield ).tableName() ;
1850 os_ << "ASDM table: " << asdmpath << LogIO::POST ;
1851 stab.rwKeywordSet().define( name, asdmpath ) ;
1852 }
1853 }
1854
1855 //double endSec = mathutil::gettimeofday_sec() ;
1856 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1857}
1858
1859void MSFiller::close()
1860{
1861 //tablesel_.closeSubTables() ;
1862 mstable_.closeSubTables() ;
1863 //tablesel_.unlock() ;
1864 mstable_.unlock() ;
1865}
1866
1867void MSFiller::fillWeather()
1868{
1869 //double startSec = mathutil::gettimeofday_sec() ;
1870 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1871
1872 if ( !isWeather_ ) {
1873 // add dummy row
1874 table_->weather().table().addRow(1,True) ;
1875 return ;
1876 }
1877
1878 Table mWeather = mstable_.weather() ;
1879 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1880 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1881 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1882 if ( mWeatherSel.nrow() == 0 ) {
1883 os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1884 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1885 if ( mWeatherSel.nrow() == 0 ) {
1886 os_ << "No rows in WEATHER table" << LogIO::POST ;
1887 }
1888 }
1889 uInt wnrow = mWeatherSel.nrow() ;
1890 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1891
1892 if ( wnrow == 0 )
1893 return ;
1894
1895 Table wtab = table_->weather().table() ;
1896 wtab.addRow( wnrow ) ;
1897
1898 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
1899 Int stationId = -1 ;
1900 if ( stationInfoExists ) {
1901 // determine which station is closer
1902 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1903 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
1904 Vector<Int> stationIds = stationCol.getColumn() ;
1905 Vector<Int> stationIdList( 0 ) ;
1906 Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
1907 uInt numStation = 0 ;
1908 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
1909 if ( !anyEQ( stationIdList, stationIds[i] ) ) {
1910 numStation++ ;
1911 stationIdList.resize( numStation, True ) ;
1912 stationIdList[numStation-1] = stationIds[i] ;
1913 stationPosList.resize( numStation, 3, True ) ;
1914 stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
1915 }
1916 }
1917 //os_ << "staionIdList = " << stationIdList << endl ;
1918 Table mAntenna = mstable_.antenna() ;
1919 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
1920 Vector<Double> antpos = antposCol( antenna_ ) ;
1921 Double minDiff = -1.0 ;
1922 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
1923 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
1924 if ( minDiff < 0.0 || minDiff > diff ) {
1925 minDiff = diff ;
1926 stationId = stationIdList[i] ;
1927 }
1928 }
1929 }
1930 //os_ << "stationId = " << stationId << endl ;
1931
1932 ScalarColumn<Float> *fCol ;
1933 ROScalarColumn<Float> *sharedFloatCol ;
1934 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1935 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1936 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1937 fCol->putColumn( *sharedFloatCol ) ;
1938 delete sharedFloatCol ;
1939 delete fCol ;
1940 }
1941 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1942 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1943 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1944 fCol->putColumn( *sharedFloatCol ) ;
1945 delete sharedFloatCol ;
1946 delete fCol ;
1947 }
1948 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1949 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1950 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1951 fCol->putColumn( *sharedFloatCol ) ;
1952 delete sharedFloatCol ;
1953 delete fCol ;
1954 }
1955 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1956 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1957 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1958 fCol->putColumn( *sharedFloatCol ) ;
1959 delete sharedFloatCol ;
1960 delete fCol ;
1961 }
1962 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1963 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1964 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1965 fCol->putColumn( *sharedFloatCol ) ;
1966 delete sharedFloatCol ;
1967 delete fCol ;
1968 }
1969 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1970 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1971 idCol.put( irow, irow ) ;
1972
1973 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1974 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1975 String tUnit = tqCol.getUnits() ;
1976 Vector<Double> mwTime = tCol.getColumn() ;
1977 if ( tUnit == "d" )
1978 mwTime *= 86400.0 ;
1979 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1980 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1981 String iUnit = tqCol.getUnits() ;
1982 Vector<Double> mwInterval = tCol.getColumn() ;
1983 if ( iUnit == "d" )
1984 mwInterval *= 86400.0 ;
1985
1986 if ( stationId > 0 ) {
1987 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1988 Vector<Int> stationVec = stationCol.getColumn() ;
1989 uInt wsnrow = ntrue( stationVec == stationId ) ;
1990 mwTime_.resize( wsnrow ) ;
1991 mwInterval_.resize( wsnrow ) ;
1992 mwIndex_.resize( wsnrow ) ;
1993 uInt wsidx = 0 ;
1994 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
1995 if ( stationId == stationVec[irow] ) {
1996 mwTime_[wsidx] = mwTime[irow] ;
1997 mwInterval_[wsidx] = mwInterval[irow] ;
1998 mwIndex_[wsidx] = irow ;
1999 wsidx++ ;
2000 }
2001 }
2002 }
2003 else {
2004 mwTime_ = mwTime ;
2005 mwInterval_ = mwInterval ;
2006 mwIndex_.resize( mwTime_.size() ) ;
2007 indgen( mwIndex_ ) ;
2008 }
2009 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
2010 //double endSec = mathutil::gettimeofday_sec() ;
2011 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2012}
2013
2014void MSFiller::fillFocus()
2015{
2016 //double startSec = mathutil::gettimeofday_sec() ;
2017 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
2018 // tentative
2019 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
2020 //double endSec = mathutil::gettimeofday_sec() ;
2021 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2022}
2023
2024void MSFiller::fillTcal()
2025{
2026 //double startSec = mathutil::gettimeofday_sec() ;
2027 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
2028
2029 if ( !isSysCal_ ) {
2030 // add dummy row
2031 os_ << "No SYSCAL rows" << LogIO::POST ;
2032 table_->tcal().table().addRow(1,True) ;
2033 Vector<Float> defaultTcal( 1, 1.0 ) ;
2034 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2035 tcalCol.put( 0, defaultTcal ) ;
2036 return ;
2037 }
2038
2039 if ( colTcal_ == "NONE" ) {
2040 // add dummy row
2041 os_ << "No TCAL column" << LogIO::POST ;
2042 table_->tcal().table().addRow(1,True) ;
2043 Vector<Float> defaultTcal( 1, 1.0 ) ;
2044 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2045 tcalCol.put( 0, defaultTcal ) ;
2046 return ;
2047 }
2048
2049 Table &sctab = mstable_.sysCal() ;
2050 if ( sctab.nrow() == 0 ) {
2051 os_ << "No SYSCAL rows" << LogIO::POST ;
2052 return ;
2053 }
2054 ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
2055 Vector<Int> ant = antCol.getColumn() ;
2056 if ( allNE( ant, antenna_ ) ) {
2057 os_ << "No SYSCAL rows" << LogIO::POST ;
2058 return ;
2059 }
2060 ROTableColumn tcalCol( sctab, colTcal_ ) ;
2061 Bool notDefined = False ;
2062 for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
2063 if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
2064 notDefined = True ;
2065 break ;
2066 }
2067 }
2068 if ( notDefined ) {
2069 os_ << "No TCAL value" << LogIO::POST ;
2070 table_->tcal().table().addRow(1,True) ;
2071 Vector<Float> defaultTcal( 1, 1.0 ) ;
2072 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2073 tcalCol.put( 0, defaultTcal ) ;
2074 return ;
2075 }
2076
2077 static const char *cols[] = {
2078 "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
2079 NULL
2080 };
2081 static const TypeManagerImpl<Int> tmInt;
2082 static const TypeManagerImpl<Double> tmDouble;
2083 static const TypeManager *const tms[] = {
2084 &tmInt, &tmInt, &tmInt, &tmDouble, NULL
2085 };
2086 Table tab = table_->tcal().table() ;
2087 TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
2088 visitor.setTcalColumn( colTcal_ ) ;
2089
2090 traverseTable(sctab, cols, tms, &visitor);
2091
2092 infillTcal();
2093
2094 //tcalrec_.print( std::cout ) ;
2095 //double endSec = mathutil::gettimeofday_sec() ;
2096 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2097}
2098
2099void MSFiller::infillTcal()
2100{
2101 uInt nfields = tcalrec_.nfields() ;
2102 set<Int> spwAvailable;
2103 for (uInt i = 0; i < nfields; i++) {
2104 String name = tcalrec_.name(i);
2105 size_t pos1 = name.find(':') + 4;
2106 size_t pos2 = name.find(':',pos1);
2107 Int spwid = String::toInt(name.substr(pos1,pos2-pos1));
2108 //cout << "spwid=" << spwid << endl;
2109 spwAvailable.insert(spwid);
2110 }
2111 Table spwtab = mstable_.spectralWindow();
2112 Table tcaltab = table_->tcal().table();
2113 ScalarColumn<uInt> idCol(tcaltab, "ID");
2114 ScalarColumn<String> timeCol(tcaltab, "TIME");
2115 ArrayColumn<Float> tcalCol(tcaltab, "TCAL");
2116 ROScalarColumn<Int> numChanCol(spwtab, "NUM_CHAN");
2117 Int numSpw = spwtab.nrow();
2118 Int dummyFeed = 0;
2119 Double dummyTime = 0.0;
2120 Vector<uInt> idminmax(2);
2121 for (Int i = 0; i < numSpw; i++) {
2122 if (spwAvailable.find(i) == spwAvailable.end()) {
2123 String key = keyTcal(dummyFeed, i, dummyTime);
2124 Vector<Float> tcal(numChanCol(i), 1.0);
2125 uInt nrow = tcaltab.nrow();
2126 tcaltab.addRow(1);
2127 idCol.put(nrow, nrow);
2128 timeCol.put(nrow, "");
2129 tcalCol.put(nrow, tcal);
2130 idminmax = nrow;
2131 tcalrec_.define(key, idminmax);
2132 }
2133 }
2134 //tcalrec_.print(cout);
2135}
2136
2137string MSFiller::getFrame()
2138{
2139 MFrequency::Types frame = MFrequency::DEFAULT ;
2140 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
2141 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
2142 uInt nrow = numChanCol.nrow() ;
2143 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
2144 uInt nref = 0 ;
2145 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2146 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
2147 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
2148 nref++ ;
2149 }
2150 }
2151 if ( nref > 0 )
2152 frame = (MFrequency::Types)measFreqRef[0] ;
2153
2154 return MFrequency::showType( frame ) ;
2155}
2156
2157void MSFiller::initHeader( STHeader &header )
2158{
2159 header.nchan = 0 ;
2160 header.npol = 0 ;
2161 header.nif = 0 ;
2162 header.nbeam = 0 ;
2163 header.observer = "" ;
2164 header.project = "" ;
2165 header.obstype = "" ;
2166 header.antennaname = "" ;
2167 header.antennaposition.resize( 3 ) ;
2168 header.equinox = 0.0 ;
2169 header.freqref = "" ;
2170 header.reffreq = -1.0 ;
2171 header.bandwidth = 0.0 ;
2172 header.utc = 0.0 ;
2173 header.fluxunit = "" ;
2174 header.epoch = "" ;
2175 header.poltype = "" ;
2176}
2177
2178string MSFiller::frameFromSpwTable()
2179{
2180 string frameString;
2181 Table tab = mstable_.spectralWindow();
2182 ROScalarColumn<Int> mfrCol(tab, "MEAS_FREQ_REF");
2183 Vector<Int> mfr = mfrCol.getColumn();
2184 if (allEQ(mfr,mfr[0])) {
2185 frameString = MFrequency::showType(mfr[0]);
2186 //cout << "all rows have same frame: " << frameString << endl;
2187 }
2188 else {
2189 mfrCol.attach(tab, "NUM_CHAN");
2190 for (uInt i = 0; i < tab.nrow(); i++) {
2191 if (mfrCol(i) != 4) {
2192 frameString = MFrequency::showType(mfr[i]);
2193 break;
2194 }
2195 }
2196 if (frameString.size() == 0) {
2197 frameString = "TOPO";
2198 }
2199 }
2200
2201 //cout << "frameString = " << frameString << endl;
2202
2203 return frameString;
2204}
2205
2206};
Note: See TracBrowser for help on using the repository browser.