- Timestamp:
- 01/20/11 12:13:52 (14 years ago)
- Location:
- trunk/src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/MSFiller.cpp
r1974 r1975 273 273 if ( !sysVels.empty() ) { 274 274 //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ; 275 // NB: assume all SYSVEL values are the same 275 276 Double sysVel = sysVels( IPosition(1,0) ) ; 276 277 } -
trunk/src/MSWriter.cpp
r1974 r1975 18 18 #include <casa/BasicSL/String.h> 19 19 20 #include <tables/Tables/ExprNode.h> 20 21 #include <tables/Tables/TableDesc.h> 21 22 #include <tables/Tables/SetupNewTab.h> … … 27 28 #include <ms/MeasurementSets/MSPolIndex.h> 28 29 #include <ms/MeasurementSets/MSDataDescIndex.h> 30 #include <ms/MeasurementSets/MSSourceIndex.h> 29 31 30 32 #include "MSWriter.h" 31 33 #include "STHeader.h" 32 34 #include "STFrequencies.h" 35 #include "STMolecules.h" 33 36 34 37 using namespace casa ; … … 96 99 // ANTENNA 97 100 fillAntenna() ; 101 102 // PROCESSOR 103 fillProcessor() ; 104 105 // SOURCE 106 fillSource() ; 98 107 99 108 // MAIN … … 109 118 while( !iter0.pastEnd() ) { 110 119 Table t0( iter0.table() ) ; 111 ROScalarColumn<String> fieldNameCol( t0, "FIELDNAME" ) ; 112 String fieldName = fieldNameCol(0) ; 120 ROScalarColumn<String> sharedStrCol( t0, "FIELDNAME" ) ; 121 String fieldName = sharedStrCol( 0 ) ; 122 sharedStrCol.attach( t0, "SRCNAME" ) ; 123 String srcName = sharedStrCol( 0 ) ; 124 ROScalarColumn<Double> timeCol( t0, "TIME" ) ; 125 Double minTime = min( timeCol.getColumn() ) ; 126 ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ; 127 Vector<Double> scanRate = scanRateCol.getColumn()[0] ; 113 128 String::size_type pos = fieldName.find( "__" ) ; 114 129 Int fieldId = -1 ; … … 133 148 Table t1( iter1.table() ) ; 134 149 ROScalarColumn<uInt> beamNoCol( t1, "BEAMNO" ) ; 135 uInt beamNo = beamNoCol( 0) ;150 uInt beamNo = beamNoCol( 0 ) ; 136 151 os_ << "beamNo = " << beamNo << LogIO::POST ; 137 152 // … … 144 159 Table t2( iter2.table() ) ; 145 160 ROScalarColumn<uInt> scanNoCol( t2, "SCANNO" ) ; 146 uInt scanNo = scanNoCol( 0) ;161 uInt scanNo = scanNoCol( 0 ) ; 147 162 os_ << "scanNo = " << scanNo << LogIO::POST ; 148 163 // … … 163 178 Table t4( iter4.table() ) ; 164 179 ROScalarColumn<uInt> ifNoCol( t4, "IFNO" ) ; 165 uInt ifNo = ifNoCol( 0) ;180 uInt ifNo = ifNoCol( 0 ) ; 166 181 os_ << "ifNo = " << ifNo << LogIO::POST ; 167 182 ROScalarColumn<uInt> freqIdCol( t4, "FREQ_ID" ) ; 168 uInt freqId = freqIdCol( 0) ;183 uInt freqId = freqIdCol( 0 ) ; 169 184 os_ << "freqId = " << freqId << LogIO::POST ; 170 185 // … … 196 211 ROScalarMeasColumn<MEpoch> timeCol( t5, "TIME" ) ; 197 212 ScalarMeasColumn<MEpoch> msTimeCol( *mstable_, "TIME" ) ; 198 msTimeCol.put( prevnr, timeCol( 0) ) ;213 msTimeCol.put( prevnr, timeCol( 0 ) ) ; 199 214 msTimeCol.attach( *mstable_, "TIME_CENTROID" ) ; 200 msTimeCol.put( prevnr, timeCol( 0) ) ;215 msTimeCol.put( prevnr, timeCol( 0 ) ) ; 201 216 202 217 // INTERVAL and EXPOSURE 203 218 ROScalarColumn<Double> intervalCol( t5, "INTERVAL" ) ; 204 219 ScalarColumn<Double> msIntervalCol( *mstable_, "INTERVAL" ) ; 205 msIntervalCol.put( prevnr, intervalCol( 0) ) ;220 msIntervalCol.put( prevnr, intervalCol( 0 ) ) ; 206 221 msIntervalCol.attach( *mstable_, "EXPOSURE" ) ; 207 msIntervalCol.put( prevnr, intervalCol( 0) ) ;222 msIntervalCol.put( prevnr, intervalCol( 0 ) ) ; 208 223 209 224 // add DATA_DESCRIPTION row … … 268 283 fieldIdCol.putColumnCells( rows1, fieldIds ) ; 269 284 285 // add FIELD row 286 addField( fieldId, fieldName, srcName, minTime, scanRate ) ; 287 270 288 added0 += added1 ; 271 289 os_ << "added0 = " << added0 << " current0 = " << current0 << LogIO::POST ; … … 288 306 sharedIntArr = 0 ; 289 307 sharedIntCol.attach( *mstable_, "ARRAY_ID" ) ; 308 sharedIntCol.putColumn( sharedIntArr ) ; 309 310 // PROCESSOR_ID is tentatively set to 0 311 sharedIntArr = 0 ; 312 sharedIntCol.attach( *mstable_, "PROCESSOR_ID" ) ; 290 313 sharedIntCol.putColumn( sharedIntArr ) ; 291 314 … … 392 415 393 416 TableDesc sourceDesc = MSSource::requiredTableDesc() ; 417 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ; 418 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ; 419 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ; 394 420 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ; 395 421 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ; … … 480 506 } 481 507 508 void MSWriter::fillProcessor() 509 { 510 os_ << "set up PROCESSOR subtable" << LogIO::POST ; 511 512 // only add empty 1 row 513 MSProcessor msProc = mstable_->processor() ; 514 msProc.addRow( 1, True ) ; 515 } 516 517 void MSWriter::fillSource() 518 { 519 os_ << "set up SOURCE subtable" << LogIO::POST ; 520 521 522 // access to MS SOURCE subtable 523 MSSource msSrc = mstable_->source() ; 524 525 // access to MOLECULE subtable 526 STMolecules stm = table_->molecules() ; 527 528 Int srcId = 0 ; 529 Vector<Double> restFreq ; 530 Vector<String> molName ; 531 Vector<String> fMolName ; 532 // 533 // ITERATION: SRCNAME 534 // 535 Int added0 = 0 ; 536 Int current0 = msSrc.nrow() ; 537 TableIterator iter0( table_->table(), "SRCNAME" ) ; 538 while( !iter0.pastEnd() ) { 539 Table t0( iter0.table() ) ; 540 541 // get necessary information 542 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ; 543 String srcName = srcNameCol( 0 ) ; 544 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ; 545 Vector<Double> srcPM = sharedDArrRCol( 0 ) ; 546 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ; 547 Vector<Double> srcDir = sharedDArrRCol( 0 ) ; 548 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ; 549 Double srcVel = srcVelCol( 0 ) ; 550 551 // 552 // ITERATION: MOLECULE_ID 553 // 554 Int added1 = 0 ; 555 Int current1 = msSrc.nrow() ; 556 TableIterator iter1( t0, "MOLECULE_ID" ) ; 557 while( !iter1.pastEnd() ) { 558 Table t1( iter1.table() ) ; 559 560 // get necessary information 561 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ; 562 uInt molId = molIdCol( 0 ) ; 563 stm.getEntry( restFreq, molName, fMolName, molId ) ; 564 565 // 566 // ITERATION: IFNO 567 // 568 Int added2 = 0 ; 569 Int current2 = msSrc.nrow() ; 570 TableIterator iter2( t1, "IFNO" ) ; 571 while( !iter2.pastEnd() ) { 572 Table t2( iter2.table() ) ; 573 uInt prevnr = msSrc.nrow() ; 574 575 // get necessary information 576 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ; 577 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID 578 MEpoch midTimeMeas ; 579 Double interval ; 580 getValidTimeRange( midTimeMeas, interval, t2 ) ; 581 582 // add row 583 msSrc.addRow( 1, True ) ; 584 585 // fill SPECTRAL_WINDOW_ID 586 ScalarColumn<Int> sSpwIdCol( msSrc, "SPECTRAL_WINDOW_ID" ) ; 587 sSpwIdCol.put( prevnr, ifno ) ; 588 589 // fill TIME, INTERVAL 590 ScalarMeasColumn<MEpoch> sTimeMeasCol( msSrc, "TIME" ) ; 591 sTimeMeasCol.put( prevnr, midTimeMeas ) ; 592 ScalarColumn<Double> sIntervalCol( msSrc, "INTERVAL" ) ; 593 sIntervalCol.put( prevnr, interval ) ; 594 595 added2++ ; 596 iter2.next() ; 597 } 598 599 // fill NUM_LINES, REST_FREQUENCY, TRANSITION, SYSVEL 600 RefRows rows2( current2, current2+added2-1 ) ; 601 uInt numFreq = restFreq.size() ; 602 Vector<Int> numLines( added2, numFreq ) ; 603 ScalarColumn<Int> numLinesCol( msSrc, "NUM_LINES" ) ; 604 numLinesCol.putColumnCells( rows2, numLines ) ; 605 if ( numFreq != 0 ) { 606 Array<Double> rf( IPosition( 2, numFreq, added2 ) ) ; 607 Array<String> trans( IPosition( 2, numFreq, added2 ) ) ; 608 Array<Double> srcVelArr( IPosition( 2, numFreq, added2 ), srcVel ) ; 609 for ( uInt ifreq = 0 ; ifreq < numFreq ; ifreq++ ) { 610 Slicer slice( Slice(ifreq),Slice(0,added2,1) ) ; 611 rf( slice ) = restFreq[ifreq] ; 612 String transStr = fMolName[ifreq] ; 613 if ( transStr.size() == 0 ) { 614 transStr = molName[ifreq] ; 615 } 616 trans( slice ) = transStr ; 617 } 618 ArrayColumn<Double> sharedDArrCol( msSrc, "REST_FREQUENCY" ) ; 619 sharedDArrCol.putColumnCells( rows2, rf ) ; 620 sharedDArrCol.attach( msSrc, "SYSVEL" ) ; 621 sharedDArrCol.putColumnCells( rows2, srcVelArr ) ; 622 ArrayColumn<String> transCol( msSrc, "TRANSITION" ) ; 623 transCol.putColumnCells( rows2, trans ) ; 624 } 625 626 added1 += added2 ; 627 iter1.next() ; 628 } 629 630 // fill NAME, SOURCE_ID 631 RefRows rows1( current1, current1+added1-1 ) ; 632 Vector<String> nameArr( added1, srcName ) ; 633 Vector<Int> srcIdArr( added1, srcId ) ; 634 ScalarColumn<String> sNameCol( msSrc, "NAME" ) ; 635 ScalarColumn<Int> srcIdCol( msSrc, "SOURCE_ID" ) ; 636 sNameCol.putColumnCells( rows1, nameArr ) ; 637 srcIdCol.putColumnCells( rows1, srcIdArr ) ; 638 639 // fill DIRECTION, PROPER_MOTION 640 Array<Double> srcPMArr( IPosition( 2, 2, added1 ) ) ; 641 Array<Double> srcDirArr( IPosition( 2, 2, added1 ) ) ; 642 for ( uInt i = 0 ; i < 2 ; i++ ) { 643 Slicer slice( Slice(i),Slice(0,added1,1) ) ; 644 srcPMArr( slice ) = srcPM[i] ; 645 srcDirArr( slice ) = srcDir[i] ; 646 } 647 ArrayColumn<Double> sharedDArrCol( msSrc, "DIRECTION" ) ; 648 sharedDArrCol.putColumnCells( rows1, srcDirArr ) ; 649 sharedDArrCol.attach( msSrc, "PROPER_MOTION" ) ; 650 sharedDArrCol.putColumnCells( rows1, srcPMArr ) ; 651 652 // fill TIME, INTERVAL 653 654 // increment srcId if SRCNAME changed 655 srcId++ ; 656 657 added0 += added1 ; 658 iter0.next() ; 659 } 660 } 661 482 662 void MSWriter::addFeed( Int id ) 483 663 { … … 547 727 } 548 728 729 void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate ) 730 { 731 os_ << "set up FIELD subtable" << LogIO::POST ; 732 733 MSField msField = mstable_->field() ; 734 while( (Int)msField.nrow() <= fid ) { 735 msField.addRow( 1, True ) ; 736 } 737 MSFieldColumns msFieldCols( msField ) ; 738 739 // Access to SOURCE table 740 MSSource msSrc = mstable_->source() ; 741 742 // fill target row 743 msFieldCols.name().put( fid, fieldname ) ; 744 msFieldCols.time().put( fid, t ) ; 745 Int numPoly = 0 ; 746 if ( anyNE( rate, 0.0 ) ) 747 numPoly = 1 ; 748 msFieldCols.numPoly().put( fid, numPoly ) ; 749 MSSourceIndex msSrcIdx( msSrc ) ; 750 Int srcId = -1 ; 751 Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ; 752 if ( srcIdArr.size() != 0 ) { 753 srcId = srcIdArr[0] ; 754 MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId ) ; 755 ROMSSourceColumns msSrcCols( msSrcSel ) ; 756 Vector<Double> srcDir = msSrcCols.direction()( 0 ) ; 757 Array<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ; 758 os_ << "srcDirA = " << srcDirA << LogIO::POST ; 759 os_ << "sliced srcDirA = " << srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) << LogIO::POST ; 760 srcDirA( Slicer( Slice(0,2,1), Slice(0) ) ) = srcDir.reform( IPosition(2,2,1) ) ; 761 os_ << "srcDirA = " << srcDirA << LogIO::POST ; 762 if ( numPoly != 0 ) 763 srcDirA( Slicer( Slice(0,2,1), Slice(1) ) ) = rate ; 764 msFieldCols.phaseDir().put( fid, srcDirA ) ; 765 msFieldCols.referenceDir().put( fid, srcDirA ) ; 766 msFieldCols.delayDir().put( fid, srcDirA ) ; 767 } 768 msFieldCols.sourceId().put( fid, srcId ) ; 769 } 770 549 771 Int MSWriter::addPolarization( Vector<Int> polnos ) 550 772 { 551 os_ << "set up SPECTRAL_WINDOWsubtable" << LogIO::POST ;773 os_ << "set up POLARIZATION subtable" << LogIO::POST ; 552 774 553 775 os_ << "polnos = " << polnos << LogIO::POST ; … … 584 806 } 585 807 else if ( npol == 2 ) { 586 corrProd.column( 0) = 0 ;587 corrProd.column( 1) = 1 ;808 corrProd.column( 0 ) = 0 ; 809 corrProd.column( 1 ) = 1 ; 588 810 } 589 811 else { 590 corrProd.column( 0) = 0 ;591 corrProd.column( 3) = 1 ;592 corrProd( 0,1) = 0 ;593 corrProd( 1,1) = 1 ;594 corrProd( 0,2) = 1 ;595 corrProd( 1,2) = 0 ;812 corrProd.column( 0 ) = 0 ; 813 corrProd.column( 3 ) = 1 ; 814 corrProd( 0,1 ) = 0 ; 815 corrProd( 1,1 ) = 1 ; 816 corrProd( 0,2 ) = 1 ; 817 corrProd( 1,2 ) = 0 ; 596 818 } 597 819 msPolCols.corrProduct().put( nrow, corrProd ) ; … … 604 826 return polid ; 605 827 } 828 829 Int MSWriter::addDataDescription( Int polid, Int spwid ) 830 { 831 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ; 832 833 MSDataDescription msDataDesc = mstable_->dataDescription() ; 834 uInt nrow = msDataDesc.nrow() ; 835 836 MSDataDescIndex msDataDescIdx( msDataDesc ) ; 837 838 Vector<Int> ddids = msDataDescIdx.matchSpwIdAndPolznId( spwid, polid ) ; 839 os_ << "ddids = " << ddids << LogIO::POST ; 840 841 Int ddid = -1 ; 842 if ( ddids.size() == 0 ) { 843 msDataDesc.addRow( 1, True ) ; 844 MSDataDescColumns msDataDescCols( msDataDesc ) ; 845 msDataDescCols.polarizationId().put( nrow, polid ) ; 846 msDataDescCols.spectralWindowId().put( nrow, spwid ) ; 847 ddid = (Int)nrow ; 848 } 849 else { 850 ddid = ddids[0] ; 851 } 852 853 return ddid ; 854 } 606 855 607 856 Vector<Int> MSWriter::toCorrType( Vector<Int> polnos ) … … 694 943 } 695 944 696 Int MSWriter::addDataDescription( Int polid, Int spwid ) 697 { 698 os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ; 699 700 MSDataDescription msDataDesc = mstable_->dataDescription() ; 701 uInt nrow = msDataDesc.nrow() ; 702 703 MSDataDescIndex msDataDescIdx( msDataDesc ) ; 704 705 Vector<Int> ddids = msDataDescIdx.matchSpwIdAndPolznId( spwid, polid ) ; 706 os_ << "ddids = " << ddids << LogIO::POST ; 707 708 Int ddid = -1 ; 709 if ( ddids.size() == 0 ) { 710 msDataDesc.addRow( 1, True ) ; 711 MSDataDescColumns msDataDescCols( msDataDesc ) ; 712 msDataDescCols.polarizationId().put( nrow, polid ) ; 713 msDataDescCols.spectralWindowId().put( nrow, spwid ) ; 714 ddid = (Int)nrow ; 715 } 716 else { 717 ddid = ddids[0] ; 718 } 719 720 return ddid ; 721 } 722 } 945 void MSWriter::getValidTimeRange( MEpoch &me, Double &interval, Table &tab ) 946 { 947 ROScalarMeasColumn<MEpoch> timeMeasCol( tab, "TIME" ) ; 948 ROScalarColumn<Double> timeCol( tab, "TIME" ) ; 949 String refStr = timeMeasCol( 0 ).getRefString() ; 950 Vector<Double> timeArr = timeCol.getColumn() ; 951 MEpoch::Types meType ; 952 MEpoch::getType( meType, refStr ) ; 953 Unit tUnit = timeMeasCol.measDesc().getUnits()( 0 ) ; 954 Double minTime ; 955 Double maxTime ; 956 minMax( minTime, maxTime, timeArr ) ; 957 Double midTime = 0.5 * ( minTime + maxTime ) ; 958 me = MEpoch( Quantity( midTime, tUnit ), meType ) ; 959 interval = Quantity( maxTime-minTime, tUnit ).getValue( "s" ) ; 960 } 961 962 } -
trunk/src/MSWriter.h
r1974 r1975 22 22 #include <casa/Logging/LogIO.h> 23 23 24 #include <tables/Tables/Table.h> 24 25 #include <tables/Tables/RefRows.h> 25 26 26 27 #include <ms/MeasurementSets/MeasurementSet.h> 27 28 #include <ms/MeasurementSets/MSColumns.h> 29 30 #include <measures/Measures/MEpoch.h> 28 31 29 32 #include "Scantable.h" … … 55 58 void fillObservation() ; 56 59 void fillAntenna() ; 60 void fillProcessor() ; 61 void fillSource() ; 62 63 // add rows to subtables 57 64 void addFeed( casa::Int id ) ; 58 65 void addSpectralWindow( casa::Int spwid, casa::Int freqid ) ; 66 void addField( casa::Int fid, casa::String fieldname, casa::String srcname, casa::Double t, casa::Vector<casa::Double> scanrate ) ; 59 67 casa::Int addPolarization( casa::Vector<casa::Int> polnos ) ; 60 68 casa::Int addDataDescription( casa::Int polid, casa::Int spwid ) ; … … 62 70 // utility 63 71 casa::Vector<casa::Int> toCorrType( casa::Vector<casa::Int> polnos ) ; 72 void getValidTimeRange( casa::MEpoch &me, casa::Double &interval, casa::Table &tab ) ; 64 73 65 74 casa::CountedPtr<Scantable> table_ ;
Note:
See TracChangeset
for help on using the changeset viewer.