source: trunk/src/MSFiller.cpp@ 2291

Last change on this file since 2291 was 2291, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: MSFillerUtils and MSWriterUtils added

Test Programs: sd regressions, test_sdsave

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote MSFiller and MSWriter based on TableTraverse.
TableTraverse is changed a bit since it doesn't work on plain table.


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