source: branches/hpc34/src/MSFiller.cpp@ 2643

Last change on this file since 2643 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

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