source: trunk/src/STCalibration.cpp@ 2955

Last change on this file since 2955 was 2955, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6585, CAS-6571

Ready for Test: Yes

Interface Changes: Yes/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...

Added FLAGTRA column to sky and Tsys caltable.
To support previous versions of caltabls that don't have FLAGTRA
column, the code automatically add FLAGTRA column and initialize
it with initial value (uChar)0 if no FLAGTRA exists.


File size: 4.3 KB
Line 
1//
2// C++ Implementation: STCalibration
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 "STCalibration.h"
14#include "STCalSkyTable.h"
15#include "RowAccumulator.h"
16#include "STIdxIter.h"
17
18using namespace casa;
19
20namespace asap {
21STCalibration::STCalibration(CountedPtr<Scantable> &s, const String target_column)
22 : scantable_(s),
23 target_column_(target_column)
24{
25}
26
27void STCalibration::calibrate()
28{
29 STSelector selOrg = scantable_->getSelection();
30 setupSelector(selOrg);
31 scantable_->setSelection(sel_);
32
33 fillCalTable();
34
35 scantable_->setSelection(selOrg);
36}
37
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 STIdxIter2 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 Record current = iter.currentValue();
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.asuInt("BEAMNO"), current.asuInt("IFNO"), current.asuInt("POLNO"),
81 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow),
82 flagCol(irow));
83 iter.next();
84 continue;
85 }
86
87 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
88 Vector<uChar> flag(nchan);
89 Vector<Bool> bflag(nchan);
90 Vector<Float> spec(nchan);
91
92 Vector<Double> timeSep(len);
93 for (uInt i = 0; i < len-1; i++) {
94 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
95 }
96 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
97 timeSep[len-1] = tMedian * 10000.0 ; // any large value
98
99 uInt irow ;
100 uInt jrow ;
101 for (uInt i = 0; i < len; i++) {
102 //os_ << "start row " << rows[i] << LogIO::POST;
103 irow = rows[i];
104 jrow = (i < len-1) ? rows[i+1] : rows[i];
105 // accumulate data
106 flagCol.get(irow, flag);
107 convertArray(bflag, flag);
108 specCol.get(irow, spec);
109 if ( !allEQ(bflag,True) )
110 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
111 timeCen += timeSec[irow];
112 elCen += elevation[irow];
113 count++;
114
115 // check time gap
116 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
117 if ( gap > 1.1 ) {
118 if ( acc.state() ) {
119 acc.replaceNaN() ;
120// const Vector<Bool> &msk = acc.getMask();
121// convertArray(flag, !msk);
122// for (uInt k = 0; k < nchan; ++k) {
123// uChar userFlag = 1 << 7;
124// if (msk[k]==True) userFlag = 0 << 7;
125// flag(k) = userFlag;
126// }
127 timeCen /= (Double)count * 86400.0; // sec->day
128 elCen /= (Float)count;
129 const Vector<Bool> &mask = acc.getMask();
130 Vector<uChar> flag(mask.shape(), (uChar)0);
131 const uChar userFlag = 1 << 7;
132 for (uInt k = 0; k < flag.nelements(); ++k) {
133 if (mask[k] == True)
134 flag[k] = userFlag;
135 }
136 appenddata(0, 0, current.asuInt("BEAMNO"), current.asuInt("IFNO"),
137 current.asuInt("POLNO"),
138 freqidCol(irow), timeCen, elCen,
139 acc.getSpectrum(), flag);
140 }
141 acc.reset() ;
142 timeCen = 0.0;
143 elCen = 0.0;
144 count = 0;
145 }
146 }
147
148 iter.next() ;
149 //os_ << "end " << current << LogIO::POST;
150 }
151}
152
153}
Note: See TracBrowser for help on using the repository browser.