Changeset 2919


Ignore:
Timestamp:
04/04/14 20:08:17 (11 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...

Refactoring.


Location:
trunk/src
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/CalibrationHelper.h

    r2917 r2919  
    1414
    1515namespace {
    16 // Iteration Helper
    17 
    18 // template<class T>
    19 // class IterationHelper
    20 // {
    21 // public:
    22 //   static void Iterate(T &calibrator)
    23 //   {
    24 //     vector<string> cols( 3 ) ;
    25 //     cols[0] = "BEAMNO" ;
    26 //     cols[1] = "POLNO" ;
    27 //     cols[2] = "IFNO" ;
    28 //     STIdxIter2 calibrator.GetIterator(cols);//iter( out, cols ) ;
    29 //     STSelector sel ;
    30 //     while ( !iter.pastEnd() ) {
    31 //       Record ids = iter.currentValue() ;
    32 //       stringstream ss ;
    33 //       ss << "SELECT FROM $1 WHERE "
    34 //          << "BEAMNO==" << ids.asuInt(cols[0]) << "&&"
    35 //          << "POLNO==" << ids.asuInt(cols[1]) << "&&"
    36 //          << "IFNO==" << ids.asuInt(cols[2]) ;
    37 //       //cout << "TaQL string: " << ss.str() << endl ;
    38 //       sel.setTaQL( ss.str() ) ;
    39 //       Vector<uInt> rows = iter.getRows( SHARE ) ;
    40 //       calibrator.Calibrate(sel, rows);
    41 //       aoff->setSelection( sel ) ;
    42 //       // out should be an exact copy of s except that SPECTRA column is empty
    43 //       calibrateALMA( out, s, aoff, rows ) ;
    44 //       aoff->unsetSelection() ;
    45 //       sel.reset() ;
    46 //       iter.next() ;
    47 //     }   
    48 //   }
    49 // };
    50 
    5116// Interpolation Helper
    5217class TcalData
    … …  
    8045  const String method_name() const {return "getTsysFromTime";}
    8146  uInt nrow() const {return tsyscolumn_.nrow();}
    82   Vector<Float> GetEntry(int idx) const
    83   {
    84     return tsyscolumn_(idx);
    85   } 
     47  Vector<Float> GetEntry(int idx) const {return tsyscolumn_(idx);}
    8648private:
    8749  ROArrayColumn<Float> tsyscolumn_;
    … …  
    9759  const String method_name() const {return "getSpectraFromTime";}
    9860  uInt nrow() const {return data_.ncolumn();}
    99   Vector<Float> GetEntry(int idx) const
    100   {
    101     return data_.column(idx);
    102   } 
     61  Vector<Float> GetEntry(int idx) const {return data_.column(idx);}
    10362private:
    10463  Matrix<Float> data_;
    … …  
    255214};
    256215 
     216// Calibration Helper
     217class CalibrationHelper
     218{
     219public:
     220  static void CalibrateALMA( CountedPtr<Scantable>& out,
     221                             const CountedPtr<Scantable>& on,
     222                             const CountedPtr<Scantable>& off,
     223                             const Vector<uInt>& rows )
     224  {
     225    // 2012/05/22 TN
     226    // Assume that out has empty SPECTRA column
     227
     228    // if rows is empty, just return
     229    if ( rows.nelements() == 0 )
     230      return ;
     231
     232    ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
     233    ROArrayColumn<Float> in_tsys_column(on->table(), "TSYS");
     234    ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
     235    ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
     236    ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
     237   
     238    Vector<Double> timeVec = GetScalarColumn<Double>(off->table(), "TIME");
     239    Vector<Double> refTimeVec = GetScalarColumn<Double>(on->table(), "TIME");
     240    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
     241    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
     242    vector<int> ids( 2 ) ;
     243    for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     244      uInt row = rows[irow];
     245      double reftime = refTimeVec[row];
     246      ids = getRowIdFromTime( reftime, timeVec ) ;
     247      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
     248      Vector<Float> spec = in_spectra_column(row);
     249      Vector<Float> tsys = in_tsys_column(row);
     250      Vector<uChar> flag = in_flagtra_column(row);
     251      // ALMA Calibration
     252      //
     253      // Ta* = Tsys * ( ON - OFF ) / OFF
     254      //
     255      // 2010/01/07 Takeshi Nakazato
     256      unsigned int tsyssize = tsys.nelements() ;
     257      for ( unsigned int j = 0 ; j < spsize ; j++ ) {
     258        if ( spoff[j] == 0.0 ) {
     259          spec[j] = 0.0 ;
     260          flag[j] = (uChar)True;
     261        }
     262        else {
     263          spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
     264        }
     265        spec[j] *= (tsyssize == spsize) ? tsys[j] : tsys[0];
     266      }
     267      out_spectra_column.put(row, spec);
     268      out_flagtra_column.put(row, flag);
     269    }
     270  }
     271  static void CalibrateChopperWheel( CountedPtr<Scantable> &out,
     272                                     const CountedPtr<Scantable>& on,
     273                                     const CountedPtr<Scantable>& off,
     274                                     const CountedPtr<Scantable>& sky,
     275                                     const CountedPtr<Scantable>& hot,
     276                                     const CountedPtr<Scantable>& cold,
     277                                     const Vector<uInt> &rows )
     278  {
     279    // 2012/05/22 TN
     280    // Assume that out has empty SPECTRA column
     281   
     282    // if rows is empty, just return
     283    if ( rows.nelements() == 0 )
     284      return ;
     285
     286    string antenna_name = out->getAntennaName();
     287    ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
     288    ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
     289    ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
     290    ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
     291    ArrayColumn<Float> out_tsys_column(out->table(), "TSYS");
     292       
     293    Vector<Double> timeOff = GetScalarColumn<Double>(off->table(), "TIME");
     294    Vector<Double> timeSky = GetScalarColumn<Double>(sky->table(), "TIME");
     295    Vector<Double> timeHot = GetScalarColumn<Double>(hot->table(), "TIME");
     296    Vector<Double> timeOn = GetScalarColumn<Double>(on->table(), "TIME");
     297    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
     298    SpectralData skyspectra(Matrix<Float>(GetArrayColumn<Float>(sky->table(), "SPECTRA")));
     299    SpectralData hotspectra(Matrix<Float>(GetArrayColumn<Float>(hot->table(), "SPECTRA")));
     300    TcalData tcaldata(sky);
     301    TsysData tsysdata(sky);
     302    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
     303    vector<int> ids( 2 ) ;
     304    for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     305      uInt row = rows[irow];
     306      double reftime = timeOn[row];
     307      ids = getRowIdFromTime( reftime, timeOff ) ;
     308      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
     309      ids = getRowIdFromTime( reftime, timeSky ) ;
     310      Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
     311      Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
     312      Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
     313      ids = getRowIdFromTime( reftime, timeHot ) ;
     314      Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
     315      Vector<Float> spec = in_spectra_column(row);
     316      Vector<uChar> flag = in_flagtra_column(row);
     317      if ( antenna_name.find( "APEX" ) != String::npos ) {
     318        // using gain array
     319        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     320          if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
     321            spec[j] = 0.0 ;
     322            flag[j] = (uChar)True;
     323          }
     324          else {
     325            spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
     326              * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     327          }
     328        }
     329      }
     330      else {
     331        // Chopper-Wheel calibration (Ulich & Haas 1976)
     332        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     333          if ( (sphot[j]-spsky[j]) == 0.0 ) {
     334            spec[j] = 0.0 ;
     335            flag[j] = (uChar)True;
     336          }
     337          else {
     338            spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     339          }
     340        }
     341      }
     342      out_spectra_column.put(row, spec);
     343      out_flagtra_column.put(row, flag);
     344      out_tsys_column.put(row, tsys);
     345    }
     346  }
     347  static void GetSelector(STSelector &sel, const vector<string> &names, const Record &values)
     348  {
     349    stringstream ss ;
     350    ss << "SELECT FROM $1 WHERE ";
     351    string separator = "";
     352    for (vector<string>::const_iterator i = names.begin(); i != names.end(); ++i) {
     353      ss << separator << *i << "==";
     354      switch (values.dataType(*i)) {
     355      case TpUInt:
     356        ss << values.asuInt(*i);
     357        break;
     358      case TpInt:
     359        ss << values.asInt(*i);
     360        break;
     361      case TpFloat:
     362        ss << values.asFloat(*i);
     363        break;
     364      case TpDouble:
     365        ss << values.asDouble(*i);
     366        break;
     367      case TpComplex:
     368        ss << values.asComplex(*i);
     369        break;
     370      case TpString:
     371        ss << values.asString(*i);
     372        break;
     373      default:
     374        break;
     375      }
     376      separator = "&&";
     377    }
     378    sel.setTaQL(ss.str());
     379  }
     380private:
     381  template<class T>
     382  static Vector<T> GetScalarColumn(const Table &table, const String &name)
     383  {
     384    ROScalarColumn<T> column(table, name);
     385    return column.getColumn();
     386  }
     387  template<class T>
     388  static Array<T> GetArrayColumn(const Table &table, const String &name)
     389  {
     390    ROArrayColumn<T> column(table, name);
     391    return column.getColumn();
     392  }
     393};
     394 
     395class AlmaCalibrator
     396{
     397public:
     398  AlmaCalibrator(CountedPtr<Scantable> &out,
     399                 const CountedPtr<Scantable> &on,
     400                 const CountedPtr<Scantable> &off)
     401    : target_(out),
     402      selector_(),
     403      on_(on),
     404      off_(off)
     405  {}
     406  ~AlmaCalibrator() {}
     407  CountedPtr<Scantable> target() const {return target_;}
     408  void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
     409    CalibrationHelper::GetSelector(selector_, cols, values);
     410    off_->setSelection(selector_);
     411    CalibrationHelper::CalibrateALMA(target_, on_, off_, rows);
     412    off_->unsetSelection();
     413  }
     414private:
     415  CountedPtr<Scantable> target_;
     416  STSelector selector_;
     417  const CountedPtr<Scantable> on_;
     418  const CountedPtr<Scantable> off_;
     419};
     420
     421class ChopperWheelCalibrator
     422{
     423public:
     424  ChopperWheelCalibrator(CountedPtr<Scantable> &out,
     425                         const CountedPtr<Scantable> &on,
     426                         const CountedPtr<Scantable> &sky,
     427                         const CountedPtr<Scantable> &off,
     428                         const CountedPtr<Scantable> &hot,
     429                         const CountedPtr<Scantable> &cold)
     430    : target_(out),
     431      selector_(),
     432      on_(on),
     433      off_(off),
     434      sky_(sky),
     435      hot_(hot),
     436      cold_(cold)
     437  {}
     438  ~ChopperWheelCalibrator() {}
     439  CountedPtr<Scantable> target() const {return target_;}
     440  void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
     441    CalibrationHelper::GetSelector(selector_, cols, values);
     442    off_->setSelection(selector_);
     443    sky_->setSelection(selector_);
     444    hot_->setSelection(selector_);
     445    CalibrationHelper::CalibrateChopperWheel(target_, on_, off_, sky_, hot_, cold_, rows);
     446    off_->unsetSelection();
     447    sky_->unsetSelection();
     448    hot_->unsetSelection();
     449  }
     450private:
     451  CountedPtr<Scantable> target_;
     452  STSelector selector_;
     453  const CountedPtr<Scantable> on_;
     454  const CountedPtr<Scantable> off_;
     455  const CountedPtr<Scantable> sky_;
     456  const CountedPtr<Scantable> hot_;
     457  const CountedPtr<Scantable> cold_;
     458};
    257459 
    258460} // anonymous namespace

  • trunk/src/STMath.cpp

    r2918 r2919  
    5858
    5959#include "CalibrationHelper.h"
     60#include "IterationHelper.h"
    6061
    6162using namespace casa;
    6263using namespace asap;
    63 
    64 // namespace {
    65 // class CalibrationIterator
    66 // {
    67 // public:
    68 //   CalibrationIterator()
    69 //     : table_(table)
    70 //   {}
    71 //   ~CalibrationIterator() {}
    72 //   void Iterate() {
    73 //     STSelector sel;
    74 //     vector<string> cols( 3 ) ;
    75 //     cols[0] = "BEAMNO" ;
    76 //     cols[1] = "POLNO" ;
    77 //     cols[2] = "IFNO" ;
    78 //     STIdxIter2 iter(out, cols);
    79 //     while ( !iter.pastEnd() ) {
    80 //       Record ids = iter.currentValue();
    81 //       stringstream ss ;
    82 //       ss << "SELECT FROM $1 WHERE "
    83 //          << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    84 //          << "POLNO==" << ids.asuInt("POLNO") << "&&"
    85 //          << "IFNO==" << ids.asuInt("IFNO") ;
    86 //       //cout << "TaQL string: " << ss.str() << endl ;
    87 //       sel.setTaQL( ss.str() ) ;
    88 //       ////
    89 //       // begin individual processing
    90 //       aoff->setSelection( sel ) ;
    91 //       ahot->setSelection( sel ) ;
    92 //       asky->setSelection( sel ) ;
    93 //       Vector<uInt> rows = iter.getRows( SHARE ) ;
    94 //       // out should be an exact copy of s except that SPECTRA column is empty
    95 //       calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
    96 //       aoff->unsetSelection() ;
    97 //       ahot->unsetSelection() ;
    98 //       asky->unsetSelection() ;
    99 //       // end individual processing
    100 //       ////
    101 //       sel.reset() ;
    102 //       iter.next() ;
    103 //     }   
    104 //   };
    105 // protected:
    106 //   CountedPtr<Scantable> table_;
    107 // };
    108 
    109 
    110 // } // anonymous namespace
    11164
    11265// 2012/02/17 TN
    … …  
    35813534    copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
    35823535   
    3583     // process each on scan
    3584     STSelector sel ;
    3585     vector<string> cols( 3 ) ;
    3586     cols[0] = "BEAMNO" ;
    3587     cols[1] = "POLNO" ;
    3588     cols[2] = "IFNO" ;
    3589     STIdxIter2 iter(out, cols);
    3590     while ( !iter.pastEnd() ) {
    3591       Record ids = iter.currentValue();
    3592       stringstream ss ;
    3593       ss << "SELECT FROM $1 WHERE "
    3594          << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3595          << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3596          << "IFNO==" << ids.asuInt("IFNO") ;
    3597       //cout << "TaQL string: " << ss.str() << endl ;
    3598       sel.setTaQL( ss.str() ) ;
    3599       aoff->setSelection( sel ) ;
    3600       ahot->setSelection( sel ) ;
    3601       asky->setSelection( sel ) ;
    3602       Vector<uInt> rows = iter.getRows( SHARE ) ;
    3603       // out should be an exact copy of s except that SPECTRA column is empty
    3604       calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
    3605       aoff->unsetSelection() ;
    3606       ahot->unsetSelection() ;
    3607       asky->unsetSelection() ;
    3608       sel.reset() ;
    3609       iter.next() ;
    3610     }
     3536    // iterate throgh IterationHelper
     3537    ChopperWheelCalibrator cal(out, s, aoff, asky, ahot, acold);
     3538    IterationHelper<ChopperWheelCalibrator>::Iterate(cal, "BEAMNO,POLNO,IFNO");
     3539
    36113540    s->table_ = torg ;
    36123541    s->attach() ;
    … …  
    36643593//     t0 = mathutil::gettimeofday_sec() ;
    36653594
    3666     // using STIdxIter2
    3667     vector<string> cols( 3 ) ;
    3668     cols[0] = "BEAMNO" ;
    3669     cols[1] = "POLNO" ;
    3670     cols[2] = "IFNO" ;
    3671     STIdxIter2 iter( out, cols ) ;
    3672     STSelector sel ;
    3673     while ( !iter.pastEnd() ) {
    3674       Record ids = iter.currentValue() ;
    3675       stringstream ss ;
    3676       ss << "SELECT FROM $1 WHERE "
    3677          << "BEAMNO==" << ids.asuInt(cols[0]) << "&&"
    3678          << "POLNO==" << ids.asuInt(cols[1]) << "&&"
    3679          << "IFNO==" << ids.asuInt(cols[2]) ;
    3680       //cout << "TaQL string: " << ss.str() << endl ;
    3681       sel.setTaQL( ss.str() ) ;
    3682       aoff->setSelection( sel ) ;
    3683       Vector<uInt> rows = iter.getRows( SHARE ) ;
    3684       // out should be an exact copy of s except that SPECTRA column is empty
    3685       calibrateALMA( out, s, aoff, rows ) ;
    3686       aoff->unsetSelection() ;
    3687       sel.reset() ;
    3688       iter.next() ;
    3689     }
     3595    // iterate throgh IterationHelper
     3596    AlmaCalibrator cal(out, s, aoff);
     3597    IterationHelper<AlmaCalibrator>::Iterate(cal, "BEAMNO,POLNO,IFNO");
     3598
     3599    // finalize
    36903600    s->table_ = torg ;
    36913601    s->attach() ;
    … …  
    38483758     
    38493759      // process each sig and ref scan
    3850 //       STSelector sel ;
    3851       vector<string> cols( 3 ) ;
    3852       cols[0] = "BEAMNO" ;
    3853       cols[1] = "POLNO" ;
    3854       cols[2] = "IFNO" ;
    3855       STIdxIter2 iter(ssig, cols);
    3856       STSelector sel ;
    3857       while ( !iter.pastEnd() ) {
    3858         Record ids = iter.currentValue() ;
    3859         stringstream ss ;
    3860         ss << "SELECT FROM $1 WHERE "
    3861            << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3862            << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3863            << "IFNO==" << ids.asuInt("IFNO") ;
    3864         //cout << "TaQL string: " << ss.str() << endl ;
    3865         sel.setTaQL( ss.str() ) ;
    3866         aofflo->setSelection( sel ) ;
    3867         ahotlo->setSelection( sel ) ;
    3868         askylo->setSelection( sel ) ;
    3869         Vector<uInt> rows = iter.getRows( SHARE ) ;
    3870         calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ;
    3871         aofflo->unsetSelection() ;
    3872         ahotlo->unsetSelection() ;
    3873         askylo->unsetSelection() ;
    3874         sel.reset() ;
    3875         iter.next() ;
    3876       }
    3877       iter = STIdxIter2(sref, cols);
    3878       while ( !iter.pastEnd() ) {
    3879         Record ids = iter.currentValue() ;
    3880         stringstream ss ;
    3881         ss << "SELECT FROM $1 WHERE "
    3882            << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3883            << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3884            << "IFNO==" << ids.asuInt("IFNO") ;
    3885         //cout << "TaQL string: " << ss.str() << endl ;
    3886         sel.setTaQL( ss.str() ) ;
    3887         aoffhi->setSelection( sel ) ;
    3888         ahothi->setSelection( sel ) ;
    3889         askyhi->setSelection( sel ) ;
    3890         Vector<uInt> rows = iter.getRows( SHARE ) ;
    3891         calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ;
    3892         aoffhi->unsetSelection() ;
    3893         ahothi->unsetSelection() ;
    3894         askyhi->unsetSelection() ;
    3895         sel.reset() ;
    3896         iter.next() ;
    3897       }
     3760      // iterate throgh IterationHelper
     3761      ChopperWheelCalibrator cal_sig(ssig, rsig, aofflo, askylo, ahotlo, acoldlo);
     3762      IterationHelper<ChopperWheelCalibrator>::Iterate(cal_sig, "BEAMNO,POLNO,IFNO");
     3763      ChopperWheelCalibrator cal_ref(sref, rref, aoffhi, askyhi, ahothi, acoldhi);
     3764      IterationHelper<ChopperWheelCalibrator>::Iterate(cal_ref, "BEAMNO,POLNO,IFNO");
    38983765    }
    38993766  }
    … …  
    40533920
    40543921  return out ;
    4055 }
    4056 
    4057 void STMath::calibrateCW( CountedPtr<Scantable> &out,
    4058                           const CountedPtr<Scantable>& on,
    4059                           const CountedPtr<Scantable>& off,
    4060                           const CountedPtr<Scantable>& sky,
    4061                           const CountedPtr<Scantable>& hot,
    4062                           const CountedPtr<Scantable>& cold,
    4063                           const Vector<uInt> &rows,
    4064                           const String &antname )
    4065 {
    4066   // 2012/05/22 TN
    4067   // Assume that out has empty SPECTRA column
    4068 
    4069   // if rows is empty, just return
    4070   if ( rows.nelements() == 0 )
    4071     return ;
    4072   ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    4073   Vector<Double> timeOff = timeCol.getColumn() ;
    4074   timeCol.attach( sky->table(), "TIME" ) ;
    4075   Vector<Double> timeSky = timeCol.getColumn() ;
    4076   timeCol.attach( hot->table(), "TIME" ) ;
    4077   Vector<Double> timeHot = timeCol.getColumn() ;
    4078   timeCol.attach( on->table(), "TIME" ) ;
    4079   ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
    4080   SpectralData offspectra(arrayFloatCol.getColumn());
    4081   arrayFloatCol.attach( sky->table(), "SPECTRA" ) ;
    4082   SpectralData skyspectra(arrayFloatCol.getColumn());
    4083   arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
    4084   SpectralData hotspectra(arrayFloatCol.getColumn());
    4085   TcalData tcaldata(sky);
    4086   TsysData tsysdata(sky);
    4087   unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    4088   // I know that the data is contiguous
    4089   const uInt *p = rows.data() ;
    4090   vector<int> ids( 2 ) ;
    4091   Block<uInt> flagchan( spsize ) ;
    4092   uInt nflag = 0 ;
    4093   for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    4094     double reftime = timeCol.asdouble(*p) ;
    4095     ids = getRowIdFromTime( reftime, timeOff ) ;
    4096     Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
    4097     ids = getRowIdFromTime( reftime, timeSky ) ;
    4098     Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
    4099     Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
    4100     Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
    4101     ids = getRowIdFromTime( reftime, timeHot ) ;
    4102     Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
    4103     Vector<Float> spec = on->specCol_( *p ) ;
    4104     if ( antname.find( "APEX" ) != String::npos ) {
    4105       // using gain array
    4106       for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4107         if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
    4108           spec[j] = 0.0 ;
    4109           flagchan[nflag++] = j ;
    4110         }
    4111         else {
    4112           spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
    4113             * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
    4114         }
    4115       }
    4116     }
    4117     else {
    4118       // Chopper-Wheel calibration (Ulich & Haas 1976)
    4119       for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4120         if ( (sphot[j]-spsky[j]) == 0.0 ) {
    4121           spec[j] = 0.0 ;
    4122           flagchan[nflag++] = j ;
    4123         }
    4124         else {
    4125           spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
    4126         }
    4127       }
    4128     }
    4129     out->specCol_.put( *p, spec ) ;
    4130     out->tsysCol_.put( *p, tsys ) ;
    4131     if ( nflag > 0 ) {
    4132       Vector<uChar> fl = out->flagsCol_( *p ) ;
    4133       for ( unsigned int j = 0 ; j < nflag ; j++ ) {
    4134         fl[flagchan[j]] = (uChar)True ;
    4135       }
    4136       out->flagsCol_.put( *p, fl ) ;
    4137     }
    4138     nflag = 0 ;
    4139     p++ ;
    4140   }
    4141 }
    4142 
    4143 void STMath::calibrateALMA( CountedPtr<Scantable>& out,
    4144                             const CountedPtr<Scantable>& on,
    4145                             const CountedPtr<Scantable>& off,
    4146                             const Vector<uInt>& rows )
    4147 {
    4148   // 2012/05/22 TN
    4149   // Assume that out has empty SPECTRA column
    4150 
    4151   // if rows is empty, just return
    4152   if ( rows.nelements() == 0 )
    4153     return ;
    4154   ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    4155   Vector<Double> timeVec = timeCol.getColumn() ;
    4156   timeCol.attach( on->table(), "TIME" ) ;
    4157   ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
    4158   SpectralData offspectra(arrayFloatCol.getColumn());
    4159   unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    4160   // I know that the data is contiguous
    4161   const uInt *p = rows.data() ;
    4162   vector<int> ids( 2 ) ;
    4163   Block<uInt> flagchan( spsize ) ;
    4164   uInt nflag = 0 ;
    4165   for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    4166     double reftime = timeCol.asdouble(*p) ;
    4167     ids = getRowIdFromTime( reftime, timeVec ) ;
    4168     Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
    4169     Vector<Float> spec = on->specCol_( *p ) ;
    4170     Vector<Float> tsys = on->tsysCol_( *p ) ;
    4171     // ALMA Calibration
    4172     //
    4173     // Ta* = Tsys * ( ON - OFF ) / OFF
    4174     //
    4175     // 2010/01/07 Takeshi Nakazato
    4176     unsigned int tsyssize = tsys.nelements() ;
    4177     for ( unsigned int j = 0 ; j < spsize ; j++ ) {
    4178       if ( spoff[j] == 0.0 ) {
    4179         spec[j] = 0.0 ;
    4180         flagchan[nflag++] = j ;
    4181       }
    4182       else {
    4183         spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
    4184       }
    4185       if ( tsyssize == spsize )
    4186         spec[j] *= tsys[j] ;
    4187       else
    4188         spec[j] *= tsys[0] ;
    4189     }
    4190     out->specCol_.put( *p, spec ) ;
    4191     if ( nflag > 0 ) {
    4192       Vector<uChar> fl = out->flagsCol_( *p ) ;
    4193       for ( unsigned int j = 0 ; j < nflag ; j++ ) {
    4194         fl[flagchan[j]] = (uChar)True ;
    4195       }
    4196       out->flagsCol_.put( *p, fl ) ;
    4197     }
    4198     nflag = 0 ;
    4199     p++ ;
    4200   }
    42013922}
    42023923

  • trunk/src/STMath.h

    r2917 r2919  
    401401  casa::Vector<casa::uChar>
    402402    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
    403 
    404   // Chopper-Wheel type calibration
    405   void calibrateCW( casa::CountedPtr<Scantable> &out,
    406                     const casa::CountedPtr<Scantable> &on,
    407                     const casa::CountedPtr<Scantable> &off,
    408                     const casa::CountedPtr<Scantable> &sky,
    409                     const casa::CountedPtr<Scantable> &hot,
    410                     const casa::CountedPtr<Scantable> &cold,
    411                     const casa::Vector<casa::uInt> &rows,
    412                     const casa::String &antname ) ;
    413 
    414   // Tsys * (ON-OFF)/OFF
    415   void calibrateALMA( casa::CountedPtr<Scantable>& out,
    416                       const casa::CountedPtr<Scantable>& on,
    417                       const casa::CountedPtr<Scantable>& off,
    418                       const casa::Vector<casa::uInt>& rows ) ;
    419403
    420404  // Frequency switching
Note: See TracChangeset for help on using the changeset viewer.