source: trunk/src/STCalTsys.cpp@ 2777

Last change on this file since 2777 was 2749, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-4770

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...

Bug fix for the case when only one row is selected in a certain iteration
loop.


File size: 4.1 KB
RevLine 
[2696]1//
2// C++ Implementation: STCalTsys
3//
4// Description:
5//
6//
[2703]7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp> (C) 2012
[2696]8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
[2703]13#include <vector>
14#include "STSelector.h"
[2696]15#include "STCalTsys.h"
[2703]16#include "RowAccumulator.h"
17#include "STIdxIter.h"
18#include "STDefs.h"
19#include <atnf/PKSIO/SrcType.h>
[2696]20
[2703]21#include <casa/Arrays/ArrayMath.h>
[2696]22
[2703]23using namespace std;
[2696]24using namespace casa;
25
26namespace asap {
[2703]27 STCalTsys::STCalTsys(CountedPtr<Scantable> &s, vector<int> &iflist)
28 : STCalibration(s),
29 iflist_(iflist)
[2696]30{
[2703]31 applytable_ = new STCalTsysTable(*s);
[2696]32}
33
[2742]34void STCalTsys::setupSelector()
[2696]35{
[2742]36 sel_.reset();
37 sel_.setIFs(iflist_);
[2696]38}
39
[2703]40void STCalTsys::fillCalTable()
[2696]41{
[2703]42 RowAccumulator acc(W_TINT);
43
44 vector<string> cols(3);
45 cols[0] = "IFNO";
46 cols[1] = "POLNO";
47 cols[2] = "BEAMNO";
48 STIdxIterAcc iter(scantable_, cols);
[2696]49
[2703]50 ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME");
51 Vector<Double> timeSec = tcol->getColumn() * 86400.0;
52 tcol->attach(scantable_->table(), "INTERVAL");
53 Vector<Double> intervalSec = tcol->getColumn();
54 delete tcol;
55 ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION");
56 Vector<Float> elevation = ecol->getColumn();
57 delete ecol;
[2696]58
[2703]59 ROArrayColumn<Float> specCol(scantable_->table(), "TSYS");
60 ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA");
[2720]61 ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID");
[2696]62
[2703]63 // dummy Tsys: the following process doesn't need Tsys but RowAccumulator
64 // requires to set it with spectral data
65 Vector<Float> tsys(1, 1.0);
[2696]66
[2703]67 Double timeCen = 0.0;
68 Float elCen = 0.0;
69 uInt count = 0;
[2696]70
[2703]71 while(!iter.pastEnd()) {
72 Vector<uInt> rows = iter.getRows(SHARE);
73 Vector<uInt> current = iter.current();
[2749]74 //os_ << "current=" << current << LogIO::POST;
[2703]75 uInt len = rows.nelements();
76 if (len == 0) {
77 iter.next();
78 continue;
79 }
[2749]80 else if (len == 1) {
81 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
82 uInt irow = rows[0];
83 p->appenddata(0, 0, current[2], current[0], current[1],
84 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow));
85 iter.next();
86 continue;
87 }
[2703]88
89 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
90 Vector<uChar> flag(nchan);
91 Vector<Bool> bflag(nchan);
92 Vector<Float> spec(nchan);
[2696]93
[2703]94 Vector<Double> timeSep(len);
95 for (uInt i = 0; i < len-1; i++) {
96 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
97 }
98 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
99 timeSep[len-1] = tMedian * 10000.0 ; // any large value
100
101 uInt irow ;
102 uInt jrow ;
103 for (uInt i = 0; i < len; i++) {
[2749]104 //os_ << "start row " << rows[i] << LogIO::POST;
[2703]105 irow = rows[i];
106 jrow = (i < len-1) ? rows[i+1] : rows[i];
107 // accumulate data
108 flagCol.get(irow, flag);
109 convertArray(bflag, flag);
110 specCol.get(irow, spec);
111 if ( !allEQ(bflag,True) )
112 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
113 timeCen += timeSec[irow];
114 elCen += elevation[irow];
115 count++;
116
117 // check time gap
118 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
119 if ( gap > 5.0 ) {
120 if ( acc.state() ) {
121 acc.replaceNaN() ;
122// const Vector<Bool> &msk = acc.getMask();
123// convertArray(flag, !msk);
124// for (uInt k = 0; k < nchan; ++k) {
125// uChar userFlag = 1 << 7;
126// if (msk[k]==True) userFlag = 0 << 7;
127// flag(k) = userFlag;
128// }
129 timeCen /= (Double)count * 86400.0; // sec->day
130 elCen /= (Float)count;
131 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_));
132 p->appenddata(0, 0, current[2], current[0], current[1],
[2720]133 freqidCol(irow), timeCen, elCen, acc.getSpectrum());
[2703]134 }
135 acc.reset() ;
136 timeCen = 0.0;
137 elCen = 0.0;
138 count = 0;
139 }
140 }
141
142 iter.next() ;
[2749]143 //os_ << "end " << current << LogIO::POST;
[2703]144 }
[2696]145}
146
147}
Note: See TracBrowser for help on using the repository browser.