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