Changeset 1789 for branches/mergetest/src
- Timestamp:
- 07/30/10 16:59:56 (14 years ago)
- Location:
- branches/mergetest/src
- Files:
-
- 30 edited
- 16 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/mergetest/src
-
Property svn:mergeinfo
set to
/branches/alma/src merged eligible
-
Property svn:mergeinfo
set to
-
branches/mergetest/src/MathUtils.cpp
r1570 r1789 38 38 #include <casa/BasicSL/String.h> 39 39 #include <scimath/Mathematics/MedianSlider.h> 40 #include <casa/Exceptions/Error.h> 40 41 41 42 #include <scimath/Fitting/LinearFit.h> … … 53 54 String str(which); 54 55 str.upcase(); 55 if (str. contains(String("MIN"))) {56 if (str.matches(String("MIN"))) { 56 57 return min(data); 57 } else if (str. contains(String("MAX"))) {58 } else if (str.matches(String("MAX"))) { 58 59 return max(data); 59 } else if (str. contains(String("SUMSQ"))) {60 } else if (str.matches(String("SUMSQ"))) { 60 61 return sumsquares(data); 61 } else if (str. contains(String("SUM"))) {62 } else if (str.matches(String("SUM"))) { 62 63 return sum(data); 63 } else if (str. contains(String("MEAN"))) {64 } else if (str.matches(String("MEAN"))) { 64 65 return mean(data); 65 } else if (str. contains(String("VAR"))) {66 } else if (str.matches(String("VAR"))) { 66 67 return variance(data); 67 } else if (str.contains(String("STDDEV"))) {68 } else if (str.matches(String("STDDEV"))) { 68 69 return stddev(data); 69 } else if (str. contains(String("AVDEV"))) {70 } else if (str.matches(String("AVDEV"))) { 70 71 return avdev(data); 71 } else if (str. contains(String("RMS"))) {72 } else if (str.matches(String("RMS"))) { 72 73 uInt n = data.nelementsValid(); 73 74 return sqrt(sumsquares(data)/n); 74 } else if (str. contains(String("MED"))) {75 } else if (str.matches(String("MEDIAN"))) { 75 76 return median(data); 76 } 77 } else { 78 String msg = str + " is not a valid type of statistics"; 79 throw(AipsError(msg)); 80 } 77 81 return 0.0; 78 82 } 79 83 84 IPosition mathutil::minMaxPos(const String& which, 85 const MaskedArray<Float>& data) 86 { 87 Float minVal, maxVal; 88 IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0); 89 minMax(minVal, maxVal, minPos, maxPos, data); 90 String str(which); 91 str.upcase(); 92 if (str.contains(String("MIN"))) { 93 return minPos; 94 } else if (str.contains(String("MAX"))) { 95 return maxPos; 96 } else { 97 String msg = str + " is not a valid type of statistics"; 98 throw(AipsError(msg)); 99 } 100 //return 0.0; 101 } 80 102 81 103 void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask) -
branches/mergetest/src/MathUtils.h
r1570 r1789 37 37 #include <casa/Arrays/Vector.h> 38 38 #include <casa/BasicSL/String.h> 39 #include <casa/Arrays/IPosition.h> 39 40 40 41 namespace mathutil { … … 79 80 float hwidth, int order); 80 81 81 82 82 // Generate specified statistic 83 83 float statistics(const casa::String& which, 84 const casa::MaskedArray<casa::Float>& data); 85 86 // Return a position of min or max value 87 casa::IPosition minMaxPos(const casa::String& which, 84 88 const casa::MaskedArray<casa::Float>& data); 85 89 -
branches/mergetest/src/RowAccumulator.cpp
r1569 r1789 152 152 userMask_ = m; 153 153 } 154 155 // Added by TT check the state of RowAccumulator 156 casa::Bool RowAccumulator::state() const 157 { 158 return initialized_; 159 } 160 -
branches/mergetest/src/RowAccumulator.h
r1569 r1789 85 85 */ 86 86 void reset(); 87 /** 88 * check the initialization state 89 */ 90 casa::Bool state() const; 87 91 88 92 private: -
branches/mergetest/src/SConscript
-
Property svn:mergeinfo
set to
/branches/alma/src/SConscript merged eligible /branches/asap-3.x/src/SConscript merged eligible
-
Property svn:mergeinfo
set to
-
branches/mergetest/src/STAsciiWriter.cpp
r1552 r1789 1 1 2 //#--------------------------------------------------------------------------- 2 3 //# STAsciiWriter.cc: ASAP class to write out single dish spectra as FITS images … … 88 89 89 90 String rootName(fileName); 90 91 91 92 Block<String> cols(4); 92 93 cols[0] = String("SCANNO"); … … 132 133 String wcs = stable.frequencies().print(rec.asuInt("FREQ_ID"), True); 133 134 addLine(of, "WCS", wcs); 134 addLine(of, "Rest Freq.", 135 stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") )); 135 std::vector<double> restfreqs= stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID")); 136 int nf = restfreqs.size(); 137 //addLine(of, "Rest Freq.", 138 // stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") )); 139 addLine(of, "Rest Freq.", restfreqs[0]); 140 for ( unsigned int i=1; i<nf; ++i) { 141 addLine(of, " ", restfreqs[i]); 142 } 143 ostringstream osflagrow; 144 for ( unsigned int i=0; i<t.nrow(); ++i) { 145 osflagrow << "Pol" << i << ":" << ((row.get(i).asuInt("FLAGROW") > 0) ? "True" : "False") << " "; 146 } 147 addLine(of, "Row_Flagged", String(osflagrow)); 136 148 of << setfill('#') << setw(70) << "" << setfill(' ') << endl; 137 149 -
branches/mergetest/src/STFiller.cpp
r1725 r1789 25 25 #include <tables/Tables/TableRow.h> 26 26 27 #include <measures/Measures/MDirection.h> 28 #include <measures/Measures/MeasConvert.h> 29 27 30 #include <atnf/PKSIO/PKSrecord.h> 28 31 #include <atnf/PKSIO/PKSreader.h> … … 30 33 #include <casa/System/ProgressMeter.h> 31 34 #endif 35 #include <casa/System/ProgressMeter.h> 36 #include <atnf/PKSIO/NROReader.h> 37 #include <casa/Logging/LogIO.h> 38 39 #include <time.h> 40 32 41 33 42 #include "STDefs.h" … … 45 54 header_(0), 46 55 table_(0), 47 refRx_(".*(e|w|_R)$") 56 refRx_(".*(e|w|_R)$"), 57 nreader_(0) 48 58 { 49 59 } … … 53 63 header_(0), 54 64 table_(stbl), 55 refRx_(".*(e|w|_R)$") 65 refRx_(".*(e|w|_R)$"), 66 nreader_(0) 56 67 { 57 68 } … … 61 72 header_(0), 62 73 table_(0), 63 refRx_(".*(e|w|_R)$") 74 refRx_(".*(e|w|_R)$"), 75 nreader_(0) 64 76 { 65 open(filename, whichIF, whichBeam);77 open(filename, "", whichIF, whichBeam); 66 78 } 67 79 … … 70 82 close(); 71 83 } 72 73 void STFiller::open( const std::string& filename, int whichIF, int whichBeam)84 85 void STFiller::open( const std::string& filename, const std::string& antenna, int whichIF, int whichBeam, casa::Bool getPt ) 74 86 { 75 87 if (table_.null()) { … … 93 105 Vector<Bool> beams, ifs; 94 106 Vector<uInt> nchans,npols; 95 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs, 107 108 // 109 // if isNRO_ is true, try NROReader 110 // 111 // 2008/11/11 Takeshi Nakazato 112 isNRO_ = fileCheck() ; 113 if ( isNRO_ ) { 114 if ( (nreader_ = getNROReader( inName, format )) == 0 ) { 115 throw(AipsError("Creation of NROReader failed")) ; 116 } 117 else { 118 openNRO( whichIF, whichBeam ) ; 119 return ; 120 } 121 } 122 // 123 124 if ( (reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs, 96 125 nchans, npols, haveXPol_,haveBase, haveSpectra 97 126 )) == 0 ) { … … 118 147 header_->npol = max(npols); 119 148 header_->nbeam = nBeam_; 120 149 121 150 Int status = reader_->getHeader(header_->observer, header_->project, 122 151 header_->antennaname, header_->antennaposition, … … 197 226 Vector<Int> start(nIF_, 1); 198 227 Vector<Int> end(nIF_, 0); 199 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0] );228 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt); 200 229 table_->setHeader(*header_); 201 230 //For MS, add the location of POINTING of the input MS so one get … … 216 245 } 217 246 } 247 String freqFrame = header_->freqref; 248 //translate frequency reference frame back to 249 //MS style (as PKSMS2reader converts the original frame 250 //in FITS standard style) 251 if (freqFrame == "TOPOCENT") { 252 freqFrame = "TOPO"; 253 } else if (freqFrame == "GEOCENER") { 254 freqFrame = "GEO"; 255 } else if (freqFrame == "BARYCENT") { 256 freqFrame = "BARY"; 257 } else if (freqFrame == "GALACTOC") { 258 freqFrame = "GALACTO"; 259 } else if (freqFrame == "LOCALGRP") { 260 freqFrame = "LGROUP"; 261 } else if (freqFrame == "CMBDIPOL") { 262 freqFrame = "CMB"; 263 } else if (freqFrame == "SOURCE") { 264 freqFrame = "REST"; 265 } 266 // set both "FRAME" and "BASEFRAME" 267 table_->frequencies().setFrame(freqFrame, false); 268 table_->frequencies().setFrame(freqFrame,true); 218 269 //table_->focus().setParallactify(true); 219 270 } … … 222 273 { 223 274 delete reader_;reader_=0; 275 delete nreader_;nreader_=0; 224 276 delete header_;header_=0; 225 277 table_ = 0; … … 229 281 { 230 282 int status = 0; 283 284 // 285 // for NRO data 286 // 287 // 2008/11/12 Takeshi Nakazato 288 if ( isNRO_ ) { 289 status = readNRO() ; 290 return status ; 291 } 292 // 293 294 /** 295 Int beamNo, IFno, refBeam, scanNo, cycleNo; 296 Float azimuth, elevation, focusAxi, focusRot, focusTan, 297 humidity, parAngle, pressure, temperature, windAz, windSpeed; 298 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel; 299 String fieldName, srcName, tcalTime, obsType; 300 Vector<Float> calFctr, sigma, tcal, tsys; 301 Matrix<Float> baseLin, baseSub; 302 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1); 303 Matrix<Float> spectra; 304 Matrix<uChar> flagtra; 305 Complex xCalFctr; 306 Vector<Complex> xPol; 307 **/ 231 308 232 309 Double min = 0.0; … … 236 313 #endif 237 314 PKSrecord pksrec; 315 pksrec.srcType=-1; 238 316 int n = 0; 317 bool isGBTFITS = false ; 318 if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) { 319 FILE *fp = fopen( filename_.c_str(), "r" ) ; 320 fseek( fp, 640, SEEK_SET ) ; 321 char buf[81] ; 322 fread( buf, 80, 1, fp ) ; 323 buf[80] = '\0' ; 324 if ( strstr( buf, "NRAO_GBT" ) != NULL ) { 325 isGBTFITS = true ; 326 } 327 fclose( fp ) ; 328 } 239 329 while ( status == 0 ) { 240 330 status = reader_->read(pksrec); … … 288 378 //*srcnCol = pksrec.srcName;//.before(rx2); 289 379 *srctCol = match; 380 if ( pksrec.srcType != -1 ) { 381 *srctCol = pksrec.srcType ; 382 } 290 383 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO"); 291 384 *beamCol = pksrec.beamNo-beamOffset_-1; … … 298 391 uInt id; 299 392 /// @todo this has to change when nchan isn't global anymore 300 id = table_->frequencies().addEntry(Double(header_->nchan/2), 301 pksrec.refFreq, pksrec.freqInc); 393 //id = table_->frequencies().addEntry(Double(header_->nchan/2), 394 // pksrec.refFreq, pksrec.freqInc); 395 if ( pksrec.nchan == 1 ) { 396 id = table_->frequencies().addEntry(Double(0), 397 pksrec.refFreq, pksrec.freqInc); 398 } 399 else { 400 id = table_->frequencies().addEntry(Double(pksrec.nchan/2), 401 pksrec.refFreq, pksrec.freqInc); 402 } 302 403 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID"); 303 404 *mfreqidCol = id; 405 //*ifCol = id; 304 406 305 407 id = table_->molecules().addEntry(pksrec.restFreq); … … 335 437 // into 2-4 rows in the scantable 336 438 Vector<Float> tsysvec(1); 337 // Why is pksrec.spectra.ncolumn() == 3 for haveXPol_ == True439 // Why is spectra.ncolumn() == 3 for haveXPol_ == True 338 440 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2); 339 441 for ( uInt i=0; i< npol; ++i ) { 340 442 tsysvec = pksrec.tsys(i); 341 443 *tsysCol = tsysvec; 342 *polnoCol = i; 444 if (isGBTFITS) 445 *polnoCol = pksrec.polNo ; 446 else 447 *polnoCol = i; 343 448 344 449 *specCol = pksrec.spectra.column(i); … … 347 452 row.put(table_->table().nrow()-1, rec); 348 453 } 454 455 RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW"); 456 *flagrowCol = pksrec.flagrow; 457 349 458 if ( haveXPol_[0] ) { 350 459 // no tsys given for xpol, so emulate it … … 381 490 } 382 491 492 /** 493 * For NRO data 494 * 495 * 2008/11/11 Takeshi Nakazato 496 **/ 497 void STFiller::openNRO( int whichIF, int whichBeam ) 498 { 499 // open file 500 // DEBUG 501 time_t t0 ; 502 time( &t0 ) ; 503 tm *ttm = localtime( &t0 ) ; 504 LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ; 505 // cout << "STFiller::openNRO() Start time = " << t0 506 // << " (" 507 // << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 508 // << " " 509 // << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 510 // << ")" << endl ; 511 os << "Start time = " << t0 512 << " (" 513 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 514 << " " 515 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 516 << ")" << LogIO::POST ; 517 518 // fill STHeader 519 header_ = new STHeader() ; 520 521 if ( nreader_->getHeaderInfo( header_->nchan, 522 header_->npol, 523 nIF_, 524 nBeam_, 525 header_->observer, 526 header_->project, 527 header_->obstype, 528 header_->antennaname, 529 header_->antennaposition, 530 header_->equinox, 531 header_->freqref, 532 header_->reffreq, 533 header_->bandwidth, 534 header_->utc, 535 header_->fluxunit, 536 header_->epoch, 537 header_->poltype ) ) { 538 // cout << "STFiller::openNRO() Failed to get header information." << endl ; 539 // return ; 540 throw( AipsError("Failed to get header information.") ) ; 541 } 542 543 // set FRAME and BASEFRAME keyword of FREQUENCIES table 544 if ( header_->freqref != "TOPO" ) { 545 table_->frequencies().setFrame( header_->freqref, false ) ; 546 table_->frequencies().setFrame( header_->freqref, true ) ; 547 } 548 549 ifOffset_ = 0; 550 vector<Bool> ifs = nreader_->getIFs() ; 551 if ( whichIF >= 0 ) { 552 if ( whichIF >= 0 && whichIF < nIF_ ) { 553 for ( int i = 0 ; i < nIF_ ; i++ ) 554 ifs[i] = False ; 555 ifs[whichIF] = True ; 556 header_->nif = 1; 557 nIF_ = 1; 558 ifOffset_ = whichIF; 559 } else { 560 delete reader_; 561 reader_ = 0; 562 delete header_; 563 header_ = 0; 564 throw(AipsError("Illegal IF selection")); 565 } 566 } 567 568 beamOffset_ = 0; 569 vector<Bool> beams = nreader_->getBeams() ; 570 if (whichBeam>=0) { 571 if (whichBeam>=0 && whichBeam<nBeam_) { 572 for ( int i = 0 ; i < nBeam_ ; i++ ) 573 beams[i] = False ; 574 beams[whichBeam] = True; 575 header_->nbeam = 1; 576 nBeam_ = 1; 577 beamOffset_ = whichBeam; 578 } else { 579 delete reader_; 580 reader_ = 0; 581 delete header_; 582 header_ = 0; 583 throw(AipsError("Illegal Beam selection")); 584 } 585 } 586 587 header_->nbeam = nBeam_ ; 588 header_->nif = nIF_ ; 589 590 // set header 591 table_->setHeader( *header_ ) ; 592 593 // DEBUG 594 time_t t1 ; 595 time( &t1 ) ; 596 ttm = localtime( &t1 ) ; 597 // cout << "STFiller::openNRO() End time = " << t1 598 // << " (" 599 // << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 600 // << " " 601 // << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 602 // << ")" << endl ; 603 // cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ; 604 os << "End time = " << t1 605 << " (" 606 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 607 << " " 608 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 609 << ")" << endl ; 610 os << "Elapsed time = " << t1 - t0 << " sec" << endl ; 611 os.post() ; 612 // 613 614 return ; 615 } 616 617 int STFiller::readNRO() 618 { 619 // DEBUG 620 time_t t0 ; 621 time( &t0 ) ; 622 tm *ttm = localtime( &t0 ) ; 623 LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ; 624 // cout << "STFiller::readNRO() Start time = " << t0 625 // << " (" 626 // << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 627 // << " " 628 // << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 629 // << ")" << endl ; 630 os << "Start time = " << t0 631 << " (" 632 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 633 << " " 634 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 635 << ")" << LogIO::POST ; 636 // 637 638 // fill row 639 uInt id ; 640 uInt imax = nreader_->getRowNum() ; 641 vector< vector<double > > freqs ; 642 uInt i = 0 ; 643 int count = 0 ; 644 uInt scanno ; 645 uInt cycleno ; 646 uInt beamno ; 647 uInt polno ; 648 vector<double> fqs ; 649 Vector<Double> restfreq ; 650 uInt refbeamno ; 651 Double scantime ; 652 Double interval ; 653 String srcname ; 654 String fieldname ; 655 Array<Float> spectra ; 656 Array<uChar> flagtra ; 657 Array<Float> tsys ; 658 Array<Double> direction ; 659 Float azimuth ; 660 Float elevation ; 661 Float parangle ; 662 Float opacity ; 663 uInt tcalid ; 664 Int fitid ; 665 uInt focusid ; 666 Float temperature ; 667 Float pressure ; 668 Float humidity ; 669 Float windvel ; 670 Float winddir ; 671 Double srcvel ; 672 Array<Double> propermotion ; 673 Vector<Double> srcdir ; 674 Array<Double> scanrate ; 675 for ( i = 0 ; i < imax ; i++ ) { 676 string scanType = nreader_->getScanType( i ) ; 677 Int srcType = -1 ; 678 if ( scanType.compare( 0, 2, "ON") == 0 ) { 679 // os << "ON srcType: " << i << LogIO::POST ; 680 srcType = 0 ; 681 } 682 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) { 683 //os << "OFF srcType: " << i << LogIO::POST ; 684 srcType = 1 ; 685 } 686 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) { 687 //os << "ZERO srcType: " << i << LogIO::POST ; 688 srcType = 2 ; 689 } 690 else { 691 //os << "Undefined srcType: " << i << LogIO::POST ; 692 srcType = 3 ; 693 } 694 695 // if srcType is 2 (ZERO scan), ignore scan 696 if ( srcType != 2 && srcType != -1 && srcType != 3 ) { 697 TableRow row( table_->table() ) ; 698 TableRecord& rec = row.record(); 699 700 if ( nreader_->getScanInfo( i, 701 scanno, 702 cycleno, 703 beamno, 704 polno, 705 fqs, 706 restfreq, 707 refbeamno, 708 scantime, 709 interval, 710 srcname, 711 fieldname, 712 spectra, 713 flagtra, 714 tsys, 715 direction, 716 azimuth, 717 elevation, 718 parangle, 719 opacity, 720 tcalid, 721 fitid, 722 focusid, 723 temperature, 724 pressure, 725 humidity, 726 windvel, 727 winddir, 728 srcvel, 729 propermotion, 730 srcdir, 731 scanrate ) ) { 732 // cerr << "STFiller::readNRO() Failed to get scan information." << endl ; 733 // return 1 ; 734 throw( AipsError("Failed to get scan information.") ) ; 735 } 736 737 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ; 738 *scannoCol = scanno ; 739 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ; 740 *cyclenoCol = cycleno ; 741 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ; 742 *beamCol = beamno ; 743 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ; 744 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ; 745 *polnoCol = polno ; 746 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ; 747 if ( freqs.size() == 0 ) { 748 id = table_->frequencies().addEntry( Double( fqs[0] ), 749 Double( fqs[1] ), 750 Double( fqs[2] ) ) ; 751 *mfreqidCol = id ; 752 *ifCol = id ; 753 freqs.push_back( fqs ) ; 754 } 755 else { 756 int iadd = -1 ; 757 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) { 758 //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ; 759 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ; 760 //os << "fdiff = " << fdiff << LogIO::POST ; 761 if ( fdiff < 1.0e-8 ) { 762 iadd = iif ; 763 break ; 764 } 765 } 766 if ( iadd == -1 ) { 767 id = table_->frequencies().addEntry( Double( fqs[0] ), 768 Double( fqs[1] ), 769 Double( fqs[2] ) ) ; 770 *mfreqidCol = id ; 771 *ifCol = id ; 772 freqs.push_back( fqs ) ; 773 } 774 else { 775 *mfreqidCol = iadd ; 776 *ifCol = iadd ; 777 } 778 } 779 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ; 780 id = table_->molecules().addEntry( restfreq ) ; 781 *molidCol = id ; 782 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ; 783 *rbCol = refbeamno ; 784 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ; 785 *mjdCol = scantime ; 786 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ; 787 *intervalCol = interval ; 788 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ; 789 *srcnCol = srcname ; 790 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ; 791 *srctCol = srcType ; 792 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME"); 793 *fieldnCol = fieldname ; 794 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ; 795 *specCol = spectra ; 796 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ; 797 *flagCol = flagtra ; 798 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ; 799 *tsysCol = tsys ; 800 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ; 801 *dirCol = direction ; 802 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ; 803 *azCol = azimuth ; 804 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ; 805 *elCol = elevation ; 806 RecordFieldPtr<Float> parCol(rec, "PARANGLE") ; 807 *parCol = parangle ; 808 RecordFieldPtr<Float> tauCol(rec, "OPACITY") ; 809 *tauCol = opacity ; 810 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ; 811 *mcalidCol = tcalid ; 812 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ; 813 *fitCol = fitid ; 814 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ; 815 *mfocusidCol = focusid ; 816 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ; 817 id = table_->weather().addEntry( temperature, 818 pressure, 819 humidity, 820 windvel, 821 winddir ) ; 822 *mweatheridCol = id ; 823 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ; 824 *svelCol = srcvel ; 825 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ; 826 *spmCol = propermotion ; 827 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ; 828 *sdirCol = srcdir ; 829 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE"); 830 *srateCol = scanrate ; 831 832 table_->table().addRow() ; 833 row.put(table_->table().nrow()-1, rec) ; 834 } 835 else { 836 count++ ; 837 } 838 // DEBUG 839 //int rownum = nreader_->getRowNum() ; 840 //os << "Finished row " << i << "/" << rownum << LogIO::POST ; 841 // 842 } 843 844 // DEBUG 845 time_t t1 ; 846 time( &t1 ) ; 847 ttm = localtime( &t1 ) ; 848 // cout << "STFiller::readNRO() Processed " << i << " rows" << endl ; 849 // cout << "STFiller::readNRO() Added " << i - count << " rows (ignored " 850 // << count << " \"ZERO\" scans)" << endl ; 851 // cout << "STFiller::readNRO() End time = " << t1 852 // << " (" 853 // << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 854 // << " " 855 // << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 856 // << ")" << endl ; 857 // cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ; 858 os << "Processed " << i << " rows" << endl ; 859 os << "Added " << i - count << " rows (ignored " 860 << count << " \"ZERO\" scans)" << endl ; 861 os.post() ; 862 os << "End time = " << t1 863 << " (" 864 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday 865 << " " 866 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec 867 << ")" << endl ; 868 os << "Elapsed time = " << t1 - t0 << " sec" << endl ; 869 os.post() ; 870 // 871 872 return 0 ; 873 } 874 875 Bool STFiller::fileCheck() 876 { 877 bool bval = false ; 878 879 // if filename_ is directory, return false 880 File inFile( filename_ ) ; 881 if ( inFile.isDirectory() ) 882 return bval ; 883 884 // if beginning of header data is "RW", return true 885 // otherwise, return false ; 886 FILE *fp = fopen( filename_.c_str(), "r" ) ; 887 char buf[9] ; 888 char buf2[80] ; 889 fread( buf, 4, 1, fp ) ; 890 buf[4] = '\0' ; 891 fseek( fp, 640, SEEK_SET ) ; 892 fread( buf2, 80, 1, fp ) ; 893 if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) { 894 bval = true ; 895 } 896 fclose( fp ) ; 897 return bval ; 898 } 899 383 900 }//namespace asap -
branches/mergetest/src/STFiller.h
r1504 r1789 27 27 28 28 class PKSreader; 29 class NROReader; 29 30 30 31 namespace asap { … … 61 62 */ 62 63 explicit STFiller( const std::string& filename, int whichIF=-1, 63 int whichBeam=-1 );64 int whichBeam=-1 ); 64 65 65 66 /** … … 75 76 * @exception AipsError Creation of PKSreader failed 76 77 */ 77 void open( const std::string& filename, int whichIF=-1, int whichBeam=-1);78 void open( const std::string& filename, const std::string& antenna, int whichIF=-1, int whichBeam=-1, casa::Bool getPt=casa::False ); 78 79 79 80 /** … … 93 94 casa::CountedPtr<Scantable> getTable() const { return table_;} 94 95 96 /** 97 * For NRO data 98 * 99 * 2008/11/11 Takeshi Nakazato 100 * 101 * openNRO : NRO version of open(), which performs to open file and 102 * read header data. 103 * 104 * readNRO : NRO version of read(), which performs to read scan 105 * records. 106 * 107 * fileCheck: Identify a type (NRO data or not) of filename_. 108 **/ 109 void openNRO( int whichIF=-1, int whichBeam=-1 ) ; 110 int readNRO() ; 111 casa::Bool fileCheck() ; 112 95 113 void setReferenceExpr(const std::string& rx) { refRx_ = rx; } 96 114 … … 105 123 casa::Vector<casa::Bool> haveXPol_; 106 124 casa::String refRx_; 125 NROReader *nreader_ ; 126 casa::Bool isNRO_ ; 107 127 }; 108 128 -
branches/mergetest/src/STFitter.cpp
r1391 r1789 32 32 #include <casa/Arrays/ArrayMath.h> 33 33 #include <casa/Arrays/ArrayLogical.h> 34 #include <casa/Logging/LogIO.h> 34 35 #include <scimath/Fitting.h> 35 36 #include <scimath/Fitting/LinearFit.h> … … 37 38 #include <scimath/Functionals/CompoundFunction.h> 38 39 #include <scimath/Functionals/Gaussian1D.h> 40 #include "Lorentzian1D.h" 39 41 #include <scimath/Functionals/Polynomial.h> 40 42 #include <scimath/Mathematics/AutoDiff.h> … … 146 148 funcs_.resize(1); 147 149 funcs_[0] = new Polynomial<Float>(ncomp); 150 } else if (expr == "lorentz") { 151 if (ncomp < 1) throw (AipsError("Need at least one lorentzian to fit.")); 152 funcs_.resize(ncomp); 153 for (Int k=0; k<ncomp; ++k) { 154 funcs_[k] = new Lorentzian1D<Float>(); 155 } 148 156 } else { 149 cerr << " compiled functions not yet implemented" << endl; 157 //cerr << " compiled functions not yet implemented" << endl; 158 LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ; 159 os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST; 150 160 //funcs_.resize(1); 151 161 //funcs_[0] = new CompiledFunction<Float>(); … … 227 237 (funcs_[0]->parameters())[i] = tmppar[i]; 228 238 } 239 } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) { 240 uInt count = 0; 241 for (uInt j=0; j < funcs_.nelements(); ++j) { 242 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) { 243 (funcs_[j]->parameters())[i] = tmppar[count]; 244 parameters_[count] = tmppar[count]; 245 ++count; 246 } 247 } 229 248 } 230 249 // reset … … 260 279 funcs_[0]->mask(i) = !fixed[i]; 261 280 } 281 } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) { 282 uInt count = 0; 283 for (uInt j=0; j < funcs_.nelements(); ++j) { 284 for (uInt i=0; i < funcs_[j]->nparameters(); ++i) { 285 funcs_[j]->mask(i) = !fixed[count]; 286 fixedpar_[count] = fixed[count]; 287 ++count; 288 } 289 } 262 290 } 263 291 return true; -
branches/mergetest/src/STFrequencies.cpp
r1694 r1789 128 128 } 129 129 130 void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id ) 131 { 132 Table t = table_(table_.col("ID") == Int(id) ); 133 if (t.nrow() == 0 ) { 134 throw(AipsError("STFrequencies::getEntry - freqID out of range")); 135 } 136 for ( uInt i = 0 ; i < table_.nrow() ; i++ ) { 137 uInt fid ; 138 idCol_.get( i, fid ) ; 139 if ( fid == id ) { 140 refpixCol_.put( i, refpix ) ; 141 refvalCol_.put( i, refval ) ; 142 incrCol_.put( i, inc ) ; 143 } 144 } 145 } 146 130 147 SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const 131 148 { … … 145 162 } 146 163 164 /** 147 165 SpectralCoordinate 148 166 STFrequencies::getSpectralCoordinate( const MDirection& md, … … 150 168 const MEpoch& me, 151 169 Double restfreq, uInt id ) const 170 **/ 171 SpectralCoordinate 172 STFrequencies::getSpectralCoordinate( const MDirection& md, 173 const MPosition& mp, 174 const MEpoch& me, 175 Vector<Double> restfreq, uInt id ) const 152 176 { 153 177 SpectralCoordinate spc = getSpectralCoordinate(id); 154 spc.setRestFrequency(restfreq, True); 178 //spc.setRestFrequency(restfreq, True); 179 // for now just use the first rest frequency 180 if (restfreq.nelements()==0 ) { 181 restfreq.resize(1); 182 restfreq[0] = 0; 183 } 184 spc.setRestFrequency(restfreq[0], True); 155 185 if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) { 156 186 throw(AipsError("Couldn't convert frequency frame.")); -
branches/mergetest/src/STFrequencies.h
r1375 r1789 59 59 casa::Double& inc, casa::uInt id ); 60 60 61 /*** 62 * Set the frequency values for a specific id via references 63 * @param refpix the reference pixel 64 * @param refval the reference value 65 * @param inc the increment 66 * @param id the identifier 67 * 68 * 17/09/2008 Takeshi Nakazato 69 ***/ 70 void setEntry( casa::Double refpix, casa::Double refval, 71 casa::Double inc, casa::uInt id ) ; 72 61 73 62 74 bool conformant(const STFrequencies& other) const; … … 69 81 casa::SpectralCoordinate getSpectralCoordinate( casa::uInt freqID ) const; 70 82 83 /** 71 84 casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md, 72 85 const casa::MPosition& mp, … … 75 88 casa::uInt freqID 76 89 ) const; 90 **/ 91 casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md, 92 const casa::MPosition& mp, 93 const casa::MEpoch& me, 94 casa::Vector<casa::Double> restfreq, 95 casa::uInt freqID 96 ) const; 77 97 78 98 /** -
branches/mergetest/src/STHeader.cpp
r1439 r1789 37 37 #include <casa/Arrays/IPosition.h> 38 38 #include <casa/Quanta/MVTime.h> 39 #include <casa/Logging/LogIO.h> 39 40 40 41 #include <sstream> 41 42 42 43 #include "STDefs.h" … … 79 80 MVTime mvt(this->utc); 80 81 mvt.setFormat(MVTime::YMD); 81 cout << "Observer: " << this->observer << endl 82 << "Project: " << this->project << endl 83 << "Obstype: " << this->obstype << endl 84 << "Antenna: " << this->antennaname << endl 85 << "Ant. Position: " << this->antennaposition << endl 86 << "Equinox: " << this->equinox << endl 87 << "Freq. ref.: " << this->freqref << endl 88 << "Ref. frequency: " << this->reffreq << endl 89 << "Bandwidth: " << this->bandwidth << endl 90 << "Time (utc): " 91 << mvt 92 << endl; 82 // cout << "Observer: " << this->observer << endl 83 // << "Project: " << this->project << endl 84 // << "Obstype: " << this->obstype << endl 85 // << "Antenna: " << this->antennaname << endl 86 // << "Ant. Position: " << this->antennaposition << endl 87 // << "Equinox: " << this->equinox << endl 88 // << "Freq. ref.: " << this->freqref << endl 89 // << "Ref. frequency: " << this->reffreq << endl 90 // << "Bandwidth: " << this->bandwidth << endl 91 // << "Time (utc): " 92 // << mvt 93 // << endl; 94 LogIO os( LogOrigin( "STHeader", "print()", WHERE ) ) ; 95 os << "Observer: " << this->observer << endl 96 << "Project: " << this->project << endl 97 << "Obstype: " << this->obstype << endl 98 << "Antenna: " << this->antennaname << endl 99 << "Ant. Position: " << this->antennaposition << endl 100 << "Equinox: " << this->equinox << endl 101 << "Freq. ref.: " << this->freqref << endl 102 << "Ref. frequency: " << this->reffreq << endl 103 << "Bandwidth: " << this->bandwidth << endl 104 << "Time (utc): " 105 << mvt 106 << LogIO::POST ; 93 107 //setprecision(10) << this->utc << endl; 94 108 } … … 129 143 { 130 144 if (n_>0) { 131 cerr << "Source ID" << endl; 132 for (uInt i=0; i<n_; i++) { 133 cout << setw(11) << source_(i) << ID_(i) << endl; 134 } 145 // cerr << "Source ID" << endl; 146 // for (uInt i=0; i<n_; i++) { 147 // cout << setw(11) << source_(i) << ID_(i) << endl; 148 LogIO os( LogOrigin( "SDDataDesc", "summary()", WHERE ) ) ; 149 ostringstream oss ; 150 oss << "Source ID" << endl; 151 for (uInt i=0; i<n_; i++) { 152 oss << setw(11) << source_(i) << ID_(i) << endl; 153 } 154 os << oss.str() << LogIO::POST ; 135 155 } 136 156 } -
branches/mergetest/src/STMath.cpp
r1689 r1789 44 44 #include <scimath/Functionals/Polynomial.h> 45 45 46 #include <atnf/PKSIO/SrcType.h> 47 48 #include <casa/Logging/LogIO.h> 49 #include <sstream> 50 46 51 #include "MathUtils.h" 47 52 #include "RowAccumulator.h" … … 53 58 54 59 using namespace asap; 60 61 // tolerance for direction comparison (rad) 62 #define TOL_OTF 1.0e-15 63 #define TOL_POINT 2.9088821e-4 // 1 arcmin 55 64 56 65 STMath::STMath(bool insitu) : … … 70 79 const std::string& avmode) 71 80 { 81 LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ; 72 82 if ( avmode == "SCAN" && in.size() != 1 ) 73 83 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n" 74 84 "Use merge first.")); 75 85 WeightType wtype = stringToWeight(weight); 86 87 // check if OTF observation 88 String obstype = in[0]->getHeader().obstype ; 89 Double tol = 0.0 ; 90 if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) { 91 tol = TOL_OTF ; 92 } 93 else { 94 tol = TOL_POINT ; 95 } 76 96 77 97 // output … … 113 133 } 114 134 if ( avmode == "SCAN" && in.size() == 1) { 115 cols.resize(4); 116 cols[3] = String("SCANNO"); 135 //cols.resize(4); 136 //cols[3] = String("SCANNO"); 137 cols.resize(5); 138 cols[3] = String("SRCNAME"); 139 cols[4] = String("SCANNO"); 117 140 } 118 141 uInt outrowCount = 0; 119 142 TableIterator iter(baset, cols); 143 // int count = 0 ; 120 144 while (!iter.pastEnd()) { 121 145 Table subt = iter.table(); 122 // copy the first row of this selection into the new table 123 tout.addRow(); 124 TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 125 // re-index to 0 126 if ( avmode != "SCAN" && avmode != "SOURCE" ) { 127 scanColOut.put(outrowCount, uInt(0)); 128 } 129 ++outrowCount; 146 // // copy the first row of this selection into the new table 147 // tout.addRow(); 148 // TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 149 // // re-index to 0 150 // if ( avmode != "SCAN" && avmode != "SOURCE" ) { 151 // scanColOut.put(outrowCount, uInt(0)); 152 // } 153 // ++outrowCount; 154 MDirection::ScalarColumn dircol ; 155 dircol.attach( subt, "DIRECTION" ) ; 156 Int length = subt.nrow() ; 157 vector< Vector<Double> > dirs ; 158 vector<int> indexes ; 159 for ( Int i = 0 ; i < length ; i++ ) { 160 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ; 161 //os << << count++ << ": " ; 162 //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ; 163 bool adddir = true ; 164 for ( uInt j = 0 ; j < dirs.size() ; j++ ) { 165 //if ( allTrue( t == dirs[j] ) ) { 166 Double dx = t[0] - dirs[j][0] ; 167 Double dy = t[1] - dirs[j][1] ; 168 Double dd = sqrt( dx * dx + dy * dy ) ; 169 //if ( allNearAbs( t, dirs[j], tol ) ) { 170 if ( dd <= tol ) { 171 adddir = false ; 172 break ; 173 } 174 } 175 if ( adddir ) { 176 dirs.push_back( t ) ; 177 indexes.push_back( i ) ; 178 } 179 } 180 uInt rowNum = dirs.size() ; 181 tout.addRow( rowNum ) ; 182 for ( uInt i = 0 ; i < rowNum ; i++ ) { 183 TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ; 184 // re-index to 0 185 if ( avmode != "SCAN" && avmode != "SOURCE" ) { 186 scanColOut.put(outrowCount+i, uInt(0)); 187 } 188 } 189 outrowCount += rowNum ; 130 190 ++iter; 131 191 } 132 133 192 RowAccumulator acc(wtype); 134 193 Vector<Bool> cmask(mask); … … 155 214 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") ); 156 215 } else if (avmode == "SCAN") { 157 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) ); 216 //subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) ); 217 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) 218 && basesubt.col("SRCNAME") == rec.asString("SRCNAME") ); 158 219 } else { 159 220 subt = basesubt; 160 221 } 222 223 vector<uInt> removeRows ; 224 uInt nrsubt = subt.nrow() ; 225 for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) { 226 //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) { 227 Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ; 228 Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ; 229 double dx = x0[0] - x1[0] ; 230 double dy = x0[0] - x1[0] ; 231 Double dd = sqrt( dx * dx + dy * dy ) ; 232 //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) { 233 if ( dd > tol ) { 234 removeRows.push_back( irow ) ; 235 } 236 } 237 if ( removeRows.size() != 0 ) { 238 subt.removeRow( removeRows ) ; 239 } 240 241 if ( nrsubt == removeRows.size() ) 242 throw(AipsError("Averaging data is empty.")) ; 243 161 244 specCol.attach(subt,"SPECTRA"); 162 245 flagCol.attach(subt,"FLAGTRA"); … … 192 275 } 193 276 //write out 194 Vector<uChar> flg(msk.shape()); 195 convertArray(flg, !msk); 196 flagColOut.put(i, flg); 197 specColOut.put(i, acc.getSpectrum()); 198 tsysColOut.put(i, acc.getTsys()); 199 intColOut.put(i, acc.getInterval()); 200 mjdColOut.put(i, acc.getTime()); 201 // we should only have one cycle now -> reset it to be 0 202 // frequency switched data has different CYCLENO for different IFNO 203 // which requires resetting this value 204 cycColOut.put(i, uInt(0)); 277 if (acc.state()) { 278 Vector<uChar> flg(msk.shape()); 279 convertArray(flg, !msk); 280 flagColOut.put(i, flg); 281 specColOut.put(i, acc.getSpectrum()); 282 tsysColOut.put(i, acc.getTsys()); 283 intColOut.put(i, acc.getInterval()); 284 mjdColOut.put(i, acc.getTime()); 285 // we should only have one cycle now -> reset it to be 0 286 // frequency switched data has different CYCLENO for different IFNO 287 // which requires resetting this value 288 cycColOut.put(i, uInt(0)); 289 } else { 290 ostringstream oss; 291 oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl; 292 pushLog(String(oss)); 293 } 205 294 acc.reset(); 206 295 } 207 296 if (rowstodelete.nelements() > 0) { 297 //cout << rowstodelete << endl; 298 os << rowstodelete << LogIO::POST ; 208 299 tout.removeRow(rowstodelete); 209 300 if (tout.nrow() == 0) { … … 219 310 const std::string& avmode ) 220 311 { 312 // check if OTF observation 313 String obstype = in->getHeader().obstype ; 314 Double tol = 0.0 ; 315 if ( obstype.find( "OTF" ) != String::npos ) { 316 tol = TOL_OTF ; 317 } 318 else { 319 tol = TOL_POINT ; 320 } 321 221 322 // clone as this is non insitu 222 323 bool insitu = insitu_; … … 250 351 flagCol.attach(subt,"FLAGTRA"); 251 352 tsysCol.attach(subt,"TSYS"); 252 tout.addRow(); 253 TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 254 if ( avmode != "SCAN") { 255 scanColOut.put(outrowCount, uInt(0)); 256 } 257 Vector<Float> tmp; 258 specCol.get(0, tmp); 259 uInt nchan = tmp.nelements(); 260 // have to do channel by channel here as MaskedArrMath 261 // doesn't have partialMedians 262 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 263 Vector<Float> outspec(nchan); 264 Vector<uChar> outflag(nchan,0); 265 Vector<Float> outtsys(1);/// @fixme when tsys is channel based 266 for (uInt i=0; i<nchan; ++i) { 267 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 268 MaskedArray<Float> ma = maskedArray(specs,flags); 269 outspec[i] = median(ma); 270 if ( allEQ(ma.getMask(), False) ) 271 outflag[i] = userflag;// flag data 272 } 273 outtsys[0] = median(tsysCol.getColumn()); 274 specColOut.put(outrowCount, outspec); 275 flagColOut.put(outrowCount, outflag); 276 tsysColOut.put(outrowCount, outtsys); 277 Double intsum = sum(intCol.getColumn()); 278 intColOut.put(outrowCount, intsum); 279 ++outrowCount; 353 // tout.addRow(); 354 // TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 355 // if ( avmode != "SCAN") { 356 // scanColOut.put(outrowCount, uInt(0)); 357 // } 358 // Vector<Float> tmp; 359 // specCol.get(0, tmp); 360 // uInt nchan = tmp.nelements(); 361 // // have to do channel by channel here as MaskedArrMath 362 // // doesn't have partialMedians 363 // Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 364 // Vector<Float> outspec(nchan); 365 // Vector<uChar> outflag(nchan,0); 366 // Vector<Float> outtsys(1);/// @fixme when tsys is channel based 367 // for (uInt i=0; i<nchan; ++i) { 368 // Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 369 // MaskedArray<Float> ma = maskedArray(specs,flags); 370 // outspec[i] = median(ma); 371 // if ( allEQ(ma.getMask(), False) ) 372 // outflag[i] = userflag;// flag data 373 // } 374 // outtsys[0] = median(tsysCol.getColumn()); 375 // specColOut.put(outrowCount, outspec); 376 // flagColOut.put(outrowCount, outflag); 377 // tsysColOut.put(outrowCount, outtsys); 378 // Double intsum = sum(intCol.getColumn()); 379 // intColOut.put(outrowCount, intsum); 380 // ++outrowCount; 381 // ++iter; 382 MDirection::ScalarColumn dircol ; 383 dircol.attach( subt, "DIRECTION" ) ; 384 Int length = subt.nrow() ; 385 vector< Vector<Double> > dirs ; 386 vector<int> indexes ; 387 for ( Int i = 0 ; i < length ; i++ ) { 388 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ; 389 bool adddir = true ; 390 for ( uInt j = 0 ; j < dirs.size() ; j++ ) { 391 //if ( allTrue( t == dirs[j] ) ) { 392 Double dx = t[0] - dirs[j][0] ; 393 Double dy = t[1] - dirs[j][1] ; 394 Double dd = sqrt( dx * dx + dy * dy ) ; 395 //if ( allNearAbs( t, dirs[j], tol ) ) { 396 if ( dd <= tol ) { 397 adddir = false ; 398 break ; 399 } 400 } 401 if ( adddir ) { 402 dirs.push_back( t ) ; 403 indexes.push_back( i ) ; 404 } 405 } 406 uInt rowNum = dirs.size() ; 407 tout.addRow( rowNum ); 408 for ( uInt i = 0 ; i < rowNum ; i++ ) { 409 TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ; 410 if ( avmode != "SCAN") { 411 //scanColOut.put(outrowCount+i, uInt(0)); 412 } 413 } 414 MDirection::ScalarColumn dircolOut ; 415 dircolOut.attach( tout, "DIRECTION" ) ; 416 for ( uInt irow = 0 ; irow < rowNum ; irow++ ) { 417 Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ; 418 Vector<Float> tmp; 419 specCol.get(0, tmp); 420 uInt nchan = tmp.nelements(); 421 // have to do channel by channel here as MaskedArrMath 422 // doesn't have partialMedians 423 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 424 // mask spectra for different DIRECTION 425 for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) { 426 Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ; 427 //if ( t[0] != direction[0] || t[1] != direction[1] ) { 428 Double dx = t[0] - direction[0] ; 429 Double dy = t[1] - direction[1] ; 430 Double dd = sqrt( dx * dx + dy * dy ) ; 431 //if ( !allNearAbs( t, direction, tol ) ) { 432 if ( dd > tol ) { 433 flags[jrow] = userflag ; 434 } 435 } 436 Vector<Float> outspec(nchan); 437 Vector<uChar> outflag(nchan,0); 438 Vector<Float> outtsys(1);/// @fixme when tsys is channel based 439 for (uInt i=0; i<nchan; ++i) { 440 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 441 MaskedArray<Float> ma = maskedArray(specs,flags); 442 outspec[i] = median(ma); 443 if ( allEQ(ma.getMask(), False) ) 444 outflag[i] = userflag;// flag data 445 } 446 outtsys[0] = median(tsysCol.getColumn()); 447 specColOut.put(outrowCount+irow, outspec); 448 flagColOut.put(outrowCount+irow, outflag); 449 tsysColOut.put(outrowCount+irow, outtsys); 450 Vector<Double> integ = intCol.getColumn() ; 451 MaskedArray<Double> mi = maskedArray( integ, flags ) ; 452 Double intsum = sum(mi); 453 intColOut.put(outrowCount+irow, intsum); 454 } 455 outrowCount += rowNum ; 280 456 ++iter; 281 457 } … … 323 499 if ( tsys ) { 324 500 ts += val; 501 tsysCol.put(i, ts); 502 } 503 } 504 } 505 return out; 506 } 507 508 CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in, 509 const std::vector<float> val, 510 const std::string& mode, 511 const std::string& opmode, 512 bool tsys ) 513 { 514 CountedPtr< Scantable > out ; 515 if ( opmode == "channel" ) { 516 out = arrayOperateChannel( in, val, mode, tsys ) ; 517 } 518 else if ( opmode == "row" ) { 519 out = arrayOperateRow( in, val, mode, tsys ) ; 520 } 521 else { 522 throw( AipsError( "Unknown array operation mode." ) ) ; 523 } 524 return out ; 525 } 526 527 CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in, 528 const std::vector<float> val, 529 const std::string& mode, 530 bool tsys ) 531 { 532 if ( val.size() == 1 ){ 533 return unaryOperate( in, val[0], mode, tsys ) ; 534 } 535 536 // conformity of SPECTRA and TSYS 537 if ( tsys ) { 538 TableIterator titer(in->table(), "IFNO"); 539 while ( !titer.pastEnd() ) { 540 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ; 541 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ; 542 Array<Float> spec = specCol.getColumn() ; 543 Array<Float> ts = tsysCol.getColumn() ; 544 if ( !spec.conform( ts ) ) { 545 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ; 546 } 547 titer.next() ; 548 } 549 } 550 551 // check if all spectra in the scantable have the same number of channel 552 vector<uInt> nchans; 553 vector<uInt> ifnos = in->getIFNos() ; 554 for ( uInt i = 0 ; i < ifnos.size() ; i++ ) { 555 nchans.push_back( in->nchan( ifnos[i] ) ) ; 556 } 557 Vector<uInt> mchans( nchans ) ; 558 if ( anyNE( mchans, mchans[0] ) ) { 559 throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ; 560 } 561 562 // check if vector size is equal to nchan 563 Vector<Float> fact( val ) ; 564 if ( fact.nelements() != mchans[0] ) { 565 throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ; 566 } 567 568 // check divided by zero 569 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) { 570 throw( AipsError("Divided by zero is not recommended." ) ) ; 571 } 572 573 CountedPtr< Scantable > out = getScantable(in, false); 574 Table& tab = out->table(); 575 ArrayColumn<Float> specCol(tab,"SPECTRA"); 576 ArrayColumn<Float> tsysCol(tab,"TSYS"); 577 for (uInt i=0; i<tab.nrow(); ++i) { 578 Vector<Float> spec; 579 Vector<Float> ts; 580 specCol.get(i, spec); 581 tsysCol.get(i, ts); 582 if (mode == "MUL" || mode == "DIV") { 583 if (mode == "DIV") fact = (float)1.0 / fact; 584 spec *= fact; 585 specCol.put(i, spec); 586 if ( tsys ) { 587 ts *= fact; 588 tsysCol.put(i, ts); 589 } 590 } else if ( mode == "ADD" || mode == "SUB") { 591 if (mode == "SUB") fact *= (float)-1.0 ; 592 spec += fact; 593 specCol.put(i, spec); 594 if ( tsys ) { 595 ts += fact; 596 tsysCol.put(i, ts); 597 } 598 } 599 } 600 return out; 601 } 602 603 CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in, 604 const std::vector<float> val, 605 const std::string& mode, 606 bool tsys ) 607 { 608 if ( val.size() == 1 ) { 609 return unaryOperate( in, val[0], mode, tsys ) ; 610 } 611 612 // conformity of SPECTRA and TSYS 613 if ( tsys ) { 614 TableIterator titer(in->table(), "IFNO"); 615 while ( !titer.pastEnd() ) { 616 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ; 617 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ; 618 Array<Float> spec = specCol.getColumn() ; 619 Array<Float> ts = tsysCol.getColumn() ; 620 if ( !spec.conform( ts ) ) { 621 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ; 622 } 623 titer.next() ; 624 } 625 } 626 627 // check if vector size is equal to nrow 628 Vector<Float> fact( val ) ; 629 if ( fact.nelements() != in->nrow() ) { 630 throw( AipsError("Vector size must be 1 or be same as number of row.") ) ; 631 } 632 633 // check divided by zero 634 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) { 635 throw( AipsError("Divided by zero is not recommended." ) ) ; 636 } 637 638 CountedPtr< Scantable > out = getScantable(in, false); 639 Table& tab = out->table(); 640 ArrayColumn<Float> specCol(tab,"SPECTRA"); 641 ArrayColumn<Float> tsysCol(tab,"TSYS"); 642 if (mode == "DIV") fact = (float)1.0 / fact; 643 if (mode == "SUB") fact *= (float)-1.0 ; 644 for (uInt i=0; i<tab.nrow(); ++i) { 645 Vector<Float> spec; 646 Vector<Float> ts; 647 specCol.get(i, spec); 648 tsysCol.get(i, ts); 649 if (mode == "MUL" || mode == "DIV") { 650 spec *= fact[i]; 651 specCol.put(i, spec); 652 if ( tsys ) { 653 ts *= fact[i]; 654 tsysCol.put(i, ts); 655 } 656 } else if ( mode == "ADD" || mode == "SUB") { 657 spec += fact[i]; 658 specCol.put(i, spec); 659 if ( tsys ) { 660 ts += fact[i]; 661 tsysCol.put(i, ts); 662 } 663 } 664 } 665 return out; 666 } 667 668 CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in, 669 const std::vector< std::vector<float> > val, 670 const std::string& mode, 671 bool tsys ) 672 { 673 // conformity of SPECTRA and TSYS 674 if ( tsys ) { 675 TableIterator titer(in->table(), "IFNO"); 676 while ( !titer.pastEnd() ) { 677 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ; 678 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ; 679 Array<Float> spec = specCol.getColumn() ; 680 Array<Float> ts = tsysCol.getColumn() ; 681 if ( !spec.conform( ts ) ) { 682 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ; 683 } 684 titer.next() ; 685 } 686 } 687 688 // some checks 689 vector<uInt> nchans; 690 for ( uInt i = 0 ; i < in->nrow() ; i++ ) { 691 nchans.push_back( (in->getSpectrum( i )).size() ) ; 692 } 693 //Vector<uInt> mchans( nchans ) ; 694 vector< Vector<Float> > facts ; 695 for ( uInt i = 0 ; i < nchans.size() ; i++ ) { 696 Vector<Float> tmp( val[i] ) ; 697 // check divided by zero 698 if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) { 699 throw( AipsError("Divided by zero is not recommended." ) ) ; 700 } 701 // conformity check 702 if ( tmp.nelements() != nchans[i] ) { 703 stringstream ss ; 704 ss << "Row " << i << ": Vector size must be same as number of channel." ; 705 throw( AipsError( ss.str() ) ) ; 706 } 707 facts.push_back( tmp ) ; 708 } 709 710 711 CountedPtr< Scantable > out = getScantable(in, false); 712 Table& tab = out->table(); 713 ArrayColumn<Float> specCol(tab,"SPECTRA"); 714 ArrayColumn<Float> tsysCol(tab,"TSYS"); 715 for (uInt i=0; i<tab.nrow(); ++i) { 716 Vector<Float> fact = facts[i] ; 717 Vector<Float> spec; 718 Vector<Float> ts; 719 specCol.get(i, spec); 720 tsysCol.get(i, ts); 721 if (mode == "MUL" || mode == "DIV") { 722 if (mode == "DIV") fact = (float)1.0 / fact; 723 spec *= fact; 724 specCol.put(i, spec); 725 if ( tsys ) { 726 ts *= fact; 727 tsysCol.put(i, ts); 728 } 729 } else if ( mode == "ADD" || mode == "SUB") { 730 if (mode == "SUB") fact *= (float)-1.0 ; 731 spec += fact; 732 specCol.put(i, spec); 733 if ( tsys ) { 734 ts += fact; 325 735 tsysCol.put(i, ts); 326 736 } … … 386 796 } 387 797 798 MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s, 799 const Vector<uChar>& f) 800 { 801 Vector<Bool> mask; 802 mask.resize(f.shape()); 803 convertArray(mask, f); 804 return MaskedArray<Double>(s,!mask); 805 } 806 388 807 Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma) 389 808 { … … 402 821 // make this operation non insitu 403 822 const Table& tin = in->table(); 404 Table ons = tin(tin.col("SRCTYPE") == Int( 0));405 Table offs = tin(tin.col("SRCTYPE") == Int( 1));823 Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON)); 824 Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF)); 406 825 if ( offs.nrow() == 0 ) 407 826 throw(AipsError("No 'off' scans present.")); … … 641 1060 //Debug 642 1061 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl; 1062 //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ; 1063 //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST; 643 1064 meanoff = sum(spoff)/noff; 644 1065 meandiff = sum(spdiff)/ndiff; … … 760 1181 //Debug 761 1182 //cerr<<"Tsys used="<<tsysrefscalar<<endl; 1183 //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ; 1184 //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST; 762 1185 // fill the result, replay signal tsys by reference tsys 763 1186 outintCol.put(i, resint); … … 782 1205 setInsitu(false); 783 1206 STSelector sel; 784 std::vector<int> scan1, scan2, beams ;1207 std::vector<int> scan1, scan2, beams, types; 785 1208 std::vector< vector<int> > scanpair; 786 std::vector<string> calstate; 1209 //std::vector<string> calstate; 1210 std::vector<int> calstate; 787 1211 String msg; 788 1212 … … 831 1255 scanpair.push_back(scan1); 832 1256 scanpair.push_back(scan2); 833 calstate.push_back("*calon"); 834 calstate.push_back("*[^calon]"); 1257 //calstate.push_back("*calon"); 1258 //calstate.push_back("*[^calon]"); 1259 calstate.push_back(SrcType::NODCAL); 1260 calstate.push_back(SrcType::NOD); 835 1261 CountedPtr< Scantable > ws = getScantable(s, false); 836 1262 uInt l=0; … … 841 1267 sel.reset(); 842 1268 sel.setScans(scanpair[i]); 843 sel.setName(calstate[k]); 1269 //sel.setName(calstate[k]); 1270 types.clear(); 1271 types.push_back(calstate[k]); 1272 sel.setTypes(types); 844 1273 beams.clear(); 845 1274 beams.push_back(j); … … 926 1355 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2); 927 1356 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl; 1357 // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ; 1358 // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST; 928 1359 tsys1[0] = sqrt(tsyssq1 + tsyssq2); 929 1360 Array<Float> avtsys = tsys1; … … 952 1383 CountedPtr< Scantable > ws = getScantable(s, false); 953 1384 CountedPtr< Scantable > sig, sigwcal, ref, refwcal; 954 CountedPtr< Scantable > calsig, calref, out; 1385 CountedPtr< Scantable > calsig, calref, out, out1, out2; 1386 Bool nofold=False; 1387 vector<int> types ; 955 1388 956 1389 //split the data 957 sel.setName("*_fs"); 1390 //sel.setName("*_fs"); 1391 types.push_back( SrcType::FSON ) ; 1392 sel.setTypes( types ) ; 958 1393 ws->setSelection(sel); 959 1394 sig = getScantable(ws,false); 960 1395 sel.reset(); 961 sel.setName("*_fs_calon"); 1396 types.clear() ; 1397 //sel.setName("*_fs_calon"); 1398 types.push_back( SrcType::FONCAL ) ; 1399 sel.setTypes( types ) ; 962 1400 ws->setSelection(sel); 963 1401 sigwcal = getScantable(ws,false); 964 1402 sel.reset(); 965 sel.setName("*_fsr"); 1403 types.clear() ; 1404 //sel.setName("*_fsr"); 1405 types.push_back( SrcType::FSOFF ) ; 1406 sel.setTypes( types ) ; 966 1407 ws->setSelection(sel); 967 1408 ref = getScantable(ws,false); 968 1409 sel.reset(); 969 sel.setName("*_fsr_calon"); 1410 types.clear() ; 1411 //sel.setName("*_fsr_calon"); 1412 types.push_back( SrcType::FOFFCAL ) ; 1413 sel.setTypes( types ) ; 970 1414 ws->setSelection(sel); 971 1415 refwcal = getScantable(ws,false); 1416 sel.reset() ; 1417 types.clear() ; 972 1418 973 1419 calsig = dototalpower(sigwcal, sig, tcal=tcal); 974 1420 calref = dototalpower(refwcal, ref, tcal=tcal); 975 1421 976 out=dosigref(calsig,calref,smoothref,tsysv,tau); 977 1422 out1=dosigref(calsig,calref,smoothref,tsysv,tau); 1423 out2=dosigref(calref,calsig,smoothref,tsysv,tau); 1424 1425 Table& tabout1=out1->table(); 1426 Table& tabout2=out2->table(); 1427 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID"); 1428 ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID"); 1429 ROArrayColumn<Float> specCol(tabout2, "SPECTRA"); 1430 Vector<Float> spec; specCol.get(0, spec); 1431 uInt nchan = spec.nelements(); 1432 uInt freqid1; freqidCol1.get(0,freqid1); 1433 uInt freqid2; freqidCol2.get(0,freqid2); 1434 Double rp1, rp2, rv1, rv2, inc1, inc2; 1435 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1); 1436 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2); 1437 //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ; 1438 //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ; 1439 //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ; 1440 if (rp1==rp2) { 1441 Double foffset = rv1 - rv2; 1442 uInt choffset = static_cast<uInt>(foffset/abs(inc2)); 1443 if (choffset >= nchan) { 1444 //cerr<<"out-band frequency switching, no folding"<<endl; 1445 LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ; 1446 os<<"out-band frequency switching, no folding"<<LogIO::POST; 1447 nofold = True; 1448 } 1449 } 1450 1451 if (nofold) { 1452 std::vector< CountedPtr< Scantable > > tabs; 1453 tabs.push_back(out1); 1454 tabs.push_back(out2); 1455 out = merge(tabs); 1456 } 1457 else { 1458 //out = out1; 1459 Double choffset = ( rv1 - rv2 ) / inc2 ; 1460 out = dofold( out1, out2, choffset ) ; 1461 } 1462 978 1463 return out; 1464 } 1465 1466 CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig, 1467 const CountedPtr<Scantable> &ref, 1468 Double choffset, 1469 Double choffset2 ) 1470 { 1471 LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ; 1472 os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ; 1473 1474 // output scantable 1475 CountedPtr<Scantable> out = getScantable( sig, false ) ; 1476 1477 // separate choffset to integer part and decimal part 1478 Int ioffset = (Int)choffset ; 1479 Double doffset = choffset - ioffset ; 1480 Int ioffset2 = (Int)choffset2 ; 1481 Double doffset2 = choffset2 - ioffset2 ; 1482 os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ; 1483 os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ; 1484 1485 // get column 1486 ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ; 1487 ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ; 1488 ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ; 1489 ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ; 1490 ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ; 1491 ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ; 1492 ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ; 1493 ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ; 1494 ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ; 1495 ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ; 1496 1497 // check 1498 if ( ioffset == 0 ) { 1499 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ; 1500 os << "channel offset is zero, no folding" << LogIO::POST ; 1501 return out ; 1502 } 1503 int nchan = ref->nchan() ; 1504 if ( abs(ioffset) >= nchan ) { 1505 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ; 1506 os << "out-band frequency switching, no folding" << LogIO::POST ; 1507 return out ; 1508 } 1509 1510 // attach column for output scantable 1511 ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ; 1512 ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ; 1513 ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ; 1514 ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ; 1515 ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ; 1516 ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ; 1517 1518 // for each row 1519 // assume that the data order are same between sig and ref 1520 RowAccumulator acc( asap::W_TINTSYS ) ; 1521 for ( int i = 0 ; i < sig->nrow() ; i++ ) { 1522 // get values 1523 Vector<Float> spsig ; 1524 specCol1.get( i, spsig ) ; 1525 Vector<Float> spref ; 1526 specCol2.get( i, spref ) ; 1527 Vector<Float> tsyssig ; 1528 tsysCol1.get( i, tsyssig ) ; 1529 Vector<Float> tsysref ; 1530 tsysCol2.get( i, tsysref ) ; 1531 Vector<uChar> flagsig ; 1532 flagCol1.get( i, flagsig ) ; 1533 Vector<uChar> flagref ; 1534 flagCol2.get( i, flagref ) ; 1535 Double timesig ; 1536 mjdCol1.get( i, timesig ) ; 1537 Double timeref ; 1538 mjdCol2.get( i, timeref ) ; 1539 Double intsig ; 1540 intervalCol1.get( i, intsig ) ; 1541 Double intref ; 1542 intervalCol2.get( i, intref ) ; 1543 1544 // shift reference spectra 1545 int refchan = spref.nelements() ; 1546 Vector<Float> sspref( spref.nelements() ) ; 1547 Vector<Float> stsysref( tsysref.nelements() ) ; 1548 Vector<uChar> sflagref( flagref.nelements() ) ; 1549 if ( ioffset > 0 ) { 1550 // SPECTRA and FLAGTRA 1551 for ( int j = 0 ; j < refchan-ioffset ; j++ ) { 1552 sspref[j] = spref[j+ioffset] ; 1553 sflagref[j] = flagref[j+ioffset] ; 1554 } 1555 for ( int j = refchan-ioffset ; j < refchan ; j++ ) { 1556 sspref[j] = spref[j-refchan+ioffset] ; 1557 sflagref[j] = flagref[j-refchan+ioffset] ; 1558 } 1559 spref = sspref.copy() ; 1560 flagref = sflagref.copy() ; 1561 for ( int j = 0 ; j < refchan - 1 ; j++ ) { 1562 sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ; 1563 sflagref[j] = flagref[j+1] + flagref[j] ; 1564 } 1565 sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ; 1566 sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ; 1567 1568 // TSYS 1569 if ( spref.nelements() == tsysref.nelements() ) { 1570 for ( int j = 0 ; j < refchan-ioffset ; j++ ) { 1571 stsysref[j] = tsysref[j+ioffset] ; 1572 } 1573 for ( int j = refchan-ioffset ; j < refchan ; j++ ) { 1574 stsysref[j] = tsysref[j-refchan+ioffset] ; 1575 } 1576 tsysref = stsysref.copy() ; 1577 for ( int j = 0 ; j < refchan - 1 ; j++ ) { 1578 stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ; 1579 } 1580 stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ; 1581 } 1582 } 1583 else { 1584 // SPECTRA and FLAGTRA 1585 for ( int j = 0 ; j < abs(ioffset) ; j++ ) { 1586 sspref[j] = spref[refchan+ioffset+j] ; 1587 sflagref[j] = flagref[refchan+ioffset+j] ; 1588 } 1589 for ( int j = abs(ioffset) ; j < refchan ; j++ ) { 1590 sspref[j] = spref[j+ioffset] ; 1591 sflagref[j] = flagref[j+ioffset] ; 1592 } 1593 spref = sspref.copy() ; 1594 flagref = sflagref.copy() ; 1595 sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ; 1596 sflagref[0] = flagref[0] + flagref[refchan-1] ; 1597 for ( int j = 1 ; j < refchan ; j++ ) { 1598 sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ; 1599 sflagref[j] = flagref[j-1] + flagref[j] ; 1600 } 1601 // TSYS 1602 if ( spref.nelements() == tsysref.nelements() ) { 1603 for ( int j = 0 ; j < abs(ioffset) ; j++ ) { 1604 stsysref[j] = tsysref[refchan+ioffset+j] ; 1605 } 1606 for ( int j = abs(ioffset) ; j < refchan ; j++ ) { 1607 stsysref[j] = tsysref[j+ioffset] ; 1608 } 1609 tsysref = stsysref.copy() ; 1610 stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ; 1611 for ( int j = 1 ; j < refchan ; j++ ) { 1612 stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ; 1613 } 1614 } 1615 } 1616 1617 // shift signal spectra if necessary (only for APEX?) 1618 if ( choffset2 != 0.0 ) { 1619 int sigchan = spsig.nelements() ; 1620 Vector<Float> sspsig( spsig.nelements() ) ; 1621 Vector<Float> stsyssig( tsyssig.nelements() ) ; 1622 Vector<uChar> sflagsig( flagsig.nelements() ) ; 1623 if ( ioffset2 > 0 ) { 1624 // SPECTRA and FLAGTRA 1625 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) { 1626 sspsig[j] = spsig[j+ioffset2] ; 1627 sflagsig[j] = flagsig[j+ioffset2] ; 1628 } 1629 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) { 1630 sspsig[j] = spsig[j-sigchan+ioffset2] ; 1631 sflagsig[j] = flagsig[j-sigchan+ioffset2] ; 1632 } 1633 spsig = sspsig.copy() ; 1634 flagsig = sflagsig.copy() ; 1635 for ( int j = 0 ; j < sigchan - 1 ; j++ ) { 1636 sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ; 1637 sflagsig[j] = flagsig[j+1] || flagsig[j] ; 1638 } 1639 sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ; 1640 sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ; 1641 // TSTS 1642 if ( spsig.nelements() == tsyssig.nelements() ) { 1643 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) { 1644 stsyssig[j] = tsyssig[j+ioffset2] ; 1645 } 1646 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) { 1647 stsyssig[j] = tsyssig[j-sigchan+ioffset2] ; 1648 } 1649 tsyssig = stsyssig.copy() ; 1650 for ( int j = 0 ; j < sigchan - 1 ; j++ ) { 1651 stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ; 1652 } 1653 stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ; 1654 } 1655 } 1656 else { 1657 // SPECTRA and FLAGTRA 1658 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) { 1659 sspsig[j] = spsig[sigchan+ioffset2+j] ; 1660 sflagsig[j] = flagsig[sigchan+ioffset2+j] ; 1661 } 1662 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) { 1663 sspsig[j] = spsig[j+ioffset2] ; 1664 sflagsig[j] = flagsig[j+ioffset2] ; 1665 } 1666 spsig = sspsig.copy() ; 1667 flagsig = sflagsig.copy() ; 1668 sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ; 1669 sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ; 1670 for ( int j = 1 ; j < sigchan ; j++ ) { 1671 sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ; 1672 sflagsig[j] = flagsig[j-1] + flagsig[j] ; 1673 } 1674 // TSYS 1675 if ( spsig.nelements() == tsyssig.nelements() ) { 1676 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) { 1677 stsyssig[j] = tsyssig[sigchan+ioffset2+j] ; 1678 } 1679 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) { 1680 stsyssig[j] = tsyssig[j+ioffset2] ; 1681 } 1682 tsyssig = stsyssig.copy() ; 1683 stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ; 1684 for ( int j = 1 ; j < sigchan ; j++ ) { 1685 stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ; 1686 } 1687 } 1688 } 1689 } 1690 1691 // folding 1692 acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ; 1693 acc.add( sspref, !sflagref, stsysref, intref, timeref ) ; 1694 1695 // put result 1696 specColOut.put( i, acc.getSpectrum() ) ; 1697 const Vector<Bool> &msk = acc.getMask() ; 1698 Vector<uChar> flg( msk.shape() ) ; 1699 convertArray( flg, !msk ) ; 1700 flagColOut.put( i, flg ) ; 1701 tsysColOut.put( i, acc.getTsys() ) ; 1702 intervalColOut.put( i, acc.getInterval() ) ; 1703 mjdColOut.put( i, acc.getTime() ) ; 1704 // change FREQ_ID to unshifted IF setting (only for APEX?) 1705 if ( choffset2 != 0.0 ) { 1706 uInt freqid = fidColOut( 0 ) ; // assume single-IF data 1707 double refpix, refval, increment ; 1708 out->frequencies().getEntry( refpix, refval, increment, freqid ) ; 1709 refval -= choffset * increment ; 1710 uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ; 1711 Vector<uInt> freqids = fidColOut.getColumn() ; 1712 for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) { 1713 if ( freqids[j] == freqid ) 1714 freqids[j] = newfreqid ; 1715 } 1716 fidColOut.putColumn( freqids ) ; 1717 } 1718 1719 acc.reset() ; 1720 } 1721 1722 return out ; 979 1723 } 980 1724 … … 1051 1795 } 1052 1796 out.push_back(outstat); 1797 } 1798 return out; 1799 } 1800 1801 std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in, 1802 const std::vector< bool > & mask, 1803 const std::string& which ) 1804 { 1805 1806 Vector<Bool> m(mask); 1807 const Table& tab = in->table(); 1808 ROArrayColumn<Float> specCol(tab, "SPECTRA"); 1809 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA"); 1810 std::vector<int> out; 1811 for (uInt i=0; i < tab.nrow(); ++i ) { 1812 Vector<Float> spec; specCol.get(i, spec); 1813 Vector<uChar> flag; flagCol.get(i, flag); 1814 MaskedArray<Float> ma = maskedArray(spec, flag); 1815 if (ma.ndim() != 1) { 1816 throw (ArrayError( 1817 "std::vector<int> STMath::minMaxChan(" 1818 "ContedPtr<Scantable> &in, std::vector<bool> &mask, " 1819 " std::string &which)" 1820 " - MaskedArray is not 1D")); 1821 } 1822 IPosition outpos(1,0); 1823 if ( spec.nelements() == m.nelements() ) { 1824 outpos = mathutil::minMaxPos(which, ma(m)); 1825 } else { 1826 outpos = mathutil::minMaxPos(which, ma); 1827 } 1828 out.push_back(outpos[0]); 1053 1829 } 1054 1830 return out; … … 1571 2347 if ( ! (*it)->conformant(*out) ) { 1572 2348 // non conformant. 1573 pushLog(String("Warning: Can't merge scantables as header info differs.")); 2349 //pushLog(String("Warning: Can't merge scantables as header info differs.")); 2350 LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ; 2351 os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ; 1574 2352 } 1575 2353 out->appendToHistoryTable((*it)->history()); … … 1593 2371 id = out->frequencies().addEntry(rp, rv, inc); 1594 2372 freqidcol.put(k,id); 1595 String name,fname;Double rf; 2373 //String name,fname;Double rf; 2374 Vector<String> name,fname;Vector<Double> rf; 1596 2375 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID")); 1597 2376 id = out->molecules().addEntry(rf, name, fname); … … 1997 2776 int fstart = -1; 1998 2777 int fend = -1; 1999 for ( int k=0; k < flag.nelements(); ++k ) {2778 for (unsigned int k=0; k < flag.nelements(); ++k ) { 2000 2779 if (flag[k] > 0) { 2001 2780 fstart = k; … … 2051 2830 return out; 2052 2831 } 2832 2833 // Averaging spectra with different channel/resolution 2834 CountedPtr<Scantable> 2835 STMath::new_average( const std::vector<CountedPtr<Scantable> >& in, 2836 const bool& compel, 2837 const std::vector<bool>& mask, 2838 const std::string& weight, 2839 const std::string& avmode ) 2840 throw ( casa::AipsError ) 2841 { 2842 LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ; 2843 if ( avmode == "SCAN" && in.size() != 1 ) 2844 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n" 2845 "Use merge first.")); 2846 2847 // check if OTF observation 2848 String obstype = in[0]->getHeader().obstype ; 2849 Double tol = 0.0 ; 2850 if ( obstype.find( "OTF" ) != String::npos ) { 2851 tol = TOL_OTF ; 2852 } 2853 else { 2854 tol = TOL_POINT ; 2855 } 2856 2857 CountedPtr<Scantable> out ; // processed result 2858 if ( compel ) { 2859 std::vector< CountedPtr<Scantable> > newin ; // input for average process 2860 uInt insize = in.size() ; // number of input scantables 2861 2862 // TEST: do normal average in each table before IF grouping 2863 os << "Do preliminary averaging" << LogIO::POST ; 2864 vector< CountedPtr<Scantable> > tmpin( insize ) ; 2865 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 2866 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ; 2867 tmpin[itable] = average( v, mask, weight, avmode ) ; 2868 } 2869 2870 // warning 2871 os << "Average spectra with different spectral resolution" << LogIO::POST ; 2872 2873 // temporarily set coordinfo 2874 vector<string> oldinfo( insize ) ; 2875 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 2876 vector<string> coordinfo = in[itable]->getCoordInfo() ; 2877 oldinfo[itable] = coordinfo[0] ; 2878 coordinfo[0] = "Hz" ; 2879 tmpin[itable]->setCoordInfo( coordinfo ) ; 2880 } 2881 2882 // columns 2883 ScalarColumn<uInt> freqIDCol ; 2884 ScalarColumn<uInt> ifnoCol ; 2885 ScalarColumn<uInt> scannoCol ; 2886 2887 2888 // check IF frequency coverage 2889 // freqid: list of FREQ_ID, which is used, in each table 2890 // iffreq: list of minimum and maximum frequency for each FREQ_ID in 2891 // each table 2892 // freqid[insize][numIF] 2893 // freqid: [[id00, id01, ...], 2894 // [id10, id11, ...], 2895 // ... 2896 // [idn0, idn1, ...]] 2897 // iffreq[insize][numIF*2] 2898 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...], 2899 // [min_id10, max_id10, min_id11, max_id11, ...], 2900 // ... 2901 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]] 2902 //os << "Check IF settings in each table" << LogIO::POST ; 2903 vector< vector<uInt> > freqid( insize ); 2904 vector< vector<double> > iffreq( insize ) ; 2905 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 2906 uInt rows = tmpin[itable]->nrow() ; 2907 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ; 2908 for ( uInt irow = 0 ; irow < rows ; irow++ ) { 2909 if ( freqid[itable].size() == freqnrows ) { 2910 break ; 2911 } 2912 else { 2913 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ; 2914 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ; 2915 uInt id = freqIDCol( irow ) ; 2916 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) { 2917 //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ; 2918 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ; 2919 freqid[itable].push_back( id ) ; 2920 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ; 2921 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ; 2922 } 2923 } 2924 } 2925 } 2926 2927 // debug 2928 //os << "IF settings summary:" << endl ; 2929 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) { 2930 //os << " Table" << i << endl ; 2931 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) { 2932 //os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ; 2933 //} 2934 //} 2935 //os << endl ; 2936 //os.post() ; 2937 2938 // IF grouping based on their frequency coverage 2939 // ifgrp: list of table index and FREQ_ID for all members in each IF group 2940 // ifgfreq: list of minimum and maximum frequency in each IF group 2941 // ifgrp[numgrp][nummember*2] 2942 // ifgrp: [[table00, freqrow00, table01, freqrow01, ...], 2943 // [table10, freqrow10, table11, freqrow11, ...], 2944 // ... 2945 // [tablen0, freqrown0, tablen1, freqrown1, ...]] 2946 // ifgfreq[numgrp*2] 2947 // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...] 2948 //os << "IF grouping based on their frequency coverage" << LogIO::POST ; 2949 vector< vector<uInt> > ifgrp ; 2950 vector<double> ifgfreq ; 2951 2952 // parameter for IF grouping 2953 // groupmode = OR retrieve all region 2954 // AND only retrieve overlaped region 2955 //string groupmode = "AND" ; 2956 string groupmode = "OR" ; 2957 uInt sizecr = 0 ; 2958 if ( groupmode == "AND" ) 2959 sizecr = 2 ; 2960 else if ( groupmode == "OR" ) 2961 sizecr = 0 ; 2962 2963 vector<double> sortedfreq ; 2964 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) { 2965 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) { 2966 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 ) 2967 sortedfreq.push_back( iffreq[i][j] ) ; 2968 } 2969 } 2970 sort( sortedfreq.begin(), sortedfreq.end() ) ; 2971 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) { 2972 ifgfreq.push_back( *i ) ; 2973 ifgfreq.push_back( *(i+1) ) ; 2974 } 2975 ifgrp.resize( ifgfreq.size()/2 ) ; 2976 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 2977 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) { 2978 double range0 = iffreq[itable][2*iif] ; 2979 double range1 = iffreq[itable][2*iif+1] ; 2980 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) { 2981 double fmin = max( range0, ifgfreq[2*j] ) ; 2982 double fmax = min( range1, ifgfreq[2*j+1] ) ; 2983 if ( fmin < fmax ) { 2984 ifgrp[j].push_back( itable ) ; 2985 ifgrp[j].push_back( freqid[itable][iif] ) ; 2986 } 2987 } 2988 } 2989 } 2990 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ; 2991 vector<double>::iterator giter = ifgfreq.begin() ; 2992 while( fiter != ifgrp.end() ) { 2993 if ( fiter->size() <= sizecr ) { 2994 fiter = ifgrp.erase( fiter ) ; 2995 giter = ifgfreq.erase( giter ) ; 2996 giter = ifgfreq.erase( giter ) ; 2997 } 2998 else { 2999 fiter++ ; 3000 advance( giter, 2 ) ; 3001 } 3002 } 3003 3004 // Grouping continuous IF groups (without frequency gap) 3005 // freqgrp: list of IF group indexes in each frequency group 3006 // freqrange: list of minimum and maximum frequency in each frequency group 3007 // freqgrp[numgrp][nummember] 3008 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...], 3009 // [ifgrp10, ifgrp11, ifgrp12, ...], 3010 // ... 3011 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]] 3012 // freqrange[numgrp*2] 3013 // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...] 3014 vector< vector<uInt> > freqgrp ; 3015 double freqrange = 0.0 ; 3016 uInt grpnum = 0 ; 3017 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) { 3018 // Assumed that ifgfreq was sorted 3019 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) { 3020 freqgrp[grpnum-1].push_back( i ) ; 3021 } 3022 else { 3023 vector<uInt> grp0( 1, i ) ; 3024 freqgrp.push_back( grp0 ) ; 3025 grpnum++ ; 3026 } 3027 freqrange = ifgfreq[2*i+1] ; 3028 } 3029 3030 3031 // print IF groups 3032 ostringstream oss ; 3033 oss << "IF Group summary: " << endl ; 3034 oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ; 3035 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) { 3036 oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ; 3037 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) { 3038 oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 3039 } 3040 oss << endl ; 3041 } 3042 oss << endl ; 3043 os << oss.str() << LogIO::POST ; 3044 3045 // print frequency group 3046 oss.str("") ; 3047 oss << "Frequency Group summary: " << endl ; 3048 oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ; 3049 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) { 3050 oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ; 3051 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) { 3052 oss << freqgrp[i][j] << " " ; 3053 } 3054 oss << endl ; 3055 } 3056 oss << endl ; 3057 os << oss.str() << LogIO::POST ; 3058 3059 // membership check 3060 // groups: list of IF group indexes whose frequency range overlaps with 3061 // that of each table and IF 3062 // groups[numtable][numIF][nummembership] 3063 // groups: [[[grp, grp,...], [grp, grp,...],...], 3064 // [[grp, grp,...], [grp, grp,...],...], 3065 // ... 3066 // [[grp, grp,...], [grp, grp,...],...]] 3067 vector< vector< vector<uInt> > > groups( insize ) ; 3068 for ( uInt i = 0 ; i < insize ; i++ ) { 3069 groups[i].resize( freqid[i].size() ) ; 3070 } 3071 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) { 3072 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) { 3073 uInt tableid = ifgrp[igrp][2*imem] ; 3074 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ; 3075 if ( iter != freqid[tableid].end() ) { 3076 uInt rowid = distance( freqid[tableid].begin(), iter ) ; 3077 groups[tableid][rowid].push_back( igrp ) ; 3078 } 3079 } 3080 } 3081 3082 // print membership 3083 //oss.str("") ; 3084 //for ( uInt i = 0 ; i < insize ; i++ ) { 3085 //oss << "Table " << i << endl ; 3086 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) { 3087 //oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ; 3088 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) { 3089 //oss << setw( 2 ) << groups[i][j][k] << " " ; 3090 //} 3091 //oss << endl ; 3092 //} 3093 //} 3094 //os << oss.str() << LogIO::POST ; 3095 3096 // set back coordinfo 3097 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3098 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ; 3099 coordinfo[0] = oldinfo[itable] ; 3100 tmpin[itable]->setCoordInfo( coordinfo ) ; 3101 } 3102 3103 // Create additional table if needed 3104 bool oldInsitu = insitu_ ; 3105 setInsitu( false ) ; 3106 vector< vector<uInt> > addrow( insize ) ; 3107 vector<uInt> addtable( insize, 0 ) ; 3108 vector<uInt> newtableids( insize ) ; 3109 vector<uInt> newifids( insize, 0 ) ; 3110 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3111 //os << "Table " << itable << ": " ; 3112 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) { 3113 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ; 3114 //os << addrow[itable][ifrow] << " " ; 3115 } 3116 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ; 3117 //os << "(" << addtable[itable] << ")" << LogIO::POST ; 3118 } 3119 newin.resize( insize ) ; 3120 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ; 3121 for ( uInt i = 0 ; i < insize ; i++ ) { 3122 newtableids[i] = i ; 3123 } 3124 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3125 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) { 3126 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ; 3127 vector<int> freqidlist ; 3128 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) { 3129 if ( groups[itable][i].size() > iadd + 1 ) { 3130 freqidlist.push_back( freqid[itable][i] ) ; 3131 } 3132 } 3133 stringstream taqlstream ; 3134 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ; 3135 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) { 3136 taqlstream << i ; 3137 if ( i < freqidlist.size() - 1 ) 3138 taqlstream << "," ; 3139 else 3140 taqlstream << "]" ; 3141 } 3142 string taql = taqlstream.str() ; 3143 //os << "taql = " << taql << LogIO::POST ; 3144 STSelector selector = STSelector() ; 3145 selector.setTaQL( taql ) ; 3146 add->setSelection( selector ) ; 3147 newin.push_back( add ) ; 3148 newtableids.push_back( itable ) ; 3149 newifids.push_back( iadd + 1 ) ; 3150 } 3151 } 3152 3153 // udpate ifgrp 3154 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3155 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) { 3156 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) { 3157 if ( groups[itable][ifrow].size() > iadd + 1 ) { 3158 uInt igrp = groups[itable][ifrow][iadd+1] ; 3159 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) { 3160 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) { 3161 ifgrp[igrp][2*imem] = insize + iadd ; 3162 } 3163 } 3164 } 3165 } 3166 } 3167 } 3168 3169 // print IF groups again for debug 3170 //oss.str( "" ) ; 3171 //oss << "IF Group summary: " << endl ; 3172 //oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ; 3173 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) { 3174 //oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ; 3175 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) { 3176 //oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 3177 //} 3178 //oss << endl ; 3179 //} 3180 //oss << endl ; 3181 //os << oss.str() << LogIO::POST ; 3182 3183 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation 3184 os << "All scan number is set to 0" << LogIO::POST ; 3185 //os << "All IF number is set to IF group index" << LogIO::POST ; 3186 insize = newin.size() ; 3187 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3188 uInt rows = newin[itable]->nrow() ; 3189 Table &tmpt = newin[itable]->table() ; 3190 freqIDCol.attach( tmpt, "FREQ_ID" ) ; 3191 scannoCol.attach( tmpt, "SCANNO" ) ; 3192 ifnoCol.attach( tmpt, "IFNO" ) ; 3193 for ( uInt irow=0 ; irow < rows ; irow++ ) { 3194 scannoCol.put( irow, 0 ) ; 3195 uInt freqID = freqIDCol( irow ) ; 3196 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ; 3197 if ( iter != freqid[newtableids[itable]].end() ) { 3198 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ; 3199 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ; 3200 } 3201 else { 3202 throw(AipsError("IF grouping was wrong in additional tables.")) ; 3203 } 3204 } 3205 } 3206 oldinfo.resize( insize ) ; 3207 setInsitu( oldInsitu ) ; 3208 3209 // temporarily set coordinfo 3210 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3211 vector<string> coordinfo = newin[itable]->getCoordInfo() ; 3212 oldinfo[itable] = coordinfo[0] ; 3213 coordinfo[0] = "Hz" ; 3214 newin[itable]->setCoordInfo( coordinfo ) ; 3215 } 3216 3217 // save column values in the vector 3218 vector< vector<uInt> > freqTableIdVec( insize ) ; 3219 vector< vector<uInt> > freqIdVec( insize ) ; 3220 vector< vector<uInt> > ifNoVec( insize ) ; 3221 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3222 ScalarColumn<uInt> freqIDs ; 3223 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ; 3224 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ; 3225 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ; 3226 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) { 3227 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ; 3228 } 3229 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) { 3230 freqIdVec[itable].push_back( freqIDCol( irow ) ) ; 3231 ifNoVec[itable].push_back( ifnoCol( irow ) ) ; 3232 } 3233 } 3234 3235 // reset spectra and flagtra: pick up common part of frequency coverage 3236 //os << "Pick common frequency range and align resolution" << LogIO::POST ; 3237 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3238 uInt rows = newin[itable]->nrow() ; 3239 int nminchan = -1 ; 3240 int nmaxchan = -1 ; 3241 vector<uInt> freqIdUpdate ; 3242 for ( uInt irow = 0 ; irow < rows ; irow++ ) { 3243 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index 3244 double minfreq = ifgfreq[2*ifno] ; 3245 double maxfreq = ifgfreq[2*ifno+1] ; 3246 //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ; 3247 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ; 3248 int nchan = abcissa.size() ; 3249 double resol = abcissa[1] - abcissa[0] ; 3250 //os << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ; 3251 if ( minfreq <= abcissa[0] ) 3252 nminchan = 0 ; 3253 else { 3254 //double cfreq = ( minfreq - abcissa[0] ) / resol ; 3255 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ; 3256 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ; 3257 } 3258 if ( maxfreq >= abcissa[abcissa.size()-1] ) 3259 nmaxchan = abcissa.size() - 1 ; 3260 else { 3261 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ; 3262 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ; 3263 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ; 3264 } 3265 //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ; 3266 if ( nmaxchan > nminchan ) { 3267 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ; 3268 int newchan = nmaxchan - nminchan + 1 ; 3269 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan ) 3270 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ; 3271 } 3272 else { 3273 throw(AipsError("Failed to pick up common part of frequency range.")) ; 3274 } 3275 } 3276 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) { 3277 uInt freqId = freqIdUpdate[i] ; 3278 Double refpix ; 3279 Double refval ; 3280 Double increment ; 3281 3282 // update row 3283 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 3284 refval = refval - ( refpix - nminchan ) * increment ; 3285 refpix = 0 ; 3286 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ; 3287 } 3288 } 3289 3290 3291 // reset spectra and flagtra: align spectral resolution 3292 //os << "Align spectral resolution" << LogIO::POST ; 3293 // gmaxdnu: the coarsest frequency resolution in the frequency group 3294 // gmemid: member index that have a resolution equal to gmaxdnu 3295 // gmaxdnu[numfreqgrp] 3296 // gmaxdnu: [dnu0, dnu1, ...] 3297 // gmemid[numfreqgrp] 3298 // gmemid: [id0, id1, ...] 3299 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ; 3300 vector<uInt> gmemid( freqgrp.size(), 0 ) ; 3301 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) { 3302 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution 3303 int minchan = INT_MAX ; // minimum channel number 3304 Double refpixref = -1 ; // reference of 'reference pixel' 3305 Double refvalref = -1 ; // reference of 'reference frequency' 3306 Double refinc = -1 ; // reference frequency resolution 3307 uInt refreqid ; 3308 uInt reftable = INT_MAX; 3309 // process only if group member > 1 3310 if ( ifgrp[igrp].size() > 2 ) { 3311 // find minchan and maxdnu in each group 3312 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) { 3313 uInt tableid = ifgrp[igrp][2*imem] ; 3314 uInt rowid = ifgrp[igrp][2*imem+1] ; 3315 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ; 3316 if ( iter != freqIdVec[tableid].end() ) { 3317 uInt index = distance( freqIdVec[tableid].begin(), iter ) ; 3318 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ; 3319 int nchan = abcissa.size() ; 3320 double dnu = abcissa[1] - abcissa[0] ; 3321 //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ; 3322 if ( nchan < minchan ) { 3323 minchan = nchan ; 3324 maxdnu = dnu ; 3325 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ; 3326 refreqid = rowid ; 3327 reftable = tableid ; 3328 } 3329 } 3330 } 3331 // regrid spectra in each group 3332 os << "GROUP " << igrp << endl ; 3333 os << " Channel number is adjusted to " << minchan << endl ; 3334 os << " Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ; 3335 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) { 3336 uInt tableid = ifgrp[igrp][2*imem] ; 3337 uInt rowid = ifgrp[igrp][2*imem+1] ; 3338 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ; 3339 //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ; 3340 //os << " regridChannel applied to " ; 3341 if ( tableid != reftable ) 3342 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ; 3343 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) { 3344 uInt tfreqid = freqIdVec[tableid][irow] ; 3345 if ( tfreqid == rowid ) { 3346 //os << irow << " " ; 3347 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ; 3348 freqIDCol.put( irow, refreqid ) ; 3349 freqIdVec[tableid][irow] = refreqid ; 3350 } 3351 } 3352 //os << LogIO::POST ; 3353 } 3354 } 3355 else { 3356 uInt tableid = ifgrp[igrp][0] ; 3357 uInt rowid = ifgrp[igrp][1] ; 3358 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ; 3359 if ( iter != freqIdVec[tableid].end() ) { 3360 uInt index = distance( freqIdVec[tableid].begin(), iter ) ; 3361 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ; 3362 minchan = abcissa.size() ; 3363 maxdnu = abcissa[1] - abcissa[0] ; 3364 } 3365 } 3366 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) { 3367 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) { 3368 if ( maxdnu > gmaxdnu[i] ) { 3369 gmaxdnu[i] = maxdnu ; 3370 gmemid[i] = igrp ; 3371 } 3372 break ; 3373 } 3374 } 3375 } 3376 3377 // set back coordinfo 3378 for ( uInt itable = 0 ; itable < insize ; itable++ ) { 3379 vector<string> coordinfo = newin[itable]->getCoordInfo() ; 3380 coordinfo[0] = oldinfo[itable] ; 3381 newin[itable]->setCoordInfo( coordinfo ) ; 3382 } 3383 3384 // accumulate all rows into the first table 3385 // NOTE: assumed in.size() = 1 3386 vector< CountedPtr<Scantable> > tmp( 1 ) ; 3387 if ( newin.size() == 1 ) 3388 tmp[0] = newin[0] ; 3389 else 3390 tmp[0] = merge( newin ) ; 3391 3392 //return tmp[0] ; 3393 3394 // average 3395 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ; 3396 3397 //return tmpout ; 3398 3399 // combine frequency group 3400 os << "Combine spectra based on frequency grouping" << LogIO::POST ; 3401 os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ; 3402 vector<string> coordinfo = tmpout->getCoordInfo() ; 3403 oldinfo[0] = coordinfo[0] ; 3404 coordinfo[0] = "Hz" ; 3405 tmpout->setCoordInfo( coordinfo ) ; 3406 // create proformas of output table 3407 stringstream taqlstream ; 3408 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ; 3409 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) { 3410 taqlstream << gmemid[i] ; 3411 if ( i < gmemid.size() - 1 ) 3412 taqlstream << "," ; 3413 else 3414 taqlstream << "]" ; 3415 } 3416 string taql = taqlstream.str() ; 3417 //os << "taql = " << taql << LogIO::POST ; 3418 STSelector selector = STSelector() ; 3419 selector.setTaQL( taql ) ; 3420 oldInsitu = insitu_ ; 3421 setInsitu( false ) ; 3422 out = getScantable( tmpout, false ) ; 3423 setInsitu( oldInsitu ) ; 3424 out->setSelection( selector ) ; 3425 // regrid rows 3426 ifnoCol.attach( tmpout->table(), "IFNO" ) ; 3427 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) { 3428 uInt ifno = ifnoCol( irow ) ; 3429 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) { 3430 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) { 3431 vector<double> abcissa = tmpout->getAbcissa( irow ) ; 3432 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ; 3433 int nchan = (int)( bw / gmaxdnu[igrp] ) ; 3434 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ; 3435 break ; 3436 } 3437 } 3438 } 3439 // combine spectra 3440 ArrayColumn<Float> specColOut ; 3441 specColOut.attach( out->table(), "SPECTRA" ) ; 3442 ArrayColumn<uChar> flagColOut ; 3443 flagColOut.attach( out->table(), "FLAGTRA" ) ; 3444 ScalarColumn<uInt> ifnoColOut ; 3445 ifnoColOut.attach( out->table(), "IFNO" ) ; 3446 ScalarColumn<uInt> polnoColOut ; 3447 polnoColOut.attach( out->table(), "POLNO" ) ; 3448 ScalarColumn<uInt> freqidColOut ; 3449 freqidColOut.attach( out->table(), "FREQ_ID" ) ; 3450 MDirection::ScalarColumn dirColOut ; 3451 dirColOut.attach( out->table(), "DIRECTION" ) ; 3452 Table &tab = tmpout->table() ; 3453 Block<String> cols(1); 3454 cols[0] = String("POLNO") ; 3455 TableIterator iter( tab, cols ) ; 3456 bool done = false ; 3457 vector< vector<uInt> > sizes( freqgrp.size() ) ; 3458 while( !iter.pastEnd() ) { 3459 vector< vector<Float> > specout( freqgrp.size() ) ; 3460 vector< vector<uChar> > flagout( freqgrp.size() ) ; 3461 ArrayColumn<Float> specCols ; 3462 specCols.attach( iter.table(), "SPECTRA" ) ; 3463 ArrayColumn<uChar> flagCols ; 3464 flagCols.attach( iter.table(), "FLAGTRA" ) ; 3465 ifnoCol.attach( iter.table(), "IFNO" ) ; 3466 ScalarColumn<uInt> polnos ; 3467 polnos.attach( iter.table(), "POLNO" ) ; 3468 MDirection::ScalarColumn dircol ; 3469 dircol.attach( iter.table(), "DIRECTION" ) ; 3470 uInt polno = polnos( 0 ) ; 3471 //os << "POLNO iteration: " << polno << LogIO::POST ; 3472 // for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) { 3473 // sizes[igrp].resize( freqgrp[igrp].size() ) ; 3474 // for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) { 3475 // for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) { 3476 // uInt ifno = ifnoCol( irow ) ; 3477 // if ( ifno == freqgrp[igrp][imem] ) { 3478 // Vector<Float> spec = specCols( irow ) ; 3479 // Vector<uChar> flag = flagCols( irow ) ; 3480 // vector<Float> svec ; 3481 // spec.tovector( svec ) ; 3482 // vector<uChar> fvec ; 3483 // flag.tovector( fvec ) ; 3484 // //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ; 3485 // specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ; 3486 // flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ; 3487 // //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ; 3488 // sizes[igrp][imem] = spec.nelements() ; 3489 // } 3490 // } 3491 // } 3492 // for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) { 3493 // uInt ifout = ifnoColOut( irow ) ; 3494 // uInt polout = polnoColOut( irow ) ; 3495 // if ( ifout == gmemid[igrp] && polout == polno ) { 3496 // // set SPECTRA and FRAGTRA 3497 // Vector<Float> newspec( specout[igrp] ) ; 3498 // Vector<uChar> newflag( flagout[igrp] ) ; 3499 // specColOut.put( irow, newspec ) ; 3500 // flagColOut.put( irow, newflag ) ; 3501 // // IFNO renumbering 3502 // ifnoColOut.put( irow, igrp ) ; 3503 // } 3504 // } 3505 // } 3506 // get a list of number of channels for each frequency group member 3507 if ( !done ) { 3508 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) { 3509 sizes[igrp].resize( freqgrp[igrp].size() ) ; 3510 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) { 3511 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) { 3512 uInt ifno = ifnoCol( irow ) ; 3513 if ( ifno == freqgrp[igrp][imem] ) { 3514 Vector<Float> spec = specCols( irow ) ; 3515 sizes[igrp][imem] = spec.nelements() ; 3516 break ; 3517 } 3518 } 3519 } 3520 } 3521 done = true ; 3522 } 3523 // combine spectra 3524 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) { 3525 uInt polout = polnoColOut( irow ) ; 3526 if ( polout == polno ) { 3527 uInt ifout = ifnoColOut( irow ) ; 3528 Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ; 3529 uInt igrp ; 3530 for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) { 3531 if ( ifout == gmemid[jgrp] ) { 3532 igrp = jgrp ; 3533 break ; 3534 } 3535 } 3536 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) { 3537 for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) { 3538 uInt ifno = ifnoCol( jrow ) ; 3539 Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ; 3540 //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) { 3541 Double dx = tdir[0] - direction[0] ; 3542 Double dy = tdir[1] - direction[1] ; 3543 Double dd = sqrt( dx * dx + dy * dy ) ; 3544 //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) { 3545 if ( ifno == freqgrp[igrp][imem] && dd <= tol ) { 3546 Vector<Float> spec = specCols( jrow ) ; 3547 Vector<uChar> flag = flagCols( jrow ) ; 3548 vector<Float> svec ; 3549 spec.tovector( svec ) ; 3550 vector<uChar> fvec ; 3551 flag.tovector( fvec ) ; 3552 //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ; 3553 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ; 3554 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ; 3555 //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ; 3556 } 3557 } 3558 } 3559 // set SPECTRA and FRAGTRA 3560 Vector<Float> newspec( specout[igrp] ) ; 3561 Vector<uChar> newflag( flagout[igrp] ) ; 3562 specColOut.put( irow, newspec ) ; 3563 flagColOut.put( irow, newflag ) ; 3564 // IFNO renumbering 3565 ifnoColOut.put( irow, igrp ) ; 3566 } 3567 } 3568 iter++ ; 3569 } 3570 // update FREQUENCIES subtable 3571 vector<bool> updated( freqgrp.size(), false ) ; 3572 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) { 3573 uInt index = 0 ; 3574 uInt pixShift = 0 ; 3575 while ( freqgrp[igrp][index] != gmemid[igrp] ) { 3576 pixShift += sizes[igrp][index++] ; 3577 } 3578 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) { 3579 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) { 3580 uInt freqidOut = freqidColOut( irow ) ; 3581 //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ; 3582 double refpix ; 3583 double refval ; 3584 double increm ; 3585 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ; 3586 refpix += pixShift ; 3587 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ; 3588 updated[igrp] = true ; 3589 } 3590 } 3591 } 3592 3593 //out = tmpout ; 3594 3595 coordinfo = tmpout->getCoordInfo() ; 3596 coordinfo[0] = oldinfo[0] ; 3597 tmpout->setCoordInfo( coordinfo ) ; 3598 } 3599 else { 3600 // simple average 3601 out = average( in, mask, weight, avmode ) ; 3602 } 3603 3604 return out ; 3605 } 3606 3607 CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s, 3608 const String calmode, 3609 const String antname ) 3610 { 3611 // frequency switch 3612 if ( calmode == "fs" ) { 3613 return cwcalfs( s, antname ) ; 3614 } 3615 else { 3616 vector<bool> masks = s->getMask( 0 ) ; 3617 vector<int> types ; 3618 3619 // sky scan 3620 STSelector sel = STSelector() ; 3621 types.push_back( SrcType::SKY ) ; 3622 sel.setTypes( types ) ; 3623 s->setSelection( sel ) ; 3624 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ; 3625 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ; 3626 s->unsetSelection() ; 3627 sel.reset() ; 3628 types.clear() ; 3629 3630 // hot scan 3631 types.push_back( SrcType::HOT ) ; 3632 sel.setTypes( types ) ; 3633 s->setSelection( sel ) ; 3634 tmp.clear() ; 3635 tmp.push_back( getScantable( s, false ) ) ; 3636 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ; 3637 s->unsetSelection() ; 3638 sel.reset() ; 3639 types.clear() ; 3640 3641 // cold scan 3642 CountedPtr<Scantable> acold ; 3643 // types.push_back( SrcType::COLD ) ; 3644 // sel.setTypes( types ) ; 3645 // s->setSelection( sel ) ; 3646 // tmp.clear() ; 3647 // tmp.push_back( getScantable( s, false ) ) ; 3648 // CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ; 3649 // s->unsetSelection() ; 3650 // sel.reset() ; 3651 // types.clear() ; 3652 3653 // off scan 3654 types.push_back( SrcType::PSOFF ) ; 3655 sel.setTypes( types ) ; 3656 s->setSelection( sel ) ; 3657 tmp.clear() ; 3658 tmp.push_back( getScantable( s, false ) ) ; 3659 CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ; 3660 s->unsetSelection() ; 3661 sel.reset() ; 3662 types.clear() ; 3663 3664 // on scan 3665 bool insitu = insitu_ ; 3666 insitu_ = false ; 3667 CountedPtr<Scantable> out = getScantable( s, true ) ; 3668 insitu_ = insitu ; 3669 types.push_back( SrcType::PSON ) ; 3670 sel.setTypes( types ) ; 3671 s->setSelection( sel ) ; 3672 TableCopy::copyRows( out->table(), s->table() ) ; 3673 s->unsetSelection() ; 3674 sel.reset() ; 3675 types.clear() ; 3676 3677 // process each on scan 3678 ArrayColumn<Float> tsysCol ; 3679 tsysCol.attach( out->table(), "TSYS" ) ; 3680 for ( int i = 0 ; i < out->nrow() ; i++ ) { 3681 vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ; 3682 out->setSpectrum( sp, i ) ; 3683 string reftime = out->getTime( i ) ; 3684 vector<int> ii( 1, out->getIF( i ) ) ; 3685 vector<int> ib( 1, out->getBeam( i ) ) ; 3686 vector<int> ip( 1, out->getPol( i ) ) ; 3687 sel.setIFs( ii ) ; 3688 sel.setBeams( ib ) ; 3689 sel.setPolarizations( ip ) ; 3690 asky->setSelection( sel ) ; 3691 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ; 3692 const Vector<Float> Vtsys( sptsys ) ; 3693 tsysCol.put( i, Vtsys ) ; 3694 asky->unsetSelection() ; 3695 sel.reset() ; 3696 } 3697 3698 // flux unit 3699 out->setFluxUnit( "K" ) ; 3700 3701 return out ; 3702 } 3703 } 3704 3705 CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s, 3706 const String calmode ) 3707 { 3708 // frequency switch 3709 if ( calmode == "fs" ) { 3710 return almacalfs( s ) ; 3711 } 3712 else { 3713 vector<bool> masks = s->getMask( 0 ) ; 3714 3715 // off scan 3716 STSelector sel = STSelector() ; 3717 vector<int> types ; 3718 types.push_back( SrcType::PSOFF ) ; 3719 sel.setTypes( types ) ; 3720 s->setSelection( sel ) ; 3721 // TODO 2010/01/08 TN 3722 // Grouping by time should be needed before averaging. 3723 // Each group must have own unique SCANNO (should be renumbered). 3724 // See PIPELINE/SDCalibration.py 3725 CountedPtr<Scantable> soff = getScantable( s, false ) ; 3726 Table ttab = soff->table() ; 3727 ROScalarColumn<Double> timeCol( ttab, "TIME" ) ; 3728 uInt nrow = timeCol.nrow() ; 3729 Vector<Double> timeSep( nrow - 1 ) ; 3730 for ( uInt i = 0 ; i < nrow - 1 ; i++ ) { 3731 timeSep[i] = timeCol(i+1) - timeCol(i) ; 3732 } 3733 ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ; 3734 Vector<Double> interval = intervalCol.getColumn() ; 3735 interval /= 86400.0 ; 3736 ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ; 3737 vector<uInt> glist ; 3738 for ( uInt i = 0 ; i < nrow - 1 ; i++ ) { 3739 double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ; 3740 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ; 3741 if ( gap > 1.1 ) { 3742 glist.push_back( i ) ; 3743 } 3744 } 3745 Vector<uInt> gaplist( glist ) ; 3746 //cout << "gaplist = " << gaplist << endl ; 3747 uInt newid = 0 ; 3748 for ( uInt i = 0 ; i < nrow ; i++ ) { 3749 scanCol.put( i, newid ) ; 3750 if ( i == gaplist[newid] ) { 3751 newid++ ; 3752 } 3753 } 3754 //cout << "new scancol = " << scanCol.getColumn() << endl ; 3755 vector< CountedPtr<Scantable> > tmp( 1, soff ) ; 3756 CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ; 3757 //cout << "aoff.nrow = " << aoff->nrow() << endl ; 3758 s->unsetSelection() ; 3759 sel.reset() ; 3760 types.clear() ; 3761 3762 // on scan 3763 bool insitu = insitu_ ; 3764 insitu_ = false ; 3765 CountedPtr<Scantable> out = getScantable( s, true ) ; 3766 insitu_ = insitu ; 3767 types.push_back( SrcType::PSON ) ; 3768 sel.setTypes( types ) ; 3769 s->setSelection( sel ) ; 3770 TableCopy::copyRows( out->table(), s->table() ) ; 3771 s->unsetSelection() ; 3772 sel.reset() ; 3773 types.clear() ; 3774 3775 // process each on scan 3776 ArrayColumn<Float> tsysCol ; 3777 tsysCol.attach( out->table(), "TSYS" ) ; 3778 for ( int i = 0 ; i < out->nrow() ; i++ ) { 3779 vector<float> sp = getCalibratedSpectra( out, aoff, i ) ; 3780 out->setSpectrum( sp, i ) ; 3781 } 3782 3783 // flux unit 3784 out->setFluxUnit( "K" ) ; 3785 3786 return out ; 3787 } 3788 } 3789 3790 CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s, 3791 const String antname ) 3792 { 3793 vector<int> types ; 3794 3795 // APEX calibration mode 3796 int apexcalmode = 1 ; 3797 3798 if ( antname.find( "APEX" ) != string::npos ) { 3799 // check if off scan exists or not 3800 STSelector sel = STSelector() ; 3801 //sel.setName( offstr1 ) ; 3802 types.push_back( SrcType::FLOOFF ) ; 3803 sel.setTypes( types ) ; 3804 try { 3805 s->setSelection( sel ) ; 3806 } 3807 catch ( AipsError &e ) { 3808 apexcalmode = 0 ; 3809 } 3810 sel.reset() ; 3811 } 3812 s->unsetSelection() ; 3813 types.clear() ; 3814 3815 vector<bool> masks = s->getMask( 0 ) ; 3816 CountedPtr<Scantable> ssig, sref ; 3817 CountedPtr<Scantable> out ; 3818 3819 if ( antname.find( "APEX" ) != string::npos ) { 3820 // APEX calibration 3821 // sky scan 3822 STSelector sel = STSelector() ; 3823 types.push_back( SrcType::FLOSKY ) ; 3824 sel.setTypes( types ) ; 3825 s->setSelection( sel ) ; 3826 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ; 3827 CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ; 3828 s->unsetSelection() ; 3829 sel.reset() ; 3830 types.clear() ; 3831 types.push_back( SrcType::FHISKY ) ; 3832 sel.setTypes( types ) ; 3833 s->setSelection( sel ) ; 3834 tmp.clear() ; 3835 tmp.push_back( getScantable( s, false ) ) ; 3836 CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ; 3837 s->unsetSelection() ; 3838 sel.reset() ; 3839 types.clear() ; 3840 3841 // hot scan 3842 types.push_back( SrcType::FLOHOT ) ; 3843 sel.setTypes( types ) ; 3844 s->setSelection( sel ) ; 3845 tmp.clear() ; 3846 tmp.push_back( getScantable( s, false ) ) ; 3847 CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ; 3848 s->unsetSelection() ; 3849 sel.reset() ; 3850 types.clear() ; 3851 types.push_back( SrcType::FHIHOT ) ; 3852 sel.setTypes( types ) ; 3853 s->setSelection( sel ) ; 3854 tmp.clear() ; 3855 tmp.push_back( getScantable( s, false ) ) ; 3856 CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ; 3857 s->unsetSelection() ; 3858 sel.reset() ; 3859 types.clear() ; 3860 3861 // cold scan 3862 CountedPtr<Scantable> acoldlo, acoldhi ; 3863 // types.push_back( SrcType::FLOCOLD ) ; 3864 // sel.setTypes( types ) ; 3865 // s->setSelection( sel ) ; 3866 // tmp.clear() ; 3867 // tmp.push_back( getScantable( s, false ) ) ; 3868 // CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ; 3869 // s->unsetSelection() ; 3870 // sel.reset() ; 3871 // types.clear() ; 3872 // types.push_back( SrcType::FHICOLD ) ; 3873 // sel.setTypes( types ) ; 3874 // s->setSelection( sel ) ; 3875 // tmp.clear() ; 3876 // tmp.push_back( getScantable( s, false ) ) ; 3877 // CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ; 3878 // s->unsetSelection() ; 3879 // sel.reset() ; 3880 // types.clear() ; 3881 3882 // ref scan 3883 bool insitu = insitu_ ; 3884 insitu_ = false ; 3885 sref = getScantable( s, true ) ; 3886 insitu_ = insitu ; 3887 types.push_back( SrcType::FSLO ) ; 3888 sel.setTypes( types ) ; 3889 s->setSelection( sel ) ; 3890 TableCopy::copyRows( sref->table(), s->table() ) ; 3891 s->unsetSelection() ; 3892 sel.reset() ; 3893 types.clear() ; 3894 3895 // sig scan 3896 insitu_ = false ; 3897 ssig = getScantable( s, true ) ; 3898 insitu_ = insitu ; 3899 types.push_back( SrcType::FSHI ) ; 3900 sel.setTypes( types ) ; 3901 s->setSelection( sel ) ; 3902 TableCopy::copyRows( ssig->table(), s->table() ) ; 3903 s->unsetSelection() ; 3904 sel.reset() ; 3905 types.clear() ; 3906 3907 if ( apexcalmode == 0 ) { 3908 // APEX fs data without off scan 3909 // process each sig and ref scan 3910 ArrayColumn<Float> tsysCollo ; 3911 tsysCollo.attach( ssig->table(), "TSYS" ) ; 3912 ArrayColumn<Float> tsysColhi ; 3913 tsysColhi.attach( sref->table(), "TSYS" ) ; 3914 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 3915 vector< CountedPtr<Scantable> > sky( 2 ) ; 3916 sky[0] = askylo ; 3917 sky[1] = askyhi ; 3918 vector< CountedPtr<Scantable> > hot( 2 ) ; 3919 hot[0] = ahotlo ; 3920 hot[1] = ahothi ; 3921 vector< CountedPtr<Scantable> > cold( 2 ) ; 3922 //cold[0] = acoldlo ; 3923 //cold[1] = acoldhi ; 3924 vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ; 3925 ssig->setSpectrum( sp, i ) ; 3926 string reftime = ssig->getTime( i ) ; 3927 vector<int> ii( 1, ssig->getIF( i ) ) ; 3928 vector<int> ib( 1, ssig->getBeam( i ) ) ; 3929 vector<int> ip( 1, ssig->getPol( i ) ) ; 3930 sel.setIFs( ii ) ; 3931 sel.setBeams( ib ) ; 3932 sel.setPolarizations( ip ) ; 3933 askylo->setSelection( sel ) ; 3934 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ; 3935 const Vector<Float> Vtsyslo( sptsys ) ; 3936 tsysCollo.put( i, Vtsyslo ) ; 3937 askylo->unsetSelection() ; 3938 sel.reset() ; 3939 sky[0] = askyhi ; 3940 sky[1] = askylo ; 3941 hot[0] = ahothi ; 3942 hot[1] = ahotlo ; 3943 cold[0] = acoldhi ; 3944 cold[1] = acoldlo ; 3945 sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ; 3946 sref->setSpectrum( sp, i ) ; 3947 reftime = sref->getTime( i ) ; 3948 ii[0] = sref->getIF( i ) ; 3949 ib[0] = sref->getBeam( i ) ; 3950 ip[0] = sref->getPol( i ) ; 3951 sel.setIFs( ii ) ; 3952 sel.setBeams( ib ) ; 3953 sel.setPolarizations( ip ) ; 3954 askyhi->setSelection( sel ) ; 3955 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ; 3956 const Vector<Float> Vtsyshi( sptsys ) ; 3957 tsysColhi.put( i, Vtsyshi ) ; 3958 askyhi->unsetSelection() ; 3959 sel.reset() ; 3960 } 3961 } 3962 else if ( apexcalmode == 1 ) { 3963 // APEX fs data with off scan 3964 // off scan 3965 types.push_back( SrcType::FLOOFF ) ; 3966 sel.setTypes( types ) ; 3967 s->setSelection( sel ) ; 3968 tmp.clear() ; 3969 tmp.push_back( getScantable( s, false ) ) ; 3970 CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ; 3971 s->unsetSelection() ; 3972 sel.reset() ; 3973 types.clear() ; 3974 types.push_back( SrcType::FHIOFF ) ; 3975 sel.setTypes( types ) ; 3976 s->setSelection( sel ) ; 3977 tmp.clear() ; 3978 tmp.push_back( getScantable( s, false ) ) ; 3979 CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ; 3980 s->unsetSelection() ; 3981 sel.reset() ; 3982 types.clear() ; 3983 3984 // process each sig and ref scan 3985 ArrayColumn<Float> tsysCollo ; 3986 tsysCollo.attach( ssig->table(), "TSYS" ) ; 3987 ArrayColumn<Float> tsysColhi ; 3988 tsysColhi.attach( sref->table(), "TSYS" ) ; 3989 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 3990 vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ; 3991 ssig->setSpectrum( sp, i ) ; 3992 sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ; 3993 string reftime = ssig->getTime( i ) ; 3994 vector<int> ii( 1, ssig->getIF( i ) ) ; 3995 vector<int> ib( 1, ssig->getBeam( i ) ) ; 3996 vector<int> ip( 1, ssig->getPol( i ) ) ; 3997 sel.setIFs( ii ) ; 3998 sel.setBeams( ib ) ; 3999 sel.setPolarizations( ip ) ; 4000 askylo->setSelection( sel ) ; 4001 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ; 4002 const Vector<Float> Vtsyslo( sptsys ) ; 4003 tsysCollo.put( i, Vtsyslo ) ; 4004 askylo->unsetSelection() ; 4005 sel.reset() ; 4006 sref->setSpectrum( sp, i ) ; 4007 reftime = sref->getTime( i ) ; 4008 ii[0] = sref->getIF( i ) ; 4009 ib[0] = sref->getBeam( i ) ; 4010 ip[0] = sref->getPol( i ) ; 4011 sel.setIFs( ii ) ; 4012 sel.setBeams( ib ) ; 4013 sel.setPolarizations( ip ) ; 4014 askyhi->setSelection( sel ) ; 4015 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ; 4016 const Vector<Float> Vtsyshi( sptsys ) ; 4017 tsysColhi.put( i, Vtsyshi ) ; 4018 askyhi->unsetSelection() ; 4019 sel.reset() ; 4020 } 4021 } 4022 } 4023 else { 4024 // non-APEX fs data 4025 // sky scan 4026 STSelector sel = STSelector() ; 4027 types.push_back( SrcType::SKY ) ; 4028 sel.setTypes( types ) ; 4029 s->setSelection( sel ) ; 4030 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ; 4031 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ; 4032 s->unsetSelection() ; 4033 sel.reset() ; 4034 types.clear() ; 4035 4036 // hot scan 4037 types.push_back( SrcType::HOT ) ; 4038 sel.setTypes( types ) ; 4039 s->setSelection( sel ) ; 4040 tmp.clear() ; 4041 tmp.push_back( getScantable( s, false ) ) ; 4042 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ; 4043 s->unsetSelection() ; 4044 sel.reset() ; 4045 types.clear() ; 4046 4047 // cold scan 4048 CountedPtr<Scantable> acold ; 4049 // types.push_back( SrcType::COLD ) ; 4050 // sel.setTypes( types ) ; 4051 // s->setSelection( sel ) ; 4052 // tmp.clear() ; 4053 // tmp.push_back( getScantable( s, false ) ) ; 4054 // CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ; 4055 // s->unsetSelection() ; 4056 // sel.reset() ; 4057 // types.clear() ; 4058 4059 // ref scan 4060 bool insitu = insitu_ ; 4061 insitu_ = false ; 4062 sref = getScantable( s, true ) ; 4063 insitu_ = insitu ; 4064 types.push_back( SrcType::FSOFF ) ; 4065 sel.setTypes( types ) ; 4066 s->setSelection( sel ) ; 4067 TableCopy::copyRows( sref->table(), s->table() ) ; 4068 s->unsetSelection() ; 4069 sel.reset() ; 4070 types.clear() ; 4071 4072 // sig scan 4073 insitu_ = false ; 4074 ssig = getScantable( s, true ) ; 4075 insitu_ = insitu ; 4076 types.push_back( SrcType::FSON ) ; 4077 sel.setTypes( types ) ; 4078 s->setSelection( sel ) ; 4079 TableCopy::copyRows( ssig->table(), s->table() ) ; 4080 s->unsetSelection() ; 4081 sel.reset() ; 4082 types.clear() ; 4083 4084 // process each sig and ref scan 4085 ArrayColumn<Float> tsysColsig ; 4086 tsysColsig.attach( ssig->table(), "TSYS" ) ; 4087 ArrayColumn<Float> tsysColref ; 4088 tsysColref.attach( ssig->table(), "TSYS" ) ; 4089 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 4090 vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ; 4091 ssig->setSpectrum( sp, i ) ; 4092 string reftime = ssig->getTime( i ) ; 4093 vector<int> ii( 1, ssig->getIF( i ) ) ; 4094 vector<int> ib( 1, ssig->getBeam( i ) ) ; 4095 vector<int> ip( 1, ssig->getPol( i ) ) ; 4096 sel.setIFs( ii ) ; 4097 sel.setBeams( ib ) ; 4098 sel.setPolarizations( ip ) ; 4099 asky->setSelection( sel ) ; 4100 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ; 4101 const Vector<Float> Vtsys( sptsys ) ; 4102 tsysColsig.put( i, Vtsys ) ; 4103 asky->unsetSelection() ; 4104 sel.reset() ; 4105 sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ; 4106 sref->setSpectrum( sp, i ) ; 4107 tsysColref.put( i, Vtsys ) ; 4108 } 4109 } 4110 4111 // do folding if necessary 4112 Table sigtab = ssig->table() ; 4113 Table reftab = sref->table() ; 4114 ScalarColumn<uInt> sigifnoCol ; 4115 ScalarColumn<uInt> refifnoCol ; 4116 ScalarColumn<uInt> sigfidCol ; 4117 ScalarColumn<uInt> reffidCol ; 4118 Int nchan = (Int)ssig->nchan() ; 4119 sigifnoCol.attach( sigtab, "IFNO" ) ; 4120 refifnoCol.attach( reftab, "IFNO" ) ; 4121 sigfidCol.attach( sigtab, "FREQ_ID" ) ; 4122 reffidCol.attach( reftab, "FREQ_ID" ) ; 4123 Vector<uInt> sfids( sigfidCol.getColumn() ) ; 4124 Vector<uInt> rfids( reffidCol.getColumn() ) ; 4125 vector<uInt> sfids_unique ; 4126 vector<uInt> rfids_unique ; 4127 vector<uInt> sifno_unique ; 4128 vector<uInt> rifno_unique ; 4129 for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) { 4130 if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) { 4131 sfids_unique.push_back( sfids[i] ) ; 4132 sifno_unique.push_back( ssig->getIF( i ) ) ; 4133 } 4134 if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) { 4135 rfids_unique.push_back( rfids[i] ) ; 4136 rifno_unique.push_back( sref->getIF( i ) ) ; 4137 } 4138 } 4139 double refpix_sig, refval_sig, increment_sig ; 4140 double refpix_ref, refval_ref, increment_ref ; 4141 vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ; 4142 for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) { 4143 ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ; 4144 sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ; 4145 if ( refpix_sig == refpix_ref ) { 4146 double foffset = refval_ref - refval_sig ; 4147 int choffset = static_cast<int>(foffset/increment_sig) ; 4148 double doffset = foffset / increment_sig ; 4149 if ( abs(choffset) >= nchan ) { 4150 LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ; 4151 os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ; 4152 os << "Just return signal data" << LogIO::POST ; 4153 //std::vector< CountedPtr<Scantable> > tabs ; 4154 //tabs.push_back( ssig ) ; 4155 //tabs.push_back( sref ) ; 4156 //out = merge( tabs ) ; 4157 tmp[i] = ssig ; 4158 } 4159 else { 4160 STSelector sel = STSelector() ; 4161 vector<int> v( 1, sifno_unique[i] ) ; 4162 sel.setIFs( v ) ; 4163 ssig->setSelection( sel ) ; 4164 sel.reset() ; 4165 v[0] = rifno_unique[i] ; 4166 sel.setIFs( v ) ; 4167 sref->setSelection( sel ) ; 4168 sel.reset() ; 4169 if ( antname.find( "APEX" ) != string::npos ) { 4170 tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ; 4171 //tmp[i] = dofold( ssig, sref, doffset ) ; 4172 } 4173 else { 4174 tmp[i] = dofold( ssig, sref, doffset ) ; 4175 } 4176 ssig->unsetSelection() ; 4177 sref->unsetSelection() ; 4178 } 4179 } 4180 } 4181 4182 if ( tmp.size() > 1 ) { 4183 out = merge( tmp ) ; 4184 } 4185 else { 4186 out = tmp[0] ; 4187 } 4188 4189 // flux unit 4190 out->setFluxUnit( "K" ) ; 4191 4192 return out ; 4193 } 4194 4195 CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s ) 4196 { 4197 CountedPtr<Scantable> out ; 4198 4199 return out ; 4200 } 4201 4202 vector<float> STMath::getSpectrumFromTime( string reftime, 4203 CountedPtr<Scantable>& s, 4204 string mode ) 4205 { 4206 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ; 4207 vector<float> sp ; 4208 4209 if ( s->nrow() == 0 ) { 4210 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ; 4211 return sp ; 4212 } 4213 else if ( s->nrow() == 1 ) { 4214 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ; 4215 return s->getSpectrum( 0 ) ; 4216 } 4217 else { 4218 vector<int> idx = getRowIdFromTime( reftime, s ) ; 4219 if ( mode == "before" ) { 4220 int id = -1 ; 4221 if ( idx[0] != -1 ) { 4222 id = idx[0] ; 4223 } 4224 else if ( idx[1] != -1 ) { 4225 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 4226 id = idx[1] ; 4227 } 4228 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4229 sp = s->getSpectrum( id ) ; 4230 } 4231 else if ( mode == "after" ) { 4232 int id = -1 ; 4233 if ( idx[1] != -1 ) { 4234 id = idx[1] ; 4235 } 4236 else if ( idx[0] != -1 ) { 4237 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 4238 id = idx[1] ; 4239 } 4240 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4241 sp = s->getSpectrum( id ) ; 4242 } 4243 else if ( mode == "nearest" ) { 4244 int id = -1 ; 4245 if ( idx[0] == -1 ) { 4246 id = idx[1] ; 4247 } 4248 else if ( idx[1] == -1 ) { 4249 id = idx[0] ; 4250 } 4251 else if ( idx[0] == idx[1] ) { 4252 id = idx[0] ; 4253 } 4254 else { 4255 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4256 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4257 double tref = getMJD( reftime ) ; 4258 if ( abs( t0 - tref ) > abs( t1 - tref ) ) { 4259 id = idx[1] ; 4260 } 4261 else { 4262 id = idx[0] ; 4263 } 4264 } 4265 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4266 sp = s->getSpectrum( id ) ; 4267 } 4268 else if ( mode == "linear" ) { 4269 if ( idx[0] == -1 ) { 4270 // use after 4271 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4272 int id = idx[1] ; 4273 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4274 sp = s->getSpectrum( id ) ; 4275 } 4276 else if ( idx[1] == -1 ) { 4277 // use before 4278 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4279 int id = idx[0] ; 4280 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4281 sp = s->getSpectrum( id ) ; 4282 } 4283 else if ( idx[0] == idx[1] ) { 4284 // use before 4285 //os << "No need to interporate." << LogIO::POST ; 4286 int id = idx[0] ; 4287 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4288 sp = s->getSpectrum( id ) ; 4289 } 4290 else { 4291 // do interpolation 4292 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 4293 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4294 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4295 double tref = getMJD( reftime ) ; 4296 vector<float> sp0 = s->getSpectrum( idx[0] ) ; 4297 vector<float> sp1 = s->getSpectrum( idx[1] ) ; 4298 for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) { 4299 float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ; 4300 sp.push_back( v ) ; 4301 } 4302 } 4303 } 4304 else { 4305 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 4306 } 4307 return sp ; 4308 } 4309 } 4310 4311 double STMath::getMJD( string strtime ) 4312 { 4313 if ( strtime.find("/") == string::npos ) { 4314 // MJD time string 4315 return atof( strtime.c_str() ) ; 4316 } 4317 else { 4318 // string in YYYY/MM/DD/HH:MM:SS format 4319 uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ; 4320 uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ; 4321 uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ; 4322 uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ; 4323 uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ; 4324 uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ; 4325 Time t( year, month, day, hour, minute, sec ) ; 4326 return t.modifiedJulianDay() ; 4327 } 4328 } 4329 4330 vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s ) 4331 { 4332 double reft = getMJD( reftime ) ; 4333 double dtmin = 1.0e100 ; 4334 double dtmax = -1.0e100 ; 4335 vector<double> dt ; 4336 int just_before = -1 ; 4337 int just_after = -1 ; 4338 for ( int i = 0 ; i < s->nrow() ; i++ ) { 4339 dt.push_back( getMJD( s->getTime( i ) ) - reft ) ; 4340 } 4341 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) { 4342 if ( dt[i] > 0.0 ) { 4343 // after reftime 4344 if ( dt[i] < dtmin ) { 4345 just_after = i ; 4346 dtmin = dt[i] ; 4347 } 4348 } 4349 else if ( dt[i] < 0.0 ) { 4350 // before reftime 4351 if ( dt[i] > dtmax ) { 4352 just_before = i ; 4353 dtmax = dt[i] ; 4354 } 4355 } 4356 else { 4357 // just a reftime 4358 just_before = i ; 4359 just_after = i ; 4360 dtmax = 0 ; 4361 dtmin = 0 ; 4362 break ; 4363 } 4364 } 4365 4366 vector<int> v ; 4367 v.push_back( just_before ) ; 4368 v.push_back( just_after ) ; 4369 4370 return v ; 4371 } 4372 4373 vector<float> STMath::getTcalFromTime( string reftime, 4374 CountedPtr<Scantable>& s, 4375 string mode ) 4376 { 4377 LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ; 4378 vector<float> tcal ; 4379 STTcal tcalTable = s->tcal() ; 4380 String time ; 4381 Vector<Float> tcalval ; 4382 if ( s->nrow() == 0 ) { 4383 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ; 4384 return tcal ; 4385 } 4386 else if ( s->nrow() == 1 ) { 4387 uInt tcalid = s->getTcalId( 0 ) ; 4388 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4389 tcalTable.getEntry( time, tcalval, tcalid ) ; 4390 tcalval.tovector( tcal ) ; 4391 return tcal ; 4392 } 4393 else { 4394 vector<int> idx = getRowIdFromTime( reftime, s ) ; 4395 if ( mode == "before" ) { 4396 int id = -1 ; 4397 if ( idx[0] != -1 ) { 4398 id = idx[0] ; 4399 } 4400 else if ( idx[1] != -1 ) { 4401 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 4402 id = idx[1] ; 4403 } 4404 uInt tcalid = s->getTcalId( id ) ; 4405 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4406 tcalTable.getEntry( time, tcalval, tcalid ) ; 4407 tcalval.tovector( tcal ) ; 4408 } 4409 else if ( mode == "after" ) { 4410 int id = -1 ; 4411 if ( idx[1] != -1 ) { 4412 id = idx[1] ; 4413 } 4414 else if ( idx[0] != -1 ) { 4415 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 4416 id = idx[1] ; 4417 } 4418 uInt tcalid = s->getTcalId( id ) ; 4419 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4420 tcalTable.getEntry( time, tcalval, tcalid ) ; 4421 tcalval.tovector( tcal ) ; 4422 } 4423 else if ( mode == "nearest" ) { 4424 int id = -1 ; 4425 if ( idx[0] == -1 ) { 4426 id = idx[1] ; 4427 } 4428 else if ( idx[1] == -1 ) { 4429 id = idx[0] ; 4430 } 4431 else if ( idx[0] == idx[1] ) { 4432 id = idx[0] ; 4433 } 4434 else { 4435 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4436 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4437 double tref = getMJD( reftime ) ; 4438 if ( abs( t0 - tref ) > abs( t1 - tref ) ) { 4439 id = idx[1] ; 4440 } 4441 else { 4442 id = idx[0] ; 4443 } 4444 } 4445 uInt tcalid = s->getTcalId( id ) ; 4446 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4447 tcalTable.getEntry( time, tcalval, tcalid ) ; 4448 tcalval.tovector( tcal ) ; 4449 } 4450 else if ( mode == "linear" ) { 4451 if ( idx[0] == -1 ) { 4452 // use after 4453 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4454 int id = idx[1] ; 4455 uInt tcalid = s->getTcalId( id ) ; 4456 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4457 tcalTable.getEntry( time, tcalval, tcalid ) ; 4458 tcalval.tovector( tcal ) ; 4459 } 4460 else if ( idx[1] == -1 ) { 4461 // use before 4462 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4463 int id = idx[0] ; 4464 uInt tcalid = s->getTcalId( id ) ; 4465 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4466 tcalTable.getEntry( time, tcalval, tcalid ) ; 4467 tcalval.tovector( tcal ) ; 4468 } 4469 else if ( idx[0] == idx[1] ) { 4470 // use before 4471 //os << "No need to interporate." << LogIO::POST ; 4472 int id = idx[0] ; 4473 uInt tcalid = s->getTcalId( id ) ; 4474 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4475 tcalTable.getEntry( time, tcalval, tcalid ) ; 4476 tcalval.tovector( tcal ) ; 4477 } 4478 else { 4479 // do interpolation 4480 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 4481 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4482 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4483 double tref = getMJD( reftime ) ; 4484 vector<float> tcal0 ; 4485 vector<float> tcal1 ; 4486 uInt tcalid0 = s->getTcalId( idx[0] ) ; 4487 uInt tcalid1 = s->getTcalId( idx[1] ) ; 4488 tcalTable.getEntry( time, tcalval, tcalid0 ) ; 4489 tcalval.tovector( tcal0 ) ; 4490 tcalTable.getEntry( time, tcalval, tcalid1 ) ; 4491 tcalval.tovector( tcal1 ) ; 4492 for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) { 4493 float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ; 4494 tcal.push_back( v ) ; 4495 } 4496 } 4497 } 4498 else { 4499 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 4500 } 4501 return tcal ; 4502 } 4503 } 4504 4505 vector<float> STMath::getTsysFromTime( string reftime, 4506 CountedPtr<Scantable>& s, 4507 string mode ) 4508 { 4509 LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ; 4510 ArrayColumn<Float> tsysCol ; 4511 tsysCol.attach( s->table(), "TSYS" ) ; 4512 vector<float> tsys ; 4513 String time ; 4514 Vector<Float> tsysval ; 4515 if ( s->nrow() == 0 ) { 4516 os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ; 4517 return tsys ; 4518 } 4519 else if ( s->nrow() == 1 ) { 4520 //os << "use row " << 0 << LogIO::POST ; 4521 tsysval = tsysCol( 0 ) ; 4522 tsysval.tovector( tsys ) ; 4523 return tsys ; 4524 } 4525 else { 4526 vector<int> idx = getRowIdFromTime( reftime, s ) ; 4527 if ( mode == "before" ) { 4528 int id = -1 ; 4529 if ( idx[0] != -1 ) { 4530 id = idx[0] ; 4531 } 4532 else if ( idx[1] != -1 ) { 4533 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 4534 id = idx[1] ; 4535 } 4536 //os << "use row " << id << LogIO::POST ; 4537 tsysval = tsysCol( id ) ; 4538 tsysval.tovector( tsys ) ; 4539 } 4540 else if ( mode == "after" ) { 4541 int id = -1 ; 4542 if ( idx[1] != -1 ) { 4543 id = idx[1] ; 4544 } 4545 else if ( idx[0] != -1 ) { 4546 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 4547 id = idx[1] ; 4548 } 4549 //os << "use row " << id << LogIO::POST ; 4550 tsysval = tsysCol( id ) ; 4551 tsysval.tovector( tsys ) ; 4552 } 4553 else if ( mode == "nearest" ) { 4554 int id = -1 ; 4555 if ( idx[0] == -1 ) { 4556 id = idx[1] ; 4557 } 4558 else if ( idx[1] == -1 ) { 4559 id = idx[0] ; 4560 } 4561 else if ( idx[0] == idx[1] ) { 4562 id = idx[0] ; 4563 } 4564 else { 4565 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4566 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4567 double tref = getMJD( reftime ) ; 4568 if ( abs( t0 - tref ) > abs( t1 - tref ) ) { 4569 id = idx[1] ; 4570 } 4571 else { 4572 id = idx[0] ; 4573 } 4574 } 4575 //os << "use row " << id << LogIO::POST ; 4576 tsysval = tsysCol( id ) ; 4577 tsysval.tovector( tsys ) ; 4578 } 4579 else if ( mode == "linear" ) { 4580 if ( idx[0] == -1 ) { 4581 // use after 4582 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4583 int id = idx[1] ; 4584 //os << "use row " << id << LogIO::POST ; 4585 tsysval = tsysCol( id ) ; 4586 tsysval.tovector( tsys ) ; 4587 } 4588 else if ( idx[1] == -1 ) { 4589 // use before 4590 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4591 int id = idx[0] ; 4592 //os << "use row " << id << LogIO::POST ; 4593 tsysval = tsysCol( id ) ; 4594 tsysval.tovector( tsys ) ; 4595 } 4596 else if ( idx[0] == idx[1] ) { 4597 // use before 4598 //os << "No need to interporate." << LogIO::POST ; 4599 int id = idx[0] ; 4600 //os << "use row " << id << LogIO::POST ; 4601 tsysval = tsysCol( id ) ; 4602 tsysval.tovector( tsys ) ; 4603 } 4604 else { 4605 // do interpolation 4606 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 4607 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4608 double t1 = getMJD( s->getTime( idx[1] ) ) ; 4609 double tref = getMJD( reftime ) ; 4610 vector<float> tsys0 ; 4611 vector<float> tsys1 ; 4612 tsysval = tsysCol( idx[0] ) ; 4613 tsysval.tovector( tsys0 ) ; 4614 tsysval = tsysCol( idx[1] ) ; 4615 tsysval.tovector( tsys1 ) ; 4616 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) { 4617 float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ; 4618 tsys.push_back( v ) ; 4619 } 4620 } 4621 } 4622 else { 4623 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 4624 } 4625 return tsys ; 4626 } 4627 } 4628 4629 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on, 4630 CountedPtr<Scantable>& off, 4631 CountedPtr<Scantable>& sky, 4632 CountedPtr<Scantable>& hot, 4633 CountedPtr<Scantable>& cold, 4634 int index, 4635 string antname ) 4636 { 4637 string reftime = on->getTime( index ) ; 4638 vector<int> ii( 1, on->getIF( index ) ) ; 4639 vector<int> ib( 1, on->getBeam( index ) ) ; 4640 vector<int> ip( 1, on->getPol( index ) ) ; 4641 vector<int> ic( 1, on->getScan( index ) ) ; 4642 STSelector sel = STSelector() ; 4643 sel.setIFs( ii ) ; 4644 sel.setBeams( ib ) ; 4645 sel.setPolarizations( ip ) ; 4646 sky->setSelection( sel ) ; 4647 hot->setSelection( sel ) ; 4648 //cold->setSelection( sel ) ; 4649 off->setSelection( sel ) ; 4650 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ; 4651 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ; 4652 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ; 4653 vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ; 4654 vector<float> spec = on->getSpectrum( index ) ; 4655 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ; 4656 vector<float> sp( tcal.size() ) ; 4657 if ( antname.find( "APEX" ) != string::npos ) { 4658 // using gain array 4659 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 4660 float v = ( ( spec[j] - spoff[j] ) / spoff[j] ) 4661 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ; 4662 sp[j] = v ; 4663 } 4664 } 4665 else { 4666 // Chopper-Wheel calibration (Ulich & Haas 1976) 4667 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 4668 float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ; 4669 sp[j] = v ; 4670 } 4671 } 4672 sel.reset() ; 4673 sky->unsetSelection() ; 4674 hot->unsetSelection() ; 4675 //cold->unsetSelection() ; 4676 off->unsetSelection() ; 4677 4678 return sp ; 4679 } 4680 4681 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on, 4682 CountedPtr<Scantable>& off, 4683 int index ) 4684 { 4685 string reftime = on->getTime( index ) ; 4686 vector<int> ii( 1, on->getIF( index ) ) ; 4687 vector<int> ib( 1, on->getBeam( index ) ) ; 4688 vector<int> ip( 1, on->getPol( index ) ) ; 4689 vector<int> ic( 1, on->getScan( index ) ) ; 4690 STSelector sel = STSelector() ; 4691 sel.setIFs( ii ) ; 4692 sel.setBeams( ib ) ; 4693 sel.setPolarizations( ip ) ; 4694 off->setSelection( sel ) ; 4695 vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ; 4696 vector<float> spec = on->getSpectrum( index ) ; 4697 //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ; 4698 //vector<float> tsys = on->getTsysVec( index ) ; 4699 ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ; 4700 Vector<Float> tsys = tsysCol( index ) ; 4701 vector<float> sp( spec.size() ) ; 4702 // ALMA Calibration 4703 // 4704 // Ta* = Tsys * ( ON - OFF ) / OFF 4705 // 4706 // 2010/01/07 Takeshi Nakazato 4707 unsigned int tsyssize = tsys.nelements() ; 4708 unsigned int spsize = sp.size() ; 4709 for ( unsigned int j = 0 ; j < sp.size() ; j++ ) { 4710 float tscale = 0.0 ; 4711 if ( tsyssize == spsize ) 4712 tscale = tsys[j] ; 4713 else 4714 tscale = tsys[0] ; 4715 float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ; 4716 sp[j] = v ; 4717 } 4718 sel.reset() ; 4719 off->unsetSelection() ; 4720 4721 return sp ; 4722 } 4723 4724 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig, 4725 CountedPtr<Scantable>& ref, 4726 CountedPtr<Scantable>& sky, 4727 CountedPtr<Scantable>& hot, 4728 CountedPtr<Scantable>& cold, 4729 int index ) 4730 { 4731 string reftime = sig->getTime( index ) ; 4732 vector<int> ii( 1, sig->getIF( index ) ) ; 4733 vector<int> ib( 1, sig->getBeam( index ) ) ; 4734 vector<int> ip( 1, sig->getPol( index ) ) ; 4735 vector<int> ic( 1, sig->getScan( index ) ) ; 4736 STSelector sel = STSelector() ; 4737 sel.setIFs( ii ) ; 4738 sel.setBeams( ib ) ; 4739 sel.setPolarizations( ip ) ; 4740 sky->setSelection( sel ) ; 4741 hot->setSelection( sel ) ; 4742 //cold->setSelection( sel ) ; 4743 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ; 4744 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ; 4745 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ; 4746 vector<float> spref = ref->getSpectrum( index ) ; 4747 vector<float> spsig = sig->getSpectrum( index ) ; 4748 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ; 4749 vector<float> sp( tcal.size() ) ; 4750 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 4751 float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ; 4752 sp[j] = v ; 4753 } 4754 sel.reset() ; 4755 sky->unsetSelection() ; 4756 hot->unsetSelection() ; 4757 //cold->unsetSelection() ; 4758 4759 return sp ; 4760 } 4761 4762 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig, 4763 CountedPtr<Scantable>& ref, 4764 vector< CountedPtr<Scantable> >& sky, 4765 vector< CountedPtr<Scantable> >& hot, 4766 vector< CountedPtr<Scantable> >& cold, 4767 int index ) 4768 { 4769 string reftime = sig->getTime( index ) ; 4770 vector<int> ii( 1, sig->getIF( index ) ) ; 4771 vector<int> ib( 1, sig->getBeam( index ) ) ; 4772 vector<int> ip( 1, sig->getPol( index ) ) ; 4773 vector<int> ic( 1, sig->getScan( index ) ) ; 4774 STSelector sel = STSelector() ; 4775 sel.setIFs( ii ) ; 4776 sel.setBeams( ib ) ; 4777 sel.setPolarizations( ip ) ; 4778 sky[0]->setSelection( sel ) ; 4779 hot[0]->setSelection( sel ) ; 4780 //cold[0]->setSelection( sel ) ; 4781 vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ; 4782 vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ; 4783 //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ; 4784 vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ; 4785 sel.reset() ; 4786 ii[0] = ref->getIF( index ) ; 4787 sel.setIFs( ii ) ; 4788 sel.setBeams( ib ) ; 4789 sel.setPolarizations( ip ) ; 4790 sky[1]->setSelection( sel ) ; 4791 hot[1]->setSelection( sel ) ; 4792 //cold[1]->setSelection( sel ) ; 4793 vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ; 4794 vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ; 4795 //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ; 4796 vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 4797 vector<float> spref = ref->getSpectrum( index ) ; 4798 vector<float> spsig = sig->getSpectrum( index ) ; 4799 vector<float> sp( tcals.size() ) ; 4800 for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) { 4801 float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ; 4802 sp[j] = v ; 4803 } 4804 sel.reset() ; 4805 sky[0]->unsetSelection() ; 4806 hot[0]->unsetSelection() ; 4807 //cold[0]->unsetSelection() ; 4808 sky[1]->unsetSelection() ; 4809 hot[1]->unsetSelection() ; 4810 //cold[1]->unsetSelection() ; 4811 4812 return sp ; 4813 } -
branches/mergetest/src/STMath.h
r1689 r1789 119 119 const std::string& mode, bool tsys=false ); 120 120 121 // array operation 122 casa::CountedPtr<Scantable> 123 arrayOperate( const casa::CountedPtr<Scantable>& in, 124 const std::vector<float> val, 125 const std::string& mode, 126 const std::string& opmode="channel", 127 bool tsys=false ); 128 129 // channel operation 130 casa::CountedPtr<Scantable> 131 arrayOperateChannel( const casa::CountedPtr<Scantable>& in, 132 const std::vector<float> val, 133 const std::string& mode, bool tsys=false ); 134 135 // row operation 136 casa::CountedPtr<Scantable> 137 arrayOperateRow( const casa::CountedPtr<Scantable>& in, 138 const std::vector<float> val, 139 const std::string& mode, bool tsys=false ); 140 141 // 2d array operation 142 casa::CountedPtr<Scantable> 143 array2dOperate( const casa::CountedPtr<Scantable>& in, 144 const std::vector< std::vector<float> > val, 145 const std::string& mode, bool tsys=false ); 146 121 147 casa::CountedPtr<Scantable> 122 148 binaryOperate( const casa::CountedPtr<Scantable>& left, … … 134 160 /** 135 161 * Calibrate total power scans (translated from GBTIDL) 136 * @param calon uncalibrated Scantable with CAL noise signal 162 * @param calon uncalibrated Scantable with CAL noise signal 137 163 * @param caloff uncalibrated Scantable with no CAL signal 138 164 * @param tcal optional scalar Tcal, CAL temperature (K) 139 * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 165 * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 140 166 * (spectrum - average of the two CAL on and off spectra; 141 167 * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2) 142 */ 168 */ 143 169 casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon, 144 170 const casa::CountedPtr<Scantable>& caloff, … … 151 177 * @param smoothref optional Boxcar smooth width of the reference scans 152 178 * default: no smoothing (=1) 153 * @param tsysv optional scalar Tsys value at the zenith, required to 154 * set tau, as well 179 * @param tsysv optional scalar Tsys value at the zenith, required to 180 * set tau, as well 155 181 * @param tau optional scalar Tau value 156 182 * @return casa::CountedPtr<Scantable> which holds combined scans … … 163 189 casa::Float tau=0.0 ); 164 190 165 /**191 /** 166 192 * Calibrate GBT Nod scan pairs (translated from GBTIDL) 167 193 * @param s Scantable which contains Nod scans … … 170 196 * @param tsysv optional scalar Tsys value at the zenith, required to 171 197 * set tau, as well 172 * @param tau optional scalar Tau value 198 * @param tau optional scalar Tau value 173 199 * @param tcal optional scalar Tcal, CAL temperature (K) 174 200 * @return casa::CountedPtr<Scantable> which holds calibrated scans … … 199 225 casa::Float tcal=0.0 ); 200 226 227 /** 228 * Calibrate data with Chopper-Wheel like calibration method 229 * which adopts position switching by antenna motion, 230 * wobbler (nutator) switching and On-The-Fly observation. 231 * 232 * The method is applicable to APEX, and other telescopes other than GBT. 233 * 234 * @param a Scantable which contains ON and OFF scans 235 * @param a string that indicates calibration mode 236 * @param a string that indicates antenna name 237 **/ 238 casa::CountedPtr<Scantable> cwcal( const casa::CountedPtr<Scantable>& s, 239 const casa::String calmode, 240 const casa::String antname ); 241 242 /** 243 * Calibrate frequency switched scans with Chopper-Wheel like 244 * calibration method. 245 * 246 * The method is applicable to APEX, and other telescopes other than GBT. 247 * 248 * @param a Scantable which contains ON and OFF scans 249 * @param a string that indicates antenna name 250 **/ 251 casa::CountedPtr<Scantable> cwcalfs( const casa::CountedPtr<Scantable>& s, 252 const casa::String antname ); 253 254 255 /** 256 * Folding frequency-switch data 257 * @param sig 258 * @param ref 259 * @param choffset 260 **/ 261 casa::CountedPtr<Scantable> dofold( const casa::CountedPtr<Scantable> &sig, 262 const casa::CountedPtr<Scantable> &ref, 263 casa::Double choffset, 264 casa::Double choffset = 0.0 ); 265 266 /** 267 * ALMA calibration 268 **/ 269 casa::CountedPtr<Scantable> almacal( const casa::CountedPtr<Scantable>& s, 270 const casa::String calmode ) ; 271 casa::CountedPtr<Scantable> almacalfs( const casa::CountedPtr<Scantable>& s ) ; 201 272 202 273 casa::CountedPtr<Scantable> … … 206 277 const std::vector<bool>& mask, 207 278 const std::string& which); 279 280 std::vector< int > minMaxChan(const casa::CountedPtr<Scantable>& in, 281 const std::vector<bool>& mask, 282 const std::string& which); 208 283 209 284 casa::CountedPtr<Scantable> bin( const casa::CountedPtr<Scantable>& in, … … 266 341 double end, const std::string& mode="frequency"); 267 342 343 // test for average spectra with different channel/resolution 344 casa::CountedPtr<Scantable> 345 new_average( const std::vector<casa::CountedPtr<Scantable> >& in, 346 const bool& compel, 347 const std::vector<bool>& mask = std::vector<bool>(), 348 const std::string& weight = "NONE", 349 const std::string& avmode = "SCAN" ) 350 throw (casa::AipsError) ; 351 268 352 private: 269 353 casa::CountedPtr<Scantable> applyToPol( const casa::CountedPtr<Scantable>& in, … … 301 385 maskedArray( const casa::Vector<casa::Float>& s, 302 386 const casa::Vector<casa::uChar>& f ); 387 casa::MaskedArray<casa::Double> 388 maskedArray( const casa::Vector<casa::Double>& s, 389 const casa::Vector<casa::uChar>& f ); 303 390 casa::Vector<casa::uChar> 304 391 flagsFromMA(const casa::MaskedArray<casa::Float>& ma); 305 392 393 vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 394 vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 395 vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 396 vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ; 397 398 // Chopper-Wheel type calibration 399 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 400 casa::CountedPtr<Scantable>& off, 401 casa::CountedPtr<Scantable>& sky, 402 casa::CountedPtr<Scantable>& hot, 403 casa::CountedPtr<Scantable>& cold, 404 int index, 405 string antname ) ; 406 // Tsys * (ON-OFF)/OFF 407 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 408 casa::CountedPtr<Scantable>& off, 409 int index ) ; 410 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 411 casa::CountedPtr<Scantable>& ref, 412 casa::CountedPtr<Scantable>& sky, 413 casa::CountedPtr<Scantable>& hot, 414 casa::CountedPtr<Scantable>& cold, 415 int index ) ; 416 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 417 casa::CountedPtr<Scantable>& ref, 418 vector< casa::CountedPtr<Scantable> >& sky, 419 vector< casa::CountedPtr<Scantable> >& hot, 420 vector< casa::CountedPtr<Scantable> >& cold, 421 int index ) ; 422 double getMJD( string strtime ) ; 423 306 424 bool insitu_; 307 425 }; -
branches/mergetest/src/STMathWrapper.h
r1689 r1789 73 73 { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); } 74 74 75 ScantableWrapper arrayOperate( const ScantableWrapper& in, 76 const std::vector<float> val, 77 const std::string& mode, 78 bool tsys=false ) 79 { return ScantableWrapper(STMath::arrayOperateChannel(in.getCP(), val, mode, tsys)); } 80 81 ScantableWrapper array2dOperate( const ScantableWrapper& in, 82 const std::vector< std::vector<float> > val, 83 const std::string& mode, bool tsys=false ) 84 { return ScantableWrapper(STMath::array2dOperate(in.getCP(), val, mode, tsys)); } 85 75 86 ScantableWrapper binaryOperate( const ScantableWrapper& left, 76 87 const ScantableWrapper& right, … … 121 132 { return STMath::statistic(in.getCP(), mask, which); } 122 133 134 std::vector<int> minMaxChan(const ScantableWrapper& in, 135 const std::vector<bool>& mask, 136 const std::string& which) 137 { return STMath::minMaxChan(in.getCP(), mask, which); } 138 123 139 ScantableWrapper bin( const ScantableWrapper& in, int width=5) 124 140 { return ScantableWrapper(STMath::bin(in.getCP(), width)); } … … 175 191 const std::string& refTime, 176 192 const std::string& method ) 177 { return ScantableWrapper(STMath::frequencyAlign(in.getCP() , refTime, method)); }193 { return ScantableWrapper(STMath::frequencyAlign(in.getCP())); } 178 194 179 195 ScantableWrapper convertPolarisation( const ScantableWrapper& in, … … 191 207 mode)); } 192 208 209 // test for average spectra with different channel/resolution 210 ScantableWrapper 211 new_average( const std::vector<ScantableWrapper>& in, 212 const bool& compel, 213 const std::vector<bool>& mask, 214 const std::string& weight, 215 const std::string& avmode ) 216 { 217 std::vector<casa::CountedPtr<Scantable> > sts; 218 for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP()); 219 return ScantableWrapper(STMath::new_average(sts, compel, mask, weight, avmode)); 220 } 221 222 // cwcal 223 ScantableWrapper cwcal( const ScantableWrapper &in, 224 const std::string calmode, 225 const std::string antname ) 226 { 227 casa::CountedPtr<Scantable> tab = in.getCP() ; 228 casa::String mode( calmode ) ; 229 casa::String name( antname ) ; 230 return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ; 231 } 232 // almacal 233 ScantableWrapper almacal( const ScantableWrapper &in, 234 const std::string calmode ) 235 { 236 casa::CountedPtr<Scantable> tab = in.getCP() ; 237 casa::String mode( calmode ) ; 238 return ScantableWrapper( STMath::almacal( tab, mode ) ) ; 239 } 193 240 }; 194 241 -
branches/mergetest/src/STMolecules.cpp
r870 r1789 14 14 #include <tables/Tables/SetupNewTab.h> 15 15 #include <tables/Tables/ScaColDesc.h> 16 #include <tables/Tables/ArrColDesc.h> 16 17 #include <tables/Tables/TableRecord.h> 17 18 #include <tables/Tables/TableParse.h> 18 19 #include <tables/Tables/TableRow.h> 19 20 #include <casa/Containers/RecordField.h> 21 22 #include <tables/Tables/TableProxy.h> 20 23 21 24 #include "STMolecules.h" … … 59 62 { 60 63 // add to base class table 61 table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY")); 62 table_.addColumn(ScalarColumnDesc<String>("NAME")); 63 table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME")); 64 //table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY")); 65 table_.addColumn(ArrayColumnDesc<Double>("RESTFREQUENCY")); 66 //table_.addColumn(ScalarColumnDesc<String>("NAME")); 67 table_.addColumn(ArrayColumnDesc<String>("NAME")); 68 //table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME")); 69 table_.addColumn(ArrayColumnDesc<String>("FORMATTEDNAME")); 64 70 table_.rwKeywordSet().define("UNIT", String("Hz")); 65 71 // new cached columns … … 69 75 } 70 76 77 /*** 71 78 uInt STMolecules::addEntry( Double restfreq, const String& name, 72 79 const String& formattedname ) … … 94 101 return resultid; 95 102 } 96 103 ***/ 104 uInt STMolecules::addEntry( Vector<Double> restfreq, const Vector<String>& name, 105 const Vector<String>& formattedname ) 106 { 107 // How to handle this...? 108 Table result = 109 table_( nelements(table_.col("RESTFREQUENCY")) == uInt (restfreq.size()) && 110 all(table_.col("RESTFREQUENCY")== restfreq) ); 111 uInt resultid = 0; 112 if ( result.nrow() > 0) { 113 ROScalarColumn<uInt> c(result, "ID"); 114 c.get(0, resultid); 115 } else { 116 uInt rno = table_.nrow(); 117 table_.addRow(); 118 // get last assigned _id and increment 119 if ( rno > 0 ) { 120 idCol_.get(rno-1, resultid); 121 resultid++; 122 } 123 restfreqCol_.put(rno, restfreq); 124 nameCol_.put(rno, name); 125 formattednameCol_.put(rno, formattedname); 126 idCol_.put(rno, resultid); 127 } 128 return resultid; 129 } 130 131 132 133 134 /*** 97 135 void STMolecules::getEntry( Double& restfreq, String& name, 98 136 String& formattedname, uInt id ) const … … 104 142 ROTableRow row(t); 105 143 // get first row - there should only be one matching id 144 106 145 const TableRecord& rec = row.get(0); 107 146 restfreq = rec.asDouble("RESTFREQUENCY"); … … 109 148 formattedname = rec.asString("FORMATTEDNAME"); 110 149 } 150 ***/ 151 void STMolecules::getEntry( Vector<Double>& restfreq, Vector<String>& name, 152 Vector<String>& formattedname, uInt id ) const 153 { 154 Table t = table_(table_.col("ID") == Int(id) ); 155 if (t.nrow() == 0 ) { 156 throw(AipsError("STMolecules::getEntry - id out of range")); 157 } 158 ROTableRow row(t); 159 // get first row - there should only be one matching id 160 161 const TableRecord& rec = row.get(0); 162 //restfreq = rec.asDouble("RESTFREQUENCY"); 163 restfreq = rec.asArrayDouble("RESTFREQUENCY"); 164 //name = rec.asString("NAME"); 165 name = rec.asArrayString("NAME"); 166 //formattedname = rec.asString("FORMATTEDNAME"); 167 formattedname = rec.asArrayString("FORMATTEDNAME"); 168 } 111 169 112 170 std::vector< double > asap::STMolecules::getRestFrequencies( ) const 113 171 { 114 172 std::vector<double> out; 173 //TableProxy itsTable(table_); 174 //Record rec; 115 175 Vector<Double> rfs = restfreqCol_.getColumn(); 116 176 rfs.tovector(out); 177 //rec = itsTable.getVarColumn("RESTFREQUENCY", 0, 1, 1); 117 178 return out; 118 179 } 119 180 120 double asap::STMolecules::getRestFrequency( uInt id ) const 121 { 181 std::vector< double > asap::STMolecules::getRestFrequency( uInt id ) const 182 { 183 std::vector<double> out; 122 184 Table t = table_(table_.col("ID") == Int(id) ); 123 185 if (t.nrow() == 0 ) { … … 126 188 ROTableRow row(t); 127 189 const TableRecord& rec = row.get(0); 128 return double(rec.asDouble("RESTFREQUENCY")); 129 } 130 190 //return double(rec.asDouble("RESTFREQUENCY")); 191 Vector<Double> rfs = rec.asArrayDouble("RESTFREQUENCY"); 192 rfs.tovector(out); 193 return out; 194 } 195 196 int asap::STMolecules::nrow() const 197 { 198 return int(table_.nrow()); 199 } 131 200 132 201 }//namespace -
branches/mergetest/src/STMolecules.h
r1353 r1789 17 17 #include <tables/Tables/Table.h> 18 18 #include <tables/Tables/ScalarColumn.h> 19 #include <tables/Tables/ArrayColumn.h> 20 #include <casa/Arrays/Array.h> 19 21 20 22 #include "STSubTable.h" … … 37 39 STMolecules& operator=(const STMolecules& other); 38 40 41 /*** 39 42 casa::uInt addEntry( casa::Double restfreq, const casa::String& name="", 40 43 const casa::String& formattedname=""); 44 ***/ 41 45 46 casa::uInt addEntry( casa::Vector<casa::Double> restfreq, const casa::Vector<casa::String>& name=casa::Vector<casa::String>(0), 47 const casa::Vector<casa::String>& formattedname=casa::Vector<casa::String>(0)); 48 49 /*** 42 50 void getEntry( casa::Double& restfreq, casa::String& name, 43 51 casa::String& formattedname, casa::uInt id) const; 52 ***/ 53 void getEntry( casa::Vector<casa::Double>& restfreq, casa::Vector<casa::String>& name, 54 casa::Vector<casa::String>& formattedname, casa::uInt id) const; 44 55 45 56 std::vector<double> getRestFrequencies() const; 46 doublegetRestFrequency( casa::uInt id ) const;57 std::vector<double> getRestFrequency( casa::uInt id ) const; 47 58 const casa::String& name() const { return name_; } 59 int nrow() const; 48 60 49 61 private: … … 52 64 //casa::Table table_; 53 65 //casa::ScalarColumn<casa::uInt> freqidCol_; 54 casa::ScalarColumn<casa::Double> restfreqCol_; 55 casa::ScalarColumn<casa::String> nameCol_; 56 casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex 66 //casa::ScalarColumn<casa::Double> restfreqCol_; 67 casa::ArrayColumn<casa::Double> restfreqCol_; 68 //casa::ScalarColumn<casa::String> nameCol_; 69 casa::ArrayColumn<casa::String> nameCol_; 70 //casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex 71 casa::ArrayColumn<casa::String> formattednameCol_; // e.g. latex 57 72 58 73 }; -
branches/mergetest/src/STSelector.cpp
r1542 r1789 87 87 } 88 88 89 void STSelector::setTypes( const std::vector< int >& types ) 90 { 91 setint("SRCTYPE", types); 92 } 93 89 94 void STSelector::setint(const std::string& key, const std::vector< int >& val) 90 95 { … … 116 121 } 117 122 123 void STSelector::setRows( const std::vector< int >& rows ) 124 { 125 rowselection_ = rows; 126 } 127 128 // Table STSelector::apply( const Table& tab ) 129 // { 130 // if ( empty() ) { 131 // return sort(tab); 132 // } 133 // TableExprNode query; 134 // intidmap::const_iterator it; 135 // for (it = intselections_.begin(); it != intselections_.end(); ++it) { 136 // TableExprNode theset(Vector<Int>( (*it).second )); 137 // if ( query.isNull() ) { 138 // query = tab.col((*it).first).in(theset); 139 // } else { 140 // query = tab.col((*it).first).in(theset) && query; 141 // } 142 // } 143 // stringidmap::const_iterator it1; 144 // for (it1 = stringselections_.begin(); it1 != stringselections_.end(); ++it1) { 145 // TableExprNode theset(mathutil::toVectorString( (*it1).second )); 146 // if ( query.isNull() ) { 147 // query = tab.col((*it1).first).in(theset); 148 // } else { 149 // query = tab.col((*it1).first).in(theset) && query; 150 // } 151 // } 152 // // add taql query 153 // if ( taql_.size() > 0 ) { 154 // Table tmpt = tab; 155 // std::string pytaql = "USING STYLE PYTHON " + taql_; 156 157 // if ( !query.isNull() ) { // taql and selection 158 // tmpt = tableCommand(pytaql, tab(query)); 159 // } else { // taql only 160 // tmpt = tableCommand(pytaql, tab); 161 // } 162 // return sort(tmpt); 163 // } else { 164 // if ( query.isNull() ) { 165 // return sort(tab); 166 // } else { 167 // return sort(tab(query)); 168 // } 169 // } 170 // } 118 171 Table STSelector::apply( const Table& tab ) 119 172 { 120 173 if ( empty() ) { 121 174 return sort(tab); 175 } 176 Table basetab = tab; 177 // Important!! Be sure to apply row selection first. 178 if (rowselection_.size() > 0){ 179 //Vector<Int> intrownrs(rowselection_); 180 Vector<uInt> rownrs( rowselection_.size() ); 181 convertArray(rownrs, Vector<Int> ( rowselection_ )); 182 basetab = tab( rownrs ); 183 ///TableExprNode theset(Vector<Int>( rowselection_ )); 184 ///query = tab.nodeRownr().in(theset); 122 185 } 123 186 TableExprNode query; … … 126 189 TableExprNode theset(Vector<Int>( (*it).second )); 127 190 if ( query.isNull() ) { 128 query = tab.col((*it).first).in(theset); 191 //query = tab.col((*it).first).in(theset); 192 query = basetab.col((*it).first).in(theset); 129 193 } else { 130 query = tab.col((*it).first).in(theset) && query; 194 //query = tab.col((*it).first).in(theset) && query; 195 query = basetab.col((*it).first).in(theset) && query; 131 196 } 132 197 } … … 135 200 TableExprNode theset(mathutil::toVectorString( (*it1).second )); 136 201 if ( query.isNull() ) { 137 query = tab.col((*it1).first).in(theset); 202 //query = tab.col((*it1).first).in(theset); 203 query = basetab.col((*it1).first).in(theset); 138 204 } else { 139 query = tab.col((*it1).first).in(theset) && query; 205 //query = tab.col((*it1).first).in(theset) && query; 206 query = basetab.col((*it1).first).in(theset) && query; 140 207 } 141 208 } 142 209 // add taql query 143 210 if ( taql_.size() > 0 ) { 144 Table tmpt = tab; 211 //Table tmpt = tab; 212 Table tmpt = basetab; 145 213 std::string pytaql = "USING STYLE PYTHON " + taql_; 146 214 147 215 if ( !query.isNull() ) { // taql and selection 148 tmpt = tableCommand(pytaql, tab(query)); 216 //tmpt = tableCommand(pytaql, tab(query)); 217 tmpt = tableCommand(pytaql, basetab(query)); 149 218 } else { // taql only 150 tmpt = tableCommand(pytaql, tab); 219 //tmpt = tableCommand(pytaql, tab); 220 tmpt = tableCommand(pytaql, basetab); 151 221 } 152 222 return sort(tmpt); 153 223 } else { 154 224 if ( query.isNull() ) { 155 return sort(tab); 225 //return sort(tab); 226 return sort(basetab); 156 227 } else { 157 return sort(tab(query)); 228 //return sort(tab(query)); 229 return sort(basetab(query)); 158 230 } 159 231 } … … 191 263 { 192 264 return getint("CYCLENO"); 265 } 266 267 std::vector< int > asap::STSelector::getTypes( ) const 268 { 269 return getint("SRCTYPE") ; 193 270 } 194 271 … … 227 304 bool asap::STSelector::empty( ) const 228 305 { 229 return (intselections_.empty() && taql_.size() == 0 ); 306 //return (intselections_.empty() && taql_.size() == 0 ); 307 return (intselections_.empty() && taql_.size() == 0 && rowselection_.size() == 0); 230 308 } 231 309 -
branches/mergetest/src/STSelector.h
r939 r1789 45 45 void setCycles(const std::vector<int>& cycs); 46 46 void setName(const std::string&); 47 void setTypes(const std::vector<int>& types); 47 48 virtual void setTaQL(const std::string& taql); 48 49 49 50 void setSortOrder(const std::vector<std::string>& order); 51 void setRows(const std::vector<int>& rows); 50 52 51 53 std::vector<int> getScans() const; … … 54 56 std::vector<int> getPols() const; 55 57 std::vector<int> getCycles() const; 58 std::vector<int> getTypes() const; 56 59 std::vector<std::string> getPolTypes() const; 57 60 std::string getTaQL() const { return taql_; } … … 86 89 casa::Block<casa::String> order_; 87 90 std::string taql_; 91 std::vector<int> rowselection_; 88 92 }; 89 93 -
branches/mergetest/src/STWriter.cpp
r1688 r1789 42 42 #include <atnf/PKSIO/PKSMS2writer.h> 43 43 #include <atnf/PKSIO/PKSSDwriter.h> 44 #include <atnf/PKSIO/SrcType.h> 44 45 45 46 #include <tables/Tables/Table.h> … … 154 155 havexpol(ifs[i]) = nPol(ifs[i]) > 2; 155 156 } 156 Vector<String> obstypes(2);157 obstypes(0) = "TR";//on158 obstypes(1) = "RF TR";//off157 // Vector<String> obstypes(2); 158 // obstypes(0) = "TR";//on 159 // obstypes(1) = "RF TR";//off 159 160 const Table table = inst->table(); 160 161 … … 182 183 while (!beamit.pastEnd() ) { 183 184 Table btable = beamit.table(); 185 //MDirection::ScalarColumn dirCol(btable, "DIRECTION"); 186 //pksrec.direction = dirCol(0).getAngle("rad").getValue(); 184 187 TableIterator cycit(btable, "CYCLENO"); 185 188 ROArrayColumn<Double> srateCol(btable, "SCANRATE"); … … 188 191 Vector<Float> srateflt(sratedbl.nelements()); 189 192 convertArray(srateflt, sratedbl); 190 pksrec.scanRate = srateflt; 193 //pksrec.scanRate = srateflt; 194 pksrec.scanRate = sratedbl; 191 195 ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION"); 192 196 spmCol.get(0, pksrec.srcPM); … … 200 204 while (!cycit.pastEnd() ) { 201 205 Table ctable = cycit.table(); 202 TableIterator ifit(ctable, "IFNO" );206 TableIterator ifit(ctable, "IFNO", TableIterator::Ascending, TableIterator::HeapSort); 203 207 MDirection::ScalarColumn dirCol(ctable, "DIRECTION"); 204 208 pksrec.direction = dirCol(0).getAngle("rad").getValue(); … … 213 217 uInt nchan = specCol(0).nelements(); 214 218 Double crval,crpix; 219 //Vector<Double> restfreq; 215 220 Float tmp0,tmp1,tmp2,tmp3,tmp4; 216 String stmp0,stmp1; 221 String tcalt; 222 Vector<String> stmp0, stmp1; 217 223 inst->frequencies().getEntry(crpix,crval, pksrec.freqInc, 218 224 rec.asuInt("FREQ_ID")); … … 239 245 pksrec.fieldName = rec.asString("FIELDNAME"); 240 246 pksrec.srcName = rec.asString("SRCNAME"); 241 pksrec.obsType = obstypes[rec.asInt("SRCTYPE")]; 247 //pksrec.obsType = obstypes[rec.asInt("SRCTYPE")]; 248 pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ; 242 249 pksrec.bandwidth = nchan * abs(pksrec.freqInc); 243 250 pksrec.azimuth = rec.asFloat("AZIMUTH"); … … 253 260 pksrec.baseSub = 0.0f; 254 261 pksrec.xCalFctr = 0.0; 262 pksrec.flagrow = rec.asuInt("FLAGROW"); 255 263 256 264 status = writer_->write(pksrec); … … 348 356 } 349 357 350 } 358 // get obsType string from SRCTYPE value 359 String STWriter::getObsTypes( Int srctype ) 360 { 361 String obsType ; 362 switch( srctype ) { 363 case Int(SrcType::PSON): 364 obsType = "PSON" ; 365 break ; 366 case Int(SrcType::PSOFF): 367 obsType = "PSOFF" ; 368 break ; 369 case Int(SrcType::NOD): 370 obsType = "NOD" ; 371 break ; 372 case Int(SrcType::FSON): 373 obsType = "FSON" ; 374 break ; 375 case Int(SrcType::FSOFF): 376 obsType = "FSOFF" ; 377 break ; 378 case Int(SrcType::SKY): 379 obsType = "SKY" ; 380 break ; 381 case Int(SrcType::HOT): 382 obsType = "HOT" ; 383 break ; 384 case Int(SrcType::WARM): 385 obsType = "WARM" ; 386 break ; 387 case Int(SrcType::COLD): 388 obsType = "COLD" ; 389 break ; 390 case Int(SrcType::PONCAL): 391 obsType = "PSON:CALON" ; 392 break ; 393 case Int(SrcType::POFFCAL): 394 obsType = "PSOFF:CALON" ; 395 break ; 396 case Int(SrcType::NODCAL): 397 obsType = "NOD:CALON" ; 398 break ; 399 case Int(SrcType::FONCAL): 400 obsType = "FSON:CALON" ; 401 break ; 402 case Int(SrcType::FOFFCAL): 403 obsType = "FSOFF:CALOFF" ; 404 break ; 405 case Int(SrcType::FSLO): 406 obsType = "FSLO" ; 407 break ; 408 case Int(SrcType::FLOOFF): 409 obsType = "FS:LOWER:OFF" ; 410 break ; 411 case Int(SrcType::FLOSKY): 412 obsType = "FS:LOWER:SKY" ; 413 break ; 414 case Int(SrcType::FLOHOT): 415 obsType = "FS:LOWER:HOT" ; 416 break ; 417 case Int(SrcType::FLOWARM): 418 obsType = "FS:LOWER:WARM" ; 419 break ; 420 case Int(SrcType::FLOCOLD): 421 obsType = "FS:LOWER:COLD" ; 422 break ; 423 case Int(SrcType::FSHI): 424 obsType = "FSHI" ; 425 break ; 426 case Int(SrcType::FHIOFF): 427 obsType = "FS:HIGHER:OFF" ; 428 break ; 429 case Int(SrcType::FHISKY): 430 obsType = "FS:HIGHER:SKY" ; 431 break ; 432 case Int(SrcType::FHIHOT): 433 obsType = "FS:HIGHER:HOT" ; 434 break ; 435 case Int(SrcType::FHIWARM): 436 obsType = "FS:HIGHER:WARM" ; 437 break ; 438 case Int(SrcType::FHICOLD): 439 obsType = "FS:HIGHER:COLD" ; 440 break ; 441 default: 442 obsType = "NOTYPE" ; 443 } 444 445 return obsType ; 446 } 447 448 } -
branches/mergetest/src/STWriter.h
r1391 r1789 36 36 #include <casa/aips.h> 37 37 #include <casa/Utilities/CountedPtr.h> 38 #include <casa/BasicSL/String.h> 38 39 39 40 #include "Logger.h" … … 84 85 void replacePtTab(const casa::Table& tab, const std::string& fname); 85 86 87 casa::String getObsTypes( casa::Int srctype ) ; 88 86 89 std::string format_; 87 90 PKSwriter* writer_; -
branches/mergetest/src/Scantable.cpp
r1743 r1789 11 11 // 12 12 #include <map> 13 #include <fstream> 13 14 14 15 #include <casa/aips.h> … … 24 25 #include <casa/Arrays/Vector.h> 25 26 #include <casa/Arrays/VectorSTLIterator.h> 27 #include <casa/Arrays/Slice.h> 26 28 #include <casa/BasicMath/Math.h> 27 29 #include <casa/BasicSL/Constants.h> … … 29 31 #include <casa/Containers/RecordField.h> 30 32 #include <casa/Utilities/GenSort.h> 33 #include <casa/Logging/LogIO.h> 31 34 32 35 #include <tables/Tables/TableParse.h> … … 106 109 { 107 110 initFactories(); 111 108 112 Table tab(name, Table::Update); 109 113 uInt version = tab.keywordSet().asuInt("VERSION"); … … 121 125 attach(); 122 126 } 127 /* 128 Scantable::Scantable(const std::string& name, Table::TableType ttype) : 129 type_(ttype) 130 { 131 initFactories(); 132 Table tab(name, Table::Update); 133 uInt version = tab.keywordSet().asuInt("VERSION"); 134 if (version != version_) { 135 throw(AipsError("Unsupported version of ASAP file.")); 136 } 137 if ( type_ == Table::Memory ) { 138 table_ = tab.copyToMemoryTable(generateName()); 139 } else { 140 table_ = tab; 141 } 142 143 attachSubtables(); 144 originalTable_ = table_; 145 attach(); 146 } 147 */ 123 148 124 149 Scantable::Scantable( const Scantable& other, bool clear ) … … 201 226 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO")); 202 227 228 td.addColumn(ScalarColumnDesc<uInt>("FLAGROW")); 229 203 230 td.addColumn(ScalarColumnDesc<Double>("TIME")); 204 231 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default … … 257 284 originalTable_ = table_; 258 285 } 259 260 286 261 287 void Scantable::attach() … … 285 311 mfocusidCol_.attach(table_, "FOCUS_ID"); 286 312 mmolidCol_.attach(table_, "MOLECULE_ID"); 313 314 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki) 315 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0); 316 317 } 318 319 template<class T, class T2> 320 void Scantable::attachAuxColumnDef(ScalarColumn<T>& col, 321 const String& colName, 322 const T2& defValue) 323 { 324 try { 325 col.attach(table_, colName); 326 } catch (TableError& err) { 327 String errMesg = err.getMesg(); 328 if (errMesg == "Table column " + colName + " is unknown") { 329 table_.addColumn(ScalarColumnDesc<T>(colName)); 330 col.attach(table_, colName); 331 col.fillColumn(static_cast<T>(defValue)); 332 } else { 333 throw; 334 } 335 } catch (...) { 336 throw; 337 } 338 } 339 340 template<class T, class T2> 341 void Scantable::attachAuxColumnDef(ArrayColumn<T>& col, 342 const String& colName, 343 const Array<T2>& defValue) 344 { 345 try { 346 col.attach(table_, colName); 347 } catch (TableError& err) { 348 String errMesg = err.getMesg(); 349 if (errMesg == "Table column " + colName + " is unknown") { 350 table_.addColumn(ArrayColumnDesc<T>(colName)); 351 col.attach(table_, colName); 352 353 int size = 0; 354 ArrayIterator<T2>& it = defValue.begin(); 355 while (it != defValue.end()) { 356 ++size; 357 ++it; 358 } 359 IPosition ip(1, size); 360 Array<T>& arr(ip); 361 for (int i = 0; i < size; ++i) 362 arr[i] = static_cast<T>(defValue[i]); 363 364 col.fillColumn(arr); 365 } else { 366 throw; 367 } 368 } catch (...) { 369 throw; 370 } 287 371 } 288 372 … … 627 711 } 628 712 713 void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag) 714 { 715 for (uInt i=0; i<table_.nrow(); ++i) { 716 Vector<uChar> flgs = flagsCol_(i); 717 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs); 718 flagsCol_.put(i, flgs); 719 } 720 } 721 722 std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag) 723 { 724 Vector<uChar> flags; 725 flagsCol_.get(uInt(whichrow), flags); 726 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags); 727 Vector<Bool> bflag(flags.shape()); 728 convertArray(bflag, flags); 729 //bflag = !bflag; 730 731 std::vector<bool> mask; 732 bflag.tovector(mask); 733 return mask; 734 } 735 736 void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag, 737 Vector<uChar> flgs) 738 { 739 Vector<Float> spcs = specCol_(whichrow); 740 uInt nchannel = nchan(); 741 if (spcs.nelements() != nchannel) { 742 throw(AipsError("Data has incorrect number of channels")); 743 } 744 uChar userflag = 1 << 7; 745 if (unflag) { 746 userflag = 0 << 7; 747 } 748 if (clipoutside) { 749 for (uInt j = 0; j < nchannel; ++j) { 750 Float spc = spcs(j); 751 if ((spc >= uthres) || (spc <= dthres)) { 752 flgs(j) = userflag; 753 } 754 } 755 } else { 756 for (uInt j = 0; j < nchannel; ++j) { 757 Float spc = spcs(j); 758 if ((spc < uthres) && (spc > dthres)) { 759 flgs(j) = userflag; 760 } 761 } 762 } 763 } 764 629 765 void Scantable::flag(const std::vector<bool>& msk, bool unflag) 630 766 { … … 673 809 flagsCol_.put(i, flgs); 674 810 } 811 } 812 813 void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag) 814 { 815 if ( selector_.empty() && (rows.size() == table_.nrow()) ) 816 throw(AipsError("Trying to flag whole scantable.")); 817 818 uInt rowflag = (unflag ? 0 : 1); 819 std::vector<uInt>::const_iterator it; 820 for (it = rows.begin(); it != rows.end(); ++it) 821 flagrowCol_.put(*it, rowflag); 675 822 } 676 823 … … 793 940 table_.keywordSet().get("FluxUnit", tmp); 794 941 oss << setw(15) << "Flux Unit:" << tmp << endl; 795 Vector<Double> vec(moleculeTable_.getRestFrequencies()); 942 //Vector<Double> vec(moleculeTable_.getRestFrequencies()); 943 int nid = moleculeTable_.nrow(); 944 Bool firstline = True; 796 945 oss << setw(15) << "Rest Freqs:"; 797 if (vec.nelements() > 0) { 798 oss << setprecision(10) << vec << " [Hz]" << endl; 799 } else { 800 oss << "none" << endl; 946 for (int i=0; i<nid; i++) { 947 Table t = table_(table_.col("MOLECULE_ID") == i); 948 if (t.nrow() > 0) { 949 Vector<Double> vec(moleculeTable_.getRestFrequency(i)); 950 if (vec.nelements() > 0) { 951 if (firstline) { 952 oss << setprecision(10) << vec << " [Hz]" << endl; 953 firstline=False; 954 } 955 else{ 956 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl; 957 } 958 } else { 959 oss << "none" << endl; 960 } 961 } 801 962 } 802 963 … … 901 1062 const MDirection& md = getDirection(whichrow); 902 1063 const MEpoch& me = timeCol_(whichrow); 903 Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1064 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 1065 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 904 1066 return freqTable_.getSpectralCoordinate(md, mp, me, rf, 905 1067 mfreqidCol_(whichrow)); … … 975 1137 const MDirection& md = getDirection(whichrow); 976 1138 const MEpoch& me = timeCol_(whichrow); 977 const Double& rf = mmolidCol_(whichrow); 1139 //const Double& rf = mmolidCol_(whichrow); 1140 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow)); 978 1141 SpectralCoordinate spc = 979 1142 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow)); … … 992 1155 } 993 1156 994 void Scantable::setRestFrequencies( double rf, const std::string& name, 1157 /** 1158 void asap::Scantable::setRestFrequencies( double rf, const std::string& name, 995 1159 const std::string& unit ) 1160 **/ 1161 void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name, 1162 const std::string& unit ) 1163 996 1164 { 997 1165 ///@todo lookup in line table to fill in name and formattedname 998 1166 Unit u(unit); 999 Quantum<Double> urf(rf, u); 1000 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, ""); 1167 //Quantum<Double> urf(rf, u); 1168 Quantum<Vector<Double> >urf(rf, u); 1169 Vector<String> formattedname(0); 1170 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl; 1171 1172 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, ""); 1173 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname); 1001 1174 TableVector<uInt> tabvec(table_, "MOLECULE_ID"); 1002 1175 tabvec = id; 1003 1176 } 1004 1177 1005 void Scantable::setRestFrequencies( const std::string& name ) 1178 /** 1179 void asap::Scantable::setRestFrequencies( const std::string& name ) 1006 1180 { 1007 1181 throw(AipsError("setRestFrequencies( const std::string& name ) NYI")); 1182 ///@todo implement 1183 } 1184 **/ 1185 void Scantable::setRestFrequencies( const vector<std::string>& name ) 1186 { 1187 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI")); 1008 1188 ///@todo implement 1009 1189 } … … 1044 1224 void Scantable::addFit( const STFitEntry& fit, int row ) 1045 1225 { 1046 cout << mfitidCol_(uInt(row)) << endl; 1226 //cout << mfitidCol_(uInt(row)) << endl; 1227 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ; 1228 os << mfitidCol_(uInt(row)) << LogIO::POST ; 1047 1229 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row))); 1048 1230 mfitidCol_.put(uInt(row), id); … … 1078 1260 Table subt = t( t.col("SCAN") == scanlist[i]+1 ); 1079 1261 if (subt.nrow()==0) { 1080 cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl; 1262 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl; 1263 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1264 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST; 1081 1265 ret = 1; 1082 1266 break; … … 1090 1274 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 ); 1091 1275 if ( subt2.nrow() == 0) { 1092 cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl; 1276 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1277 1278 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl; 1279 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST; 1093 1280 ret = 1; 1094 1281 break; … … 1100 1287 if (scan1seqn == 1 && scan2seqn == 2) { 1101 1288 if (laston1 == laston2) { 1102 cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1289 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1290 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1291 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST; 1103 1292 i +=1; 1104 1293 } 1105 1294 else { 1106 cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1295 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1296 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl; 1297 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST; 1107 1298 } 1108 1299 } 1109 1300 else if (scan1seqn==2 && scan2seqn == 1) { 1110 1301 if (laston1 == laston2) { 1111 cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl; 1302 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1303 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl; 1304 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST; 1112 1305 ret = 1; 1113 1306 break; … … 1115 1308 } 1116 1309 else { 1117 cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl; 1310 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1311 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl; 1312 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST; 1118 1313 ret = 1; 1119 1314 break; … … 1122 1317 } 1123 1318 else { 1124 cerr<<"The scan does not appear to be standard obsevation."<<endl; 1319 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1320 //cerr<<"The scan does not appear to be standard obsevation."<<endl; 1321 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST; 1125 1322 } 1126 1323 //if ( i >= nscan ) break; … … 1128 1325 } 1129 1326 else { 1130 cerr<<"No reference to GBT_GO table."<<endl; 1327 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ; 1328 //cerr<<"No reference to GBT_GO table."<<endl; 1329 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST; 1131 1330 ret = 1; 1132 1331 } … … 1142 1341 } 1143 1342 1343 void asap::Scantable::reshapeSpectrum( int nmin, int nmax ) 1344 throw( casa::AipsError ) 1345 { 1346 // assumed that all rows have same nChan 1347 Vector<Float> arr = specCol_( 0 ) ; 1348 int nChan = arr.nelements() ; 1349 1350 // if nmin < 0 or nmax < 0, nothing to do 1351 if ( nmin < 0 ) { 1352 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ; 1353 } 1354 if ( nmax < 0 ) { 1355 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ; 1356 } 1357 1358 // if nmin > nmax, exchange values 1359 if ( nmin > nmax ) { 1360 int tmp = nmax ; 1361 nmax = nmin ; 1362 nmin = tmp ; 1363 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1364 os << "Swap values. Applied range is [" 1365 << nmin << ", " << nmax << "]" << LogIO::POST ; 1366 } 1367 1368 // if nmin exceeds nChan, nothing to do 1369 if ( nmin >= nChan ) { 1370 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ; 1371 } 1372 1373 // if nmax exceeds nChan, reset nmax to nChan 1374 if ( nmax >= nChan ) { 1375 if ( nmin == 0 ) { 1376 // nothing to do 1377 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1378 os << "Whole range is selected. Nothing to do." << LogIO::POST ; 1379 return ; 1380 } 1381 else { 1382 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ; 1383 os << "Specified maximum exceeds nChan. Applied range is [" 1384 << nmin << ", " << nChan-1 << "]." << LogIO::POST ; 1385 nmax = nChan - 1 ; 1386 } 1387 } 1388 1389 // reshape specCol_ and flagCol_ 1390 for ( int irow = 0 ; irow < nrow() ; irow++ ) { 1391 reshapeSpectrum( nmin, nmax, irow ) ; 1392 } 1393 1394 // update FREQUENCIES subtable 1395 Double refpix ; 1396 Double refval ; 1397 Double increment ; 1398 int freqnrow = freqTable_.table().nrow() ; 1399 Vector<uInt> oldId( freqnrow ) ; 1400 Vector<uInt> newId( freqnrow ) ; 1401 for ( int irow = 0 ; irow < freqnrow ; irow++ ) { 1402 freqTable_.getEntry( refpix, refval, increment, irow ) ; 1403 /*** 1404 * need to shift refpix to nmin 1405 * note that channel nmin in old index will be channel 0 in new one 1406 ***/ 1407 refval = refval - ( refpix - nmin ) * increment ; 1408 refpix = 0 ; 1409 freqTable_.setEntry( refpix, refval, increment, irow ) ; 1410 } 1411 1412 // update nchan 1413 int newsize = nmax - nmin + 1 ; 1414 table_.rwKeywordSet().define( "nChan", newsize ) ; 1415 1416 // update bandwidth 1417 // assumed all spectra in the scantable have same bandwidth 1418 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ; 1419 1420 return ; 1421 } 1422 1423 void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow ) 1424 { 1425 // reshape specCol_ and flagCol_ 1426 Vector<Float> oldspec = specCol_( irow ) ; 1427 Vector<uChar> oldflag = flagsCol_( irow ) ; 1428 uInt newsize = nmax - nmin + 1 ; 1429 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ; 1430 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ; 1431 1432 return ; 1433 } 1434 1435 void asap::Scantable::regridChannel( int nChan, double dnu ) 1436 { 1437 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ; 1438 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ; 1439 // assumed that all rows have same nChan 1440 Vector<Float> arr = specCol_( 0 ) ; 1441 int oldsize = arr.nelements() ; 1442 1443 // if oldsize == nChan, nothing to do 1444 if ( oldsize == nChan ) { 1445 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ; 1446 return ; 1447 } 1448 1449 // if oldChan < nChan, unphysical operation 1450 if ( oldsize < nChan ) { 1451 os << "Unphysical operation. Nothing to do." << LogIO::POST ; 1452 return ; 1453 } 1454 1455 // change channel number for specCol_ and flagCol_ 1456 Vector<Float> newspec( nChan, 0 ) ; 1457 Vector<uChar> newflag( nChan, false ) ; 1458 vector<string> coordinfo = getCoordInfo() ; 1459 string oldinfo = coordinfo[0] ; 1460 coordinfo[0] = "Hz" ; 1461 setCoordInfo( coordinfo ) ; 1462 for ( int irow = 0 ; irow < nrow() ; irow++ ) { 1463 regridChannel( nChan, dnu, irow ) ; 1464 } 1465 coordinfo[0] = oldinfo ; 1466 setCoordInfo( coordinfo ) ; 1467 1468 1469 // NOTE: this method does not update metadata such as 1470 // FREQUENCIES subtable, nChan, Bandwidth, etc. 1471 1472 return ; 1473 } 1474 1475 void asap::Scantable::regridChannel( int nChan, double dnu, int irow ) 1476 { 1477 // logging 1478 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ; 1479 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ; 1480 1481 Vector<Float> oldspec = specCol_( irow ) ; 1482 Vector<uChar> oldflag = flagsCol_( irow ) ; 1483 Vector<Float> newspec( nChan, 0 ) ; 1484 Vector<uChar> newflag( nChan, false ) ; 1485 1486 // regrid 1487 vector<double> abcissa = getAbcissa( irow ) ; 1488 int oldsize = abcissa.size() ; 1489 double olddnu = abcissa[1] - abcissa[0] ; 1490 //int refChan = 0 ; 1491 //double frac = 0.0 ; 1492 //double wedge = 0.0 ; 1493 //double pile = 0.0 ; 1494 int ichan = 0 ; 1495 double wsum = 0.0 ; 1496 Vector<Float> z( nChan ) ; 1497 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ; 1498 for ( int ii = 1 ; ii < nChan ; ii++ ) 1499 z[ii] = z[ii-1] + dnu ; 1500 Vector<Float> zi( nChan+1 ) ; 1501 Vector<Float> yi( oldsize + 1 ) ; 1502 zi[0] = z[0] - 0.5 * dnu ; 1503 zi[1] = z[0] + 0.5 * dnu ; 1504 for ( int ii = 2 ; ii < nChan ; ii++ ) 1505 zi[ii] = zi[ii-1] + dnu ; 1506 zi[nChan] = z[nChan-1] + 0.5 * dnu ; 1507 yi[0] = abcissa[0] - 0.5 * olddnu ; 1508 yi[1] = abcissa[1] + 0.5 * olddnu ; 1509 for ( int ii = 2 ; ii < oldsize ; ii++ ) 1510 yi[ii] = abcissa[ii-1] + olddnu ; 1511 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ; 1512 if ( dnu > 0.0 ) { 1513 for ( int ii = 0 ; ii < nChan ; ii++ ) { 1514 double zl = zi[ii] ; 1515 double zr = zi[ii+1] ; 1516 for ( int j = ichan ; j < oldsize ; j++ ) { 1517 double yl = yi[j] ; 1518 double yr = yi[j+1] ; 1519 if ( yl <= zl ) { 1520 if ( yr <= zl ) { 1521 continue ; 1522 } 1523 else if ( yr <= zr ) { 1524 newspec[ii] += oldspec[j] * ( yr - zl ) ; 1525 newflag[ii] = newflag[ii] || oldflag[j] ; 1526 wsum += ( yr - zl ) ; 1527 } 1528 else { 1529 newspec[ii] += oldspec[j] * dnu ; 1530 newflag[ii] = newflag[ii] || oldflag[j] ; 1531 wsum += dnu ; 1532 ichan = j ; 1533 break ; 1534 } 1535 } 1536 else if ( yl < zr ) { 1537 if ( yr <= zr ) { 1538 newspec[ii] += oldspec[j] * ( yr - yl ) ; 1539 newflag[ii] = newflag[ii] || oldflag[j] ; 1540 wsum += ( yr - yl ) ; 1541 } 1542 else { 1543 newspec[ii] += oldspec[j] * ( zr - yl ) ; 1544 newflag[ii] = newflag[ii] || oldflag[j] ; 1545 wsum += ( zr - yl ) ; 1546 ichan = j ; 1547 break ; 1548 } 1549 } 1550 else { 1551 ichan = j - 1 ; 1552 break ; 1553 } 1554 } 1555 newspec[ii] /= wsum ; 1556 wsum = 0.0 ; 1557 } 1558 } 1559 else if ( dnu < 0.0 ) { 1560 for ( int ii = 0 ; ii < nChan ; ii++ ) { 1561 double zl = zi[ii] ; 1562 double zr = zi[ii+1] ; 1563 for ( int j = ichan ; j < oldsize ; j++ ) { 1564 double yl = yi[j] ; 1565 double yr = yi[j+1] ; 1566 if ( yl >= zl ) { 1567 if ( yr >= zl ) { 1568 continue ; 1569 } 1570 else if ( yr >= zr ) { 1571 newspec[ii] += oldspec[j] * abs( yr - zl ) ; 1572 newflag[ii] = newflag[ii] || oldflag[j] ; 1573 wsum += abs( yr - zl ) ; 1574 } 1575 else { 1576 newspec[ii] += oldspec[j] * abs( dnu ) ; 1577 newflag[ii] = newflag[ii] || oldflag[j] ; 1578 wsum += abs( dnu ) ; 1579 ichan = j ; 1580 break ; 1581 } 1582 } 1583 else if ( yl > zr ) { 1584 if ( yr >= zr ) { 1585 newspec[ii] += oldspec[j] * abs( yr - yl ) ; 1586 newflag[ii] = newflag[ii] || oldflag[j] ; 1587 wsum += abs( yr - yl ) ; 1588 } 1589 else { 1590 newspec[ii] += oldspec[j] * abs( zr - yl ) ; 1591 newflag[ii] = newflag[ii] || oldflag[j] ; 1592 wsum += abs( zr - yl ) ; 1593 ichan = j ; 1594 break ; 1595 } 1596 } 1597 else { 1598 ichan = j - 1 ; 1599 break ; 1600 } 1601 } 1602 newspec[ii] /= wsum ; 1603 wsum = 0.0 ; 1604 } 1605 } 1606 // * ichan = 0 1607 // ***/ 1608 // //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ; 1609 // pile += dnu ; 1610 // wedge = olddnu * ( refChan + 1 ) ; 1611 // while ( wedge < pile ) { 1612 // newspec[0] += olddnu * oldspec[refChan] ; 1613 // newflag[0] = newflag[0] || oldflag[refChan] ; 1614 // //ofs << "channel " << refChan << " is included in new channel 0" << endl ; 1615 // refChan++ ; 1616 // wedge += olddnu ; 1617 // wsum += olddnu ; 1618 // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ; 1619 // } 1620 // frac = ( wedge - pile ) / olddnu ; 1621 // wsum += ( 1.0 - frac ) * olddnu ; 1622 // newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ; 1623 // newflag[0] = newflag[0] || oldflag[refChan] ; 1624 // //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ; 1625 // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ; 1626 // newspec[0] /= wsum ; 1627 // //ofs << "newspec[0] = " << newspec[0] << endl ; 1628 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1629 1630 // /*** 1631 // * ichan = 1 - nChan-2 1632 // ***/ 1633 // for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) { 1634 // pile += dnu ; 1635 // newspec[ichan] += frac * olddnu * oldspec[refChan] ; 1636 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1637 // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ; 1638 // refChan++ ; 1639 // wedge += olddnu ; 1640 // wsum = frac * olddnu ; 1641 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1642 // while ( wedge < pile ) { 1643 // newspec[ichan] += olddnu * oldspec[refChan] ; 1644 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1645 // //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ; 1646 // refChan++ ; 1647 // wedge += olddnu ; 1648 // wsum += olddnu ; 1649 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1650 // } 1651 // frac = ( wedge - pile ) / olddnu ; 1652 // wsum += ( 1.0 - frac ) * olddnu ; 1653 // newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ; 1654 // newflag[ichan] = newflag[ichan] || oldflag[refChan] ; 1655 // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ; 1656 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1657 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ; 1658 // newspec[ichan] /= wsum ; 1659 // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ; 1660 // } 1661 1662 // /*** 1663 // * ichan = nChan-1 1664 // ***/ 1665 // // NOTE: Assumed that all spectra have the same bandwidth 1666 // pile += dnu ; 1667 // newspec[nChan-1] += frac * olddnu * oldspec[refChan] ; 1668 // newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ; 1669 // //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ; 1670 // refChan++ ; 1671 // wedge += olddnu ; 1672 // wsum = frac * olddnu ; 1673 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1674 // for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) { 1675 // newspec[nChan-1] += olddnu * oldspec[jchan] ; 1676 // newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ; 1677 // wsum += olddnu ; 1678 // //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ; 1679 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1680 // } 1681 // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ; 1682 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ; 1683 // newspec[nChan-1] /= wsum ; 1684 // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ; 1685 1686 // specCol_.put( irow, newspec ) ; 1687 // flagsCol_.put( irow, newflag ) ; 1688 1689 // // ofs.close() ; 1690 1691 1692 return ; 1693 } 1694 1144 1695 std::vector<float> Scantable::getWeather(int whichrow) const 1145 1696 { … … 1154 1705 1155 1706 } 1156 1707 //namespace asap -
branches/mergetest/src/Scantable.h
r1730 r1789 29 29 30 30 #include <coordinates/Coordinates/SpectralCoordinate.h> 31 32 #include <casa/Arrays/Vector.h> 33 #include <casa/Quanta/Quantum.h> 34 35 #include <casa/Exceptions/Error.h> 31 36 32 37 #include "Logger.h" … … 226 231 227 232 /** 233 * Flag the data in a row-based manner. (CAS-1433 Wataru Kawasaki) 234 * param[in] rows list of row numbers to be flagged 235 */ 236 void flagRow( const std::vector<casa::uInt>& rows = std::vector<casa::uInt>(), bool unflag=false); 237 238 /** 239 * Get flagRow info at the specified row. If true, the whole data 240 * at the row should be flagged. 241 */ 242 bool getFlagRow(int whichrow) const 243 { return (flagrowCol_(whichrow) > 0); } 244 245 /** 246 * Flag the data outside a specified range (in a channel-based manner). 247 * (CAS-1807 Wataru Kawasaki) 248 */ 249 void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag); 250 251 /** 252 * Return a list of booleans with the size of nchan for a specified row, to get info 253 * about which channel is clipped. 254 */ 255 std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag); 256 void srchChannelsToClip(casa::uInt whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag, 257 casa::Vector<casa::uChar> flgs); 258 259 /** 228 260 * Return a list of row numbers with respect to the original table. 229 261 * @return a list of unsigned ints … … 276 308 std::vector<uint> getScanNos() const { return getNumbers(scanCol_); } 277 309 int getScan(int whichrow) const { return scanCol_(whichrow); } 310 311 //TT addition 312 std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); } 278 313 279 314 /** … … 300 335 { return azCol_(whichrow); } 301 336 float getParAngle(int whichrow) const 302 { return focus().getParAngle(mfocusidCol_(whichrow)); } 337 { return focus().getParAngle(mfocusidCol_(whichrow)); } 338 int getTcalId(int whichrow) const 339 { return mtcalidCol_(whichrow); } 303 340 304 341 std::string getSourceName(int whichrow) const … … 353 390 std::vector<double> getRestFrequencies() const 354 391 { return moleculeTable_.getRestFrequencies(); } 355 392 std::vector<double> getRestFrequency(int id) const 393 { return moleculeTable_.getRestFrequency(id); } 394 395 /** 356 396 void setRestFrequencies(double rf, const std::string& name = "", 357 397 const std::string& = "Hz"); 358 void setRestFrequencies(const std::string& name); 398 **/ 399 // Modified by Takeshi Nakazato 05/09/2008 400 /*** 401 void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "", 402 const std::string& = "Hz"); 403 ***/ 404 void setRestFrequencies(vector<double> rf, 405 const vector<std::string>& name = vector<std::string>(1,""), 406 const std::string& = "Hz"); 407 408 //void setRestFrequencies(const std::string& name); 409 void setRestFrequencies(const vector<std::string>& name); 359 410 360 411 void shift(int npix); … … 415 466 { focus().setParallactify(flag); } 416 467 468 /** 469 * Reshape spectrum 470 * @param[in] nmin, nmax minimum and maximum channel 471 * @param[in] irow row number 472 * 473 * 30/07/2008 Takeshi Nakazato 474 **/ 475 void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError ); 476 void reshapeSpectrum( int nmin, int nmax, int irow ) ; 477 478 /** 479 * Change channel number under fixed bandwidth 480 * @param[in] nchan, dnu new channel number and spectral resolution 481 * @param[in] irow row number 482 * 483 * 27/08/2008 Takeshi Nakazato 484 **/ 485 void regridChannel( int nchan, double dnu ) ; 486 void regridChannel( int nchan, double dnu, int irow ) ; 487 488 417 489 private: 418 490 … … 489 561 casa::ScalarColumn<casa::Float> elCol_; 490 562 casa::ScalarColumn<casa::String> srcnCol_, fldnCol_; 491 casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_ ;563 casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_; 492 564 casa::ScalarColumn<casa::Int> rbeamCol_, srctCol_; 493 565 casa::ArrayColumn<casa::Float> specCol_, tsysCol_; … … 510 582 void initFactories(); 511 583 584 /** 585 * Add an auxiliary column to the main table and attach it to a 586 * cached column. Use for adding new columns that the original asap2 587 * tables do not have. 588 * @param[in] col reference to the cached column to be attached 589 * @param[in] colName column name in asap table 590 * @param[in] defValue default value to fill in the column 591 * 592 * 25/10/2009 Wataru Kawasaki 593 */ 594 template<class T, class T2> void attachAuxColumnDef(casa::ScalarColumn<T>&, 595 const casa::String&, 596 const T2&); 597 template<class T, class T2> void attachAuxColumnDef(casa::ArrayColumn<T>&, 598 const casa::String&, 599 const casa::Array<T2>&); 512 600 }; 513 601 -
branches/mergetest/src/ScantableWrapper.h
r1730 r1789 109 109 { table_->flag(msk, unflag); } 110 110 111 void flagRow(const std::vector<casa::uInt>& rows=std::vector<casa::uInt>(), bool unflag=false) 112 { table_->flagRow(rows, unflag); } 113 114 bool getFlagRow(int whichrow=0) const 115 { return table_->getFlagRow(whichrow); } 116 117 void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside=true, bool unflag=false) 118 { table_->clip(uthres, dthres, clipoutside, unflag); } 119 120 std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag) const 121 { return table_->getClipMask(whichrow, uthres, dthres, clipoutside, unflag); } 122 111 123 std::string getSourceName(int whichrow=0) const 112 124 { return table_->getSourceName(whichrow); } … … 134 146 std::vector<uint> getScanNos() { return table_->getScanNos(); } 135 147 int getScan(int whichrow) const {return table_->getScan(whichrow);} 148 std::vector<uint> getMolNos() { return table_->getMolNos();} 136 149 137 150 STSelector getSelection() const { return table_->getSelection(); } … … 158 171 { table_->shift(npix); } 159 172 173 /** 174 commented out by TT 160 175 void setRestFrequencies(double rf, const std::string& name, 161 176 const std::string& unit) 162 177 { table_->setRestFrequencies(rf, name, unit); } 178 **/ 179 void setRestFrequencies(vector<double> rf, const vector<std::string>& name, 180 const std::string& unit) 181 { table_->setRestFrequencies(rf, name, unit); } 182 163 183 /* 164 184 void setRestFrequencies(const std::string& name) { … … 167 187 */ 168 188 189 /* 169 190 std::vector<double> getRestFrequencies() const 170 191 { return table_->getRestFrequencies(); } 192 */ 193 std::vector<double> getRestFrequency(int id) const 194 { return table_->getRestFrequency(id); } 171 195 172 196 void setCoordInfo(std::vector<string> theinfo) { … … 209 233 int checkScanInfo(const vector<int>& scanlist) const 210 234 { return table_->checkScanInfo(scanlist); } 211 235 212 236 std::vector<double> getDirectionVector(int whichrow) const 213 237 { return table_->getDirectionVector(whichrow); } … … 222 246 std::vector<float> getWeather(int whichrow) const 223 247 { return table_->getWeather(whichrow); } 248 249 void reshapeSpectrum( int nmin, int nmax ) 250 { table_->reshapeSpectrum( nmin, nmax ); } 224 251 225 252 private: … … 229 256 } // namespace 230 257 #endif 258 -
branches/mergetest/src/python_STMath.cpp
r1391 r1789 49 49 .def("_averagebeams", &STMathWrapper::averageBeams) 50 50 .def("_unaryop", &STMathWrapper::unaryOperate) 51 .def("_arrayop", &STMathWrapper::arrayOperate) 52 //.def("_array2dop", &STMathWrapper::array2dOperate) 51 53 .def("_binaryop", &STMathWrapper::binaryOperate) 52 54 .def("_auto_quotient", &STMathWrapper::autoQuotient) … … 57 59 .def("_dofs", &STMathWrapper::dofs) 58 60 .def("_stats", &STMathWrapper::statistic) 61 .def("_minmaxchan", &STMathWrapper::minMaxChan) 59 62 .def("_freqswitch", &STMathWrapper::freqSwitch) 60 63 .def("_bin", &STMathWrapper::bin) … … 73 76 .def("_mx_extract", &STMathWrapper::mxExtract) 74 77 .def("_lag_flag", &STMathWrapper::lagFlag) 78 // testing average spectra with different channel/resolution 79 .def("_new_average", &STMathWrapper::new_average) 80 // cwcal 81 .def("cwcal", &STMathWrapper::cwcal) 82 .def("almacal", &STMathWrapper::almacal) 75 83 ; 76 84 }; -
branches/mergetest/src/python_STSelector.cpp
r939 r1789 29 29 .def("_getscans", &STSelector::getScans) 30 30 .def("_getcycles", &STSelector::getCycles) 31 .def("_gettypes", &STSelector::getTypes) 31 32 .def("_gettaql", &STSelector::getTaQL) 32 33 .def("_getorder", &STSelector::getSortOrder) … … 41 42 .def("_settaql", &STSelector::setTaQL) 42 43 .def("_setorder", &STSelector::setSortOrder) 44 .def("_setrows", &STSelector::setRows) 45 .def("_settypes", &STSelector::setTypes) 43 46 .def("_empty", &STSelector::empty) 44 47 ; -
branches/mergetest/src/python_Scantable.cpp
r1730 r1789 58 58 .def("getscannos", &ScantableWrapper::getScanNos) 59 59 .def("getcycle", &ScantableWrapper::getCycle) 60 .def("getmolnos", &ScantableWrapper::getMolNos) 60 61 .def("nif", &ScantableWrapper::nif, 61 62 (boost::python::arg("scanno")=-1) ) … … 87 88 .def("_getmask", &ScantableWrapper::getMask, 88 89 (boost::python::arg("whichrow")=0) ) 90 .def("_getclipmask", &ScantableWrapper::getClipMask, 91 (boost::python::arg("whichrow")=0) ) 89 92 .def("_gettsys", &ScantableWrapper::getTsys) 90 93 .def("_getsourcename", &ScantableWrapper::getSourceName, … … 103 106 (boost::python::arg("whichrow")=0) ) 104 107 .def("get_antennaname", &ScantableWrapper::getAntennaName) 105 .def("_flag", &ScantableWrapper::flag, 106 (boost::python::arg("unflag") = false) ) 108 .def("_flag", &ScantableWrapper::flag) 109 .def("_flag_row", &ScantableWrapper::flagRow) 110 .def("_getflagrow", &ScantableWrapper::getFlagRow, 111 (boost::python::arg("whichrow")=0) ) 112 .def("_clip", &ScantableWrapper::clip, 113 (boost::python::arg("clipoutside")=true, 114 boost::python::arg("unflag")=false) ) 107 115 .def("_save", &ScantableWrapper::makePersistent) 108 116 .def("_summary", &ScantableWrapper::summary, 109 117 (boost::python::arg("verbose")=true) ) 110 .def("_getrestfreqs", &ScantableWrapper::getRestFrequencies) 118 //.def("_getrestfreqs", &ScantableWrapper::getRestFrequencies) 119 .def("_getrestfreqs", &ScantableWrapper::getRestFrequency) 111 120 .def("_setrestfreqs", &ScantableWrapper::setRestFrequencies) 112 121 .def("shift_refpix", &ScantableWrapper::shift) … … 127 136 .def("get_coordinate", &ScantableWrapper::getCoordinate) 128 137 .def("_get_weather", &ScantableWrapper::getWeather) 138 .def("_reshape", &ScantableWrapper::reshapeSpectrum, 139 (boost::python::arg("nmin")=-1, 140 boost::python::arg("nmax")=-1) ) 129 141 ; 130 142 }; -
branches/mergetest/src/python_asap.cpp
r1715 r1789 77 77 asap::python::python_STCoordinate(); 78 78 asap::python::python_STAtmosphere(); 79 asap::python::python_SrcType(); 79 80 80 81 #ifndef HAVE_LIBPYRAP -
branches/mergetest/src/python_asap.h
r1715 r1789 49 49 void python_STCoordinate(); 50 50 void python_STAtmosphere(); 51 void python_SrcType(); 51 52 52 53 } // python
Note:
See TracChangeset
for help on using the changeset viewer.