Changeset 2915 for trunk


Ignore:
Timestamp:
04/02/14 18:38:27 (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 STCalibration and its derived classes.
Fix for threshold of time gap analysis.


Location:
trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STCalSkyPSAlma.cpp

    r2786 r2915  
    1717#include "STSelector.h"
    1818#include "STCalSkyPSAlma.h"
    19 #include "RowAccumulator.h"
    20 #include "STIdxIter.h"
    2119#include "STDefs.h"
    2220#include <atnf/PKSIO/SrcType.h>
    … …  
    2725namespace asap {
    2826STCalSkyPSAlma::STCalSkyPSAlma(CountedPtr<Scantable> &s)
    29   : STCalibration(s)
     27  : STCalibration(s, "SPECTRA")
    3028{
    3129  applytable_ = new STCalSkyTable(*s, "PSALMA");
    … …  
    5149  }
    5250}
    53 
    54 void STCalSkyPSAlma::fillCalTable()
     51 
     52void STCalSkyPSAlma::appenddata(uInt scanno, uInt cycleno,
     53                                uInt beamno, uInt ifno, uInt polno,
     54                                uInt freqid, Double time, Float elevation,
     55                                Vector<Float> any_data)
    5556{
    56   RowAccumulator acc(W_TINT);
    57  
    58   vector<string> cols(3);
    59   cols[0] = "IFNO";
    60   cols[1] = "POLNO";
    61   cols[2] = "BEAMNO";
    62   STIdxIterAcc iter(scantable_, cols);
    63 
    64   ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME");
    65   Vector<Double> timeSec = tcol->getColumn() * 86400.0;
    66   tcol->attach(scantable_->table(), "INTERVAL");
    67   Vector<Double> intervalSec = tcol->getColumn();
    68   delete tcol;
    69   ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION");
    70   Vector<Float> elevation = ecol->getColumn();
    71   delete ecol;
    72 
    73   ROArrayColumn<Float> specCol(scantable_->table(), "SPECTRA");
    74   ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA");
    75   ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID");
    76 
    77   // dummy Tsys: the following process doesn't need Tsys but RowAccumulator
    78   //             requires to set it with spectral data
    79   Vector<Float> tsys(1, 1.0);
    80 
    81   Double timeCen = 0.0;
    82   Float elCen = 0.0;
    83   uInt count = 0;
    84 
    85   while(!iter.pastEnd()) {
    86     Vector<uInt> rows = iter.getRows(SHARE);
    87     Vector<uInt> current = iter.current();
    88     uInt len = rows.nelements();
    89     if (len == 0) {
    90       iter.next();
    91       continue;
    92     }
    93     else if (len == 1) {
    94       STCalSkyTable *p = dynamic_cast<STCalSkyTable *>(&(*applytable_));
    95       uInt irow = rows[0];
    96       p->appenddata(0, 0, current[2], current[0], current[1],
    97                     freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow));
    98       iter.next();
    99       continue;
    100     }
    101    
    102     uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
    103     Vector<uChar> flag(nchan);
    104     Vector<Bool> bflag(nchan);
    105     Vector<Float> spec(nchan);
    106 
    107     Vector<Double> timeSep(len);
    108     for (uInt i = 0; i < len-1; i++) {
    109       timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
    110     }
    111     Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
    112     timeSep[len-1] = tMedian * 10000.0 ; // any large value
    113 
    114     uInt irow ;
    115     uInt jrow ;
    116     for (uInt i = 0; i < len; i++) {
    117       irow = rows[i];
    118       jrow = (i < len-1) ? rows[i+1] : rows[i];
    119       // accumulate data
    120       flagCol.get(irow, flag);
    121       convertArray(bflag, flag);
    122       specCol.get(irow, spec);
    123       if ( !allEQ(bflag,True) )
    124         acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
    125       timeCen += timeSec[irow];
    126       elCen += elevation[irow];
    127       count++;
    128 
    129       // check time gap
    130       double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
    131       if ( gap > 1.1 ) {
    132         if ( acc.state() ) {
    133           acc.replaceNaN() ;
    134 //           const Vector<Bool> &msk = acc.getMask();
    135 //           convertArray(flag, !msk);
    136 //           for (uInt k = 0; k < nchan; ++k) {
    137 //             uChar userFlag = 1 << 7;
    138 //             if (msk[k]==True) userFlag = 0 << 7;
    139 //             flag(k) = userFlag;
    140 //           }
    141           timeCen /= (Double)count * 86400.0; // sec->day
    142           elCen /= (Float)count;
    143           STCalSkyTable *p = dynamic_cast<STCalSkyTable *>(&(*applytable_));
    144           p->appenddata(0, 0, current[2], current[0], current[1],
    145                         freqidCol(irow), timeCen, elCen, acc.getSpectrum());
    146         }
    147         acc.reset() ;
    148         timeCen = 0.0;
    149         elCen = 0.0;
    150         count = 0;
    151       }
    152     }
    153 
    154     iter.next() ;
    155   }
     57  STCalSkyTable *p = dynamic_cast<STCalSkyTable *>(&(*applytable_));
     58  p->appenddata(scanno, cycleno, beamno, ifno, polno,
     59                freqid, time, elevation, any_data);
    15660}
    15761

  • trunk/src/STCalSkyPSAlma.h

    r2786 r2915  
    4343protected:
    4444  virtual void setupSelector(const STSelector &sel);
    45   virtual void fillCalTable();
     45  virtual void appenddata(casa::uInt scanno, casa::uInt cycleno,
     46                          casa::uInt beamno, casa::uInt ifno, casa::uInt polno,
     47                          casa::uInt freqid, casa::Double time, casa::Float elevation,
     48                          casa::Vector<casa::Float> any_data);
    4649};
    4750

  • trunk/src/STCalTsys.cpp

    r2786 r2915  
    1818#include "STSelector.h"
    1919#include "STCalTsys.h"
    20 #include "RowAccumulator.h"
    21 #include "STIdxIter.h"
    2220#include "STDefs.h"
    2321#include <atnf/PKSIO/SrcType.h>
    … …  
    2826namespace asap {
    2927  STCalTsys::STCalTsys(CountedPtr<Scantable> &s, vector<int> &iflist)
    30     : STCalibration(s),
     28    : STCalibration(s, "TSYS"),
    3129      iflist_(iflist)
    3230{
    … …  
    5957}
    6058
    61 void STCalTsys::fillCalTable()
     59void STCalTsys::appenddata(uInt scanno, uInt cycleno,
     60                           uInt beamno, uInt ifno, uInt polno,
     61                           uInt freqid, Double time, Float elevation,
     62                           Vector<Float> any_data)
    6263{
    63   RowAccumulator acc(W_TINT);
    64  
    65   vector<string> cols(3);
    66   cols[0] = "IFNO";
    67   cols[1] = "POLNO";
    68   cols[2] = "BEAMNO";
    69   STIdxIterAcc iter(scantable_, cols);
    70 
    71   ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME");
    72   Vector<Double> timeSec = tcol->getColumn() * 86400.0;
    73   tcol->attach(scantable_->table(), "INTERVAL");
    74   Vector<Double> intervalSec = tcol->getColumn();
    75   delete tcol;
    76   ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION");
    77   Vector<Float> elevation = ecol->getColumn();
    78   delete ecol;
    79 
    80   ROArrayColumn<Float> specCol(scantable_->table(), "TSYS");
    81   ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA");
    82   ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID");
    83 
    84   // dummy Tsys: the following process doesn't need Tsys but RowAccumulator
    85   //             requires to set it with spectral data
    86   Vector<Float> tsys(1, 1.0);
    87 
    88   Double timeCen = 0.0;
    89   Float elCen = 0.0;
    90   uInt count = 0;
    91 
    92   while(!iter.pastEnd()) {
    93     Vector<uInt> rows = iter.getRows(SHARE);
    94     Vector<uInt> current = iter.current();
    95     //os_ << "current=" << current << LogIO::POST;
    96     uInt len = rows.nelements();
    97     if (len == 0) {
    98       iter.next();
    99       continue;
    100     }
    101     else if (len == 1) {
    102       STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
    103       uInt irow = rows[0];
    104       p->appenddata(0, 0, current[2], current[0], current[1],
    105                     freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow));
    106       iter.next();
    107       continue;
    108     }
    109    
    110     uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
    111     Vector<uChar> flag(nchan);
    112     Vector<Bool> bflag(nchan);
    113     Vector<Float> spec(nchan);
    114 
    115     Vector<Double> timeSep(len);
    116     for (uInt i = 0; i < len-1; i++) {
    117       timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
    118     }
    119     Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
    120     timeSep[len-1] = tMedian * 10000.0 ; // any large value
    121 
    122     uInt irow ;
    123     uInt jrow ;
    124     for (uInt i = 0; i < len; i++) {
    125       //os_ << "start row " << rows[i] << LogIO::POST;
    126       irow = rows[i];
    127       jrow = (i < len-1) ? rows[i+1] : rows[i];
    128       // accumulate data
    129       flagCol.get(irow, flag);
    130       convertArray(bflag, flag);
    131       specCol.get(irow, spec);
    132       if ( !allEQ(bflag,True) )
    133         acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
    134       timeCen += timeSec[irow];
    135       elCen += elevation[irow];
    136       count++;
    137 
    138       // check time gap
    139       double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
    140       if ( gap > 5.0 ) {
    141         if ( acc.state() ) {
    142           acc.replaceNaN() ;
    143 //           const Vector<Bool> &msk = acc.getMask();
    144 //           convertArray(flag, !msk);
    145 //           for (uInt k = 0; k < nchan; ++k) {
    146 //             uChar userFlag = 1 << 7;
    147 //             if (msk[k]==True) userFlag = 0 << 7;
    148 //             flag(k) = userFlag;
    149 //           }
    150           timeCen /= (Double)count * 86400.0; // sec->day
    151           elCen /= (Float)count;
    152           STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
    153           p->appenddata(0, 0, current[2], current[0], current[1],
    154                         freqidCol(irow), timeCen, elCen, acc.getSpectrum());
    155         }
    156         acc.reset() ;
    157         timeCen = 0.0;
    158         elCen = 0.0;
    159         count = 0;
    160       }
    161     }
    162 
    163     iter.next() ;
    164     //os_ << "end " << current << LogIO::POST;
    165   }
     64  STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
     65  p->appenddata(scanno, cycleno, beamno, ifno, polno,
     66                freqid, time, elevation, any_data);
    16667}
    16768

  • trunk/src/STCalTsys.h

    r2786 r2915  
    4444private:
    4545  void setupSelector(const STSelector &sel);
    46   void fillCalTable();
     46  virtual void appenddata(casa::uInt scanno, casa::uInt cycleno,
     47                          casa::uInt beamno, casa::uInt ifno, casa::uInt polno,
     48                          casa::uInt freqid, casa::Double time, casa::Float elevation,
     49                          casa::Vector<casa::Float> any_data);
    4750
    4851  vector<int> iflist_;

  • trunk/src/STCalibration.cpp

    r2786 r2915  
    1313#include "STCalibration.h"
    1414#include "STCalSkyTable.h"
     15#include "RowAccumulator.h"
     16#include "STIdxIter.h"
    1517
    1618using namespace casa;
    1719
    1820namespace asap {
    19 STCalibration::STCalibration(CountedPtr<Scantable> &s)
    20   : scantable_(s)
     21STCalibration::STCalibration(CountedPtr<Scantable> &s, const String target_column)
     22  : scantable_(s),
     23    target_column_(target_column)
    2124{
    2225}
    … …  
    3336}
    3437
     38void STCalibration::fillCalTable()
     39{
     40  RowAccumulator acc(W_TINT);
     41 
     42  vector<string> cols(3);
     43  cols[0] = "IFNO";
     44  cols[1] = "POLNO";
     45  cols[2] = "BEAMNO";
     46  STIdxIterAcc iter(scantable_, cols);
     47
     48  ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME");
     49  Vector<Double> timeSec = tcol->getColumn() * 86400.0;
     50  tcol->attach(scantable_->table(), "INTERVAL");
     51  Vector<Double> intervalSec = tcol->getColumn();
     52  delete tcol;
     53  ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION");
     54  Vector<Float> elevation = ecol->getColumn();
     55  delete ecol;
     56
     57  ROArrayColumn<Float> specCol(scantable_->table(), target_column_);
     58  ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA");
     59  ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID");
     60
     61  // dummy Tsys: the following process doesn't need Tsys but RowAccumulator
     62  //             requires to set it with spectral data
     63  Vector<Float> tsys(1, 1.0);
     64
     65  Double timeCen = 0.0;
     66  Float elCen = 0.0;
     67  uInt count = 0;
     68
     69  while(!iter.pastEnd()) {
     70    Vector<uInt> rows = iter.getRows(SHARE);
     71    Vector<uInt> current = iter.current();
     72    //os_ << "current=" << current << LogIO::POST;
     73    uInt len = rows.nelements();
     74    if (len == 0) {
     75      iter.next();
     76      continue;
     77    }
     78    else if (len == 1) {
     79      uInt irow = rows[0];
     80      appenddata(0, 0, current[2], current[0], current[1],
     81                 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow));
     82      iter.next();
     83      continue;
     84    }
     85   
     86    uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
     87    Vector<uChar> flag(nchan);
     88    Vector<Bool> bflag(nchan);
     89    Vector<Float> spec(nchan);
     90
     91    Vector<Double> timeSep(len);
     92    for (uInt i = 0; i < len-1; i++) {
     93      timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
     94    }
     95    Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
     96    timeSep[len-1] = tMedian * 10000.0 ; // any large value
     97
     98    uInt irow ;
     99    uInt jrow ;
     100    for (uInt i = 0; i < len; i++) {
     101      //os_ << "start row " << rows[i] << LogIO::POST;
     102      irow = rows[i];
     103      jrow = (i < len-1) ? rows[i+1] : rows[i];
     104      // accumulate data
     105      flagCol.get(irow, flag);
     106      convertArray(bflag, flag);
     107      specCol.get(irow, spec);
     108      if ( !allEQ(bflag,True) )
     109        acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
     110      timeCen += timeSec[irow];
     111      elCen += elevation[irow];
     112      count++;
     113
     114      // check time gap
     115      double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
     116      if ( gap > 1.1 ) {
     117        if ( acc.state() ) {
     118          acc.replaceNaN() ;
     119//           const Vector<Bool> &msk = acc.getMask();
     120//           convertArray(flag, !msk);
     121//           for (uInt k = 0; k < nchan; ++k) {
     122//             uChar userFlag = 1 << 7;
     123//             if (msk[k]==True) userFlag = 0 << 7;
     124//             flag(k) = userFlag;
     125//           }
     126          timeCen /= (Double)count * 86400.0; // sec->day
     127          elCen /= (Float)count;
     128          appenddata(0, 0, current[2], current[0], current[1],
     129                     freqidCol(irow), timeCen, elCen, acc.getSpectrum());
     130        }
     131        acc.reset() ;
     132        timeCen = 0.0;
     133        elCen = 0.0;
     134        count = 0;
     135      }
     136    }
     137
     138    iter.next() ;
     139    //os_ << "end " << current << LogIO::POST;
     140  } 
    35141}
     142 
     143}

  • trunk/src/STCalibration.h

    r2786 r2915  
    3535class STCalibration {
    3636public:
    37   STCalibration(casa::CountedPtr<Scantable> &s);
     37  STCalibration(casa::CountedPtr<Scantable> &s, const casa::String target_column);
    3838
    3939  void calibrate();
    … …  
    4949protected:
    5050  virtual void setupSelector(const STSelector &sel) = 0;
    51   virtual void fillCalTable() = 0;
     51  virtual void fillCalTable();
     52  virtual void appenddata(casa::uInt scanno, casa::uInt cycleno,
     53                          casa::uInt beamno, casa::uInt ifno, casa::uInt polno,
     54                          casa::uInt freqid, casa::Double time, casa::Float elevation,
     55                          casa::Vector<casa::Float> any_data) = 0;
    5256
    5357  STSelector sel_;
    … …  
    5660  casa::LogIO os_;
    5761  casa::Record options_;
     62  const casa::String target_column_;
    5863};
    59 
     64 
    6065}
    6166#endif
Note: See TracChangeset for help on using the changeset viewer.