source: trunk/src/STCalSkyPSAlma.cpp@ 2835

Last change on this file since 2835 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.4 KB
Line 
1//
2// C++ Implementation: STCalSkyPSAlma
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/Logging/LogIO.h>
16
17#include "STSelector.h"
18#include "STCalSkyPSAlma.h"
19#include "RowAccumulator.h"
20#include "STIdxIter.h"
21#include "STDefs.h"
22#include <atnf/PKSIO/SrcType.h>
23
24using namespace std;
25using namespace casa;
26
27namespace asap {
28STCalSkyPSAlma::STCalSkyPSAlma(CountedPtr<Scantable> &s)
29 : STCalibration(s)
30{
31 applytable_ = new STCalSkyTable(*s, "PSALMA");
32}
33
34void STCalSkyPSAlma::setupSelector(const STSelector &sel)
35{
36 sel_ = sel;
37 vector<int> types = sel_.getTypes();
38 if (types.size() == 0) {
39 types.resize(1);
40 types[0] = SrcType::PSOFF;
41 sel_.setTypes(types);
42 }
43 else if (find(types.begin(), types.end(), SrcType::PSOFF) == types.end()) {
44 LogIO os(LogOrigin("STCalSkyPSAlma", "setupSelector", WHERE));
45 os << LogIO::SEVERE << "Selection contains no data." << LogIO::EXCEPTION;
46 }
47 else {
48 types.resize(1);
49 types[0] = SrcType::PSOFF;
50 sel_.setTypes(types);
51 }
52}
53
54void STCalSkyPSAlma::fillCalTable()
55{
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 }
156}
157
158}
Note: See TracBrowser for help on using the repository browser.