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.


File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.