source: branches/hpc33/src/MSFiller.cpp@ 2550

Last change on this file since 2550 was 2550, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Clean up code.


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