Changeset 2310
- Timestamp:
- 09/27/11 14:20:58 (13 years ago)
- Location:
- branches/casa-prerelease/pre-asap
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/casa-prerelease/pre-asap
-
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
- Property svn:mergeinfo changed
/trunk/src merged: 2309
- Property svn:mergeinfo changed
-
branches/casa-prerelease/pre-asap/src/MSWriter.cpp
r2293 r2310 180 180 const String *lastFieldName; 181 181 uInt lastRecordNo; 182 uInt lastBeamNo, lastScanNo, lastIfNo ;182 uInt lastBeamNo, lastScanNo, lastIfNo, lastPolNo; 183 183 Int lastSrcType; 184 184 uInt lastCycleNo; 185 185 Double lastTime; 186 Int lastPolNo;187 186 protected: 188 187 const Table &table; … … 233 232 uInt cycleNo; 234 233 Double time; 235 Int polNo;234 uInt polNo; 236 235 { // prologue 237 236 uInt i = 0; … … 423 422 makePolMap() ; 424 423 initFrequencies() ; 425 initTcal() ;426 424 427 425 // … … 590 588 holder.reset() ; 591 589 } 592 if ( tcalKey != -1 ) {593 tcalNotYet[tcalKey] = False ;594 tcalKey = -1 ;595 }596 590 } 597 591 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { 598 592 //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 }631 593 } 632 594 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { … … 690 652 srcRec = r ; 691 653 } 692 map< Int,Vector<uInt> > &getTcalIdRecord() { return tcalIdRec ; }693 map< Int,Vector<uInt> > &getTcalRowRecord() { return tcalRowRec ; }694 654 private: 695 655 void addField( Int &fid, String &fname, String &srcName, … … 1250 1210 } 1251 1211 } 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 }1263 1212 1264 1213 Table &ms; … … 1325 1274 MFrequency::Types freqframe; 1326 1275 Record srcRec; 1327 map< Int,Vector<uInt> > tcalIdRec; 1328 map< Int,Vector<uInt> > tcalRowRec; 1329 Int tcalKey; 1330 Vector<Bool> tcalNotYet; 1276 }; 1277 1278 class BaseMSSysCalVisitor: public TableVisitor { 1279 uInt lastRecordNo; 1280 uInt lastBeamNo, lastIfNo, lastPolNo; 1281 Double lastTime; 1282 protected: 1283 const Table &table; 1284 uInt count; 1285 public: 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 1391 class BaseTsysHolder 1392 { 1393 public: 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 } 1439 protected: 1440 ROArrayColumn<Float> col ; 1441 Matrix<Float> tsys ; 1442 uInt nchan ; 1443 }; 1444 1445 class TsysHolder : public BaseTsysHolder 1446 { 1447 public: 1448 TsysHolder( ROArrayColumn<Float> &tsysCol ) 1449 : BaseTsysHolder( tsysCol ) 1450 {} 1451 virtual Array<Float> getTsys() 1452 { 1453 return tsys.column( 0 ) ; 1454 } 1455 }; 1456 1457 class TsysSpectrumHolder : public BaseTsysHolder 1458 { 1459 public: 1460 TsysSpectrumHolder( ROArrayColumn<Float> &tsysCol ) 1461 : BaseTsysHolder( tsysCol ) 1462 {} 1463 virtual Array<Float> getTsys() 1464 { 1465 return tsys ; 1466 } 1467 }; 1468 1469 class BaseTcalProcessor 1470 { 1471 public: 1472 BaseTcalProcessor( ROArrayColumn<Float> &tcalCol ) 1473 : col( tcalCol ) 1474 {} 1475 void setTcalId( Vector<uInt> &tcalId ) { id.assign( tcalId ) ; } 1476 virtual Array<Float> getTcal() = 0 ; 1477 protected: 1478 ROArrayColumn<Float> col ; 1479 Vector<uInt> id ; 1480 }; 1481 1482 class TcalProcessor : public BaseTcalProcessor 1483 { 1484 public: 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 1499 class TcalSpectrumProcessor : public BaseTcalProcessor 1500 { 1501 public: 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 1518 class MSSysCalVisitor : public BaseMSSysCalVisitor 1519 { 1520 public: 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 1655 private: 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; 1331 1720 }; 1332 1721 1333 1722 MSWriter::MSWriter(CountedPtr<Scantable> stable) 1334 1723 : table_(stable), 1335 isTcal_(False),1336 1724 isWeather_(False), 1337 1725 tcalSpec_(False), … … 1409 1797 if ( isWeather_ ) 1410 1798 fillWeather() ; 1799 1800 // SYSCAL 1801 fillSysCal() ; 1411 1802 1412 1803 /*** … … 1424 1815 static const TypeManagerImpl<String> tmString; 1425 1816 static const TypeManager *const tms[] = { 1426 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tm Int, NULL1817 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL 1427 1818 }; 1428 1819 //double t0 = mathutil::gettimeofday_sec() ; … … 1440 1831 //double t3 = mathutil::gettimeofday_sec() ; 1441 1832 //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 // SYSCAL1446 if ( isTcal_ )1447 fillSysCal( idRec, rowRec ) ;1448 1833 } 1449 1834 /*** … … 1518 1903 1519 1904 // Check if some subtables are exists 1905 Bool isTcal = False ; 1520 1906 if ( table_->tcal().table().nrow() != 0 ) { 1521 1907 ROTableColumn col( table_->tcal().table(), "TCAL" ) ; 1522 1908 if ( col.isDefined( 0 ) ) { 1523 1909 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ; 1524 isTcal _= True ;1910 isTcal = True ; 1525 1911 } 1526 1912 else { … … 1546 1932 1547 1933 // 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 } 1555 1943 } 1556 1944 // examine spectral data … … 1671 2059 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ; 1672 2060 1673 // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary1674 2061 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ; 1675 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 2 ) ;1676 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 2 ) ;1677 2062 if ( tcalSpec_ ) 1678 2063 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ; 2064 else 2065 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ; 1679 2066 if ( tsysSpec_ ) 1680 2067 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ; 2068 else 2069 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ; 1681 2070 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ; 1682 2071 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ; … … 2035 2424 } 2036 2425 2037 void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec, 2038 map< Int,Vector<uInt> > &rowRec ) 2426 void MSWriter::fillSysCal() 2039 2427 { 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 2246 2449 } 2247 2450 -
branches/casa-prerelease/pre-asap/src/MSWriter.h
r2293 r2310 84 84 void fillSource() ; 85 85 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() ; 88 87 89 88 // utility … … 92 91 void antennaProperty( casa::String &name, casa::String &mount, casa::String &type, casa::Double &diameter ) ; 93 92 94 // tool for HPC95 // double gettimeofday_sec() ;96 97 93 casa::CountedPtr<Scantable> table_ ; 98 94 STHeader header_ ; 99 95 casa::MeasurementSet *mstable_ ; 100 96 101 casa::Bool isTcal_ ;102 97 casa::Bool isWeather_ ; 103 98 -
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.