Changeset 1975 for trunk/src


Ignore:
Timestamp:
01/20/11 12:13:52 (14 years ago)
Author:
Takeshi Nakazato
Message:

New Development: Yes

JIRA Issue: Yes CAS-2718

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...

Writer improvement.


Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/MSFiller.cpp

    r1974 r1975  
    273273          if ( !sysVels.empty() ) {
    274274            //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ;
     275            // NB: assume all SYSVEL values are the same
    275276            Double sysVel = sysVels( IPosition(1,0) ) ;
    276277          }

  • trunk/src/MSWriter.cpp

    r1974 r1975  
    1818#include <casa/BasicSL/String.h>
    1919
     20#include <tables/Tables/ExprNode.h>
    2021#include <tables/Tables/TableDesc.h>
    2122#include <tables/Tables/SetupNewTab.h>
    … …  
    2728#include <ms/MeasurementSets/MSPolIndex.h>
    2829#include <ms/MeasurementSets/MSDataDescIndex.h>
     30#include <ms/MeasurementSets/MSSourceIndex.h>
    2931
    3032#include "MSWriter.h"
    3133#include "STHeader.h"
    3234#include "STFrequencies.h"
     35#include "STMolecules.h"
    3336
    3437using namespace casa ;
    … …  
    9699  // ANTENNA
    97100  fillAntenna() ;
     101
     102  // PROCESSOR
     103  fillProcessor() ;
     104
     105  // SOURCE
     106  fillSource() ;
    98107
    99108  // MAIN
    … …  
    109118  while( !iter0.pastEnd() ) {
    110119    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] ;
    113128    String::size_type pos = fieldName.find( "__" ) ;
    114129    Int fieldId = -1 ;
    … …  
    133148      Table t1( iter1.table() ) ;
    134149      ROScalarColumn<uInt> beamNoCol( t1, "BEAMNO" ) ;
    135       uInt beamNo = beamNoCol(0) ;
     150      uInt beamNo = beamNoCol( 0 ) ;
    136151      os_ << "beamNo = " << beamNo << LogIO::POST ;
    137152      //
    … …  
    144159        Table t2( iter2.table() ) ;
    145160        ROScalarColumn<uInt> scanNoCol( t2, "SCANNO" ) ;
    146         uInt scanNo = scanNoCol(0) ;
     161        uInt scanNo = scanNoCol( 0 ) ;
    147162        os_ << "scanNo = " << scanNo << LogIO::POST ;
    148163        //
    … …  
    163178            Table t4( iter4.table() ) ;
    164179            ROScalarColumn<uInt> ifNoCol( t4, "IFNO" ) ;
    165             uInt ifNo = ifNoCol(0) ;
     180            uInt ifNo = ifNoCol( 0 ) ;
    166181            os_ << "ifNo = " << ifNo << LogIO::POST ;
    167182            ROScalarColumn<uInt> freqIdCol( t4, "FREQ_ID" ) ;
    168             uInt freqId = freqIdCol(0) ;
     183            uInt freqId = freqIdCol( 0 ) ;
    169184            os_ << "freqId = " << freqId << LogIO::POST ;
    170185            //
    … …  
    196211              ROScalarMeasColumn<MEpoch> timeCol( t5, "TIME" ) ;
    197212              ScalarMeasColumn<MEpoch> msTimeCol( *mstable_, "TIME" ) ;
    198               msTimeCol.put( prevnr, timeCol(0) ) ;
     213              msTimeCol.put( prevnr, timeCol( 0 ) ) ;
    199214              msTimeCol.attach( *mstable_, "TIME_CENTROID" ) ;
    200               msTimeCol.put( prevnr, timeCol(0) ) ;
     215              msTimeCol.put( prevnr, timeCol( 0 ) ) ;
    201216             
    202217              // INTERVAL and EXPOSURE
    203218              ROScalarColumn<Double> intervalCol( t5, "INTERVAL" ) ;
    204219              ScalarColumn<Double> msIntervalCol( *mstable_, "INTERVAL" ) ;
    205               msIntervalCol.put( prevnr, intervalCol(0) ) ;
     220              msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
    206221              msIntervalCol.attach( *mstable_, "EXPOSURE" ) ;
    207               msIntervalCol.put( prevnr, intervalCol(0) ) ;
     222              msIntervalCol.put( prevnr, intervalCol( 0 ) ) ;
    208223             
    209224              // add DATA_DESCRIPTION row
    … …  
    268283    fieldIdCol.putColumnCells( rows1, fieldIds ) ;
    269284
     285    // add FIELD row
     286    addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
     287
    270288    added0 += added1 ;
    271289    os_ << "added0 = " << added0 << " current0 = " << current0 << LogIO::POST ;
    … …  
    288306  sharedIntArr = 0 ;
    289307  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" ) ;
    290313  sharedIntCol.putColumn( sharedIntArr ) ;
    291314
    … …  
    392415
    393416  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 ) ;
    394420  SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
    395421  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
    … …  
    480506}
    481507
     508void 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
     517void 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
    482662void MSWriter::addFeed( Int id )
    483663{
    … …  
    547727}
    548728
     729void 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
    549771Int MSWriter::addPolarization( Vector<Int> polnos )
    550772{
    551   os_ << "set up SPECTRAL_WINDOW subtable" << LogIO::POST ;
     773  os_ << "set up POLARIZATION subtable" << LogIO::POST ;
    552774
    553775  os_ << "polnos = " << polnos << LogIO::POST ;
    … …  
    584806    }
    585807    else if ( npol == 2 ) {
    586       corrProd.column(0) = 0 ;
    587       corrProd.column(1) = 1 ;
     808      corrProd.column( 0 ) = 0 ;
     809      corrProd.column( 1 ) = 1 ;
    588810    }
    589811    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 ;
    596818    }
    597819    msPolCols.corrProduct().put( nrow, corrProd ) ;
    … …  
    604826  return polid ;
    605827}
     828
     829Int 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}
    606855
    607856Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
    … …  
    694943}
    695944
    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 }
     945void 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  
    2222#include <casa/Logging/LogIO.h>
    2323
     24#include <tables/Tables/Table.h>
    2425#include <tables/Tables/RefRows.h>
    2526
    2627#include <ms/MeasurementSets/MeasurementSet.h>
    2728#include <ms/MeasurementSets/MSColumns.h>
     29
     30#include <measures/Measures/MEpoch.h>
    2831
    2932#include "Scantable.h"
    … …  
    5558  void fillObservation() ;
    5659  void fillAntenna() ;
     60  void fillProcessor() ;
     61  void fillSource() ;
     62
     63  // add rows to subtables
    5764  void addFeed( casa::Int id ) ;
    5865  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 ) ;
    5967  casa::Int addPolarization( casa::Vector<casa::Int> polnos ) ;
    6068  casa::Int addDataDescription( casa::Int polid, casa::Int spwid ) ;
    … …  
    6270  // utility
    6371  casa::Vector<casa::Int> toCorrType( casa::Vector<casa::Int> polnos ) ;
     72  void getValidTimeRange( casa::MEpoch &me, casa::Double &interval, casa::Table &tab ) ;
    6473
    6574  casa::CountedPtr<Scantable> table_ ;
Note: See TracChangeset for help on using the changeset viewer.