Changeset 2919 for trunk/src/CalibrationHelper.h

Timestamp:

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

File:

• trunk/src/CalibrationHelper.h (modified) (4 diffs)

trunk/src/CalibrationHelper.h

@@ -14 +14 @@
 
 namespace {
-// Iteration Helper
-
-// template<class T>
-// class IterationHelper
-// {
-// public:
-//   static void Iterate(T &calibrator)
-//   {
-//     vector<string> cols( 3 ) ;
-//     cols[0] = "BEAMNO" ;
-//     cols[1] = "POLNO" ;
-//     cols[2] = "IFNO" ;
-//     STIdxIter2 calibrator.GetIterator(cols);//iter( out, cols ) ;
-//     STSelector sel ;
-//     while ( !iter.pastEnd() ) {
-//       Record ids = iter.currentValue() ;
-//       stringstream ss ;
-//       ss << "SELECT FROM $1 WHERE "
-//          << "BEAMNO==" << ids.asuInt(cols[0]) << "&&"
-//          << "POLNO==" << ids.asuInt(cols[1]) << "&&"
-//          << "IFNO==" << ids.asuInt(cols[2]) ;
-//       //cout << "TaQL string: " << ss.str() << endl ;
-//       sel.setTaQL( ss.str() ) ;
-//       Vector<uInt> rows = iter.getRows( SHARE ) ;
-//       calibrator.Calibrate(sel, rows);
-//       aoff->setSelection( sel ) ;
-//       // out should be an exact copy of s except that SPECTRA column is empty
-//       calibrateALMA( out, s, aoff, rows ) ;
-//       aoff->unsetSelection() ;
-//       sel.reset() ;
-//       iter.next() ;
-//     }    
-//   }
-// };
-
 // Interpolation Helper
 class TcalData
@@ -80 +45 @@
   const String method_name() const {return "getTsysFromTime";}
   uInt nrow() const {return tsyscolumn_.nrow();}
-  Vector<Float> GetEntry(int idx) const
-  {
-    return tsyscolumn_(idx);
-  }  
+  Vector<Float> GetEntry(int idx) const {return tsyscolumn_(idx);}
 private:
   ROArrayColumn<Float> tsyscolumn_;
@@ -97 +59 @@
   const String method_name() const {return "getSpectraFromTime";}
   uInt nrow() const {return data_.ncolumn();}
-  Vector<Float> GetEntry(int idx) const
-  {
-    return data_.column(idx);
-  }  
+  Vector<Float> GetEntry(int idx) const {return data_.column(idx);}
 private:
   Matrix<Float> data_;
@@ -255 +214 @@
 };
   
+// Calibration Helper
+class CalibrationHelper
+{
+public:
+  static void CalibrateALMA( CountedPtr<Scantable>& out,
+                             const CountedPtr<Scantable>& on,
+                             const CountedPtr<Scantable>& off,
+                             const Vector<uInt>& rows )
+  {
+    // 2012/05/22 TN
+    // Assume that out has empty SPECTRA column
+
+    // if rows is empty, just return
+    if ( rows.nelements() == 0 )
+      return ;
+
+    ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
+    ROArrayColumn<Float> in_tsys_column(on->table(), "TSYS");
+    ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
+    ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
+    ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
+    
+    Vector<Double> timeVec = GetScalarColumn<Double>(off->table(), "TIME");
+    Vector<Double> refTimeVec = GetScalarColumn<Double>(on->table(), "TIME");
+    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
+    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
+    vector<int> ids( 2 ) ;
+    for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
+      uInt row = rows[irow];
+      double reftime = refTimeVec[row];
+      ids = getRowIdFromTime( reftime, timeVec ) ;
+      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
+      Vector<Float> spec = in_spectra_column(row);
+      Vector<Float> tsys = in_tsys_column(row);
+      Vector<uChar> flag = in_flagtra_column(row);
+      // ALMA Calibration
+      // 
+      // Ta* = Tsys * ( ON - OFF ) / OFF
+      //
+      // 2010/01/07 Takeshi Nakazato
+      unsigned int tsyssize = tsys.nelements() ;
+      for ( unsigned int j = 0 ; j < spsize ; j++ ) {
+        if ( spoff[j] == 0.0 ) {
+          spec[j] = 0.0 ;
+          flag[j] = (uChar)True;
+        }
+        else {
+          spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
+        }
+        spec[j] *= (tsyssize == spsize) ? tsys[j] : tsys[0];
+      }
+      out_spectra_column.put(row, spec);
+      out_flagtra_column.put(row, flag);
+    }
+  }
+  static void CalibrateChopperWheel( CountedPtr<Scantable> &out,
+                                     const CountedPtr<Scantable>& on,
+                                     const CountedPtr<Scantable>& off,
+                                     const CountedPtr<Scantable>& sky,
+                                     const CountedPtr<Scantable>& hot,
+                                     const CountedPtr<Scantable>& cold,
+                                     const Vector<uInt> &rows )
+  {
+    // 2012/05/22 TN
+    // Assume that out has empty SPECTRA column
+    
+    // if rows is empty, just return
+    if ( rows.nelements() == 0 )
+      return ;
+
+    string antenna_name = out->getAntennaName();
+    ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
+    ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
+    ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
+    ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
+    ArrayColumn<Float> out_tsys_column(out->table(), "TSYS");
+        
+    Vector<Double> timeOff = GetScalarColumn<Double>(off->table(), "TIME");
+    Vector<Double> timeSky = GetScalarColumn<Double>(sky->table(), "TIME");
+    Vector<Double> timeHot = GetScalarColumn<Double>(hot->table(), "TIME");
+    Vector<Double> timeOn = GetScalarColumn<Double>(on->table(), "TIME");
+    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
+    SpectralData skyspectra(Matrix<Float>(GetArrayColumn<Float>(sky->table(), "SPECTRA")));
+    SpectralData hotspectra(Matrix<Float>(GetArrayColumn<Float>(hot->table(), "SPECTRA")));
+    TcalData tcaldata(sky);
+    TsysData tsysdata(sky);
+    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
+    vector<int> ids( 2 ) ;
+    for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
+      uInt row = rows[irow];
+      double reftime = timeOn[row];
+      ids = getRowIdFromTime( reftime, timeOff ) ;
+      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
+      ids = getRowIdFromTime( reftime, timeSky ) ; 
+      Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
+      Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
+      Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
+      ids = getRowIdFromTime( reftime, timeHot ) ;
+      Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
+      Vector<Float> spec = in_spectra_column(row);
+      Vector<uChar> flag = in_flagtra_column(row);
+      if ( antenna_name.find( "APEX" ) != String::npos ) {
+        // using gain array
+        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
+          if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
+            spec[j] = 0.0 ;
+            flag[j] = (uChar)True;
+          }
+          else {
+            spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
+              * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
+          }
+        }
+      }
+      else {
+        // Chopper-Wheel calibration (Ulich & Haas 1976)
+        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
+          if ( (sphot[j]-spsky[j]) == 0.0 ) {
+            spec[j] = 0.0 ;
+            flag[j] = (uChar)True;
+          }
+          else {
+            spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
+          }
+        }
+      }
+      out_spectra_column.put(row, spec);
+      out_flagtra_column.put(row, flag);
+      out_tsys_column.put(row, tsys);
+    }
+  }
+  static void GetSelector(STSelector &sel, const vector<string> &names, const Record &values)
+  {
+    stringstream ss ;
+    ss << "SELECT FROM $1 WHERE ";
+    string separator = "";
+    for (vector<string>::const_iterator i = names.begin(); i != names.end(); ++i) {
+      ss << separator << *i << "==";
+      switch (values.dataType(*i)) {
+      case TpUInt:
+        ss << values.asuInt(*i);
+        break;
+      case TpInt:
+        ss << values.asInt(*i);
+        break;
+      case TpFloat:
+        ss << values.asFloat(*i);
+        break;
+      case TpDouble:
+        ss << values.asDouble(*i);
+        break;
+      case TpComplex:
+        ss << values.asComplex(*i);
+        break;
+      case TpString:
+        ss << values.asString(*i);
+        break;
+      default:
+        break;
+      }
+      separator = "&&";
+    }
+    sel.setTaQL(ss.str());
+  }
+private:
+  template<class T>
+  static Vector<T> GetScalarColumn(const Table &table, const String &name)
+  {
+    ROScalarColumn<T> column(table, name);
+    return column.getColumn();
+  }
+  template<class T>
+  static Array<T> GetArrayColumn(const Table &table, const String &name)
+  {
+    ROArrayColumn<T> column(table, name);
+    return column.getColumn();
+  }
+};
+  
+class AlmaCalibrator
+{
+public:
+  AlmaCalibrator(CountedPtr<Scantable> &out,
+                 const CountedPtr<Scantable> &on,
+                 const CountedPtr<Scantable> &off)
+    : target_(out),
+      selector_(),
+      on_(on),
+      off_(off)
+  {}
+  ~AlmaCalibrator() {}
+  CountedPtr<Scantable> target() const {return target_;}
+  void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
+    CalibrationHelper::GetSelector(selector_, cols, values);
+    off_->setSelection(selector_);
+    CalibrationHelper::CalibrateALMA(target_, on_, off_, rows);
+    off_->unsetSelection();
+  }
+private:
+  CountedPtr<Scantable> target_;
+  STSelector selector_;
+  const CountedPtr<Scantable> on_;
+  const CountedPtr<Scantable> off_;
+};
+
+class ChopperWheelCalibrator
+{
+public:
+  ChopperWheelCalibrator(CountedPtr<Scantable> &out,
+                         const CountedPtr<Scantable> &on,
+                         const CountedPtr<Scantable> &sky,
+                         const CountedPtr<Scantable> &off,
+                         const CountedPtr<Scantable> &hot,
+                         const CountedPtr<Scantable> &cold)
+    : target_(out),
+      selector_(),
+      on_(on),
+      off_(off),
+      sky_(sky),
+      hot_(hot),
+      cold_(cold)
+  {}
+  ~ChopperWheelCalibrator() {}
+  CountedPtr<Scantable> target() const {return target_;}
+  void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
+    CalibrationHelper::GetSelector(selector_, cols, values);
+    off_->setSelection(selector_);
+    sky_->setSelection(selector_);
+    hot_->setSelection(selector_);
+    CalibrationHelper::CalibrateChopperWheel(target_, on_, off_, sky_, hot_, cold_, rows);
+    off_->unsetSelection();
+    sky_->unsetSelection();
+    hot_->unsetSelection();
+  }
+private:
+  CountedPtr<Scantable> target_;
+  STSelector selector_;
+  const CountedPtr<Scantable> on_;
+  const CountedPtr<Scantable> off_;
+  const CountedPtr<Scantable> sky_;
+  const CountedPtr<Scantable> hot_;
+  const CountedPtr<Scantable> cold_;
+};
   
 } // anonymous namespace