source: branches/casa-release-4_3/src/STCalibration.cpp@ 3120

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

New Development: No

JIRA Issue: Yes CAS-6582, 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...

Take into account row and channel flag when generating caltables.


File size: 4.4 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 ROScalarColumn<uInt> flagrowCol(scantable_->table(), "FLAGROW");
61 Vector<uInt> flagrow = flagrowCol.getColumn();
62
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);
66
67 Double timeCen = 0.0;
68 Float elCen = 0.0;
69 uInt count = 0;
70
71 while(!iter.pastEnd()) {
72 Vector<uInt> rows = iter.getRows(SHARE);
73 Record current = iter.currentValue();
74 //os_ << "current=" << current << LogIO::POST;
75 uInt len = rows.nelements();
76 if (len == 0) {
77 iter.next();
78 continue;
79 }
80 else if (len == 1) {
81 uInt irow = rows[0];
82 if (flagrow[irow] == 0) {
83 appenddata(0, 0, current.asuInt("BEAMNO"), current.asuInt("IFNO"), current.asuInt("POLNO"),
84 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow),
85 flagCol(irow));
86 }
87 iter.next();
88 continue;
89 }
90
91 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0]));
92 Vector<uChar> flag(nchan);
93 Vector<Bool> bflag(nchan);
94 Vector<Float> spec(nchan);
95
96 Vector<Double> timeSep(len);
97 for (uInt i = 0; i < len-1; i++) {
98 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ;
99 }
100 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2)));
101 timeSep[len-1] = tMedian * 10000.0 ; // any large value
102
103 uInt irow ;
104 uInt jrow ;
105 for (uInt i = 0; i < len; i++) {
106 //os_ << "start row " << rows[i] << LogIO::POST;
107 irow = rows[i];
108 jrow = (i < len-1) ? rows[i+1] : rows[i];
109 // accumulate data
110 flagCol.get(irow, flag);
111 convertArray(bflag, flag);
112 specCol.get(irow, spec);
113 if ( !allEQ(bflag,True) && flagrow[irow] == 0 )
114 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ;
115 timeCen += timeSec[irow];
116 elCen += elevation[irow];
117 count++;
118
119 // check time gap
120 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]);
121 if ( gap > 1.1 ) {
122 if ( acc.state() ) {
123 acc.replaceNaN() ;
124// const Vector<Bool> &msk = acc.getMask();
125// convertArray(flag, !msk);
126// for (uInt k = 0; k < nchan; ++k) {
127// uChar userFlag = 1 << 7;
128// if (msk[k]==True) userFlag = 0 << 7;
129// flag(k) = userFlag;
130// }
131 timeCen /= (Double)count * 86400.0; // sec->day
132 elCen /= (Float)count;
133 const Vector<Bool> &mask = acc.getMask();
134 Vector<uChar> flag(mask.shape(), (uChar)0);
135 const uChar userFlag = 1 << 7;
136 for (uInt k = 0; k < flag.nelements(); ++k) {
137 if (mask[k] == False)
138 flag[k] = userFlag;
139 }
140 appenddata(0, 0, current.asuInt("BEAMNO"), current.asuInt("IFNO"),
141 current.asuInt("POLNO"),
142 freqidCol(irow), timeCen, elCen,
143 acc.getSpectrum(), flag);
144 }
145 acc.reset() ;
146 timeCen = 0.0;
147 elCen = 0.0;
148 count = 0;
149 }
150 }
151
152 iter.next() ;
153 //os_ << "end " << current << LogIO::POST;
154 }
155}
156
157}
Note: See TracBrowser for help on using the repository browser.