Changeset 2580 for trunk/src


Ignore:
Timestamp:
06/28/12 14:22:10 (12 years ago)
Author:
ShinnosukeKawakami
Message:

hpc33 merged asap-trunk

Location:
trunk
Files:
15 edited
4 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/src

  • trunk/src/CMakeLists.txt

    r2465 r2580  
    6060     ${SRCDIR}/STUpgrade.cpp
    6161     ${SRCDIR}/STGrid.cpp
     62     ${SRCDIR}/STIdxIter.cpp
    6263     ${SRCDIR}/Templates.cpp )
    6364
    … …  
    8182     ${SRCDIR}/python_LogSink.cpp
    8283     ${SRCDIR}/python_STGrid.cpp
     84     ${SRCDIR}/python_Iterator.cpp
    8385     ${SRCDIR}/python_asap.cpp )
    8486

  • trunk/src/FillerWrapper.h

  • trunk/src/MSFiller.cpp

    r2554 r2580  
    304304    // set dummy epoch
    305305    mf.set( currentTime ) ;
    306 
    307     // initialize dirType
    308     dirType = MDirection::N_Types ;
    309306
    310307    //
    … …  
    677674    MDirection::getType( dirType, pointingRef ) ;
    678675    getpt = True ;
     676
     677    // initialize toj2000 and toazel
     678    initConvert() ;
    679679  }
    680680  void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it )
    … …  
    710710  uInt getFilledRowNum() { return rowidx ; }
    711711private:
     712  void initConvert()
     713  {
     714    toj2000 = MDirection::Convert( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     715    toazel = MDirection::Convert( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
     716  }
     717
    712718  void fluxUnit( String &u )
    713719  {
    … …  
    757763      String ref = dir.getRefString() ;
    758764      MDirection::getType( dirType, ref ) ;
     765     
     766      // initialize toj2000 and toazel
     767      initConvert() ;
    759768    }
    760769
    … …  
    954963                    pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
    955964    }
    956     mf.resetEpoch( currentTime ) ;
     965    mf.set( currentTime ) ;
    957966    Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
    958967    if ( dirType != MDirection::J2000 ) {
    959       MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
    960968      dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
    961969    }
    … …  
    964972    }
    965973    if ( dirType != MDirection::AZELGEO ) {
    966       MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
    967       //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
    968974      azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
    969975    }
    … …  
    977983  {
    978984    dir = sourceDir.getAngle( "rad" ).getValue() ;
    979     mf.resetEpoch( currentTime ) ;
    980     MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
     985    mf.set( currentTime ) ;
    981986    azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
    982987    if ( dirType != MDirection::J2000 ) {
    983       MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
    984988      dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
    985989    }
    … …  
    12611265  MPosition antpos;
    12621266  MEpoch currentTime;
     1267  MeasFrame mf;
     1268  MDirection::Convert toj2000;
     1269  MDirection::Convert toazel;
    12631270  map<Int,uInt> ifmap;
    12641271  Block<uInt> polnos;
    … …  
    12931300  Matrix<Float> sp;
    12941301  Matrix<uChar> fl;
    1295   MeasFrame mf;
    12961302
    12971303  // MS MAIN columns

  • trunk/src/RowAccumulator.cpp

    r2141 r2580  
    8181                                 const Double time)
    8282{
    83   doAddSpectrum(v, m, tsys, interval, time, False);
    84   doAddSpectrum(v, m, tsys, interval, time, True);  // CAS-2776
     83//   doAddSpectrum(v, m, tsys, interval, time, False);
     84//   doAddSpectrum(v, m, tsys, interval, time, True);  // CAS-2776
     85  doAddSpectrum2( v, m, tsys, interval, time ) ;
    8586}
    8687
    … …  
    112113}
    113114
     115void RowAccumulator::doAddSpectrum2(const Vector<Float>& v,
     116                                    const Vector<Bool>& m,
     117                                    const Vector<Float>& tsys,
     118                                    const Double interval,
     119                                    const Double time)
     120{
     121  const MaskedArray<Float> vadd(v, m);
     122  const MaskedArray<Float> vaddN(v, !m);
     123  Float totalWeight = getTotalWeight( vadd, tsys, interval, time, False ) ;
     124  Float totalWeightNoMask = getTotalWeight( vaddN, tsys, interval, time, True ) ;
     125  // process spectrum with input mask and spectrum with inverted mask
     126  // together
     127  // prepare pointers
     128  Bool vD, mD ;
     129  const Float *v_p = v.getStorage( vD ) ;
     130  const Float *v_wp = v_p ;
     131  const Bool *m_p = m.getStorage( mD ) ;
     132  const Bool *m_wp = m_p ;
     133
     134  Bool spD, spND, wgtD, wgtND, nD, nND ;
     135  Float *sp_p = spectrum_.getRWArrayStorage( spD ) ;
     136  Float *sp_wp = sp_p ;
     137  Float *wgt_p = weightSum_.getRWArrayStorage( wgtD ) ;
     138  Float *wgt_wp = wgt_p ;
     139  Float *spN_p = spectrumNoMask_.getRWArrayStorage( spND ) ;
     140  Float *spN_wp = spN_p ;
     141  Float *wgtN_p = weightSumNoMask_.getRWArrayStorage( wgtND ) ;
     142  Float *wgtN_wp = wgtN_p ;
     143  uInt *n_p = n_.getRWArrayStorage( nD ) ;
     144  uInt *n_wp = n_p ;
     145  uInt *nN_p = nNoMask_.getRWArrayStorage( nND ) ;
     146  uInt *nN_wp = nN_p ;
     147
     148  // actual processing
     149  uInt len = m.nelements() ;
     150  // loop over channels
     151  for ( uInt i = 0 ; i < len ; i++ ) {
     152    // masks for spectrum_ and spectrumNoMask_ are not needed since
     153    // it is initialized as True for all channels and those are constant
     154    if ( *m_wp ) {
     155      // mask is True
     156      // add v * totalWeight to spectrum_
     157      // add totalWeight to weightSum_
     158      // increment n_
     159      *sp_wp += *v_wp * totalWeight ;
     160      *wgt_wp += totalWeight ;
     161      *n_wp += 1 ;
     162    }
     163    else {
     164      // mask is False
     165      // add v * totalWeightNoMask to spectrumNoMask_
     166      // add totalWeightNoMask to weightSumNoMask_
     167      // increment nNoMask_
     168      *spN_wp += *v_wp * totalWeightNoMask ;
     169      *wgtN_wp += totalWeightNoMask ;
     170      *nN_wp += 1 ;
     171    }
     172    sp_wp++ ;
     173    wgt_wp++ ;
     174    n_wp++ ;
     175    spN_wp++ ;
     176    wgtN_wp++ ;
     177    nN_wp++ ;
     178    v_wp++ ;
     179    m_wp++ ;
     180  }
     181
     182  // free storage
     183  spectrum_.putArrayStorage( sp_p, spD ) ;
     184  weightSum_.putArrayStorage( wgt_p, wgtD ) ;
     185  spectrumNoMask_.putArrayStorage( spN_p, spND ) ;
     186  weightSumNoMask_.putArrayStorage( wgtN_p, wgtND ) ;
     187  n_.putArrayStorage( n_p, nD ) ;
     188  nNoMask_.putArrayStorage( nN_p, nND ) ;
     189
     190  v.freeStorage( v_p, vD ) ;
     191  m.freeStorage( m_p, mD ) ;
     192}
     193
    114194Float RowAccumulator::getTotalWeight(const MaskedArray<Float>& data,
    115195                                     const Vector<Float>& tsys,

  • trunk/src/RowAccumulator.h

    r2141 r2580  
    109109                     const casa::Double time,
    110110                     const casa::Bool inverseMask);
     111  void doAddSpectrum2(const casa::Vector<casa::Float>& v,
     112                      const casa::Vector<casa::Bool>& m,
     113                      const casa::Vector<casa::Float>& tsys,
     114                      const casa::Double interval,
     115                      const casa::Double time);
    111116  casa::Float getTotalWeight(const casa::MaskedArray<casa::Float>& data,
    112117                             const casa::Vector<casa::Float>& tsys,

  • trunk/src/SConscript

  • trunk/src/STFitter.cpp

    r2455 r2580  
    379379
    380380  // Fit
    381   Vector<Float> sigma(x_.nelements());
    382   sigma = 1.0;
     381//   Vector<Float> sigma(x_.nelements());
     382//   sigma = 1.0;
    383383
    384384  parameters_.resize();
    385   parameters_ = fitter.fit(x_, y_, sigma, &m_);
     385//   parameters_ = fitter.fit(x_, y_, sigma, &m_);
     386  parameters_ = fitter.fit(x_, y_, &m_); 
    386387  if ( !fitter.converged() ) {
    387388     return false;
    … …  
    396397  chisquared_ = fitter.getChi2();
    397398
    398   residual_.resize();
    399   residual_ =  y_;
    400   fitter.residual(residual_,x_);
     399//   residual_.resize();
     400//   residual_ =  y_;
     401//   fitter.residual(residual_,x_);
    401402  // use fitter.residual(model=True) to get the model
    402403  thefit_.resize(x_.nelements());
    403404  fitter.residual(thefit_,x_,True);
     405  // residual = data - model
     406  residual_.resize(x_.nelements());
     407  residual_ = y_ - thefit_ ;
    404408  return true;
    405409}
    … …  
    420424
    421425  // Fit
    422   Vector<Float> sigma(x_.nelements());
    423   sigma = 1.0;
     426//   Vector<Float> sigma(x_.nelements());
     427//   sigma = 1.0;
    424428
    425429  parameters_.resize();
    426   parameters_ = fitter.fit(x_, y_, sigma, &m_);
     430//   parameters_ = fitter.fit(x_, y_, sigma, &m_);
     431  parameters_ = fitter.fit(x_, y_, &m_);
    427432  std::vector<float> ps;
    428433  parameters_.tovector(ps);
    … …  
    434439  chisquared_ = fitter.getChi2();
    435440
    436   residual_.resize();
    437   residual_ =  y_;
    438   fitter.residual(residual_,x_);
     441//   residual_.resize();
     442//   residual_ =  y_;
     443//   fitter.residual(residual_,x_);
    439444  // use fitter.residual(model=True) to get the model
    440445  thefit_.resize(x_.nelements());
    441446  fitter.residual(thefit_,x_,True);
     447  // residual = data - model
     448  residual_.resize(x_.nelements());
     449  residual_ = y_ - thefit_ ;
    442450  return true;
    443451}

  • trunk/src/STLineFinder.cpp

    r2425 r2580  
    884884//
    885885
    886 STLineFinder::STLineFinder() throw() : edge(0,0)
     886STLineFinder::STLineFinder() throw() : edge(0,0), err("spurious")
    887887{
    888888  useScantable = true;
    … …  
    10501050                     // because all previous values may be corrupted by the
    10511051                     // presence of spectral lines
     1052
    10521053  while (true) {
    10531054     // a buffer for new lines found at this iteration
    … …  
    10621063         first_pass=False;
    10631064         if (!new_lines.size())
    1064               throw AipsError("spurious"); // nothing new - use the same
    1065                                            // code as for a real exception
     1065//               throw AipsError("spurious"); // nothing new - use the same
     1066//                                            // code as for a real exception
     1067              throw err; // nothing new - use the same
     1068                         // code as for a real exception
    10661069     }
    10671070     catch(const AipsError &ae) {

  • trunk/src/STLineFinder.h

    r2012 r2580  
    4848#include "ScantableWrapper.h"
    4949#include "Scantable.h"
     50#include "STFitter.h"
    5051
    5152namespace asap {
    … …  
    161162                   const casa::Float &in_noise_box=-1.,
    162163                   const casa::Bool &in_median = casa::False) throw();
     164
     165   void setDetailedOptions( const casa::Int &order=9 ) ;
    163166
    164167   // set the scan to work with (in_scan parameter)
    … …  
    266269   // scantable (default = true)
    267270   bool useScantable;
     271
     272   // shared object for nominal throw
     273   casa::AipsError err ;
    268274};
    269275

  • trunk/src/STMath.cpp

    r2479 r2580  
    5353#include "STMath.h"
    5454#include "STSelector.h"
     55#include "Accelerator.h"
     56#include "STIdxIter.h"
    5557
    5658using namespace casa;
    … …  
    8082                 const std::string& avmode)
    8183{
     84//    double t0, t1 ;
     85//    t0 = mathutil::gettimeofday_sec() ;
     86
    8287  LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
    8388  if ( avmode == "SCAN" && in.size() != 1 )
    … …  
    128133  const Table& baset = in[0]->table();
    129134
    130   Block<String> cols(3);
     135  RowAccumulator acc(wtype);
     136  Vector<Bool> cmask(mask);
     137  acc.setUserMask(cmask);
     138//   ROTableRow row(tout);
     139  ROArrayColumn<Float> specCol, tsysCol;
     140  ROArrayColumn<uChar> flagCol;
     141  ROScalarColumn<Double> mjdCol, intCol;
     142  ROScalarColumn<Int> scanIDCol;
     143
     144  //Vector<uInt> rowstodelete;
     145  Block<uInt> rowstodelB( in[0]->nrow() ) ;
     146  uInt nrowdel = 0 ;
     147
     148//   Block<String> cols(3);
     149  vector<string> cols(3) ;
    131150  cols[0] = String("BEAMNO");
    132151  cols[1] = String("IFNO");
    … …  
    144163  }
    145164  uInt outrowCount = 0;
    146   TableIterator iter(baset, cols);
     165  // use STIdxIterExAcc instead of TableIterator
     166  STIdxIterExAcc iter( in[0], cols ) ;
     167//   double t2 = 0 ;
     168//   double t3 = 0 ;
     169//   double t4 = 0 ;
     170//   double t5 = 0 ;
     171//   TableIterator iter(baset, cols);
    147172//   int count = 0 ;
    148173  while (!iter.pastEnd()) {
    149     Table subt = iter.table();
     174    Vector<uInt> rows = iter.getRows( SHARE ) ;
     175    if ( rows.nelements() == 0 ) {
     176      iter.next() ;
     177      continue ;
     178    }
     179    Vector<uInt> current = iter.current() ;
     180    String srcname = iter.getSrcName() ;
     181    //Table subt = iter.table();
    150182    // copy the first row of this selection into the new table
    151183    tout.addRow();
    152     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     184//     t4 = mathutil::gettimeofday_sec() ;
     185    // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
     186    // overwritten in the following process
     187    copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
     188//     t5 += mathutil::gettimeofday_sec() - t4 ;
    153189    // re-index to 0
    154190    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    155191      scanColOut.put(outrowCount, uInt(0));
    156192    }
    157     ++outrowCount;
     193
    158194    // 2012/02/17 TN
    159195    // Since STGrid is implemented, average doesn't consider direction
    … …  
    195231//     }
    196232//     outrowCount += rowNum ;
    197     ++iter;
    198   }
    199   RowAccumulator acc(wtype);
    200   Vector<Bool> cmask(mask);
    201   acc.setUserMask(cmask);
    202   ROTableRow row(tout);
    203   ROArrayColumn<Float> specCol, tsysCol;
    204   ROArrayColumn<uChar> flagCol;
    205   ROScalarColumn<Double> mjdCol, intCol;
    206   ROScalarColumn<Int> scanIDCol;
    207 
    208   Vector<uInt> rowstodelete;
    209 
    210   for (uInt i=0; i < tout.nrow(); ++i) {
    211     for ( int j=0; j < int(in.size()); ++j ) {
     233
     234    // merge loop
     235    uInt i = outrowCount ;
     236    // in[0] is already selected by iterator
     237    specCol.attach(baset,"SPECTRA");
     238    flagCol.attach(baset,"FLAGTRA");
     239    tsysCol.attach(baset,"TSYS");
     240    intCol.attach(baset,"INTERVAL");
     241    mjdCol.attach(baset,"TIME");
     242    Vector<Float> spec,tsys;
     243    Vector<uChar> flag;
     244    Double inter,time;
     245
     246    for (uInt l = 0; l < rows.nelements(); ++l ) {
     247      uInt k = rows[l] ;
     248      flagCol.get(k, flag);
     249      Vector<Bool> bflag(flag.shape());
     250      convertArray(bflag, flag);
     251      /*                                                                                                   
     252        if ( allEQ(bflag, True) ) {                                                                         
     253        continue;//don't accumulate                                                                         
     254        }                                                                                                   
     255      */
     256      specCol.get(k, spec);
     257      tsysCol.get(k, tsys);
     258      intCol.get(k, inter);
     259      mjdCol.get(k, time);
     260      // spectrum has to be added last to enable weighting by the other values                             
     261//       t2 = mathutil::gettimeofday_sec() ;
     262      acc.add(spec, !bflag, tsys, inter, time);
     263//       t3 += mathutil::gettimeofday_sec() - t2 ;
     264     
     265    }
     266
     267
     268    // in[0] is already selected by TableIterator so that index is
     269    // started from 1
     270    for ( int j=1; j < int(in.size()); ++j ) {
    212271      const Table& tin = in[j]->table();
    213       const TableRecord& rec = row.get(i);
     272      //const TableRecord& rec = row.get(i);
    214273      ROScalarColumn<Double> tmp(tin, "TIME");
    215274      Double td;tmp.get(0,td);
    216       Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
    217                          && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
    218                          && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
     275
     276#if 1
     277      static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
     278      //uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
     279      uInt const values1[] = { current[1], current[0], current[2] };
     280      SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
     281      CustomTableExprNodeRep myNodeRep(tin, myPred);
     282      myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
     283      CustomTableExprNode myExpr(myNodeRep);
     284      Table basesubt = tin(myExpr);
     285#else
     286//       Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
     287//                          && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
     288//                          && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
     289      Table basesubt = tin( tin.col("BEAMNO") == current[0]
     290                         && tin.col("IFNO") == current[1]
     291                         && tin.col("POLNO") == current[2] );
     292#endif
    219293      Table subt;
    220294      if ( avmode == "SOURCE") {
    221         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
     295//         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
     296        subt = basesubt( basesubt.col("SRCNAME") == srcname );
     297
    222298      } else if (avmode == "SCAN") {
    223         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
    224                       && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
     299//         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
     300//                    && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
     301        subt = basesubt( basesubt.col("SRCNAME") == srcname
     302                      && basesubt.col("SCANNO") == current[4] );
    225303      } else {
    226304        subt = basesubt;
    … …  
    256334      intCol.attach(subt,"INTERVAL");
    257335      mjdCol.attach(subt,"TIME");
    258       Vector<Float> spec,tsys;
    259       Vector<uChar> flag;
    260       Double inter,time;
    261336      for (uInt k = 0; k < subt.nrow(); ++k ) {
    262337        flagCol.get(k, flag);
    … …  
    274349        mjdCol.get(k, time);
    275350        // spectrum has to be added last to enable weighting by the other values
     351//         t2 = mathutil::gettimeofday_sec() ;
    276352        acc.add(spec, !bflag, tsys, inter, time);
    277       }
    278 
     353//         t3 += mathutil::gettimeofday_sec() - t2 ;
     354      }
     355
     356    }
     357    const Vector<Bool>& msk = acc.getMask();
     358    if ( allEQ(msk, False) ) {
     359      rowstodelB[nrowdel] = i ;
     360      nrowdel++ ;
     361      continue;
     362    }
     363    //write out
     364    if (acc.state()) {
    279365      // If there exists a channel at which all the input spectra are masked,
    280366      // spec has 'nan' values for that channel and it may affect the following
    … …  
    283369      // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
    284370      acc.replaceNaN();
    285     }
    286     const Vector<Bool>& msk = acc.getMask();
    287     if ( allEQ(msk, False) ) {
    288       uint n = rowstodelete.nelements();
    289       rowstodelete.resize(n+1, True);
    290       rowstodelete[n] = i;
    291       continue;
    292     }
    293     //write out
    294     if (acc.state()) {
     371
    295372      Vector<uChar> flg(msk.shape());
    296373      convertArray(flg, !msk);
    … …  
    316393    }
    317394    acc.reset();
    318   }
    319 
    320   if (rowstodelete.nelements() > 0) {
     395
     396    // merge with while loop for preparing out table
     397    ++outrowCount;
     398//     ++iter ;
     399    iter.next() ;
     400  }
     401
     402  if ( nrowdel > 0 ) {
     403    Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ;
    321404    os << rowstodelete << LogIO::POST ;
    322405    tout.removeRow(rowstodelete);
    … …  
    325408    }
    326409  }
     410
     411//    t1 = mathutil::gettimeofday_sec() ;
     412//    cout << "elapsed time for average(): " << t1-t0 << " sec" << endl ;
     413//    cout << "   elapsed time for acc.add(): " << t3 << " sec" << endl ;
     414//    cout << "   elapsed time for copyRows(): " << t5 << " sec" << endl ;
     415
    327416  return out;
    328417}
    … …  
    38163905    vector<int> types ;
    38173906
     3907    // save original table selection
     3908    Table torg  = s->table_ ;
     3909
    38183910    // sky scan
    3819     STSelector sel = STSelector() ;
    3820     types.push_back( SrcType::SKY ) ;
    3821     sel.setTypes( types ) ;
    3822     s->setSelection( sel ) ;
    3823     vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
    3824     CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
    3825     s->unsetSelection() ;
    3826     sel.reset() ;
    3827     types.clear() ;
    3828 
     3911    bool insitu = insitu_ ;
     3912    insitu_ = false ;
     3913    // share calibration scans before average with out
     3914    CountedPtr<Scantable> out = getScantable( s, true ) ;
     3915    insitu_ = insitu ;
     3916    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
     3917    out->attach() ;
     3918    CountedPtr<Scantable> asky = averageWithinSession( out,
     3919                                                       masks,
     3920                                                       "TINT" ) ;
    38293921    // hot scan
    3830     types.push_back( SrcType::HOT ) ;
    3831     sel.setTypes( types ) ;
    3832     s->setSelection( sel ) ;
    3833     tmp.clear() ;
    3834     tmp.push_back( getScantable( s, false ) ) ;
    3835     CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
    3836     s->unsetSelection() ;
    3837     sel.reset() ;
    3838     types.clear() ;
    3839    
     3922    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
     3923    out->attach() ;
     3924    CountedPtr<Scantable> ahot = averageWithinSession( out,
     3925                                                       masks,
     3926                                                       "TINT" ) ;
    38403927    // cold scan
    38413928    CountedPtr<Scantable> acold ;
    3842 //     types.push_back( SrcType::COLD ) ;
    3843 //     sel.setTypes( types ) ;
    3844 //     s->setSelection( sel ) ;
    3845 //     tmp.clear() ;
    3846 //     tmp.push_back( getScantable( s, false ) ) ;
    3847 //     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
    3848 //     s->unsetSelection() ;
    3849 //     sel.reset() ;
    3850 //     types.clear() ;
     3929//     out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
     3930//     out->attach() ;
     3931//     CountedPtr<Scantable> acold = averageWithinSession( out,
     3932//                                                         masks,
     3933//                                                         "TINT" ) ;
    38513934
    38523935    // off scan
    3853     types.push_back( SrcType::PSOFF ) ;
    3854     sel.setTypes( types ) ;
    3855     s->setSelection( sel ) ;
    3856     tmp.clear() ;
    3857     tmp.push_back( getScantable( s, false ) ) ;
    3858     CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
    3859     s->unsetSelection() ;
    3860     sel.reset() ;
    3861     types.clear() ;
     3936    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
     3937    out->attach() ;
     3938    CountedPtr<Scantable> aoff = averageWithinSession( out,
     3939                                                       masks,
     3940                                                       "TINT" ) ;
    38623941   
    38633942    // on scan
    3864     bool insitu = insitu_ ;
    3865     insitu_ = false ;
    3866     CountedPtr<Scantable> out = getScantable( s, true ) ;
    3867     insitu_ = insitu ;
    3868     types.push_back( SrcType::PSON ) ;
    3869     sel.setTypes( types ) ;
    3870     s->setSelection( sel ) ;
    3871     TableCopy::copyRows( out->table(), s->table() ) ;
    3872     s->unsetSelection() ;
    3873     sel.reset() ;
    3874     types.clear() ;
     3943    s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
     3944    s->attach() ;
     3945    out->table_ = out->originalTable_ ;
     3946    out->attach() ;
     3947    out->table().addRow( s->nrow() ) ;
     3948    copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
    38753949   
    38763950    // process each on scan
    3877     ArrayColumn<Float> tsysCol ;
    3878     tsysCol.attach( out->table(), "TSYS" ) ;
    3879     for ( int i = 0 ; i < out->nrow() ; i++ ) {
    3880       vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
    3881       out->setSpectrum( sp, i ) ;
    3882       string reftime = out->getTime( i ) ;
    3883       vector<int> ii( 1, out->getIF( i ) ) ;
    3884       vector<int> ib( 1, out->getBeam( i ) ) ;
    3885       vector<int> ip( 1, out->getPol( i ) ) ;
    3886       sel.setIFs( ii ) ;
    3887       sel.setBeams( ib ) ;
    3888       sel.setPolarizations( ip ) ;
    3889       asky->setSelection( sel ) ;   
    3890       vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
    3891       const Vector<Float> Vtsys( sptsys ) ;
    3892       tsysCol.put( i, Vtsys ) ;
     3951    STSelector sel ;
     3952    vector<string> cols( 3 ) ;
     3953    cols[0] = "BEAMNO" ;
     3954    cols[1] = "POLNO" ;
     3955    cols[2] = "IFNO" ;
     3956    STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
     3957    while ( !iter->pastEnd() ) {
     3958      Vector<uInt> ids = iter->current() ;
     3959      stringstream ss ;
     3960      ss << "SELECT FROM $1 WHERE "
     3961         << "BEAMNO==" << ids[0] << "&&"
     3962         << "POLNO==" << ids[1] << "&&"
     3963         << "IFNO==" << ids[2] ;
     3964      //cout << "TaQL string: " << ss.str() << endl ;
     3965      sel.setTaQL( ss.str() ) ;
     3966      aoff->setSelection( sel ) ;
     3967      ahot->setSelection( sel ) ;
     3968      asky->setSelection( sel ) ;
     3969      Vector<uInt> rows = iter->getRows( SHARE ) ;
     3970      // out should be an exact copy of s except that SPECTRA column is empty
     3971      calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
     3972      aoff->unsetSelection() ;
     3973      ahot->unsetSelection() ;
    38933974      asky->unsetSelection() ;
    38943975      sel.reset() ;
    3895     }
     3976      iter->next() ;
     3977    }
     3978    delete iter ;
     3979    s->table_ = torg ;
     3980    s->attach() ;
    38963981
    38973982    // flux unit
    … …  
    39103995  }
    39113996  else {
     3997//     double t0, t1 ;
     3998//     t0 = mathutil::gettimeofday_sec() ;
    39123999    vector<bool> masks = s->getMask( 0 ) ;
    3913    
     4000
     4001    // save original table selection
     4002    Table torg = s->table_ ;
     4003
    39144004    // off scan
    3915     STSelector sel = STSelector() ;
    3916     vector<int> types ;
    3917     types.push_back( SrcType::PSOFF ) ;
    3918     sel.setTypes( types ) ;
    3919     s->setSelection( sel ) ;
    39204005    // TODO 2010/01/08 TN
    39214006    // Grouping by time should be needed before averaging.
    39224007    // Each group must have own unique SCANNO (should be renumbered).
    39234008    // See PIPELINE/SDCalibration.py
    3924     CountedPtr<Scantable> soff = getScantable( s, false ) ;
    3925     Table ttab = soff->table() ;
    3926     ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
    3927     uInt nrow = timeCol.nrow() ;
    3928     Vector<Double> timeSep( nrow - 1 ) ;
    3929     for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
    3930       timeSep[i] = timeCol(i+1) - timeCol(i) ;
    3931     }
    3932     ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
    3933     Vector<Double> interval = intervalCol.getColumn() ;
    3934     interval /= 86400.0 ;
    3935     ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
    3936     vector<uInt> glist ;
    3937     for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
    3938       double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
    3939       //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
    3940       if ( gap > 1.1 ) {
    3941         glist.push_back( i ) ;
    3942       }
    3943     }
    3944     Vector<uInt> gaplist( glist ) ;
    3945     //cout << "gaplist = " << gaplist << endl ;
    3946     uInt newid = 0 ;
    3947     for ( uInt i = 0 ; i < nrow ; i++ ) {
    3948       scanCol.put( i, newid ) ;
    3949       if ( i == gaplist[newid] ) {
    3950         newid++ ;
    3951       }
    3952     }
    3953     //cout << "new scancol = " << scanCol.getColumn() << endl ;
    3954     vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
    3955     CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
    3956     //cout << "aoff.nrow = " << aoff->nrow() << endl ;
    3957     s->unsetSelection() ;
    3958     sel.reset() ;
    3959     types.clear() ;
    3960    
    3961     // on scan
    39624009    bool insitu = insitu_ ;
    39634010    insitu_ = false ;
     4011    // share off scan before average with out
    39644012    CountedPtr<Scantable> out = getScantable( s, true ) ;
     4013    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
     4014    out->attach() ;
    39654015    insitu_ = insitu ;
    3966     types.push_back( SrcType::PSON ) ;
    3967     sel.setTypes( types ) ;
    3968     s->setSelection( sel ) ;
    3969     TableCopy::copyRows( out->table(), s->table() ) ;
    3970     s->unsetSelection() ;
    3971     sel.reset() ;
    3972     types.clear() ;
    3973    
     4016    CountedPtr<Scantable> aoff = averageWithinSession( out,
     4017                                                       masks,
     4018                                                       "TINT" ) ;
     4019
     4020    // on scan
     4021//     t0 = mathutil::gettimeofday_sec() ;
     4022    s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
     4023    s->attach() ;
     4024    out->table_ = out->originalTable_ ;
     4025    out->attach() ;
     4026    out->table().addRow( s->nrow() ) ;
     4027    copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
     4028//     t1 = mathutil::gettimeofday_sec() ;
     4029//     cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
     4030
    39744031    // process each on scan
    3975     ArrayColumn<Float> tsysCol ;
    3976     tsysCol.attach( out->table(), "TSYS" ) ;
    3977     for ( int i = 0 ; i < out->nrow() ; i++ ) {
    3978       vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
    3979       out->setSpectrum( sp, i ) ;
    3980     }
     4032//     t0 = mathutil::gettimeofday_sec() ;
     4033
     4034    // using STIdxIterAcc
     4035    vector<string> cols( 3 ) ;
     4036    cols[0] = "BEAMNO" ;
     4037    cols[1] = "POLNO" ;
     4038    cols[2] = "IFNO" ;
     4039    STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
     4040    STSelector sel ;
     4041    while ( !iter->pastEnd() ) {
     4042      Vector<uInt> ids = iter->current() ;
     4043      stringstream ss ;
     4044      ss << "SELECT FROM $1 WHERE "
     4045         << "BEAMNO==" << ids[0] << "&&"
     4046         << "POLNO==" << ids[1] << "&&"
     4047         << "IFNO==" << ids[2] ;
     4048      //cout << "TaQL string: " << ss.str() << endl ;
     4049      sel.setTaQL( ss.str() ) ;
     4050      aoff->setSelection( sel ) ;
     4051      Vector<uInt> rows = iter->getRows( SHARE ) ;
     4052      // out should be an exact copy of s except that SPECTRA column is empty
     4053      calibrateALMA( out, s, aoff, rows ) ;
     4054      aoff->unsetSelection() ;
     4055      sel.reset() ;
     4056      iter->next() ;
     4057    }
     4058    delete iter ;
     4059    s->table_ = torg ;
     4060    s->attach() ;
     4061
     4062//     t1 = mathutil::gettimeofday_sec() ;
     4063//     cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
    39814064
    39824065    // flux unit
    … …  
    40144097  vector<bool> masks = s->getMask( 0 ) ;
    40154098  CountedPtr<Scantable> ssig, sref ;
    4016   CountedPtr<Scantable> out ;
     4099  //CountedPtr<Scantable> out ;
     4100  bool insitu = insitu_ ;
     4101  insitu_ = False ;
     4102  CountedPtr<Scantable> out = getScantable( s, true ) ;
     4103  insitu_ = insitu ;
    40174104
    40184105  if ( antname.find( "APEX" ) != string::npos ) {
    40194106    // APEX calibration
    40204107    // sky scan
    4021     STSelector sel = STSelector() ;
    4022     types.push_back( SrcType::FLOSKY ) ;
    4023     sel.setTypes( types ) ;
    4024     s->setSelection( sel ) ;
    4025     vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
    4026     CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
    4027     s->unsetSelection() ;
    4028     sel.reset() ;
    4029     types.clear() ;
    4030     types.push_back( SrcType::FHISKY ) ;
    4031     sel.setTypes( types ) ;
    4032     s->setSelection( sel ) ;
    4033     tmp.clear() ;
    4034     tmp.push_back( getScantable( s, false ) ) ;
    4035     CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
    4036     s->unsetSelection() ;
    4037     sel.reset() ;
    4038     types.clear() ;
     4108    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
     4109    out->attach() ;
     4110    CountedPtr<Scantable> askylo = averageWithinSession( out,
     4111                                                         masks,
     4112                                                         "TINT" ) ;
     4113    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
     4114    out->attach() ;
     4115    CountedPtr<Scantable> askyhi = averageWithinSession( out,
     4116                                                         masks,
     4117                                                         "TINT" ) ;
    40394118   
    40404119    // hot scan
    4041     types.push_back( SrcType::FLOHOT ) ;
    4042     sel.setTypes( types ) ;
    4043     s->setSelection( sel ) ;
    4044     tmp.clear() ;
    4045     tmp.push_back( getScantable( s, false ) ) ;
    4046     CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
    4047     s->unsetSelection() ;
    4048     sel.reset() ;
    4049     types.clear() ;
    4050     types.push_back( SrcType::FHIHOT ) ;
    4051     sel.setTypes( types ) ;
    4052     s->setSelection( sel ) ;
    4053     tmp.clear() ;
    4054     tmp.push_back( getScantable( s, false ) ) ;
    4055     CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
    4056     s->unsetSelection() ;
    4057     sel.reset() ;
    4058     types.clear() ;
     4120    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
     4121    out->attach() ;
     4122    CountedPtr<Scantable> ahotlo = averageWithinSession( out,
     4123                                                         masks,
     4124                                                         "TINT" ) ;
     4125    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
     4126    out->attach() ;
     4127    CountedPtr<Scantable> ahothi = averageWithinSession( out,
     4128                                                         masks,
     4129                                                         "TINT" ) ;
    40594130   
    40604131    // cold scan
    40614132    CountedPtr<Scantable> acoldlo, acoldhi ;
    4062 //     types.push_back( SrcType::FLOCOLD ) ;
    4063 //     sel.setTypes( types ) ;
    4064 //     s->setSelection( sel ) ;
    4065 //     tmp.clear() ;
    4066 //     tmp.push_back( getScantable( s, false ) ) ;
    4067 //     CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
    4068 //     s->unsetSelection() ;
    4069 //     sel.reset() ;
    4070 //     types.clear() ;
    4071 //     types.push_back( SrcType::FHICOLD ) ;
    4072 //     sel.setTypes( types ) ;
    4073 //     s->setSelection( sel ) ;
    4074 //     tmp.clear() ;
    4075 //     tmp.push_back( getScantable( s, false ) ) ;
    4076 //     CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
    4077 //     s->unsetSelection() ;
    4078 //     sel.reset() ;
    4079 //     types.clear() ;
     4133//     out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
     4134//     out->attach() ;
     4135//     CountedPtr<Scantable> acoldlo = averageWithinSession( out,
     4136//                                                           masks,
     4137//                                                           "TINT" ) ;
     4138//     out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
     4139//     out->attach() ;
     4140//     CountedPtr<Scantable> acoldhi = averageWithinSession( out,
     4141//                                                           masks,
     4142//                                                           "TINT" ) ;
    40804143
    40814144    // ref scan
    4082     bool insitu = insitu_ ;
    40834145    insitu_ = false ;
    40844146    sref = getScantable( s, true ) ;
     4147    CountedPtr<Scantable> rref = getScantable( s, true ) ;
    40854148    insitu_ = insitu ;
    4086     types.push_back( SrcType::FSLO ) ;
    4087     sel.setTypes( types ) ;
    4088     s->setSelection( sel ) ;
    4089     TableCopy::copyRows( sref->table(), s->table() ) ;
    4090     s->unsetSelection() ;
    4091     sel.reset() ;
    4092     types.clear() ;
     4149    rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
     4150    rref->attach() ;
     4151    copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
    40934152   
    40944153    // sig scan
    40954154    insitu_ = false ;
    40964155    ssig = getScantable( s, true ) ;
     4156    CountedPtr<Scantable> rsig = getScantable( s, true ) ;
    40974157    insitu_ = insitu ;
    4098     types.push_back( SrcType::FSHI ) ;
    4099     sel.setTypes( types ) ;
    4100     s->setSelection( sel ) ;
    4101     TableCopy::copyRows( ssig->table(), s->table() ) ;
    4102     s->unsetSelection() ;
    4103     sel.reset() ; 
    4104     types.clear() ;
     4158    rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
     4159    rsig->attach() ;
     4160    copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
    41054161         
    41064162    if ( apexcalmode == 0 ) {
    4107       // APEX fs data without off scan
    4108       // process each sig and ref scan
    4109       ArrayColumn<Float> tsysCollo ;
    4110       tsysCollo.attach( ssig->table(), "TSYS" ) ;
    4111       ArrayColumn<Float> tsysColhi ;
    4112       tsysColhi.attach( sref->table(), "TSYS" ) ;
    4113       for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    4114         vector< CountedPtr<Scantable> > sky( 2 ) ;
    4115         sky[0] = askylo ;
    4116         sky[1] = askyhi ;
    4117         vector< CountedPtr<Scantable> > hot( 2 ) ;
    4118         hot[0] = ahotlo ;
    4119         hot[1] = ahothi ;
    4120         vector< CountedPtr<Scantable> > cold( 2 ) ;
    4121         //cold[0] = acoldlo ;
    4122         //cold[1] = acoldhi ;
    4123         vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
    4124         ssig->setSpectrum( sp, i ) ;
    4125         string reftime = ssig->getTime( i ) ;
    4126         vector<int> ii( 1, ssig->getIF( i ) ) ;
    4127         vector<int> ib( 1, ssig->getBeam( i ) ) ;
    4128         vector<int> ip( 1, ssig->getPol( i ) ) ;
    4129         sel.setIFs( ii ) ;
    4130         sel.setBeams( ib ) ;
    4131         sel.setPolarizations( ip ) ;
    4132         askylo->setSelection( sel ) ;
    4133         vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
    4134         const Vector<Float> Vtsyslo( sptsys ) ;
    4135         tsysCollo.put( i, Vtsyslo ) ;
    4136         askylo->unsetSelection() ;
     4163      // using STIdxIterAcc
     4164      vector<string> cols( 3 ) ;
     4165      cols[0] = "BEAMNO" ;
     4166      cols[1] = "POLNO" ;
     4167      cols[2] = "IFNO" ;
     4168      STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
     4169      STSelector sel ;
     4170      vector< CountedPtr<Scantable> > on( 2 ) ;
     4171      on[0] = rsig ;
     4172      on[1] = rref ;
     4173      vector< CountedPtr<Scantable> > sky( 2 ) ;
     4174      sky[0] = askylo ;
     4175      sky[1] = askyhi ;
     4176      vector< CountedPtr<Scantable> > hot( 2 ) ;
     4177      hot[0] = ahotlo ;
     4178      hot[1] = ahothi ;
     4179      vector< CountedPtr<Scantable> > cold( 2 ) ;
     4180      while ( !iter->pastEnd() ) {
     4181        Vector<uInt> ids = iter->current() ;
     4182        stringstream ss ;
     4183        ss << "SELECT FROM $1 WHERE "
     4184           << "BEAMNO==" << ids[0] << "&&"
     4185           << "POLNO==" << ids[1] << "&&"
     4186           << "IFNO==" << ids[2] ;
     4187        //cout << "TaQL string: " << ss.str() << endl ;
     4188        sel.setTaQL( ss.str() ) ;
     4189        sky[0]->setSelection( sel ) ;
     4190        sky[1]->setSelection( sel ) ;
     4191        hot[0]->setSelection( sel ) ;
     4192        hot[1]->setSelection( sel ) ;
     4193        Vector<uInt> rows = iter->getRows( SHARE ) ;
     4194        calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
     4195        sky[0]->unsetSelection() ;
     4196        sky[1]->unsetSelection() ;
     4197        hot[0]->unsetSelection() ;
     4198        hot[1]->unsetSelection() ;
    41374199        sel.reset() ;
    4138         sky[0] = askyhi ;
    4139         sky[1] = askylo ;
    4140         hot[0] = ahothi ;
    4141         hot[1] = ahotlo ;
    4142         cold[0] = acoldhi ;
    4143         cold[1] = acoldlo ;
    4144         sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
    4145         sref->setSpectrum( sp, i ) ;
    4146         reftime = sref->getTime( i ) ;
    4147         ii[0] = sref->getIF( i )  ;
    4148         ib[0] = sref->getBeam( i ) ;
    4149         ip[0] = sref->getPol( i ) ;
    4150         sel.setIFs( ii ) ;
    4151         sel.setBeams( ib ) ;
    4152         sel.setPolarizations( ip ) ;
    4153         askyhi->setSelection( sel ) ;   
    4154         sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
    4155         const Vector<Float> Vtsyshi( sptsys ) ;
    4156         tsysColhi.put( i, Vtsyshi ) ;
    4157         askyhi->unsetSelection() ;
    4158         sel.reset() ;
    4159       }
     4200        iter->next() ;
     4201      }
     4202      delete iter ;
     4203
    41604204    }
    41614205    else if ( apexcalmode == 1 ) {
    41624206      // APEX fs data with off scan
    41634207      // off scan
    4164       types.push_back( SrcType::FLOOFF ) ;
    4165       sel.setTypes( types ) ;
    4166       s->setSelection( sel ) ;
    4167       tmp.clear() ;
    4168       tmp.push_back( getScantable( s, false ) ) ;
    4169       CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
    4170       s->unsetSelection() ;
    4171       sel.reset() ;
    4172       types.clear() ;
    4173       types.push_back( SrcType::FHIOFF ) ;
    4174       sel.setTypes( types ) ;
    4175       s->setSelection( sel ) ;
    4176       tmp.clear() ;
    4177       tmp.push_back( getScantable( s, false ) ) ;
    4178       CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
    4179       s->unsetSelection() ;
    4180       sel.reset() ;
    4181       types.clear() ;
     4208      out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
     4209      out->attach() ;
     4210      CountedPtr<Scantable> aofflo = averageWithinSession( out,
     4211                                                           masks,
     4212                                                           "TINT" ) ;
     4213      out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
     4214      out->attach() ;
     4215      CountedPtr<Scantable> aoffhi = averageWithinSession( out,
     4216                                                           masks,
     4217                                                           "TINT" ) ;
    41824218     
    41834219      // process each sig and ref scan
    4184       ArrayColumn<Float> tsysCollo ;
    4185       tsysCollo.attach( ssig->table(), "TSYS" ) ;
    4186       ArrayColumn<Float> tsysColhi ;
    4187       tsysColhi.attach( sref->table(), "TSYS" ) ;
    4188       for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    4189         vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
    4190         ssig->setSpectrum( sp, i ) ;
    4191         sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
    4192         string reftime = ssig->getTime( i ) ;
    4193         vector<int> ii( 1, ssig->getIF( i ) ) ;
    4194         vector<int> ib( 1, ssig->getBeam( i ) ) ;
    4195         vector<int> ip( 1, ssig->getPol( i ) ) ;
    4196         sel.setIFs( ii ) ;
    4197         sel.setBeams( ib ) ;
    4198         sel.setPolarizations( ip ) ;
     4220//       STSelector sel ;
     4221      vector<string> cols( 3 ) ;
     4222      cols[0] = "BEAMNO" ;
     4223      cols[1] = "POLNO" ;
     4224      cols[2] = "IFNO" ;
     4225      STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
     4226      STSelector sel ;
     4227      while ( !iter->pastEnd() ) {
     4228        Vector<uInt> ids = iter->current() ;
     4229        stringstream ss ;
     4230        ss << "SELECT FROM $1 WHERE "
     4231           << "BEAMNO==" << ids[0] << "&&"
     4232           << "POLNO==" << ids[1] << "&&"
     4233           << "IFNO==" << ids[2] ;
     4234        //cout << "TaQL string: " << ss.str() << endl ;
     4235        sel.setTaQL( ss.str() ) ;
     4236        aofflo->setSelection( sel ) ;
     4237        ahotlo->setSelection( sel ) ;
    41994238        askylo->setSelection( sel ) ;
    4200         vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
    4201         const Vector<Float> Vtsyslo( sptsys ) ;
    4202         tsysCollo.put( i, Vtsyslo ) ;
     4239        Vector<uInt> rows = iter->getRows( SHARE ) ;
     4240        calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ;
     4241        aofflo->unsetSelection() ;
     4242        ahotlo->unsetSelection() ;
    42034243        askylo->unsetSelection() ;
    42044244        sel.reset() ;
    4205         sref->setSpectrum( sp, i ) ;
    4206         reftime = sref->getTime( i ) ;
    4207         ii[0] = sref->getIF( i )  ;
    4208         ib[0] = sref->getBeam( i ) ;
    4209         ip[0] = sref->getPol( i ) ;
    4210         sel.setIFs( ii ) ;
    4211         sel.setBeams( ib ) ;
    4212         sel.setPolarizations( ip ) ;
    4213         askyhi->setSelection( sel ) ;   
    4214         sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
    4215         const Vector<Float> Vtsyshi( sptsys ) ;
    4216         tsysColhi.put( i, Vtsyshi ) ;
     4245        iter->next() ;
     4246      }
     4247      delete iter ;
     4248      iter = new STIdxIterAcc( sref, cols ) ;
     4249      while ( !iter->pastEnd() ) {
     4250        Vector<uInt> ids = iter->current() ;
     4251        stringstream ss ;
     4252        ss << "SELECT FROM $1 WHERE "
     4253           << "BEAMNO==" << ids[0] << "&&"
     4254           << "POLNO==" << ids[1] << "&&"
     4255           << "IFNO==" << ids[2] ;
     4256        //cout << "TaQL string: " << ss.str() << endl ;
     4257        sel.setTaQL( ss.str() ) ;
     4258        aoffhi->setSelection( sel ) ;
     4259        ahothi->setSelection( sel ) ;
     4260        askyhi->setSelection( sel ) ;
     4261        Vector<uInt> rows = iter->getRows( SHARE ) ;
     4262        calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ;
     4263        aoffhi->unsetSelection() ;
     4264        ahothi->unsetSelection() ;
    42174265        askyhi->unsetSelection() ;
    42184266        sel.reset() ;
    4219       }
     4267        iter->next() ;
     4268      }
     4269      delete iter ;
    42204270    }
    42214271  }
    … …  
    42234273    // non-APEX fs data
    42244274    // sky scan
     4275    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
     4276    out->attach() ;
     4277    CountedPtr<Scantable> asky = averageWithinSession( out,
     4278                                                       masks,
     4279                                                       "TINT" ) ;
    42254280    STSelector sel = STSelector() ;
    4226     types.push_back( SrcType::SKY ) ;
    4227     sel.setTypes( types ) ;
    4228     s->setSelection( sel ) ;
    4229     vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
    4230     CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
    4231     s->unsetSelection() ;
    4232     sel.reset() ;
    4233     types.clear() ;
    4234    
     4281
    42354282    // hot scan
    4236     types.push_back( SrcType::HOT ) ;
    4237     sel.setTypes( types ) ;
    4238     s->setSelection( sel ) ;
    4239     tmp.clear() ;
    4240     tmp.push_back( getScantable( s, false ) ) ;
    4241     CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
    4242     s->unsetSelection() ;
    4243     sel.reset() ;
    4244     types.clear() ;
     4283    out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
     4284    out->attach() ;
     4285    CountedPtr<Scantable> ahot = averageWithinSession( out,
     4286                                                       masks,
     4287                                                       "TINT" ) ;
    42454288
    42464289    // cold scan
    42474290    CountedPtr<Scantable> acold ;
    4248 //     types.push_back( SrcType::COLD ) ;
    4249 //     sel.setTypes( types ) ;
    4250 //     s->setSelection( sel ) ;
    4251 //     tmp.clear() ;
    4252 //     tmp.push_back( getScantable( s, false ) ) ;
    4253 //     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
    4254 //     s->unsetSelection() ;
    4255 //     sel.reset() ;
    4256 //     types.clear() ;
     4291//     out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
     4292//     out->attach() ;
     4293//     CountedPtr<Scantable> acold = averageWithinSession( out,
     4294//                                                         masks,
     4295//                                                         "TINT" ) ;
    42574296   
    42584297    // ref scan
    … …  
    42604299    insitu_ = false ;
    42614300    sref = getScantable( s, true ) ;
     4301    CountedPtr<Scantable> rref = getScantable( s, true ) ;
    42624302    insitu_ = insitu ;
    4263     types.push_back( SrcType::FSOFF ) ;
    4264     sel.setTypes( types ) ;
    4265     s->setSelection( sel ) ;
    4266     TableCopy::copyRows( sref->table(), s->table() ) ;
    4267     s->unsetSelection() ;
    4268     sel.reset() ;
    4269     types.clear() ;
     4303    rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
     4304    rref->attach() ;
     4305    copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
    42704306   
    42714307    // sig scan
    42724308    insitu_ = false ;
    42734309    ssig = getScantable( s, true ) ;
     4310    CountedPtr<Scantable> rsig = getScantable( s, true ) ;
    42744311    insitu_ = insitu ;
    4275     types.push_back( SrcType::FSON ) ;
    4276     sel.setTypes( types ) ;
    4277     s->setSelection( sel ) ;
    4278     TableCopy::copyRows( ssig->table(), s->table() ) ;
    4279     s->unsetSelection() ;
    4280     sel.reset() ;
    4281     types.clear() ;
     4312    rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
     4313    rsig->attach() ;
     4314    copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
    42824315
    42834316    // process each sig and ref scan
    4284     ArrayColumn<Float> tsysColsig ;
    4285     tsysColsig.attach( ssig->table(), "TSYS" ) ;
    4286     ArrayColumn<Float> tsysColref ;
    4287     tsysColref.attach( ssig->table(), "TSYS" ) ;
    4288     for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    4289       vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
    4290       ssig->setSpectrum( sp, i ) ;
    4291       string reftime = ssig->getTime( i ) ;
    4292       vector<int> ii( 1, ssig->getIF( i ) ) ;
    4293       vector<int> ib( 1, ssig->getBeam( i ) ) ;
    4294       vector<int> ip( 1, ssig->getPol( i ) ) ;
    4295       sel.setIFs( ii ) ;
    4296       sel.setBeams( ib ) ;
    4297       sel.setPolarizations( ip ) ;
     4317    vector<string> cols( 3 ) ;
     4318    cols[0] = "BEAMNO" ;
     4319    cols[1] = "POLNO" ;
     4320    cols[2] = "IFNO" ;
     4321    STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
     4322    while ( !iter->pastEnd() ) {
     4323      Vector<uInt> ids = iter->current() ;
     4324      stringstream ss ;
     4325      ss << "SELECT FROM $1 WHERE "
     4326         << "BEAMNO==" << ids[0] << "&&"
     4327         << "POLNO==" << ids[1] << "&&"
     4328         << "IFNO==" << ids[2] ;
     4329      //cout << "TaQL string: " << ss.str() << endl ;
     4330      sel.setTaQL( ss.str() ) ;
     4331      ahot->setSelection( sel ) ;
    42984332      asky->setSelection( sel ) ;
    4299       vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
    4300       const Vector<Float> Vtsys( sptsys ) ;
    4301       tsysColsig.put( i, Vtsys ) ;
     4333      Vector<uInt> rows = iter->getRows( SHARE ) ;
     4334      // out should be an exact copy of s except that SPECTRA column is empty
     4335      calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
     4336      ahot->unsetSelection() ;
    43024337      asky->unsetSelection() ;
    43034338      sel.reset() ;
    4304       sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
    4305       sref->setSpectrum( sp, i ) ;
    4306       tsysColref.put( i, Vtsys ) ;
    4307     }
     4339      iter->next() ;
     4340    }
     4341    delete iter ;
    43084342  }
    43094343
    … …  
    43114345  Table sigtab = ssig->table() ;
    43124346  Table reftab = sref->table() ;
    4313   ScalarColumn<uInt> sigifnoCol ;
    4314   ScalarColumn<uInt> refifnoCol ;
    4315   ScalarColumn<uInt> sigfidCol ;
    43164347  ScalarColumn<uInt> reffidCol ;
    43174348  Int nchan = (Int)ssig->nchan() ;
    4318   sigifnoCol.attach( sigtab, "IFNO" ) ;
    4319   refifnoCol.attach( reftab, "IFNO" ) ;
    4320   sigfidCol.attach( sigtab, "FREQ_ID" ) ;
    43214349  reffidCol.attach( reftab, "FREQ_ID" ) ;
    4322   Vector<uInt> sfids( sigfidCol.getColumn() ) ;
    4323   Vector<uInt> rfids( reffidCol.getColumn() ) ;
     4350  Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
     4351  Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
    43244352  vector<uInt> sfids_unique ;
    43254353  vector<uInt> rfids_unique ;
    … …  
    44004428}
    44014429
    4402 vector<float> STMath::getSpectrumFromTime( string reftime,
    4403                                            CountedPtr<Scantable>& s,
     4430Vector<Float> STMath::getSpectrumFromTime( double reftime,
     4431                                           const Vector<Double> &timeVec,
     4432                                           const vector<int> &idx,
     4433                                           const Matrix<Float>& spectra,
    44044434                                           string mode )
    44054435{
    44064436  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
    4407   vector<float> sp ;
    4408 
    4409   if ( s->nrow() == 0 ) {
     4437  Vector<Float> sp ;
     4438  uInt ncol = spectra.ncolumn() ;
     4439
     4440  if ( ncol == 0 ) {
    44104441    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
    44114442    return sp ;
    44124443  }
    4413   else if ( s->nrow() == 1 ) {
     4444  else if ( ncol == 1 ) {
    44144445    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
    4415     return s->getSpectrum( 0 ) ;
     4446    sp.reference( spectra.column( 0 ) ) ;
     4447    return sp ;
    44164448  }
    44174449  else {
    4418     vector<int> idx = getRowIdFromTime( reftime, s ) ;
    44194450    if ( mode == "before" ) {
    44204451      int id = -1 ;
    … …  
    44274458      }
    44284459      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4429       sp = s->getSpectrum( id ) ;
     4460      sp.reference( spectra.column( id ) ) ;
    44304461    }
    44314462    else if ( mode == "after" ) {
    … …  
    44394470      }
    44404471      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4441       sp = s->getSpectrum( id ) ;
     4472      sp.reference( spectra.column( id ) ) ;
    44424473    }
    44434474    else if ( mode == "nearest" ) {
    … …  
    44534484      }
    44544485      else {
    4455         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4456         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4457         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4458         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4459         double tref = getMJD( reftime ) ;
     4486        double t0 = timeVec[idx[0]] ;
     4487        double t1 = timeVec[idx[1]] ;
     4488        double tref = reftime ;
    44604489        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
    44614490          id = idx[1] ;
    … …  
    44664495      }
    44674496      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4468       sp = s->getSpectrum( id ) ;     
     4497      sp.reference( spectra.column( id ) ) ;
    44694498    }
    44704499    else if ( mode == "linear" ) {
    … …  
    44744503        int id = idx[1] ;
    44754504        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4476         sp = s->getSpectrum( id ) ;
     4505        sp.reference( spectra.column( id ) ) ;
    44774506      }
    44784507      else if ( idx[1] == -1 ) {
    … …  
    44814510        int id = idx[0] ;
    44824511        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4483         sp = s->getSpectrum( id ) ;
     4512        sp.reference( spectra.column( id ) ) ;
    44844513      }
    44854514      else if ( idx[0] == idx[1] ) {
    … …  
    44884517        int id = idx[0] ;
    44894518        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4490         sp = s->getSpectrum( id ) ;
     4519        sp.reference( spectra.column( id ) ) ;
    44914520      }
    44924521      else {
    44934522        // do interpolation
    44944523        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    4495         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4496         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4497         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4498         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4499         double tref = getMJD( reftime ) ;
    4500         vector<float> sp0 = s->getSpectrum( idx[0] ) ;
    4501         vector<float> sp1 = s->getSpectrum( idx[1] ) ;
    4502         for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
    4503           float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
    4504           sp.push_back( v ) ;
     4524        double t0 = timeVec[idx[0]] ;
     4525        double t1 = timeVec[idx[1]] ;
     4526        double tref = reftime ;
     4527        sp = spectra.column( idx[0] ).copy() ;
     4528        Vector<Float> sp1( spectra.column( idx[1] ) ) ;
     4529        double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
     4530        for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
     4531          sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
    45054532        }
    45064533      }
    … …  
    45134540}
    45144541
    4515 double STMath::getMJD( string strtime )
    4516 {
    4517   if ( strtime.find("/") == string::npos ) {
    4518     // MJD time string
    4519     return atof( strtime.c_str() ) ;
    4520   }
    4521   else {
    4522     // string in YYYY/MM/DD/HH:MM:SS format
    4523     uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
    4524     uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
    4525     uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
    4526     uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
    4527     uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
    4528     uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
    4529     Time t( year, month, day, hour, minute, sec ) ;
    4530     return t.modifiedJulianDay() ;
    4531   }
    4532 }
    4533 
    4534 vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
    4535 {
    4536   double reft = getMJD( reftime ) ;
     4542vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t )
     4543{
     4544//   double reft = reftime ;
    45374545  double dtmin = 1.0e100 ;
    45384546  double dtmax = -1.0e100 ;
    4539   vector<double> dt ;
     4547//   vector<double> dt ;
    45404548  int just_before = -1 ;
    45414549  int just_after = -1 ;
    4542   for ( int i = 0 ; i < s->nrow() ; i++ ) {
    4543     dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
    4544   }
     4550  Vector<Double> dt = t - reftime ;
    45454551  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
    45464552    if ( dt[i] > 0.0 ) {
    … …  
    45684574  }
    45694575
    4570   vector<int> v ;
    4571   v.push_back( just_before ) ;
    4572   v.push_back( just_after ) ;
     4576  vector<int> v(2) ;
     4577  v[0] = just_before ;
     4578  v[1] = just_after ;
    45734579
    45744580  return v ;
    45754581}
    45764582
    4577 vector<float> STMath::getTcalFromTime( string reftime,
    4578                                        CountedPtr<Scantable>& s,
     4583Vector<Float> STMath::getTcalFromTime( double reftime,
     4584                                       const Vector<Double> &timeVec,
     4585                                       const vector<int> &idx,
     4586                                       const CountedPtr<Scantable>& s,
    45794587                                       string mode )
    45804588{
    45814589  LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
    4582   vector<float> tcal ;
    45834590  STTcal tcalTable = s->tcal() ;
    45844591  String time ;
    … …  
    45864593  if ( s->nrow() == 0 ) {
    45874594    os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
    4588     return tcal ;
     4595    return tcalval ;
    45894596  }
    45904597  else if ( s->nrow() == 1 ) {
    … …  
    45924599    //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    45934600    tcalTable.getEntry( time, tcalval, tcalid ) ;
    4594     tcalval.tovector( tcal ) ;
    4595     return tcal ;
     4601    return tcalval ;
    45964602  }
    45974603  else {
    4598     vector<int> idx = getRowIdFromTime( reftime, s ) ;
    45994604    if ( mode == "before" ) {
    46004605      int id = -1 ;
    … …  
    46094614      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46104615      tcalTable.getEntry( time, tcalval, tcalid ) ;
    4611       tcalval.tovector( tcal ) ;
    46124616    }
    46134617    else if ( mode == "after" ) {
    … …  
    46234627      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46244628      tcalTable.getEntry( time, tcalval, tcalid ) ;
    4625       tcalval.tovector( tcal ) ;
    46264629    }
    46274630    else if ( mode == "nearest" ) {
    … …  
    46374640      }
    46384641      else {
    4639         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4640         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4641         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4642         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4643         double tref = getMJD( reftime ) ;
    4644         if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4642        double t0 = timeVec[idx[0]] ;
     4643        double t1 = timeVec[idx[1]] ;
     4644        if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
    46454645          id = idx[1] ;
    46464646        }
    … …  
    46524652      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46534653      tcalTable.getEntry( time, tcalval, tcalid ) ;
    4654       tcalval.tovector( tcal ) ;
    46554654    }
    46564655    else if ( mode == "linear" ) {
    … …  
    46624661        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46634662        tcalTable.getEntry( time, tcalval, tcalid ) ;
    4664         tcalval.tovector( tcal ) ;
    46654663      }
    46664664      else if ( idx[1] == -1 ) {
    … …  
    46714669        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46724670        tcalTable.getEntry( time, tcalval, tcalid ) ;
    4673         tcalval.tovector( tcal ) ;
    46744671      }
    46754672      else if ( idx[0] == idx[1] ) {
    … …  
    46804677        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    46814678        tcalTable.getEntry( time, tcalval, tcalid ) ;
    4682         tcalval.tovector( tcal ) ;
    46834679      }
    46844680      else {
    46854681        // do interpolation
    46864682        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    4687         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4688         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4689         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4690         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4691         double tref = getMJD( reftime ) ;
    4692         vector<float> tcal0 ;
    4693         vector<float> tcal1 ;
     4683        double t0 = timeVec[idx[0]] ;
     4684        double t1 = timeVec[idx[1]] ;
     4685        Vector<Float> tcal0 ;
    46944686        uInt tcalid0 = s->getTcalId( idx[0] ) ;
    46954687        uInt tcalid1 = s->getTcalId( idx[1] ) ;
    4696         tcalTable.getEntry( time, tcalval, tcalid0 ) ;
    4697         tcalval.tovector( tcal0 ) ;
     4688        tcalTable.getEntry( time, tcal0, tcalid0 ) ;
    46984689        tcalTable.getEntry( time, tcalval, tcalid1 ) ;
    4699         tcalval.tovector( tcal1 ) ;       
     4690        double tfactor = (reftime - t0) / (t1 - t0) ;
    47004691        for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
    4701           float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
    4702           tcal.push_back( v ) ;
     4692          tcalval[i] = ( tcalval[i] - tcal0[i] ) * tfactor + tcal0[i] ;
    47034693        }
    47044694      }
    … …  
    47074697      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
    47084698    }
    4709     return tcal ;
    4710   }
    4711 }
    4712 
    4713 vector<float> STMath::getTsysFromTime( string reftime,
    4714                                        CountedPtr<Scantable>& s,
     4699    return tcalval ;
     4700  }
     4701}
     4702
     4703Vector<Float> STMath::getTsysFromTime( double reftime,
     4704                                       const Vector<Double> &timeVec,
     4705                                       const vector<int> &idx,
     4706                                       const CountedPtr<Scantable> &s,
    47154707                                       string mode )
    47164708{
    … …  
    47184710  ArrayColumn<Float> tsysCol ;
    47194711  tsysCol.attach( s->table(), "TSYS" ) ;
    4720   vector<float> tsys ;
    4721   String time ;
    47224712  Vector<Float> tsysval ;
    47234713  if ( s->nrow() == 0 ) {
    47244714    os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
    4725     return tsys ;
     4715    return tsysval ;
    47264716  }
    47274717  else if ( s->nrow() == 1 ) {
    47284718    //os << "use row " << 0 << LogIO::POST ;
    47294719    tsysval = tsysCol( 0 ) ;
    4730     tsysval.tovector( tsys ) ;
    4731     return tsys ;
     4720    return tsysval ;
    47324721  }
    47334722  else {
    4734     vector<int> idx = getRowIdFromTime( reftime, s ) ;
    47354723    if ( mode == "before" ) {
    47364724      int id = -1 ;
    … …  
    47444732      //os << "use row " << id << LogIO::POST ;
    47454733      tsysval = tsysCol( id ) ;
    4746       tsysval.tovector( tsys ) ;
    47474734    }
    47484735    else if ( mode == "after" ) {
    … …  
    47574744      //os << "use row " << id << LogIO::POST ;
    47584745      tsysval = tsysCol( id ) ;
    4759       tsysval.tovector( tsys ) ;
    47604746    }
    47614747    else if ( mode == "nearest" ) {
    … …  
    47714757      }
    47724758      else {
    4773         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4774         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4775         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4776         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4777         double tref = getMJD( reftime ) ;
    4778         if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4759        double t0 = timeVec[idx[0]] ;
     4760        double t1 = timeVec[idx[1]] ;
     4761        if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
    47794762          id = idx[1] ;
    47804763        }
    … …  
    47854768      //os << "use row " << id << LogIO::POST ;
    47864769      tsysval = tsysCol( id ) ;
    4787       tsysval.tovector( tsys ) ;
    47884770    }
    47894771    else if ( mode == "linear" ) {
    … …  
    47944776        //os << "use row " << id << LogIO::POST ;
    47954777        tsysval = tsysCol( id ) ;
    4796         tsysval.tovector( tsys ) ;
    47974778      }
    47984779      else if ( idx[1] == -1 ) {
    … …  
    48024783        //os << "use row " << id << LogIO::POST ;
    48034784        tsysval = tsysCol( id ) ;
    4804         tsysval.tovector( tsys ) ;
    48054785      }
    48064786      else if ( idx[0] == idx[1] ) {
    … …  
    48104790        //os << "use row " << id << LogIO::POST ;
    48114791        tsysval = tsysCol( id ) ;
    4812         tsysval.tovector( tsys ) ;
    48134792      }
    48144793      else {
    48154794        // do interpolation
    48164795        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    4817         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4818         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4819         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4820         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4821         double tref = getMJD( reftime ) ;
    4822         vector<float> tsys0 ;
    4823         vector<float> tsys1 ;
    4824         tsysval = tsysCol( idx[0] ) ;
    4825         tsysval.tovector( tsys0 ) ;
     4796        double t0 = timeVec[idx[0]] ;
     4797        double t1 = timeVec[idx[1]] ;
     4798        Vector<Float> tsys0 ;
     4799        tsys0 = tsysCol( idx[0] ) ;
    48264800        tsysval = tsysCol( idx[1] ) ;
    4827         tsysval.tovector( tsys1 ) ;       
     4801        double tfactor = (reftime - t0) / (t1 - t0) ;
    48284802        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
    4829           float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
    4830           tsys.push_back( v ) ;
     4803          tsysval[i] = ( tsysval[i] - tsys0[i] ) * tfactor + tsys0[i] ;
    48314804        }
    48324805      }
    … …  
    48354808      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
    48364809    }
    4837     return tsys ;
    4838   }
    4839 }
    4840 
    4841 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
    4842                                             CountedPtr<Scantable>& off,
    4843                                             CountedPtr<Scantable>& sky,
    4844                                             CountedPtr<Scantable>& hot,
    4845                                             CountedPtr<Scantable>& cold,
    4846                                             int index,
    4847                                             string antname )
    4848 {
    4849   (void) cold; //currently unused
    4850   string reftime = on->getTime( index ) ;
    4851   vector<int> ii( 1, on->getIF( index ) ) ;
    4852   vector<int> ib( 1, on->getBeam( index ) ) ;
    4853   vector<int> ip( 1, on->getPol( index ) ) ;
    4854   vector<int> ic( 1, on->getScan( index ) ) ;
    4855   STSelector sel = STSelector() ;
    4856   sel.setIFs( ii ) ;
    4857   sel.setBeams( ib ) ;
    4858   sel.setPolarizations( ip ) ;
    4859   sky->setSelection( sel ) ;
    4860   hot->setSelection( sel ) ;
    4861   //cold->setSelection( sel ) ;
    4862   off->setSelection( sel ) ;
    4863   vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
    4864   vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
    4865   //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
    4866   vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
    4867   vector<float> spec = on->getSpectrum( index ) ;
    4868   vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
    4869   vector<float> sp( tcal.size() ) ;
    4870   if ( antname.find( "APEX" ) != string::npos ) {
     4810    return tsysval ;
     4811  }
     4812}
     4813
     4814void STMath::calibrateCW( CountedPtr<Scantable> &out,
     4815                          const CountedPtr<Scantable>& on,
     4816                          const CountedPtr<Scantable>& off,
     4817                          const CountedPtr<Scantable>& sky,
     4818                          const CountedPtr<Scantable>& hot,
     4819                          const CountedPtr<Scantable>& cold,
     4820                          const Vector<uInt> &rows,
     4821                          const String &antname )
     4822{
     4823  // 2012/05/22 TN
     4824  // Assume that out has empty SPECTRA column
     4825
     4826  // if rows is empty, just return
     4827  if ( rows.nelements() == 0 )
     4828    return ;
     4829  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
     4830  Vector<Double> timeOff = timeCol.getColumn() ;
     4831  timeCol.attach( sky->table(), "TIME" ) ;
     4832  Vector<Double> timeSky = timeCol.getColumn() ;
     4833  timeCol.attach( hot->table(), "TIME" ) ;
     4834  Vector<Double> timeHot = timeCol.getColumn() ;
     4835  timeCol.attach( on->table(), "TIME" ) ;
     4836  ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
     4837  Matrix<Float> offspectra = arrayFloatCol.getColumn() ;
     4838  arrayFloatCol.attach( sky->table(), "SPECTRA" ) ;
     4839  Matrix<Float> skyspectra = arrayFloatCol.getColumn() ;
     4840  arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
     4841  Matrix<Float> hotspectra = arrayFloatCol.getColumn() ;
     4842  unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
     4843  // I know that the data is contiguous
     4844  const uInt *p = rows.data() ;
     4845  vector<int> ids( 2 ) ;
     4846  Block<uInt> flagchan( spsize ) ;
     4847  uInt nflag = 0 ;
     4848  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     4849    double reftime = timeCol.asdouble(*p) ;
     4850    ids = getRowIdFromTime( reftime, timeOff ) ;
     4851    Vector<Float> spoff = getSpectrumFromTime( reftime, timeOff, ids, offspectra, "linear" ) ;
     4852    ids = getRowIdFromTime( reftime, timeSky ) ;
     4853    Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
     4854    Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     4855    Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     4856    ids = getRowIdFromTime( reftime, timeHot ) ;
     4857    Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra, "linear" ) ;
     4858    Vector<Float> spec = on->specCol_( *p ) ;
     4859    if ( antname.find( "APEX" ) != String::npos ) {
     4860      // using gain array
     4861      for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4862        if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
     4863          spec[j] = 0.0 ;
     4864          flagchan[nflag++] = j ;
     4865        }
     4866        else {
     4867          spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
     4868            * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     4869        }
     4870      }
     4871    }
     4872    else {
     4873      // Chopper-Wheel calibration (Ulich & Haas 1976)
     4874      for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4875        if ( (sphot[j]-spsky[j]) == 0.0 ) {
     4876          spec[j] = 0.0 ;
     4877          flagchan[nflag++] = j ;
     4878        }
     4879        else {
     4880          spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     4881        }
     4882      }
     4883    }
     4884    out->specCol_.put( *p, spec ) ;
     4885    out->tsysCol_.put( *p, tsys ) ;
     4886    if ( nflag > 0 ) {
     4887      Vector<uChar> fl = out->flagsCol_( *p ) ;
     4888      for ( unsigned int j = 0 ; j < nflag ; j++ ) {
     4889        fl[flagchan[j]] = (uChar)True ;
     4890      }
     4891      out->flagsCol_.put( *p, fl ) ;
     4892    }
     4893    nflag = 0 ;
     4894    p++ ;
     4895  }
     4896}
     4897
     4898void STMath::calibrateALMA( CountedPtr<Scantable>& out,
     4899                            const CountedPtr<Scantable>& on,
     4900                            const CountedPtr<Scantable>& off,
     4901                            const Vector<uInt>& rows )
     4902{
     4903  // 2012/05/22 TN
     4904  // Assume that out has empty SPECTRA column
     4905
     4906  // if rows is empty, just return
     4907  if ( rows.nelements() == 0 )
     4908    return ;
     4909  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
     4910  Vector<Double> timeVec = timeCol.getColumn() ;
     4911  timeCol.attach( on->table(), "TIME" ) ;
     4912  ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
     4913  Matrix<Float> offspectra = arrayFloatCol.getColumn() ;
     4914  unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
     4915  // I know that the data is contiguous
     4916  const uInt *p = rows.data() ;
     4917  vector<int> ids( 2 ) ;
     4918  Block<uInt> flagchan( spsize ) ;
     4919  uInt nflag = 0 ;
     4920  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     4921    double reftime = timeCol.asdouble(*p) ;
     4922    ids = getRowIdFromTime( reftime, timeVec ) ;
     4923    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, ids, offspectra, "linear" ) ;
     4924    //Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
     4925    Vector<Float> spec = on->specCol_( *p ) ;
     4926    Vector<Float> tsys = on->tsysCol_( *p ) ;
     4927    // ALMA Calibration
     4928    //
     4929    // Ta* = Tsys * ( ON - OFF ) / OFF
     4930    //
     4931    // 2010/01/07 Takeshi Nakazato
     4932    unsigned int tsyssize = tsys.nelements() ;
     4933    for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     4934      if ( spoff[j] == 0.0 ) {
     4935        spec[j] = 0.0 ;
     4936        flagchan[nflag++] = j ;
     4937      }
     4938      else {
     4939        spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
     4940      }
     4941      if ( tsyssize == spsize )
     4942        spec[j] *= tsys[j] ;
     4943      else
     4944        spec[j] *= tsys[0] ;
     4945    }
     4946    out->specCol_.put( *p, spec ) ;
     4947    if ( nflag > 0 ) {
     4948      Vector<uChar> fl = out->flagsCol_( *p ) ;
     4949      for ( unsigned int j = 0 ; j < nflag ; j++ ) {
     4950        fl[flagchan[j]] = (uChar)True ;
     4951      }
     4952      out->flagsCol_.put( *p, fl ) ;
     4953    }
     4954    nflag = 0 ;
     4955    p++ ;
     4956  }
     4957}
     4958
     4959void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
     4960                              CountedPtr<Scantable> &ref,
     4961                              const vector< CountedPtr<Scantable> >& on,
     4962                              const vector< CountedPtr<Scantable> >& sky,
     4963                              const vector< CountedPtr<Scantable> >& hot,
     4964                              const vector< CountedPtr<Scantable> >& cold,
     4965                              const Vector<uInt> &rows )
     4966{
     4967  // if rows is empty, just return
     4968  if ( rows.nelements() == 0 )
     4969    return ;
     4970  ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
     4971  Vector<Double> timeSkyS = timeCol.getColumn() ;
     4972  timeCol.attach( sky[1]->table(), "TIME" ) ;
     4973  Vector<Double> timeSkyR = timeCol.getColumn() ;
     4974  timeCol.attach( hot[0]->table(), "TIME" ) ;
     4975  Vector<Double> timeHotS = timeCol.getColumn() ;
     4976  timeCol.attach( hot[1]->table(), "TIME" ) ;
     4977  Vector<Double> timeHotR = timeCol.getColumn() ;
     4978  timeCol.attach( sig->table(), "TIME" ) ;
     4979  ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
     4980  ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
     4981  Matrix<Float> skyspectraS = arrayFloatCol.getColumn() ;
     4982  arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
     4983  Matrix<Float> skyspectraR = arrayFloatCol.getColumn() ;
     4984  arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
     4985  Matrix<Float> hotspectraS = arrayFloatCol.getColumn() ;
     4986  arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
     4987  Matrix<Float> hotspectraR = arrayFloatCol.getColumn() ;
     4988  unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
     4989  Vector<Float> spec( spsize ) ;
     4990  // I know that the data is contiguous
     4991  const uInt *p = rows.data() ;
     4992  vector<int> ids( 2 ) ;
     4993  Block<uInt> flagchan( spsize ) ;
     4994  uInt nflag = 0 ;
     4995  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     4996    double reftime = timeCol.asdouble(*p) ;
     4997    ids = getRowIdFromTime( reftime, timeSkyS ) ;
     4998    Vector<Float> spskyS = getSpectrumFromTime( reftime, timeSkyS, ids, skyspectraS, "linear" ) ;
     4999    Vector<Float> tcalS = getTcalFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ;
     5000    Vector<Float> tsysS = getTsysFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ;
     5001    ids = getRowIdFromTime( reftime, timeHotS ) ;
     5002    Vector<Float> sphotS = getSpectrumFromTime( reftime, timeHotS, ids, hotspectraS ) ;
     5003    reftime = timeCol2.asdouble(*p) ;
     5004    ids = getRowIdFromTime( reftime, timeSkyR ) ;
     5005    Vector<Float> spskyR = getSpectrumFromTime( reftime, timeSkyR, ids, skyspectraR, "linear" ) ;
     5006    Vector<Float> tcalR = getTcalFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ;
     5007    Vector<Float> tsysR = getTsysFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ;
     5008    ids = getRowIdFromTime( reftime, timeHotR ) ;
     5009    Vector<Float> sphotR = getSpectrumFromTime( reftime, timeHotR, ids, hotspectraR ) ;
     5010    Vector<Float> spsig = on[0]->specCol_( *p ) ;
     5011    Vector<Float> spref = on[1]->specCol_( *p ) ;
     5012    for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     5013      if ( (sphotS[j]-spskyS[j]) == 0.0 || (sphotR[j]-spskyR[j]) == 0.0 ) {
     5014        spec[j] = 0.0 ;
     5015        flagchan[nflag++] = j ;
     5016      }
     5017      else {
     5018        spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
     5019          - tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
     5020      }
     5021    }
     5022    sig->specCol_.put( *p, spec ) ;
     5023    sig->tsysCol_.put( *p, tsysS ) ;
     5024    spec *= (Float)-1.0 ;
     5025    ref->specCol_.put( *p, spec ) ;
     5026    ref->tsysCol_.put( *p, tsysR ) ;   
     5027    if ( nflag > 0 ) {
     5028      Vector<uChar> flsig = sig->flagsCol_( *p ) ;
     5029      Vector<uChar> flref = ref->flagsCol_( *p ) ;
     5030      for ( unsigned int j = 0 ; j < nflag ; j++ ) {
     5031        flsig[flagchan[j]] = (uChar)True ;
     5032        flref[flagchan[j]] = (uChar)True ;
     5033      }
     5034      sig->flagsCol_.put( *p, flsig ) ;
     5035      ref->flagsCol_.put( *p, flref ) ;
     5036    }
     5037    nflag = 0 ;
     5038    p++ ;
     5039  }
     5040}
     5041
     5042void STMath::calibrateFS( CountedPtr<Scantable> &sig,
     5043                          CountedPtr<Scantable> &ref,
     5044                          const CountedPtr<Scantable>& rsig,
     5045                          const CountedPtr<Scantable>& rref,
     5046                          const CountedPtr<Scantable>& sky,
     5047                          const CountedPtr<Scantable>& hot,
     5048                          const CountedPtr<Scantable>& cold,
     5049                          const Vector<uInt> &rows )
     5050{
     5051  // if rows is empty, just return
     5052  if ( rows.nelements() == 0 )
     5053    return ;
     5054  ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
     5055  Vector<Double> timeSky = timeCol.getColumn() ;
     5056  timeCol.attach( hot->table(), "TIME" ) ;
     5057  Vector<Double> timeHot = timeCol.getColumn() ;
     5058  timeCol.attach( sig->table(), "TIME" ) ;
     5059  ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
     5060  ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
     5061  Matrix<Float> skyspectra = arrayFloatCol.getColumn() ;
     5062  arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
     5063  Matrix<Float> hotspectra = arrayFloatCol.getColumn() ;
     5064  unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
     5065  Vector<Float> spec( spsize ) ;
     5066  // I know that the data is contiguous
     5067  const uInt *p = rows.data() ;
     5068  vector<int> ids( 2 ) ;
     5069  Block<uInt> flagchan( spsize ) ;
     5070  uInt nflag = 0 ;
     5071  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     5072    double reftime = timeCol.asdouble(*p) ;
     5073    ids = getRowIdFromTime( reftime, timeSky ) ;
     5074    Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
     5075    Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     5076    Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     5077    ids = getRowIdFromTime( reftime, timeHot ) ;
     5078    Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ;
     5079    Vector<Float> spsig = rsig->specCol_( *p ) ;
     5080    Vector<Float> spref = rref->specCol_( *p ) ;
    48715081    // using gain array
    4872     for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4873       float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
    4874         * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
    4875       sp[j] = v ;
    4876     }
    4877   }
    4878   else {
    4879     // Chopper-Wheel calibration (Ulich & Haas 1976)
    4880     for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4881       float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
    4882       sp[j] = v ;
    4883     }
    4884   }
    4885   sel.reset() ;
    4886   sky->unsetSelection() ;
    4887   hot->unsetSelection() ;
    4888   //cold->unsetSelection() ;
    4889   off->unsetSelection() ;
    4890 
    4891   return sp ;
    4892 }
    4893 
    4894 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
    4895                                             CountedPtr<Scantable>& off,
    4896                                             int index )
    4897 {
    4898   string reftime = on->getTime( index ) ;
    4899   vector<int> ii( 1, on->getIF( index ) ) ;
    4900   vector<int> ib( 1, on->getBeam( index ) ) ;
    4901   vector<int> ip( 1, on->getPol( index ) ) ;
    4902   vector<int> ic( 1, on->getScan( index ) ) ;
    4903   STSelector sel = STSelector() ;
    4904   sel.setIFs( ii ) ;
    4905   sel.setBeams( ib ) ;
    4906   sel.setPolarizations( ip ) ;
    4907   off->setSelection( sel ) ;
    4908   vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
    4909   vector<float> spec = on->getSpectrum( index ) ;
    4910   //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
    4911   //vector<float> tsys = on->getTsysVec( index ) ;
    4912   ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
    4913   Vector<Float> tsys = tsysCol( index ) ;
    4914   vector<float> sp( spec.size() ) ;
    4915   // ALMA Calibration
    4916   //
    4917   // Ta* = Tsys * ( ON - OFF ) / OFF
    4918   //
    4919   // 2010/01/07 Takeshi Nakazato
    4920   unsigned int tsyssize = tsys.nelements() ;
    4921   unsigned int spsize = sp.size() ;
    4922   for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
    4923     float tscale = 0.0 ;
    4924     if ( tsyssize == spsize )
    4925       tscale = tsys[j] ;
    4926     else
    4927       tscale = tsys[0] ;
    4928     float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
    4929     sp[j] = v ;
    4930   }
    4931   sel.reset() ;
    4932   off->unsetSelection() ;
    4933 
    4934   return sp ;
    4935 }
    4936 
    4937 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
    4938                                               CountedPtr<Scantable>& ref,
    4939                                               CountedPtr<Scantable>& sky,
    4940                                               CountedPtr<Scantable>& hot,
    4941                                               CountedPtr<Scantable>& cold,
    4942                                               int index )
    4943 {
    4944   (void) cold; //currently unused
    4945   string reftime = sig->getTime( index ) ;
    4946   vector<int> ii( 1, sig->getIF( index ) ) ;
    4947   vector<int> ib( 1, sig->getBeam( index ) ) ;
    4948   vector<int> ip( 1, sig->getPol( index ) ) ;
    4949   vector<int> ic( 1, sig->getScan( index ) ) ;
    4950   STSelector sel = STSelector() ;
    4951   sel.setIFs( ii ) ;
    4952   sel.setBeams( ib ) ;
    4953   sel.setPolarizations( ip ) ;
    4954   sky->setSelection( sel ) ;
    4955   hot->setSelection( sel ) ;
    4956   //cold->setSelection( sel ) ;
    4957   vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
    4958   vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
    4959   //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
    4960   vector<float> spref = ref->getSpectrum( index ) ;
    4961   vector<float> spsig = sig->getSpectrum( index ) ;
    4962   vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
    4963   vector<float> sp( tcal.size() ) ;
    4964   for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4965     float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
    4966     sp[j] = v ;
    4967   }
    4968   sel.reset() ;
    4969   sky->unsetSelection() ;
    4970   hot->unsetSelection() ;
    4971   //cold->unsetSelection() ;
    4972 
    4973   return sp ;
    4974 }
    4975 
    4976 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
    4977                                               CountedPtr<Scantable>& ref,
    4978                                               vector< CountedPtr<Scantable> >& sky,
    4979                                               vector< CountedPtr<Scantable> >& hot,
    4980                                               vector< CountedPtr<Scantable> >& cold,
    4981                                               int index )
    4982 {
    4983   (void) cold; //currently unused
    4984   string reftime = sig->getTime( index ) ;
    4985   vector<int> ii( 1, sig->getIF( index ) ) ;
    4986   vector<int> ib( 1, sig->getBeam( index ) ) ;
    4987   vector<int> ip( 1, sig->getPol( index ) ) ;
    4988   vector<int> ic( 1, sig->getScan( index ) ) ;
    4989   STSelector sel = STSelector() ;
    4990   sel.setIFs( ii ) ;
    4991   sel.setBeams( ib ) ;
    4992   sel.setPolarizations( ip ) ;
    4993   sky[0]->setSelection( sel ) ;
    4994   hot[0]->setSelection( sel ) ;
    4995   //cold[0]->setSelection( sel ) ;
    4996   vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
    4997   vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
    4998   //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
    4999   vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
    5000   sel.reset() ;
    5001   ii[0] = ref->getIF( index ) ;
    5002   sel.setIFs( ii ) ;
    5003   sel.setBeams( ib ) ;
    5004   sel.setPolarizations( ip ) ;
    5005   sky[1]->setSelection( sel ) ;
    5006   hot[1]->setSelection( sel ) ;
    5007   //cold[1]->setSelection( sel ) ;
    5008   vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
    5009   vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
    5010   //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
    5011   vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 
    5012   vector<float> spref = ref->getSpectrum( index ) ;
    5013   vector<float> spsig = sig->getSpectrum( index ) ;
    5014   vector<float> sp( tcals.size() ) ;
    5015   for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
    5016     float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
    5017     sp[j] = v ;
    5018   }
    5019   sel.reset() ;
    5020   sky[0]->unsetSelection() ;
    5021   hot[0]->unsetSelection() ;
    5022   //cold[0]->unsetSelection() ;
    5023   sky[1]->unsetSelection() ;
    5024   hot[1]->unsetSelection() ;
    5025   //cold[1]->unsetSelection() ;
    5026 
    5027   return sp ;
    5028 }
     5082    for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     5083      if ( spref[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
     5084        spec[j] = 0.0 ;
     5085        flagchan[nflag++] = j ;
     5086      }
     5087      else {
     5088        spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
     5089          * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     5090      }
     5091    }
     5092    sig->specCol_.put( *p, spec ) ;
     5093    sig->tsysCol_.put( *p, tsys ) ;
     5094    if ( nflag > 0 ) {
     5095      Vector<uChar> fl = sig->flagsCol_( *p ) ;
     5096      for ( unsigned int j = 0 ; j < nflag ; j++ ) {
     5097        fl[flagchan[j]] = (uChar)True ;
     5098      }
     5099      sig->flagsCol_.put( *p, fl ) ;
     5100    }
     5101    nflag = 0 ;
     5102
     5103    reftime = timeCol2.asdouble(*p) ;
     5104    spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
     5105    tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     5106    tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
     5107    ids = getRowIdFromTime( reftime, timeHot ) ;
     5108    sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ;
     5109    // using gain array
     5110    for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     5111      if ( spsig[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
     5112        spec[j] = 0.0 ;
     5113        flagchan[nflag++] = j ;
     5114      }
     5115      else {
     5116        spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
     5117          * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     5118      }
     5119    }
     5120    ref->specCol_.put( *p, spec ) ;
     5121    ref->tsysCol_.put( *p, tsys ) ;   
     5122    if ( nflag > 0 ) {
     5123      Vector<uChar> fl = ref->flagsCol_( *p ) ;
     5124      for ( unsigned int j = 0 ; j < nflag ; j++ ) {
     5125        fl[flagchan[j]] = (uChar)True ;
     5126      }
     5127      ref->flagsCol_.put( *p, fl ) ;
     5128    }
     5129    nflag = 0 ;
     5130    p++ ;
     5131  }
     5132}
     5133
     5134void STMath::copyRows( Table &out,
     5135                       const Table &in,
     5136                       uInt startout,
     5137                       uInt startin,
     5138                       uInt nrow,
     5139                       Bool copySpectra,
     5140                       Bool copyFlagtra,
     5141                       Bool copyTsys )
     5142{
     5143  uInt nexclude = 0 ;
     5144  Block<String> excludeColsBlock( 3 ) ;
     5145  if ( !copySpectra ) {
     5146    excludeColsBlock[nexclude] = "SPECTRA" ;
     5147    nexclude++ ;
     5148  }
     5149  if ( !copyFlagtra ) {
     5150    excludeColsBlock[nexclude] = "FLAGTRA" ;
     5151    nexclude++ ;
     5152  }
     5153  if ( !copyTsys ) {
     5154    excludeColsBlock[nexclude] = "TSYS" ;
     5155    nexclude++ ;
     5156  }
     5157  //  if ( nexclude < 3 ) {
     5158  //    excludeCols.resize( nexclude, True ) ;
     5159  //  }
     5160  Vector<String> excludeCols( IPosition(1,nexclude),
     5161                              excludeColsBlock.storage(),
     5162                              SHARE ) ;
     5163//   cout << "excludeCols=" << excludeCols << endl ;
     5164  TableRow rowout( out, excludeCols, True ) ;
     5165  ROTableRow rowin( in, excludeCols, True ) ;
     5166  uInt rin = startin ;
     5167  uInt rout = startout ;
     5168  for ( uInt i = 0 ; i < nrow ; i++ ) {
     5169    rowin.get( rin ) ;
     5170    rowout.putMatchingFields( rout, rowin.record() ) ;
     5171    rin++ ;
     5172    rout++ ;
     5173  }
     5174}
     5175
     5176CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
     5177                                                    vector<bool> &mask,
     5178                                                    string weight )
     5179{
     5180  // prepare output table
     5181  bool insitu = insitu_ ;
     5182  insitu_ = false ;
     5183  CountedPtr<Scantable> a = getScantable( s, true ) ;
     5184  insitu_ = insitu ;
     5185  Table &atab = a->table() ;
     5186  ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
     5187
     5188  if ( s->nrow() == 0 )
     5189    return a ;
     5190
     5191  // setup RowAccumulator
     5192  WeightType wtype = stringToWeight( weight ) ;
     5193  RowAccumulator acc( wtype ) ;
     5194  Vector<Bool> cmask( mask ) ;
     5195  acc.setUserMask( cmask ) ;
     5196
     5197  vector<string> cols( 3 ) ;
     5198  cols[0] = "IFNO" ;
     5199  cols[1] = "POLNO" ;
     5200  cols[2] = "BEAMNO" ;
     5201  STIdxIterAcc iter( s, cols ) ;
     5202
     5203  Table ttab = s->table() ;
     5204  ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
     5205  Vector<Double> timeVec = timeCol->getColumn() ;
     5206  delete timeCol ;
     5207  Vector<Double> interval = s->integrCol_.getColumn() ;
     5208  uInt nrow = timeVec.nelements() ;
     5209  uInt outrow = 0 ;
     5210
     5211  while( !iter.pastEnd() ) {
     5212
     5213    Vector<uInt> rows = iter.getRows( SHARE ) ;
     5214
     5215    uInt nchan = s->nchan(s->getIF(rows[0])) ;
     5216    Vector<uChar> flag( nchan ) ;
     5217    Vector<Bool> bflag( nchan ) ;
     5218    Vector<Float> spec( nchan ) ;
     5219    Vector<Float> tsys( nchan ) ;
     5220
     5221    uInt len = rows.nelements() ;
     5222
     5223    Vector<Double> timeSep( len-1 ) ;
     5224    for ( uInt i = 0 ; i < len-1 ; i++ ) {
     5225      timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
     5226    }
     5227
     5228    uInt irow ;
     5229    uInt jrow ;
     5230    for ( uInt i = 0 ; i < len-1 ; i++ ) {
     5231      irow = rows[i] ;
     5232      jrow = rows[i+1] ;
     5233      // accumulate data
     5234      s->flagsCol_.get( irow, flag ) ;
     5235      convertArray( bflag, flag ) ;
     5236      s->specCol_.get( irow, spec ) ;
     5237      tsys.assign( s->tsysCol_( irow ) ) ;
     5238      if ( !allEQ(bflag,True) )
     5239        acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
     5240      double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
     5241      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
     5242      if ( gap > 1.1 ) {
     5243        //cout << "detected gap between " << i << " and " << i+1 << endl ;
     5244        // put data to output table
     5245        // reset RowAccumulator
     5246        if ( acc.state() ) {
     5247          atab.addRow() ;
     5248          copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
     5249          acc.replaceNaN() ;
     5250          const Vector<Bool> &msk = acc.getMask() ;
     5251          convertArray( flag, !msk ) ;
     5252          for (uInt k = 0; k < nchan; ++k) {
     5253            uChar userFlag = 1 << 7;
     5254            if (msk[k]==True) userFlag = 0 << 7;
     5255            flag(k) = userFlag;
     5256          }
     5257          a->flagsCol_.put( outrow, flag ) ;
     5258          a->specCol_.put( outrow, acc.getSpectrum() ) ;
     5259          a->tsysCol_.put( outrow, acc.getTsys() ) ;
     5260          a->integrCol_.put( outrow, acc.getInterval() ) ;
     5261          timeColOut.put( outrow, acc.getTime() ) ;
     5262          a->cycleCol_.put( outrow, 0 ) ;
     5263        }
     5264        acc.reset() ;
     5265        outrow++ ;
     5266      }
     5267    }
     5268
     5269    // accumulate and add last data
     5270    irow = rows[len-1] ;
     5271    s->flagsCol_.get( irow, flag ) ;
     5272    convertArray( bflag, flag ) ;
     5273    s->specCol_.get( irow, spec ) ;
     5274    tsys.assign( s->tsysCol_( irow ) ) ;
     5275    if (!allEQ(bflag,True) )
     5276      acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
     5277    if ( acc.state() ) {
     5278      atab.addRow() ;
     5279      copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
     5280      acc.replaceNaN() ;
     5281      const Vector<Bool> &msk = acc.getMask() ;
     5282      convertArray( flag, !msk ) ;
     5283      for (uInt k = 0; k < nchan; ++k) {
     5284        uChar userFlag = 1 << 7;
     5285        if (msk[k]==True) userFlag = 0 << 7;
     5286        flag(k) = userFlag;
     5287      }
     5288      a->flagsCol_.put( outrow, flag ) ;
     5289      a->specCol_.put( outrow, acc.getSpectrum() ) ;
     5290      a->tsysCol_.put( outrow, acc.getTsys() ) ;
     5291      a->integrCol_.put( outrow, acc.getInterval() ) ;
     5292      timeColOut.put( outrow, acc.getTime() ) ;
     5293      a->cycleCol_.put( outrow, 0 ) ;
     5294    }
     5295    acc.reset() ;
     5296    outrow++ ;
     5297
     5298    iter.next() ;
     5299  }
     5300
     5301  return a ;
     5302}

  • trunk/src/STMath.h

    r2186 r2580  
    402402    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
    403403
    404   vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
    405   vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
    406   vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
    407   vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
     404  casa::Vector<casa::Float> getSpectrumFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::Matrix<casa::Float>& sp, string mode = "before" ) ;
     405  casa::Vector<casa::Float> getTcalFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     406  casa::Vector<casa::Float> getTsysFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     407  vector<int> getRowIdFromTime( double reftime, const casa::Vector<casa::Double>& t ) ;
    408408
    409409  // Chopper-Wheel type calibration
    410   vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
    411                                       casa::CountedPtr<Scantable>& off,
    412                                       casa::CountedPtr<Scantable>& sky,
    413                                       casa::CountedPtr<Scantable>& hot,
    414                                       casa::CountedPtr<Scantable>& cold,
    415                                       int index,
    416                                       string antname ) ;
     410  void calibrateCW( casa::CountedPtr<Scantable> &out,
     411                    const casa::CountedPtr<Scantable> &on,
     412                    const casa::CountedPtr<Scantable> &off,
     413                    const casa::CountedPtr<Scantable> &sky,
     414                    const casa::CountedPtr<Scantable> &hot,
     415                    const casa::CountedPtr<Scantable> &cold,
     416                    const casa::Vector<casa::uInt> &rows,
     417                    const casa::String &antname ) ;
     418
    417419  // Tsys * (ON-OFF)/OFF
    418   vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
    419                                       casa::CountedPtr<Scantable>& off,
    420                                       int index ) ;
    421   vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
    422                                         casa::CountedPtr<Scantable>& ref,
    423                                         casa::CountedPtr<Scantable>& sky,
    424                                         casa::CountedPtr<Scantable>& hot,
    425                                         casa::CountedPtr<Scantable>& cold,
    426                                         int index ) ;
    427   vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
    428                                         casa::CountedPtr<Scantable>& ref,
    429                                         vector< casa::CountedPtr<Scantable> >& sky,
    430                                         vector< casa::CountedPtr<Scantable> >& hot,
    431                                         vector< casa::CountedPtr<Scantable> >& cold,
    432                                         int index ) ;
    433   double getMJD( string strtime ) ;
     420  void calibrateALMA( casa::CountedPtr<Scantable>& out,
     421                      const casa::CountedPtr<Scantable>& on,
     422                      const casa::CountedPtr<Scantable>& off,
     423                      const casa::Vector<casa::uInt>& rows ) ;
     424
     425  // Frequency switching
     426  void calibrateFS( casa::CountedPtr<Scantable> &sig,
     427                    casa::CountedPtr<Scantable> &ref,
     428                    const casa::CountedPtr<Scantable> &rsig,
     429                    const casa::CountedPtr<Scantable> &rref,
     430                    const casa::CountedPtr<Scantable> &sky,
     431                    const casa::CountedPtr<Scantable> &hot,
     432                    const casa::CountedPtr<Scantable> &cold,
     433                    const casa::Vector<casa::uInt> &rows ) ;
     434  void calibrateAPEXFS( casa::CountedPtr<Scantable> &sig,
     435                        casa::CountedPtr<Scantable> &ref,
     436                        const vector< casa::CountedPtr<Scantable> > &on,
     437                        const vector< casa::CountedPtr<Scantable> > &sky,
     438                        const vector< casa::CountedPtr<Scantable> > &hot,
     439                        const vector< casa::CountedPtr<Scantable> > &cold,
     440                        const casa::Vector<casa::uInt> &rows ) ;
     441  void copyRows( casa::Table &out,
     442                 const casa::Table &in,
     443                 casa::uInt startout,
     444                 casa::uInt startin,
     445                 casa::uInt nrow,
     446                 casa::Bool copySpectra=true,
     447                 casa::Bool copyFlagtra=true,
     448                 casa::Bool copyTsys=true ) ;
     449  casa::CountedPtr<Scantable> averageWithinSession( casa::CountedPtr<Scantable> &s,
     450                                                    vector<bool> &mask,
     451                                                    string weight ) ;
    434452
    435453  bool insitu_;

  • trunk/src/python_asap.cpp

    r2526 r2580  
    8686  asap::python::python_SrcType();
    8787  asap::python::python_STGrid();
     88  asap::python::python_Iterator();
    8889
    8990#ifndef HAVE_LIBPYRAP

  • trunk/src/python_asap.h

    r2356 r2580  
    5353    void python_SrcType();
    5454    void python_STGrid();
     55    void python_Iterator();
    5556
    5657  } // python
Note: See TracChangeset for help on using the changeset viewer.