source: trunk/src/STCalTsys.cpp@ 2706

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

New Development: Yes

JIRA Issue: Yes CAS-4770, 4771, 4772

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

First version of calibration classes for sky calibration (STCalSkyPSAlma)
and Tsys calibration (STCalTsys), which are derived from base class
(STCalibration).

The caltable classes, STCalSky and STCalTsys, are renamed as STCalSkyTable
and STCalTsysTable, respectively.


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