Changeset 2546


Ignore:
Timestamp:
05/23/12 17:51:58 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Use direct access to Scantable.specCol_ instead to call Scantable::getSpectrum.


Location:
branches/hpc33/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/hpc33/src/STMath.cpp

    r2545 r2546  
    46014601}
    46024602
    4603 vector<float> STMath::getSpectrumFromTime( double reftime,
     4603Vector<Float> STMath::getSpectrumFromTime( double reftime,
    46044604                                           Vector<Double> &timeVec,
    46054605                                           const CountedPtr<Scantable>& s,
    … …  
    46074607{
    46084608  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
    4609   vector<float> sp ;
     4609  Vector<Float> sp ;
    46104610
    46114611  if ( s->nrow() == 0 ) {
    … …  
    46154615  else if ( s->nrow() == 1 ) {
    46164616    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
    4617     return s->getSpectrum( 0 ) ;
     4617    s->specCol_.get( 0, sp ) ;
     4618    return sp ;
    46184619  }
    46194620  else {
    4620     //ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
    4621     //Vector<Double> timeVec = timeCol.getColumn() ;
    46224621    vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
    46234622    if ( mode == "before" ) {
    … …  
    46314630      }
    46324631      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4633       sp = s->getSpectrum( id ) ;
     4632      s->specCol_.get( id, sp ) ;
     4633      //sp = s->getSpectrum( id ) ;
    46344634    }
    46354635    else if ( mode == "after" ) {
    … …  
    46434643      }
    46444644      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4645       sp = s->getSpectrum( id ) ;
     4645      s->specCol_.get( id, sp ) ;
     4646      //sp = s->getSpectrum( id ) ;
    46464647    }
    46474648    else if ( mode == "nearest" ) {
    … …  
    46574658      }
    46584659      else {
    4659         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4660         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4661 //         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4662 //         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    46634660        double t0 = timeVec[idx[0]] ;
    46644661        double t1 = timeVec[idx[1]] ;
    4665 //         cout << "t0-t0c=" << t0-t0c << endl ;
    4666 //         cout << "t1-t1c=" << t1-t1c << endl ;
    4667 //         double tref = getMJD( reftime ) ;
    46684662        double tref = reftime ;
    46694663        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
    … …  
    46754669      }
    46764670      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4677       sp = s->getSpectrum( id ) ;     
     4671      s->specCol_.get( id, sp ) ;
     4672      //sp = s->getSpectrum( id ) ;     
    46784673    }
    46794674    else if ( mode == "linear" ) {
    … …  
    46834678        int id = idx[1] ;
    46844679        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4685         sp = s->getSpectrum( id ) ;
     4680        s->specCol_.get( id, sp ) ;
     4681        //sp = s->getSpectrum( id ) ;
    46864682      }
    46874683      else if ( idx[1] == -1 ) {
    … …  
    46904686        int id = idx[0] ;
    46914687        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4692         sp = s->getSpectrum( id ) ;
     4688        s->specCol_.get( id, sp ) ;
     4689        //sp = s->getSpectrum( id ) ;
    46934690      }
    46944691      else if ( idx[0] == idx[1] ) {
    … …  
    46974694        int id = idx[0] ;
    46984695        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4699         sp = s->getSpectrum( id ) ;
     4696        s->specCol_.get( id, sp ) ;
     4697        //sp = s->getSpectrum( id ) ;
    47004698      }
    47014699      else {
    47024700        // do interpolation
    47034701        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    4704         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4705         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4706 //         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4707 //         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    47084702        double t0 = timeVec[idx[0]] ;
    47094703        double t1 = timeVec[idx[1]] ;
    4710 //         cout << "t0-t0c=" << t0-t0c << endl ;
    4711 //         cout << "t1-t1c=" << t1-t1c << endl ;
    4712 //         double tref = getMJD( reftime ) ;
    47134704        double tref = reftime ;
    4714         sp = s->getSpectrum( idx[0] ) ;
    4715         vector<float> sp1 = s->getSpectrum( idx[1] ) ;
     4705        s->specCol_.get( idx[0], sp ) ;
     4706        //sp = s->getSpectrum( idx[0] ) ;
     4707        //vector<float> sp1 = s->getSpectrum( idx[1] ) ;
     4708        Vector<Float> sp1 = s->specCol_( idx[1] ) ;
    47164709        double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
    47174710        for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
    … …  
    47244717    }
    47254718    return sp ;
     4719  }
     4720}
     4721
     4722void STMath::getSpectrumFromTime2( Vector<Float> &sp,
     4723                                   double reftime,
     4724                                   Vector<Double> &timeVec,
     4725                                   const CountedPtr<Scantable>& s,
     4726                                   string mode )
     4727{
     4728  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
     4729  //Vector<Float> sp ;
     4730
     4731  if ( s->nrow() == 0 ) {
     4732    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
     4733    //return sp ;
     4734  }
     4735  else if ( s->nrow() == 1 ) {
     4736    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
     4737    //return s->getSpectrum( 0 ) ;
     4738    s->specCol_.get( 0, sp ) ;
     4739    //return sp ;
     4740  }
     4741  else {
     4742    vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
     4743    if ( mode == "before" ) {
     4744      int id = -1 ;
     4745      if ( idx[0] != -1 ) {
     4746        id = idx[0] ;
     4747      }
     4748      else if ( idx[1] != -1 ) {
     4749        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     4750        id = idx[1] ;
     4751      }
     4752      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4753      //sp = s->getSpectrum( id ) ;
     4754      s->specCol_.get( id, sp ) ;
     4755    }
     4756    else if ( mode == "after" ) {
     4757      int id = -1 ;
     4758      if ( idx[1] != -1 ) {
     4759        id = idx[1] ;
     4760      }
     4761      else if ( idx[0] != -1 ) {
     4762        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     4763        id = idx[1] ;
     4764      }
     4765      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4766      //sp = s->getSpectrum( id ) ;
     4767      s->specCol_.get( id, sp ) ;
     4768    }
     4769    else if ( mode == "nearest" ) {
     4770      int id = -1 ;
     4771      if ( idx[0] == -1 ) {
     4772        id = idx[1] ;
     4773      }
     4774      else if ( idx[1] == -1 ) {
     4775        id = idx[0] ;
     4776      }
     4777      else if ( idx[0] == idx[1] ) {
     4778        id = idx[0] ;
     4779      }
     4780      else {
     4781        double t0 = timeVec[idx[0]] ;
     4782        double t1 = timeVec[idx[1]] ;
     4783        double tref = reftime ;
     4784        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4785          id = idx[1] ;
     4786        }
     4787        else {
     4788          id = idx[0] ;
     4789        }
     4790      }
     4791      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4792      //sp = s->getSpectrum( id ) ;     
     4793      s->specCol_.get( id, sp ) ;
     4794    }
     4795    else if ( mode == "linear" ) {
     4796      if ( idx[0] == -1 ) {
     4797        // use after
     4798        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4799        int id = idx[1] ;
     4800        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4801        //sp = s->getSpectrum( id ) ;
     4802        s->specCol_.get( id, sp ) ;
     4803      }
     4804      else if ( idx[1] == -1 ) {
     4805        // use before
     4806        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4807        int id = idx[0] ;
     4808        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4809        //sp = s->getSpectrum( id ) ;
     4810        s->specCol_.get( id, sp ) ;
     4811      }
     4812      else if ( idx[0] == idx[1] ) {
     4813        // use before
     4814        //os << "No need to interporate." << LogIO::POST ;
     4815        int id = idx[0] ;
     4816        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4817        //sp = s->getSpectrum( id ) ;
     4818        s->specCol_.get( id, sp ) ;
     4819      }
     4820      else {
     4821        // do interpolation
     4822        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4823        //double t0 = timeVec[idx[0]] ;
     4824        //double t1 = timeVec[idx[1]] ;
     4825        //double tref = reftime ;
     4826        //sp = s->getSpectrum( idx[0] ) ;
     4827        s->specCol_.get( idx[0], sp ) ;
     4828        //vector<float> sp1 = s->getSpectrum( idx[1] ) ;
     4829        Vector<Float> sp1 = s->specCol_( idx[1] ) ;
     4830        //double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
     4831        double tfactor = ( reftime - timeVec[idx[0]] ) / ( timeVec[idx[1]] - timeVec[idx[0]] ) ;
     4832        for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
     4833          sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
     4834        }
     4835      }
     4836    }
     4837    else {
     4838      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4839    }
     4840    //return sp ;
    47264841  }
    47274842}
    … …  
    51715286  sel.setPolarizations( ip ) ;
    51725287  off->setSelection( sel ) ;
    5173   vector<float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
     5288  Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    51745289  vector<float> spec = on->getSpectrum( index ) ;
    51755290  //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
    … …  
    52155330  for ( int index = 0 ; index < on->nrow() ; index++ ) {
    52165331    double reftime = timeCol.asdouble(index) ;
    5217     vector<float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
     5332    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    52185333    vector<float> spec = on->getSpectrum( index ) ;
    52195334    Vector<Float> tsys = tsysCol( index ) ;
    … …  
    52475362  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    52485363  Vector<Double> timeVec = timeCol.getColumn() ;
    5249   //ROTableColumn timeCol( on->table(), "TIME" ) ;
    52505364  timeCol.attach( on->table(), "TIME" ) ;
    52515365  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
    … …  
    52555369  const uInt *p = rows.data() ;
    52565370  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    5257     //int index = rows[irow] ;
    5258     //double reftime = timeCol.asdouble(index) ;
    52595371    double reftime = timeCol.asdouble(*p) ;
    5260     vector<float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    5261     //vector<float> spec = on->getSpectrum( index ) ;
    5262     //Vector<Float> tsys = tsysCol( index ) ;
     5372    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    52635373    vector<float> spec = on->getSpectrum( *p ) ;
    52645374    Vector<Float> tsys = tsysCol( *p ) ;
    … …  
    52785388      sp[j] = v ;
    52795389    }
    5280     //on->setSpectrum( sp, index ) ;
    52815390    on->setSpectrum( sp, *p ) ;
    52825391    p++ ;
    … …  
    52955404  if ( rows.nelements() == 0 )
    52965405    return ;
    5297 //   string reftime = on->getTime( index ) ;
    52985406  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    52995407  Vector<Double> timeVec = timeCol.getColumn() ;
    5300   //ROTableColumn timeCol( on->table(), "TIME" ) ;
    53015408  timeCol.attach( on->table(), "TIME" ) ;
    53025409  ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
    5303   //vector<float> sp( on->nchan( on->getIF(rows[0]) ) ) ;
    5304   Vector<Float> sp( on->nchan( on->getIF(rows[0]) ) ) ;
    5305   unsigned int spsize = sp.size() ;
     5410  unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
     5411  //Vector<Float> spoff( spsize ) ;
     5412  //Vector<Float> spec( spsize ) ;
    53065413  // I know that the data is contiguous
    53075414  const uInt *p = rows.data() ;
    53085415  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    5309     //int index = rows[irow] ;
    5310     //double reftime = timeCol.asdouble(index) ;
    53115416    double reftime = timeCol.asdouble(*p) ;
    5312     vector<float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    5313     //vector<float> spec = on->getSpectrum( index ) ;
    5314     //Vector<Float> tsys = tsysCol( index ) ;
    5315     vector<float> spec = on->getSpectrum( *p ) ;
     5417    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
     5418    Vector<Float> spec = on->specCol_( *p ) ;
     5419    //getSpectrumFromTime2( spoff, reftime, timeVec, off, "linear" ) ;
     5420    //on->specCol_.get( *p, spec ) ;
    53165421    Vector<Float> tsys = tsysCol( *p ) ;
    53175422    // ALMA Calibration
    … …  
    53215426    // 2010/01/07 Takeshi Nakazato
    53225427    unsigned int tsyssize = tsys.nelements() ;
    5323     for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
    5324       float tscale = 0.0 ;
     5428    for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     5429      spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
    53255430      if ( tsyssize == spsize )
    5326         tscale = tsys[j] ;
     5431        spec[j] *= tsys[j] ;
    53275432      else
    5328         tscale = tsys[0] ;
    5329       float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
    5330       sp[j] = v ;
    5331     }
    5332     //on->setSpectrum( sp, index ) ;
    5333     //on->setSpectrum( sp, *p ) ;
    5334     // no check for nchan
    5335     //Vector<Float> spv( sp ) ;
    5336     out->specCol_.put( *p, sp ) ;
     5433        spec[j] *= tsys[0] ;
     5434    }
     5435    out->specCol_.put( *p, spec ) ;
    53375436    p++ ;
    53385437  }

  • branches/hpc33/src/STMath.h

    r2544 r2546  
    403403
    404404  vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
    405   vector<float> getSpectrumFromTime( double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
     405  casa::Vector<casa::Float> getSpectrumFromTime( double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
     406  void getSpectrumFromTime2( casa::Vector<casa::Float> &sp, double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
    406407  vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
    407408  vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
Note: See TracChangeset for help on using the changeset viewer.