source: trunk/src/STApplyCal.cpp@ 2747

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

New Development: No

JIRA Issue: Yes CAS-4770

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

Defined python interface for calibration that supports both
on-the-fly and interferometry-style (generate caltable and apply)
calibration.

File size: 15.3 KB
Line 
1//
2// C++ Implementation: STApplyCal
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#include <assert.h>
13
14#include <casa/Arrays/Array.h>
15#include <casa/Arrays/Vector.h>
16#include <casa/Arrays/Matrix.h>
17#include <casa/Arrays/ArrayIO.h>
18#include <casa/Arrays/ArrayMath.h>
19#include <casa/BasicSL/String.h>
20#include <casa/Logging/LogIO.h>
21#include <casa/Exceptions/Error.h>
22#include <casa/Utilities/CountedPtr.h>
23#include <casa/Utilities/Sort.h>
24#include <tables/Tables/Table.h>
25
26#include "Scantable.h"
27#include "STApplyCal.h"
28#include "STApplyTable.h"
29#include "STCalTsysTable.h"
30#include "STCalSkyTable.h"
31#include "STCalEnum.h"
32#include "STIdxIter.h"
33#include "Calibrator.h"
34#include "PSAlmaCalibrator.h"
35#include "Interpolator1D.h"
36#include "NearestInterpolator1D.h"
37#include "BufferedLinearInterpolator1D.h"
38#include "PolynomialInterpolator1D.h"
39#include "CubicSplineInterpolator1D.h"
40#include <atnf/PKSIO/SrcType.h>
41
42
43using namespace casa;
44using namespace std;
45
46namespace asap {
47
48STApplyCal::STApplyCal()
49{
50 init();
51}
52
53STApplyCal::STApplyCal(CountedPtr<Scantable> target)
54 : target_(target)
55{
56 init();
57}
58
59STApplyCal::~STApplyCal()
60{
61}
62
63void STApplyCal::init()
64{
65 caltype_ = STCalEnum::NoType;
66 doTsys_ = False;
67 iTime_ = STCalEnum::DefaultInterpolation;
68 iFreq_ = STCalEnum::DefaultInterpolation;
69}
70
71void STApplyCal::reset()
72{
73 // call init
74 init();
75
76 // clear apply tables
77 // do not delete object here
78 skytable_.resize(0);
79 tsystable_.resize(0);
80
81 // clear mapping for Tsys transfer
82 spwmap_.clear();
83
84 // reset selector
85 sel_.reset();
86
87 // delete interpolators
88 interpolatorT_ = 0;
89 interpolatorS_ = 0;
90 interpolatorF_ = 0;
91
92 // clear working scantable
93 work_ = 0;
94
95 // clear calibrator
96 calibrator_ = 0;
97}
98
99void STApplyCal::completeReset()
100{
101 reset();
102 target_ = 0;
103}
104
105void STApplyCal::setTarget(CountedPtr<Scantable> target)
106{
107 target_ = target;
108}
109
110void STApplyCal::setTarget(const String &name)
111{
112 // always create PlainTable
113 target_ = new Scantable(name, Table::Plain);
114}
115
116void STApplyCal::push(STCalSkyTable *table)
117{
118 os_.origin(LogOrigin("STApplyCal","push",WHERE));
119 skytable_.push_back(table);
120 STCalEnum::CalType caltype = STApplyTable::getCalType(table);
121 os_ << "caltype=" << caltype << LogIO::POST;
122 if (caltype_ == STCalEnum::NoType ||
123 caltype_ == STCalEnum::DefaultType ||
124 caltype_ == STCalEnum::CalTsys) {
125 caltype_ = caltype;
126 }
127 os_ << "caltype_=" << caltype_ << LogIO::POST;
128}
129
130void STApplyCal::push(STCalTsysTable *table)
131{
132 tsystable_.push_back(table);
133 doTsys_ = True;
134}
135
136void STApplyCal::setTimeInterpolation(STCalEnum::InterpolationType itype, Int order)
137{
138 iTime_ = itype;
139 order_ = order;
140}
141
142void STApplyCal::setFrequencyInterpolation(STCalEnum::InterpolationType itype, Int order)
143{
144 iFreq_ = itype;
145 order_ = order;
146}
147
148void STApplyCal::setTsysTransfer(uInt from, Vector<uInt> to)
149{
150 os_.origin(LogOrigin("STApplyCal","setTsysTransfer",WHERE));
151 os_ << "from=" << from << ", to=" << to << LogIO::POST;
152 map<uInt, Vector<uInt> >::iterator i = spwmap_.find(from);
153 if (i == spwmap_.end()) {
154 spwmap_.insert(pair<uInt, Vector<uInt> >(from, to));
155 }
156 else {
157 Vector<uInt> toNew = i->second;
158 spwmap_.erase(i);
159 uInt k = toNew.nelements();
160 toNew.resize(k+to.nelements(), True);
161 for (uInt i = 0; i < to.nelements(); i++)
162 toNew[i+k] = to[i];
163 spwmap_.insert(pair<uInt, Vector<uInt> >(from, toNew));
164 }
165}
166
167void STApplyCal::apply(Bool insitu, Bool filltsys)
168{
169 os_.origin(LogOrigin("STApplyCal","apply",WHERE));
170 // calibrator
171 if (caltype_ == STCalEnum::CalPSAlma)
172 calibrator_ = new PSAlmaCalibrator();
173
174 // interpolator
175 initInterpolator();
176
177 // select data
178 sel_.reset();
179 if (caltype_ == STCalEnum::CalPSAlma ||
180 caltype_ == STCalEnum::CalPS) {
181 sel_.setTypes(vector<int>(1,(int)SrcType::PSON));
182 }
183 target_->setSelection(sel_);
184
185 //os_ << "sel_.print()=" << sel_.print() << LogIO::POST;
186
187 // working data
188 if (insitu)
189 work_ = target_;
190 else
191 work_ = new Scantable(*target_, false);
192
193 //os_ << "work_->nrow()=" << work_->nrow() << LogIO::POST;
194
195 // list of apply tables for sky calibration
196 Vector<uInt> skycalList;
197 uInt numSkyCal = 0;
198 uInt nrowSky = 0;
199
200 // list of apply tables for Tsys calibration
201 for (uInt i = 0 ; i < skytable_.size(); i++) {
202 STCalEnum::CalType caltype = STApplyTable::getCalType(skytable_[i]);
203 if (caltype == caltype_) {
204 skycalList.resize(numSkyCal+1, True);
205 skycalList[numSkyCal] = i;
206 numSkyCal++;
207 nrowSky += skytable_[i]->nrow();
208 }
209 }
210
211
212 vector<string> cols( 3 ) ;
213 cols[0] = "BEAMNO" ;
214 cols[1] = "POLNO" ;
215 cols[2] = "IFNO" ;
216 CountedPtr<STIdxIter> iter = new STIdxIterAcc(work_, cols) ;
217 while (!iter->pastEnd()) {
218 Vector<uInt> ids = iter->current();
219 Vector<uInt> rows = iter->getRows(SHARE);
220 if (rows.nelements() > 0)
221 doapply(ids[0], ids[2], ids[1], rows, skycalList);
222 iter->next();
223 }
224
225 target_->unsetSelection();
226}
227
228void STApplyCal::doapply(uInt beamno, uInt ifno, uInt polno,
229 Vector<uInt> &rows,
230 Vector<uInt> &skylist,
231 Bool filltsys)
232{
233 os_.origin(LogOrigin("STApplyCal","doapply",WHERE));
234 Bool doTsys = doTsys_;
235
236 STSelector sel;
237 vector<int> id(1);
238 id[0] = beamno;
239 sel.setBeams(id);
240 id[0] = ifno;
241 sel.setIFs(id);
242 id[0] = polno;
243 sel.setPolarizations(id);
244
245 // apply selection to apply tables
246 uInt nrowSky = 0;
247 uInt nrowTsys = 0;
248 for (uInt i = 0; i < skylist.nelements(); i++) {
249 skytable_[skylist[i]]->setSelection(sel);
250 nrowSky += skytable_[skylist[i]]->nrow();
251 os_ << "nrowSky=" << nrowSky << LogIO::POST;
252 }
253 uInt nchanTsys = 0;
254 Vector<Double> ftsys;
255 uInt tsysifno = getIFForTsys(ifno);
256 os_ << "tsysifno=" << tsysifno << LogIO::POST;
257 if (tsystable_.size() == 0) {
258 os_.origin(LogOrigin("STApplyTable", "doapply", WHERE));
259 os_ << "No Tsys tables are given. Skip Tsys calibratoin." << LogIO::POST;
260 doTsys = False;
261 }
262 else if (tsysifno == (uInt)-1) {
263 os_.origin(LogOrigin("STApplyTable", "doapply", WHERE));
264 os_ << "No corresponding Tsys for IFNO " << ifno << ". Skip Tsys calibration" << LogIO::POST;
265 doTsys = False;
266 }
267 else {
268 nchanTsys = tsystable_[0]->nchan(tsysifno);
269 ftsys = tsystable_[0]->getBaseFrequency(0);
270 interpolatorF_->setX(ftsys.data(), nchanTsys);
271 id[0] = (int)tsysifno;
272 sel.setIFs(id);
273 for (uInt i = 0; i < tsystable_.size() ; i++) {
274 tsystable_[i]->setSelection(sel);
275 nrowTsys += tsystable_[i]->nrow();
276 }
277 }
278
279 uInt nchanSp = skytable_[skylist[0]]->nchan(ifno);
280 Vector<Double> timeSky(nrowSky);
281 Matrix<Float> spoff(nchanSp, nrowSky);
282 Vector<Float> iOff(nchanSp);
283 nrowSky = 0;
284 for (uInt i = 0 ; i < skylist.nelements(); i++) {
285 STCalSkyTable *p = skytable_[skylist[i]];
286 Vector<Double> t = p->getTime();
287 Matrix<Float> sp = p->getSpectra();
288 for (uInt j = 0; j < t.nelements(); j++) {
289 timeSky[nrowSky] = t[j];
290 spoff.column(nrowSky) = sp.column(j);
291 nrowSky++;
292 }
293 }
294
295 Vector<uInt> skyIdx = timeSort(timeSky);
296
297 Double *xa = new Double[skyIdx.nelements()];
298 Float *ya = new Float[skyIdx.nelements()];
299 IPosition ipos(1, skyIdx.nelements());
300 Vector<Double> timeSkySorted(ipos, xa, TAKE_OVER);
301 Vector<Float> tmpOff(ipos, ya, TAKE_OVER);
302 for (uInt i = 0 ; i < skyIdx.nelements(); i++) {
303 timeSkySorted[i] = timeSky[skyIdx[i]];
304 }
305
306 interpolatorS_->setX(xa, skyIdx.nelements());
307
308 Vector<uInt> tsysIdx;
309 Vector<Double> timeTsys(nrowTsys);
310 Matrix<Float> tsys;
311 Vector<Double> timeTsysSorted;
312 Vector<Float> tmpTsys;
313 if (doTsys) {
314 os_ << "doTsys" << LogIO::POST;
315 timeTsys.resize(nrowTsys);
316 tsys.resize(nchanTsys, nrowTsys);
317 nrowTsys = 0;
318 for (uInt i = 0 ; i < tsystable_.size(); i++) {
319 STCalTsysTable *p = tsystable_[i];
320 Vector<Double> t = p->getTime();
321 Matrix<Float> ts = p->getTsys();
322 for (uInt j = 0; j < t.nelements(); j++) {
323 timeTsys[nrowTsys] = t[j];
324 tsys.column(nrowTsys) = ts.column(j);
325 nrowTsys++;
326 }
327 }
328 tsysIdx = timeSort(timeTsys);
329
330 Double *xb = new Double[tsysIdx.nelements()];
331 Float *yb = new Float[tsysIdx.nelements()];
332 IPosition ipos(1, tsysIdx.nelements());
333 timeTsysSorted.takeStorage(ipos, xb, TAKE_OVER);
334 tmpTsys.takeStorage(ipos, yb, TAKE_OVER);
335 for (uInt i = 0 ; i < tsysIdx.nelements(); i++) {
336 timeTsysSorted[i] = timeTsys[tsysIdx[i]];
337 }
338 interpolatorT_->setX(xb, tsysIdx.nelements());
339 }
340
341 Table tab = work_->table();
342 ArrayColumn<Float> spCol(tab, "SPECTRA");
343 ArrayColumn<Float> tsysCol(tab, "TSYS");
344 ScalarColumn<Double> timeCol(tab, "TIME");
345 Vector<Float> on;
346 for (uInt i = 0; i < rows.nelements(); i++) {
347 os_ << "start i = " << i << " (row = " << rows[i] << ")" << LogIO::POST;
348 uInt irow = rows[i];
349
350 // target spectral data
351 on = spCol(irow);
352 calibrator_->setSource(on);
353
354 // interpolation
355 Double t0 = timeCol(irow);
356 for (uInt ichan = 0; ichan < nchanSp; ichan++) {
357 Vector<Float> spOffSlice = spoff.row(ichan);
358 //os_ << "spOffSlice = " << spOffSlice << LogIO::POST;
359 for (uInt j = 0; j < skyIdx.nelements(); j++) {
360 tmpOff[j] = spOffSlice[skyIdx[j]];
361 }
362 interpolatorS_->setY(ya, skyIdx.nelements());
363 iOff[ichan] = interpolatorS_->interpolate(t0);
364 }
365 //os_ << "iOff=" << iOff << LogIO::POST;
366 calibrator_->setReference(iOff);
367
368 Float *Y = new Float[nchanSp];
369 Vector<Float> iTsys(IPosition(1,nchanSp), Y, TAKE_OVER);
370 if (doTsys) {
371 // Tsys correction
372 Float *yt = new Float[nchanTsys];
373 Vector<Float> iTsysT(IPosition(1,nchanTsys), yt, TAKE_OVER);
374 Float *yb = tmpTsys.data();
375 for (uInt ichan = 0; ichan < nchanTsys; ichan++) {
376 Vector<Float> tsysSlice = tsys.row(ichan);
377 for (uInt j = 0; j < tsysIdx.nelements(); j++) {
378 tmpTsys[j] = tsysSlice[tsysIdx[j]];
379 }
380 interpolatorT_->setY(yb, tsysIdx.nelements());
381 iTsysT[ichan] = interpolatorT_->interpolate(t0);
382 }
383 if (nchanSp == 1) {
384 // take average
385 iTsys[0] = mean(iTsysT);
386 }
387 else {
388 // interpolation on frequency axis
389 Vector<Double> fsp = getBaseFrequency(rows[i]);
390 interpolatorF_->setY(yt, nchanTsys);
391 for (uInt ichan = 0; ichan < nchanSp; ichan++) {
392 iTsys[ichan] = interpolatorF_->interpolate(fsp[ichan]);
393 }
394 }
395 }
396 else {
397 iTsys = 1.0;
398 }
399 //os_ << "iTsys=" << iTsys << LogIO::POST;
400 calibrator_->setScaler(iTsys);
401
402 // do calibration
403 calibrator_->calibrate();
404
405 // update table
406 //os_ << "calibrated=" << calibrator_->getCalibrated() << LogIO::POST;
407 spCol.put(irow, calibrator_->getCalibrated());
408 if (filltsys)
409 tsysCol.put(irow, iTsys);
410 }
411
412
413 // reset selection on apply tables
414 for (uInt i = 0; i < skylist.nelements(); i++)
415 skytable_[i]->unsetSelection();
416 for (uInt i = 0; i < tsystable_.size(); i++)
417 tsystable_[i]->unsetSelection();
418
419
420 // reset interpolator
421 interpolatorS_->reset();
422 interpolatorF_->reset();
423 interpolatorT_->reset();
424}
425
426Vector<uInt> STApplyCal::timeSort(Vector<Double> &t)
427{
428 Sort sort;
429 sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
430 Vector<uInt> idx;
431 sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
432 return idx;
433}
434
435uInt STApplyCal::getIFForTsys(uInt to)
436{
437 for (map<casa::uInt, Vector<uInt> >::iterator i = spwmap_.begin();
438 i != spwmap_.end(); i++) {
439 Vector<uInt> tolist = i->second;
440 os_ << "from=" << i->first << ": tolist=" << tolist << LogIO::POST;
441 for (uInt j = 0; j < tolist.nelements(); j++) {
442 if (tolist[j] == to)
443 return i->first;
444 }
445 }
446 return (uInt)-1;
447}
448
449void STApplyCal::save(const String &name)
450{
451 assert(!work_.null());
452
453 work_->setSelection(sel_);
454 work_->makePersistent(name);
455 work_->unsetSelection();
456}
457
458Vector<Double> STApplyCal::getBaseFrequency(uInt whichrow)
459{
460 assert(whichrow <= (uInt)work_->nrow());
461 ROTableColumn col(work_->table(), "IFNO");
462 uInt ifno = col.asuInt(whichrow);
463 col.attach(work_->table(), "FREQ_ID");
464 uInt freqid = col.asuInt(whichrow);
465 uInt nc = work_->nchan(ifno);
466 STFrequencies ftab = work_->frequencies();
467 Double rp, rf, inc;
468 ftab.getEntry(rp, rf, inc, freqid);
469 Vector<Double> r(nc);
470 indgen(r, rf-rp*inc, inc);
471 return r;
472}
473
474void STApplyCal::initInterpolator()
475{
476 os_.origin(LogOrigin("STApplyCal","initInterpolator",WHERE));
477 int order = (order_ > 0) ? order_ : 1;
478 switch (iTime_) {
479 case STCalEnum::NearestInterpolation:
480 {
481 os_ << "use NearestInterpolator in time axis" << LogIO::POST;
482 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
483 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
484 break;
485 }
486 case STCalEnum::LinearInterpolation:
487 {
488 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
489 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
490 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
491 break;
492 }
493 case STCalEnum::CubicSplineInterpolation:
494 {
495 os_ << "use CubicSplineInterpolator in time axis" << LogIO::POST;
496 interpolatorS_ = new CubicSplineInterpolator1D<Double, Float>();
497 interpolatorT_ = new CubicSplineInterpolator1D<Double, Float>();
498 break;
499 }
500 case STCalEnum::PolynomialInterpolation:
501 {
502 os_ << "use PolynomialInterpolator in time axis" << LogIO::POST;
503 if (order == 0) {
504 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
505 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
506 }
507 else {
508 interpolatorS_ = new PolynomialInterpolator1D<Double, Float>();
509 interpolatorT_ = new PolynomialInterpolator1D<Double, Float>();
510 interpolatorS_->setOrder(order);
511 interpolatorT_->setOrder(order);
512 }
513 break;
514 }
515 default:
516 {
517 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
518 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
519 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
520 break;
521 }
522 }
523
524 switch (iFreq_) {
525 case STCalEnum::NearestInterpolation:
526 {
527 os_ << "use NearestInterpolator in frequency axis" << LogIO::POST;
528 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
529 break;
530 }
531 case STCalEnum::LinearInterpolation:
532 {
533 os_ << "use BufferedLinearInterpolator in frequency axis" << LogIO::POST;
534 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
535 break;
536 }
537 case STCalEnum::CubicSplineInterpolation:
538 {
539 os_ << "use CubicSplineInterpolator in frequency axis" << LogIO::POST;
540 interpolatorF_ = new CubicSplineInterpolator1D<Double, Float>();
541 break;
542 }
543 case STCalEnum::PolynomialInterpolation:
544 {
545 os_ << "use PolynomialInterpolator in frequency axis" << LogIO::POST;
546 if (order == 0) {
547 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
548 }
549 else {
550 interpolatorF_ = new PolynomialInterpolator1D<Double, Float>();
551 interpolatorF_->setOrder(order);
552 }
553 break;
554 }
555 default:
556 {
557 os_ << "use LinearInterpolator in frequency axis" << LogIO::POST;
558 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
559 break;
560 }
561 }
562}
563}
Note: See TracBrowser for help on using the repository browser.