source: trunk/src/MSFiller.cpp@ 2844

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

New Development: No

JIRA Issue: Yes CAS-5535

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: Removed freq_tolsr option from scantable constructor

Test Programs: sdsave unit test

Put in Release Notes: No

Module(s): sd

Description: Describe your changes here...

The parameter freq_tolsr is removed from scantable constructor.
It is removed from python layer as well as C++ layer.
Now the code always behaves like 'freq_tolsr=False'.

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