source: trunk/src/STApplyCal.cpp@ 2896

Last change on this file since 2896 was 2862, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-4771

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: Describe your changes here...

Bug fix on sdcal2 sky calibration.
Support when Tsys is filled by filltsys module.


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