Changeset 2310 for branches


Ignore:
Timestamp:
09/27/11 14:20:58 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave, regressions/fls3a_hi_regression

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Merge bug fix on trunk (r2309).


Location:
branches/casa-prerelease/pre-asap
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • branches/casa-prerelease/pre-asap

    • Property svn:mergeinfo changed
      /trunkmerged: 2309
  • branches/casa-prerelease/pre-asap/Makefile

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/SConstruct

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/apps

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/external-alma

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/external-alma/atnf/pks/pks_maths.cc

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/getsvnrev.sh

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/python

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/casa-prerelease/pre-asap/src

  • branches/casa-prerelease/pre-asap/src/MSWriter.cpp

    r2293 r2310  
    180180  const String *lastFieldName;
    181181  uInt lastRecordNo;
    182   uInt lastBeamNo, lastScanNo, lastIfNo;
     182  uInt lastBeamNo, lastScanNo, lastIfNo, lastPolNo;
    183183  Int lastSrcType;
    184184  uInt lastCycleNo;
    185185  Double lastTime;
    186   Int lastPolNo;
    187186protected:
    188187  const Table &table;
     
    233232    uInt cycleNo;
    234233    Double time;
    235     Int polNo;
     234    uInt polNo;
    236235    { // prologue
    237236      uInt i = 0;
     
    423422    makePolMap() ;
    424423    initFrequencies() ;
    425     initTcal() ;
    426424
    427425    //
     
    590588      holder.reset() ;
    591589    }
    592     if ( tcalKey != -1 ) {
    593       tcalNotYet[tcalKey] = False ;
    594       tcalKey = -1 ;
    595     }
    596590  }
    597591  virtual void enterPolNo(const uInt recordNo, uInt columnValue) {
    598592    //printf("%u: PolNo: %d\n", recordNo, columnValue);
    599     uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
    600     if ( tcalKey == -1 ) {
    601       tcalKey = tcalId ;
    602     }
    603     if ( tcalNotYet[tcalKey] ) {
    604       map< Int,Vector<uInt> >::iterator itr = tcalIdRec.find( tcalKey ) ;
    605       if ( itr != tcalIdRec.end() ) {
    606         Vector<uInt> ids = itr->second ;
    607         uInt nrow = ids.nelements() ;
    608         ids.resize( nrow+1, True ) ;
    609         ids[nrow] = tcalId ;
    610         tcalIdRec.erase( tcalKey ) ;
    611         tcalIdRec[tcalKey] = ids ;
    612       }
    613       else {
    614         Vector<uInt> rows( 1, tcalId ) ;
    615         tcalIdRec[tcalKey] = rows ;
    616       }
    617     }
    618     map< Int,Vector<uInt> >::iterator itr = tcalRowRec.find( tcalKey ) ;
    619     if ( itr != tcalRowRec.end() ) {
    620       Vector<uInt> rows = itr->second ;
    621       uInt nrow = rows.nelements() ;
    622       rows.resize( nrow+1, True ) ;
    623       rows[nrow] = recordNo ;
    624       tcalRowRec.erase( tcalKey ) ;
    625       tcalRowRec[tcalKey] = rows ;
    626     }
    627     else {
    628       Vector<uInt> rows( 1, recordNo ) ;
    629       tcalRowRec[tcalKey] = rows ;
    630     }
    631593  }
    632594  virtual void leavePolNo(const uInt recordNo, uInt columnValue) {
     
    690652    srcRec = r ;
    691653  }
    692   map< Int,Vector<uInt> > &getTcalIdRecord() { return tcalIdRec ; }
    693   map< Int,Vector<uInt> > &getTcalRowRecord() { return tcalRowRec ; }
    694654private:
    695655  void addField( Int &fid, String &fname, String &srcName,
     
    12501210    }
    12511211  }
    1252   void initTcal()
    1253   {
    1254     const TableRecord &rec = table.keywordSet() ;
    1255     Table tcalTable = rec.asTable( "TCAL" ) ;
    1256     ROScalarColumn<uInt> idCol( tcalTable, "ID" ) ;
    1257     Vector<uInt> id = idCol.getColumn() ;
    1258     uInt maxId = max( id ) ;
    1259     tcalNotYet.resize( maxId+1 ) ;
    1260     tcalNotYet = True ;
    1261     tcalKey = -1 ;
    1262   }
    12631212
    12641213  Table &ms;
     
    13251274  MFrequency::Types freqframe;
    13261275  Record srcRec;
    1327   map< Int,Vector<uInt> > tcalIdRec;
    1328   map< Int,Vector<uInt> > tcalRowRec;
    1329   Int tcalKey;
    1330   Vector<Bool> tcalNotYet;
     1276};
     1277
     1278class BaseMSSysCalVisitor: public TableVisitor {
     1279  uInt lastRecordNo;
     1280  uInt lastBeamNo, lastIfNo, lastPolNo;
     1281  Double lastTime;
     1282protected:
     1283  const Table &table;
     1284  uInt count;
     1285public:
     1286  BaseMSSysCalVisitor(const Table &table)
     1287    : table(table)
     1288  {
     1289    count = 0;
     1290  }
     1291 
     1292  virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { }
     1293  virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { }
     1294  virtual void enterIfNo(const uInt recordNo, uInt columnValue) { }
     1295  virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { }
     1296  virtual void enterPolNo(const uInt recordNo, uInt columnValue) { }
     1297  virtual void leavePolNo(const uInt recordNo, uInt columnValue) { }
     1298  virtual void enterTime(const uInt recordNo, Double columnValue) { }
     1299  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     1300
     1301  virtual Bool visitRecord(const uInt recordNo,
     1302                           const uInt beamNo,
     1303                           const uInt ifNo,
     1304                           const uInt polNo,
     1305                           const Double time) { return True ;}
     1306
     1307  virtual Bool visit(Bool isFirst, const uInt recordNo,
     1308                     const uInt nCols, void const *const colValues[]) {
     1309    uInt beamNo, ifNo, polNo;
     1310    Double time;
     1311    { // prologue
     1312      uInt i = 0;
     1313      {
     1314        const uInt *col = (const uInt *)colValues[i++];
     1315        beamNo = col[recordNo];
     1316      }
     1317      {
     1318        const uInt *col = (const uInt *)colValues[i++];
     1319        ifNo = col[recordNo];
     1320      }
     1321      {
     1322        const Double *col = (const Double *)colValues[i++];
     1323        time = col[recordNo];
     1324      }
     1325      {
     1326        const uInt *col = (const uInt *)colValues[i++];
     1327        polNo = col[recordNo];
     1328      }
     1329      assert(nCols == i);
     1330    }
     1331
     1332    if (isFirst) {
     1333      enterBeamNo(recordNo, beamNo);
     1334      enterIfNo(recordNo, ifNo);
     1335      enterTime(recordNo, time);
     1336      enterPolNo(recordNo, polNo);
     1337    } else {
     1338      if (lastBeamNo != beamNo) {
     1339        leavePolNo(lastRecordNo, lastPolNo);
     1340        leaveTime(lastRecordNo, lastTime);
     1341        leaveIfNo(lastRecordNo, lastIfNo);
     1342        leaveBeamNo(lastRecordNo, lastBeamNo);
     1343
     1344        enterBeamNo(recordNo, beamNo);
     1345        enterIfNo(recordNo, ifNo);
     1346        enterTime(recordNo, time);
     1347        enterPolNo(recordNo, polNo);
     1348      } else if (lastIfNo != ifNo) {
     1349        leavePolNo(lastRecordNo, lastPolNo);
     1350        leaveTime(lastRecordNo, lastTime);
     1351        leaveIfNo(lastRecordNo, lastIfNo);
     1352       
     1353        enterIfNo(recordNo, ifNo);
     1354        enterTime(recordNo, time);
     1355        enterPolNo(recordNo, polNo);
     1356      } else if (lastTime != time) {
     1357        leavePolNo(lastRecordNo, lastPolNo);
     1358        leaveTime(lastRecordNo, lastTime);
     1359
     1360        enterTime(recordNo, time);
     1361        enterPolNo(recordNo, polNo);
     1362      } else if (lastPolNo != polNo) {
     1363        leavePolNo(lastRecordNo, lastPolNo);
     1364        enterPolNo(recordNo, polNo);
     1365      }
     1366    }
     1367    count++;
     1368    Bool result = visitRecord(recordNo, beamNo, ifNo, polNo, time);
     1369
     1370    { // epilogue
     1371      lastRecordNo = recordNo;
     1372
     1373      lastBeamNo = beamNo;
     1374      lastIfNo = ifNo;
     1375      lastPolNo = polNo;
     1376      lastTime = time;
     1377    }
     1378    return result ;
     1379  }
     1380
     1381  virtual void finish() {
     1382    if (count > 0) {
     1383      leavePolNo(lastRecordNo, lastPolNo);
     1384      leaveTime(lastRecordNo, lastTime);
     1385      leaveIfNo(lastRecordNo, lastIfNo);
     1386      leaveBeamNo(lastRecordNo, lastBeamNo);
     1387    }
     1388  }
     1389};
     1390
     1391class BaseTsysHolder
     1392{
     1393public:
     1394  BaseTsysHolder( ROArrayColumn<Float> &tsysCol )
     1395    : col( tsysCol ),
     1396      nchan(0)
     1397  {
     1398    reset() ;
     1399  }
     1400  virtual Array<Float> getTsys() = 0 ;
     1401  void setNchan( uInt n ) { nchan = n ; }
     1402  void appendTsys( uInt row )
     1403  {
     1404    Vector<Float> v = col( row ) ;
     1405    uInt len = tsys.nrow() ;
     1406    tsys.resize( len+1, nchan, True ) ;
     1407    if ( v.nelements() == nchan )
     1408      tsys.row( len ) = v ;
     1409    else
     1410      tsys.row( len ) = v[0] ;
     1411  }
     1412  void setTsys( uInt row, uInt idx )
     1413  {
     1414    if ( idx >= nrow() )
     1415      appendTsys( row ) ;
     1416    else {
     1417      Vector<Float> v = col( row ) ;
     1418      if ( v.nelements() == nchan )
     1419        tsys.row( idx ) = v ;
     1420      else
     1421        tsys.row( idx ) = v[0] ;
     1422    }
     1423  }
     1424  void reset()
     1425  {
     1426    tsys.resize() ;
     1427  }
     1428  uInt nrow() { return tsys.nrow() ; }
     1429  Bool isEffective()
     1430  {
     1431    return ( !(tsys.empty()) && anyNE( tsys, (Float)1.0 ) ) ;
     1432  }
     1433  BaseTsysHolder &operator= ( const BaseTsysHolder &v )
     1434  {
     1435    if ( this != &v )
     1436      tsys.assign( v.tsys ) ;
     1437    return *this ;
     1438  }
     1439protected:
     1440  ROArrayColumn<Float> col ;
     1441  Matrix<Float> tsys ;
     1442  uInt nchan ;
     1443};
     1444
     1445class TsysHolder : public BaseTsysHolder
     1446{
     1447public:
     1448  TsysHolder( ROArrayColumn<Float> &tsysCol )
     1449    : BaseTsysHolder( tsysCol )
     1450  {}
     1451  virtual Array<Float> getTsys()
     1452  {
     1453    return tsys.column( 0 ) ;
     1454  }
     1455};
     1456
     1457class TsysSpectrumHolder : public BaseTsysHolder
     1458{
     1459public:
     1460  TsysSpectrumHolder( ROArrayColumn<Float> &tsysCol )
     1461    : BaseTsysHolder( tsysCol )
     1462  {}
     1463  virtual Array<Float> getTsys()
     1464  {
     1465    return tsys ;
     1466  }
     1467};
     1468
     1469class BaseTcalProcessor
     1470{
     1471public:
     1472  BaseTcalProcessor( ROArrayColumn<Float> &tcalCol )
     1473    : col( tcalCol )
     1474  {}
     1475  void setTcalId( Vector<uInt> &tcalId ) { id.assign( tcalId ) ; }
     1476  virtual Array<Float> getTcal() = 0 ;
     1477protected:
     1478  ROArrayColumn<Float> col ;
     1479  Vector<uInt> id ;
     1480};
     1481
     1482class TcalProcessor : public BaseTcalProcessor
     1483{
     1484public:
     1485  TcalProcessor( ROArrayColumn<Float> &tcalCol )
     1486    : BaseTcalProcessor( tcalCol )
     1487  {}
     1488  virtual Array<Float> getTcal()
     1489  {
     1490    uInt npol = id.nelements() ;
     1491    Vector<Float> tcal( npol ) ;
     1492    for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
     1493      tcal[ipol] = col( id[ipol] ).data()[0] ;
     1494    //cout << "TcalProcessor: tcal = " << tcal << endl ;
     1495    return tcal ;
     1496  }
     1497};
     1498
     1499class TcalSpectrumProcessor : public BaseTcalProcessor
     1500{
     1501public:
     1502  TcalSpectrumProcessor( ROArrayColumn<Float> &tcalCol )
     1503    : BaseTcalProcessor( tcalCol )
     1504  {}
     1505  virtual Array<Float> getTcal()
     1506  {
     1507    uInt npol = id.nelements() ;
     1508    Vector<Float> tcal0 = col( 0 ) ;
     1509    uInt nchan = tcal0.nelements() ;
     1510    Matrix<Float> tcal( npol, nchan ) ;
     1511    tcal.row( 0 ) = tcal0 ;
     1512    for ( uInt ipol = 1 ; ipol < npol ; ipol++ )
     1513      tcal.row( ipol ) = col( id[ipol] ) ;
     1514    return tcal ;
     1515  }
     1516};
     1517
     1518class MSSysCalVisitor : public BaseMSSysCalVisitor
     1519{
     1520public:
     1521  MSSysCalVisitor( const Table &from, Table &to )
     1522    : BaseMSSysCalVisitor( from ),
     1523      sctab( to ),
     1524      rowidx( 0 )
     1525  {
     1526    scrow = TableRow( sctab ) ;
     1527
     1528    lastTcalId.resize() ;
     1529    theTcalId.resize() ;
     1530    startTime = 0.0 ;
     1531    endTime = 0.0 ;
     1532
     1533    const TableRecord &keys = table.keywordSet() ;
     1534    Table tcalTable = keys.asTable( "TCAL" ) ;
     1535    tcalCol.attach( tcalTable, "TCAL" ) ;
     1536    tsysCol.attach( table, "TSYS" ) ;
     1537    tcalIdCol.attach( table, "TCAL_ID" ) ;
     1538    intervalCol.attach( table, "INTERVAL" ) ;
     1539    effectiveTcal.resize( tcalTable.nrow() ) ;
     1540    for ( uInt irow = 0 ; irow < tcalTable.nrow() ; irow++ ) {
     1541      if ( allEQ( tcalCol( irow ), (Float)1.0 ) )
     1542        effectiveTcal[irow] = False ;
     1543      else
     1544        effectiveTcal[irow] = True ;
     1545    }
     1546   
     1547    TableRecord &r = scrow.record() ;
     1548    RecordFieldPtr<Int> antennaIdRF( r, "ANTENNA_ID" ) ;
     1549    *antennaIdRF = 0 ;
     1550    feedIdRF.attachToRecord( r, "FEED_ID" ) ;
     1551    specWinIdRF.attachToRecord( r, "SPECTRAL_WINDOW_ID" ) ;
     1552    timeRF.attachToRecord( r, "TIME" ) ;
     1553    intervalRF.attachToRecord( r, "INTERVAL" ) ;
     1554    if ( r.isDefined( "TCAL" ) ) {
     1555      tcalRF.attachToRecord( r, "TCAL" ) ;
     1556      tcalProcessor = new TcalProcessor( tcalCol ) ;
     1557    }
     1558    else if ( r.isDefined( "TCAL_SPECTRUM" ) ) {
     1559      tcalRF.attachToRecord( r, "TCAL_SPECTRUM" ) ;
     1560      tcalProcessor = new TcalSpectrumProcessor( tcalCol ) ;
     1561    }
     1562    if ( r.isDefined( "TSYS" ) ) {
     1563      tsysRF.attachToRecord( r, "TSYS" ) ;
     1564      theTsys = new TsysHolder( tsysCol ) ;
     1565      lastTsys = new TsysHolder( tsysCol ) ;
     1566    }
     1567    else {
     1568      tsysRF.attachToRecord( r, "TSYS_SPECTRUM" ) ;
     1569      theTsys = new TsysSpectrumHolder( tsysCol ) ;
     1570      lastTsys = new TsysSpectrumHolder( tsysCol ) ;
     1571    }
     1572
     1573  }
     1574
     1575  virtual void enterBeamNo(const uInt recordNo, uInt columnValue)
     1576  {
     1577    *feedIdRF = (Int)columnValue ;
     1578  }
     1579  virtual void leaveBeamNo(const uInt recordNo, uInt columnValue)
     1580  {
     1581  }
     1582  virtual void enterIfNo(const uInt recordNo, uInt columnValue)
     1583  {
     1584    //cout << "enterIfNo" << endl ;
     1585    ROArrayColumn<Float> sp( table, "SPECTRA" ) ;
     1586    uInt nchan = sp( recordNo ).nelements() ;
     1587    theTsys->setNchan( nchan ) ;
     1588    lastTsys->setNchan( nchan ) ;
     1589
     1590    *specWinIdRF = (Int)columnValue ;
     1591  }
     1592  virtual void leaveIfNo(const uInt recordNo, uInt columnValue)
     1593  {
     1594    //cout << "leaveIfNo" << endl ;
     1595    post() ;
     1596    lastTsys->reset() ;
     1597    lastTcalId.resize() ;
     1598    theTsys->reset() ;
     1599    theTcalId.resize() ;
     1600    startTime = 0.0 ;
     1601    endTime = 0.0 ;
     1602  }
     1603  virtual void enterTime(const uInt recordNo, Double columnValue)
     1604  {
     1605    //cout << "enterTime" << endl ;
     1606    interval = intervalCol.asdouble( recordNo ) ;
     1607    // start time and end time
     1608    if ( startTime == 0.0 ) {
     1609      startTime = columnValue * 86400.0 - 0.5 * interval ;
     1610      endTime = columnValue * 86400.0 + 0.5 * interval ;
     1611    }
     1612  }
     1613  virtual void leaveTime(const uInt recordNo, Double columnValue)
     1614  {
     1615    //cout << "leaveTime" << endl ;
     1616    if ( isUpdated() ) {
     1617      post() ;
     1618      *lastTsys = *theTsys ;
     1619      lastTcalId = theTcalId ;
     1620      theTsys->reset() ;
     1621      theTcalId.resize() ;
     1622      startTime = columnValue * 86400.0 - 0.5 * interval ;
     1623      endTime = columnValue * 86400.0 + 0.5 * interval ;
     1624    }
     1625    else {
     1626      endTime = columnValue * 86400.0 + 0.5 * interval ;
     1627    }
     1628  }
     1629  virtual void enterPolNo(const uInt recordNo, uInt columnValue)
     1630  {
     1631    //cout << "enterPolNo" << endl ;
     1632    Vector<Float> tsys = tsysCol( recordNo ) ;
     1633    uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
     1634    // lastTsys.nrow() must be npol
     1635    if ( lastTsys->nrow() == columnValue )
     1636      lastTsys->appendTsys( recordNo ) ;
     1637    // lastTcalId.nelements() must be npol
     1638    if ( lastTcalId.nelements() == columnValue )
     1639      appendTcalId( lastTcalId, tcalId, columnValue ) ;
     1640    // theTsys.nrow() must be npol
     1641    if ( theTsys->nrow() == columnValue )
     1642      theTsys->appendTsys( recordNo ) ;
     1643    else {
     1644      theTsys->setTsys( recordNo, columnValue ) ;
     1645    }
     1646    if ( theTcalId.nelements() == columnValue )
     1647      appendTcalId( theTcalId, tcalId, columnValue ) ;
     1648    else
     1649      setTcalId( theTcalId, tcalId, columnValue ) ;
     1650  }
     1651  virtual void leavePolNo( const uInt recordNo, uInt columnValue )
     1652  {
     1653  }
     1654   
     1655private:
     1656  void appendTcalId( Vector<uInt> &v, uInt &elem, uInt &polId )
     1657  {
     1658    v.resize( polId+1, True ) ;
     1659    v[polId] = elem ;
     1660  }
     1661  void setTcalId( Vector<uInt> &v, uInt &elem, uInt &polId )
     1662  {
     1663    v[polId] = elem ;
     1664  }
     1665  void post()
     1666  {
     1667    // check if given Tcal and Tsys is effective
     1668    Bool isEffective = False ;
     1669    for ( uInt ipol = 0 ; ipol < lastTcalId.nelements() ; ipol++ ) {
     1670      if ( effectiveTcal[lastTcalId[ipol]] ) {
     1671        isEffective = True ;
     1672        break ;
     1673      }
     1674    }
     1675    if ( !isEffective ) {
     1676      if ( !(lastTsys->isEffective()) )
     1677        return ;
     1678    }
     1679
     1680    //cout << " interval: " << (endTime-startTime) << " lastTcalId = " << lastTcalId << endl ;
     1681    Double midTime = 0.5 * ( startTime + endTime ) ;
     1682    Double interval = endTime - startTime ;
     1683    *timeRF = midTime ;
     1684    *intervalRF = interval ;
     1685    tcalProcessor->setTcalId( lastTcalId ) ;
     1686    Array<Float> tcal = tcalProcessor->getTcal() ;
     1687    tcalRF.define( tcal ) ;
     1688    tsysRF.define( lastTsys->getTsys() ) ;
     1689    sctab.addRow( 1, True ) ;
     1690    scrow.put( rowidx ) ;
     1691    rowidx++ ;
     1692  }
     1693 
     1694  Bool isUpdated()
     1695  {
     1696    Bool ret = False ;
     1697    ret = anyNE( theTcalId, lastTcalId ) ;
     1698    if ( !ret )
     1699      ret = anyNE( theTsys->getTsys(), lastTsys->getTsys() ) ;
     1700    return ret ;
     1701  }
     1702
     1703  Table &sctab;
     1704  TableRow scrow;
     1705  uInt rowidx;
     1706
     1707  Double startTime,endTime,interval;
     1708 
     1709  CountedPtr<BaseTsysHolder> lastTsys,theTsys;
     1710  Vector<uInt> lastTcalId,theTcalId;
     1711  CountedPtr<BaseTcalProcessor> tcalProcessor ;
     1712  Vector<Bool> effectiveTcal;
     1713
     1714  RecordFieldPtr<Int> feedIdRF,specWinIdRF;
     1715  RecordFieldPtr<Double> timeRF,intervalRF;
     1716  RecordFieldPtr< Array<Float> > tcalRF,tsysRF;
     1717
     1718  ROArrayColumn<Float> tsysCol,tcalCol;
     1719  ROTableColumn tcalIdCol,intervalCol;
    13311720};
    13321721
    13331722MSWriter::MSWriter(CountedPtr<Scantable> stable)
    13341723  : table_(stable),
    1335     isTcal_(False),
    13361724    isWeather_(False),
    13371725    tcalSpec_(False),
     
    14091797  if ( isWeather_ )
    14101798    fillWeather() ;
     1799
     1800  // SYSCAL
     1801  fillSysCal() ;
    14111802
    14121803  /***
     
    14241815    static const TypeManagerImpl<String> tmString;
    14251816    static const TypeManager *const tms[] = {
    1426       &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmInt, NULL
     1817      &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL
    14271818    };
    14281819    //double t0 = mathutil::gettimeofday_sec() ;
     
    14401831    //double t3 = mathutil::gettimeofday_sec() ;
    14411832    //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
    1442     map< Int,Vector<uInt> > &idRec = myVisitor.getTcalIdRecord() ;
    1443     map< Int,Vector<uInt> > &rowRec = myVisitor.getTcalRowRecord() ;
    1444 
    1445     // SYSCAL
    1446     if ( isTcal_ )
    1447       fillSysCal( idRec, rowRec ) ;
    14481833  }
    14491834  /***
     
    15181903
    15191904  // Check if some subtables are exists
     1905  Bool isTcal = False ;
    15201906  if ( table_->tcal().table().nrow() != 0 ) {
    15211907    ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
    15221908    if ( col.isDefined( 0 ) ) {
    15231909      os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
    1524       isTcal_ = True ;
     1910      isTcal = True ;
    15251911    }
    15261912    else {
     
    15461932
    15471933  // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
    1548   if ( isTcal_ && header_.nchan != 1 ) {
    1549     // examine TCAL subtable
    1550     Table tcaltab = table_->tcal().table() ;
    1551     ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
    1552     for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
    1553       if ( tcalCol( irow ).size() != 1 )
    1554         tcalSpec_ = True ;
     1934  if ( header_.nchan != 1 ) {
     1935    if ( isTcal ) {
     1936      // examine TCAL subtable
     1937      Table tcaltab = table_->tcal().table() ;
     1938      ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
     1939      for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
     1940        if ( tcalCol( irow ).size() != 1 )
     1941          tcalSpec_ = True ;
     1942      }
    15551943    }
    15561944    // examine spectral data
     
    16712059  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
    16722060
    1673   // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
    16742061  TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
    1675   MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 2 ) ;
    1676   MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 2 ) ;
    16772062  if ( tcalSpec_ )
    16782063    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
     2064  else
     2065    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ;
    16792066  if ( tsysSpec_ )
    16802067    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
     2068  else
     2069    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ;
    16812070  SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
    16822071  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
     
    20352424}
    20362425
    2037 void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec,
    2038                            map< Int,Vector<uInt> > &rowRec )
     2426void MSWriter::fillSysCal()
    20392427{
    2040   //double startSec = mathutil::gettimeofday_sec() ;
    2041   //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
    2042 
    2043   //idRec.print( cout ) ;
    2044 
    2045   // access to MS SYSCAL subtable
    2046   MSSysCal mssc = mstable_->sysCal() ;
    2047 
    2048   // access to TCAL subtable
    2049   Table stt = table_->tcal().table() ;
    2050   uInt nrow = stt.nrow() ;
    2051 
    2052   // access to MAIN table
    2053   Block<String> cols( 6 ) ;
    2054   cols[0] = "TIME" ;
    2055   cols[1] = "TCAL_ID" ;
    2056   cols[2] = "TSYS" ;
    2057   cols[3] = "BEAMNO" ;
    2058   cols[4] = "IFNO" ;
    2059   cols[5] = "INTERVAL" ;
    2060   Table tab = table_->table().project( cols ) ;
    2061 
    2062   if ( nrow == 0 )
    2063     return ;
    2064 
    2065   //nrow = idRec.nfields() ;
    2066   nrow = idRec.size() ;
    2067 
    2068   Double midTime ;
    2069   Double interval ;
    2070   String timeStr ;
    2071 
    2072   // row base
    2073   TableRow row( mssc ) ;
    2074   TableRecord &trec = row.record() ;
    2075   RecordFieldPtr<Int> antennaRF( trec, "ANTENNA_ID" ) ;
    2076   RecordFieldPtr<Int> feedRF( trec, "FEED_ID" ) ;
    2077   RecordFieldPtr<Int> spwRF( trec, "SPECTRAL_WINDOW_ID" ) ;
    2078   RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
    2079   RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
    2080   RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
    2081   RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
    2082   RecordFieldPtr< Array<Float> > tsysspRF ;
    2083   RecordFieldPtr< Array<Float> > tcalspRF ;
    2084   if ( tsysSpec_ )
    2085     tsysspRF.attachToRecord( trec, "TSYS_SPECTRUM" ) ;
    2086   if ( tcalSpec_ )
    2087     tcalspRF.attachToRecord( trec, "TCAL_SPECTRUM" ) ;
    2088 
    2089   // ANTENNA_ID is always 0
    2090   *antennaRF = 0 ;
    2091 
    2092   Table sortedstt = stt.sort( "ID" ) ;
    2093   ROArrayColumn<Float> tcalCol( sortedstt, "TCAL" ) ;
    2094   ROTableColumn idCol( sortedstt, "ID" ) ;
    2095   ROArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
    2096   ROTableColumn tcalidCol( tab, "TCAL_ID" ) ;
    2097   ROTableColumn timeCol( tab, "TIME" ) ;
    2098   ROTableColumn intervalCol( tab, "INTERVAL" ) ;
    2099   ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
    2100   ROTableColumn ifnoCol( tab, "IFNO" ) ;
    2101   map< Int,Vector<uInt> >::iterator itr0 = idRec.begin() ;
    2102   map< Int,Vector<uInt> >::iterator itr1 = rowRec.begin() ;
    2103   for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
    2104 //     double t1 = mathutil::gettimeofday_sec() ;
    2105     Vector<uInt> ids = itr0->second ;
    2106     itr0++ ;
    2107 //     os_ << "ids = " << ids << LogIO::POST ;
    2108     uInt npol = ids.size() ;
    2109     Vector<uInt> rows = itr1->second ;
    2110     itr1++ ;
    2111 //     os_ << "rows = " << rows << LogIO::POST ;
    2112     Vector<Double> atime( rows.nelements() ) ;
    2113     Vector<Double> ainterval( rows.nelements() ) ;
    2114     Vector<uInt> atcalid( rows.nelements() ) ;
    2115     for( uInt jrow = 0 ; jrow < rows.nelements() ; jrow++ ) {
    2116       atime[jrow] = (Double)timeCol.asdouble( rows[jrow] ) ;
    2117       ainterval[jrow] = (Double)intervalCol.asdouble( rows[jrow] ) ;
    2118       atcalid[jrow] = tcalidCol.asuInt( rows[jrow] ) ;
    2119     }
    2120     Vector<Float> dummy = tsysCol( rows[0] ) ;
    2121     Matrix<Float> tsys( npol,dummy.nelements() ) ;
    2122     tsys.row( 0 ) = dummy ;
    2123     for ( uInt jrow = 1 ; jrow < npol ; jrow++ )
    2124       tsys.row( jrow ) = tsysCol( rows[jrow] ) ;
    2125 
    2126     // FEED_ID
    2127     *feedRF = beamnoCol.asuInt( rows[0] ) ;
    2128 
    2129     // SPECTRAL_WINDOW_ID
    2130     *spwRF = ifnoCol.asuInt( rows[0] ) ;
    2131 
    2132     // TIME and INTERVAL
    2133     getValidTimeRange( midTime, interval, atime, ainterval ) ;
    2134     *timeRF = midTime ;
    2135     *intervalRF = interval ;
    2136 
    2137     // TCAL and TSYS
    2138     Matrix<Float> tcal ;
    2139     Table t ;
    2140     if ( idCol.asuInt( ids[0] ) == ids[0] ) {
    2141 //       os_ << "sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
    2142       Vector<Float> dummyC = tcalCol( ids[0] ) ;
    2143       tcal.resize( npol, dummyC.size() ) ;
    2144       tcal.row( 0 ) = dummyC ;
    2145     }
    2146     else {
    2147 //       os_ << "NOT sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
    2148       t = stt( stt.col("ID") == ids[0], 1 ) ;
    2149       Vector<Float> dummyC = tcalCol( 0 ) ;
    2150       tcal.resize( npol, dummyC.size(), True ) ;
    2151       tcal.row( 0 ) = dummyC ;
    2152     }
    2153     if ( npol == 2 ) {
    2154       if ( idCol.asuInt( ids[1] ) == ids[1] ) {
    2155 //         os_ << "sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
    2156         tcal.row( 1 ) = tcalCol( ids[1] ) ;
    2157       }
    2158       else {
    2159 //         os_ << "NOT sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
    2160         t = stt( stt.col("ID") == ids[1], 1 ) ;
    2161         tcalCol.attach( t, "TCAL" ) ;
    2162         tcal.row( 1 ) = tcalCol( 0 ) ;
    2163       }
    2164     }
    2165     else if ( npol == 3 ) {
    2166       if ( idCol.asuInt( ids[2] ) == ids[2] )
    2167         tcal.row( 1 ) = tcalCol( ids[2] ) ;
    2168       else {
    2169         t = stt( stt.col("ID") == ids[2], 1 ) ;
    2170         tcalCol.attach( t, "TCAL" ) ;
    2171         tcal.row( 1 ) = tcalCol( 0 ) ;
    2172       }
    2173       if ( idCol.asuInt( ids[1] ) == ids[1] )
    2174         tcal.row( 2 ) = tcalCol( ids[1] ) ;
    2175       else {
    2176         t = stt( stt.col("ID") == ids[1], 1 ) ;
    2177         tcalCol.attach( t, "TCAL" ) ;
    2178         tcal.row( 2 ) = tcalCol( 0 ) ;
    2179       }
    2180     }
    2181     else if ( npol == 4 ) {
    2182       if ( idCol.asuInt( ids[2] ) == ids[2] )
    2183         tcal.row( 1 ) = tcalCol( ids[2] ) ;
    2184       else {
    2185         t = stt( stt.col("ID") == ids[2], 1 ) ;
    2186         tcalCol.attach( t, "TCAL" ) ;
    2187         tcal.row( 1 ) = tcalCol( 0 ) ;
    2188       }
    2189       if ( idCol.asuInt( ids[3] ) == ids[3] )
    2190         tcal.row( 2 ) = tcalCol( ids[3] ) ;
    2191       else {
    2192         t = stt( stt.col("ID") == ids[3], 1 ) ;
    2193         tcalCol.attach( t, "TCAL" ) ;
    2194         tcal.row( 2 ) = tcalCol( 0 ) ;
    2195       }
    2196       if ( idCol.asuInt( ids[1] ) == ids[1] )
    2197         tcal.row( 2 ) = tcalCol( ids[1] ) ;
    2198       else {
    2199         t = stt( stt.col("ID") == ids[1], 1 ) ;
    2200         tcalCol.attach( t, "TCAL" ) ;
    2201         tcal.row( 3 ) = tcalCol( 0 ) ;
    2202       }
    2203     }
    2204     if ( tcalSpec_ ) {
    2205       // put TCAL_SPECTRUM
    2206       tcalspRF.define( tcal ) ;
    2207       // set TCAL (mean of TCAL_SPECTRUM)
    2208       Matrix<Float> tcalMean( npol, 1 ) ;
    2209       for ( uInt iid = 0 ; iid < npol ; iid++ ) {
    2210         tcalMean( iid, 0 ) = mean( tcal.row(iid) ) ;
    2211       }
    2212       // put TCAL
    2213       tcalRF.define( tcalMean ) ;
    2214     }
    2215     else {
    2216       // put TCAL
    2217       tcalRF.define( tcal ) ;
    2218     }
    2219    
    2220     if ( tsysSpec_ ) {
    2221       // put TSYS_SPECTRUM
    2222       tsysspRF.define( tsys ) ;
    2223       // set TSYS (mean of TSYS_SPECTRUM)
    2224       Matrix<Float> tsysMean( npol, 1 ) ;
    2225       for ( uInt iid = 0 ; iid < npol ; iid++ ) {
    2226         tsysMean( iid, 0 ) = mean( tsys.row(iid) ) ;
    2227       }
    2228       // put TSYS
    2229       tsysRF.define( tsysMean ) ;
    2230     }
    2231     else {
    2232       // put TSYS
    2233       tsysRF.define( tsys ) ;
    2234     }
    2235 
    2236     // add row
    2237     mssc.addRow( 1, True ) ;
    2238     row.put( mssc.nrow()-1 ) ;
    2239 
    2240 //     double t2 = mathutil::gettimeofday_sec() ;
    2241 //     os_ << irow << "th loop elapsed time = " << t2-t1 << "sec" << LogIO::POST ;
    2242   }
    2243  
    2244   //double endSec = mathutil::gettimeofday_sec() ;
    2245   //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     2428  Table mssc = mstable_->sysCal() ;
     2429
     2430  {
     2431    static const char *cols[] = {
     2432      "BEAMNO", "IFNO", "TIME", "POLNO",
     2433      NULL
     2434    };
     2435    static const TypeManagerImpl<uInt> tmUInt;
     2436    static const TypeManagerImpl<Double> tmDouble;
     2437    static const TypeManager *const tms[] = {
     2438      &tmUInt, &tmUInt, &tmDouble, &tmUInt, NULL
     2439    };
     2440    //double t0 = mathutil::gettimeofday_sec() ;
     2441    MSSysCalVisitor myVisitor(table_->table(),mssc);
     2442    //double t1 = mathutil::gettimeofday_sec() ;
     2443    //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
     2444    traverseTable(table_->table(), cols, tms, &myVisitor);
     2445    //double t3 = mathutil::gettimeofday_sec() ;
     2446    //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
     2447  }
     2448 
    22462449}
    22472450
  • branches/casa-prerelease/pre-asap/src/MSWriter.h

    r2293 r2310  
    8484  void fillSource() ;
    8585  void fillWeather() ;
    86   void fillSysCal( std::map< casa::Int,casa::Vector<casa::uInt> > &idrec,
    87                    std::map< casa::Int,casa::Vector<casa::uInt> > &rowrec ) ;
     86  void fillSysCal() ;
    8887
    8988  // utility
     
    9291  void antennaProperty( casa::String &name, casa::String &mount, casa::String &type, casa::Double &diameter ) ;
    9392
    94   // tool for HPC
    95 //   double gettimeofday_sec() ;
    96 
    9793  casa::CountedPtr<Scantable> table_ ;
    9894  STHeader header_ ;
    9995  casa::MeasurementSet *mstable_ ;
    10096
    101   casa::Bool isTcal_ ;
    10297  casa::Bool isWeather_ ;
    10398
  • branches/casa-prerelease/pre-asap/src/SConscript

    • Property svn:mergeinfo changed (with no actual effect on merging)
Note: See TracChangeset for help on using the changeset viewer.