source: branches/plotter2/src/STCalTsys.cpp@ 2841

Last change on this file since 2841 was 2786, 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: test_sdcal2

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on handling data selection of input Scantable object.


File size: 4.7 KB
Line 
1//
2// C++ Implementation: STCalTsys
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp> (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <vector>
14
15#include <casa/Arrays/ArrayMath.h>
16#include <casa/Logging/LogIO.h>
17
18#include "STSelector.h"
19#include "STCalTsys.h"
20#include "RowAccumulator.h"
21#include "STIdxIter.h"
22#include "STDefs.h"
23#include <atnf/PKSIO/SrcType.h>
24
25using namespace std;
26using namespace casa;
27
28namespace asap {
29 STCalTsys::STCalTsys(CountedPtr<Scantable> &s, vector<int> &iflist)
30 : STCalibration(s),
31 iflist_(iflist)
32{
33 applytable_ = new STCalTsysTable(*s);
34}
35
36void STCalTsys::setupSelector(const STSelector &sel)
37{
38 sel_ = sel;
39 vector<int> ifnos = sel_.getIFs();
40 if (ifnos.size() > 0) {
41 int nif = 0;
42 int nifOrg = iflist_.size();
43 vector<int> iflistNew(iflist_);
44 for (int i = 0; i < nifOrg; i++) {
45 if (find(ifnos.begin(), ifnos.end(), iflist_[i]) != ifnos.end()) {
46 iflistNew[nif] = iflist_[i];
47 ++nif;
48 }
49 }
50 if (nif == 0) {
51 LogIO os(LogOrigin("STCalTsys", "setupSelector", WHERE));
52 os << LogIO::SEVERE << "Selection contains no data." << LogIO::EXCEPTION;
53 }
54 sel_.setIFs(iflistNew);
55 }
56 else {
57 sel_.setIFs(iflist_);
58 }
59}
60
61void STCalTsys::fillCalTable()
62{
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 }
166}
167
168}
Note: See TracBrowser for help on using the repository browser.