source: branches/plotter2/src/STApplyCal.cpp@ 3121

Last change on this file since 3121 was 2806, 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: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on data selection of targetdata in STApplyCal.


File size: 16.2 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 nchanTsys = tsystable_[0]->nchan(tsysifno);
286 ftsys = tsystable_[0]->getBaseFrequency(0);
287 interpolatorF_->setX(ftsys.data(), nchanTsys);
288 id[0] = (int)tsysifno;
289 sel.setIFs(id);
290 for (uInt i = 0; i < tsystable_.size() ; i++) {
291 tsystable_[i]->setSelection(sel);
292 nrowTsys += tsystable_[i]->nrow();
293 }
294 }
295
296 uInt nchanSp = skytable_[skylist[0]]->nchan(ifno);
297 Vector<Double> timeSky(nrowSky);
298 Matrix<Float> spoff(nchanSp, nrowSky);
299 Vector<Float> iOff(nchanSp);
300 nrowSky = 0;
301 for (uInt i = 0 ; i < skylist.nelements(); i++) {
302 STCalSkyTable *p = skytable_[skylist[i]];
303 Vector<Double> t = p->getTime();
304 Matrix<Float> sp = p->getSpectra();
305 for (uInt j = 0; j < t.nelements(); j++) {
306 timeSky[nrowSky] = t[j];
307 spoff.column(nrowSky) = sp.column(j);
308 nrowSky++;
309 }
310 }
311
312 Vector<uInt> skyIdx = timeSort(timeSky);
313
314 Double *xa = new Double[skyIdx.nelements()];
315 Float *ya = new Float[skyIdx.nelements()];
316 IPosition ipos(1, skyIdx.nelements());
317 Vector<Double> timeSkySorted(ipos, xa, TAKE_OVER);
318 Vector<Float> tmpOff(ipos, ya, TAKE_OVER);
319 for (uInt i = 0 ; i < skyIdx.nelements(); i++) {
320 timeSkySorted[i] = timeSky[skyIdx[i]];
321 }
322
323 interpolatorS_->setX(xa, skyIdx.nelements());
324
325 Vector<uInt> tsysIdx;
326 Vector<Double> timeTsys(nrowTsys);
327 Matrix<Float> tsys;
328 Vector<Double> timeTsysSorted;
329 Vector<Float> tmpTsys;
330 if (doTsys) {
331 //os_ << "doTsys" << LogIO::POST;
332 timeTsys.resize(nrowTsys);
333 tsys.resize(nchanTsys, nrowTsys);
334 nrowTsys = 0;
335 for (uInt i = 0 ; i < tsystable_.size(); i++) {
336 STCalTsysTable *p = tsystable_[i];
337 Vector<Double> t = p->getTime();
338 Matrix<Float> ts = p->getTsys();
339 for (uInt j = 0; j < t.nelements(); j++) {
340 timeTsys[nrowTsys] = t[j];
341 tsys.column(nrowTsys) = ts.column(j);
342 nrowTsys++;
343 }
344 }
345 tsysIdx = timeSort(timeTsys);
346
347 Double *xb = new Double[tsysIdx.nelements()];
348 Float *yb = new Float[tsysIdx.nelements()];
349 IPosition ipos(1, tsysIdx.nelements());
350 timeTsysSorted.takeStorage(ipos, xb, TAKE_OVER);
351 tmpTsys.takeStorage(ipos, yb, TAKE_OVER);
352 for (uInt i = 0 ; i < tsysIdx.nelements(); i++) {
353 timeTsysSorted[i] = timeTsys[tsysIdx[i]];
354 }
355 interpolatorT_->setX(xb, tsysIdx.nelements());
356 }
357
358 Table tab = work_->table();
359 ArrayColumn<Float> spCol(tab, "SPECTRA");
360 ArrayColumn<Float> tsysCol(tab, "TSYS");
361 ScalarColumn<Double> timeCol(tab, "TIME");
362 Vector<Float> on;
363 for (uInt i = 0; i < rows.nelements(); i++) {
364 //os_ << "start i = " << i << " (row = " << rows[i] << ")" << LogIO::POST;
365 uInt irow = rows[i];
366
367 // target spectral data
368 on = spCol(irow);
369 calibrator_->setSource(on);
370
371 // interpolation
372 Double t0 = timeCol(irow);
373 for (uInt ichan = 0; ichan < nchanSp; ichan++) {
374 Vector<Float> spOffSlice = spoff.row(ichan);
375 //os_ << "spOffSlice = " << spOffSlice << LogIO::POST;
376 for (uInt j = 0; j < skyIdx.nelements(); j++) {
377 tmpOff[j] = spOffSlice[skyIdx[j]];
378 }
379 interpolatorS_->setY(ya, skyIdx.nelements());
380 iOff[ichan] = interpolatorS_->interpolate(t0);
381 }
382 //os_ << "iOff=" << iOff << LogIO::POST;
383 calibrator_->setReference(iOff);
384
385 Float *Y = new Float[nchanSp];
386 Vector<Float> iTsys(IPosition(1,nchanSp), Y, TAKE_OVER);
387 if (doTsys) {
388 // Tsys correction
389 Float *yt = new Float[nchanTsys];
390 Vector<Float> iTsysT(IPosition(1,nchanTsys), yt, TAKE_OVER);
391 Float *yb = tmpTsys.data();
392 for (uInt ichan = 0; ichan < nchanTsys; ichan++) {
393 Vector<Float> tsysSlice = tsys.row(ichan);
394 for (uInt j = 0; j < tsysIdx.nelements(); j++) {
395 tmpTsys[j] = tsysSlice[tsysIdx[j]];
396 }
397 interpolatorT_->setY(yb, tsysIdx.nelements());
398 iTsysT[ichan] = interpolatorT_->interpolate(t0);
399 }
400 if (nchanSp == 1) {
401 // take average
402 iTsys[0] = mean(iTsysT);
403 }
404 else {
405 // interpolation on frequency axis
406 Vector<Double> fsp = getBaseFrequency(rows[i]);
407 interpolatorF_->setY(yt, nchanTsys);
408 for (uInt ichan = 0; ichan < nchanSp; ichan++) {
409 iTsys[ichan] = interpolatorF_->interpolate(fsp[ichan]);
410 }
411 }
412 }
413 else {
414 Vector<Float> tsysInRow = tsysCol(irow);
415 if (tsysInRow.nelements() == 1) {
416 iTsys = tsysInRow[0];
417 }
418 else {
419 for (uInt ichan = 0; ichan < nchanTsys; ++ichan)
420 iTsys[ichan] = tsysInRow[ichan];
421 }
422 }
423 //os_ << "iTsys=" << iTsys << LogIO::POST;
424 calibrator_->setScaler(iTsys);
425
426 // do calibration
427 calibrator_->calibrate();
428
429 // update table
430 //os_ << "calibrated=" << calibrator_->getCalibrated() << LogIO::POST;
431 spCol.put(irow, calibrator_->getCalibrated());
432 if (filltsys)
433 tsysCol.put(irow, iTsys);
434 }
435
436
437 // reset selection on apply tables
438 for (uInt i = 0; i < skylist.nelements(); i++)
439 skytable_[i]->unsetSelection();
440 for (uInt i = 0; i < tsystable_.size(); i++)
441 tsystable_[i]->unsetSelection();
442
443
444 // reset interpolator
445 interpolatorS_->reset();
446 interpolatorF_->reset();
447 interpolatorT_->reset();
448}
449
450Vector<uInt> STApplyCal::timeSort(Vector<Double> &t)
451{
452 Sort sort;
453 sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
454 Vector<uInt> idx;
455 sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
456 return idx;
457}
458
459uInt STApplyCal::getIFForTsys(uInt to)
460{
461 for (map<casa::uInt, Vector<uInt> >::iterator i = spwmap_.begin();
462 i != spwmap_.end(); i++) {
463 Vector<uInt> tolist = i->second;
464 os_ << "from=" << i->first << ": tolist=" << tolist << LogIO::POST;
465 for (uInt j = 0; j < tolist.nelements(); j++) {
466 if (tolist[j] == to)
467 return i->first;
468 }
469 }
470 return (uInt)-1;
471}
472
473void STApplyCal::save(const String &name)
474{
475 //assert(!work_.null());
476 assert_<AipsError>(!work_.null(),"You have to execute apply method first.");
477
478 work_->setSelection(sel_);
479 work_->makePersistent(name);
480 work_->unsetSelection();
481}
482
483Vector<Double> STApplyCal::getBaseFrequency(uInt whichrow)
484{
485 //assert(whichrow <= (uInt)work_->nrow());
486 assert_<AipsError>(whichrow <= (uInt)work_->nrow(),"row index out of range.");
487 ROTableColumn col(work_->table(), "IFNO");
488 uInt ifno = col.asuInt(whichrow);
489 col.attach(work_->table(), "FREQ_ID");
490 uInt freqid = col.asuInt(whichrow);
491 uInt nc = work_->nchan(ifno);
492 STFrequencies ftab = work_->frequencies();
493 Double rp, rf, inc;
494 ftab.getEntry(rp, rf, inc, freqid);
495 Vector<Double> r(nc);
496 indgen(r, rf-rp*inc, inc);
497 return r;
498}
499
500void STApplyCal::initInterpolator()
501{
502 os_.origin(LogOrigin("STApplyCal","initInterpolator",WHERE));
503 int order = (order_ > 0) ? order_ : 1;
504 switch (iTime_) {
505 case STCalEnum::NearestInterpolation:
506 {
507 os_ << "use NearestInterpolator in time axis" << LogIO::POST;
508 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
509 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
510 break;
511 }
512 case STCalEnum::LinearInterpolation:
513 {
514 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
515 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
516 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
517 break;
518 }
519 case STCalEnum::CubicSplineInterpolation:
520 {
521 os_ << "use CubicSplineInterpolator in time axis" << LogIO::POST;
522 interpolatorS_ = new CubicSplineInterpolator1D<Double, Float>();
523 interpolatorT_ = new CubicSplineInterpolator1D<Double, Float>();
524 break;
525 }
526 case STCalEnum::PolynomialInterpolation:
527 {
528 os_ << "use PolynomialInterpolator in time axis" << LogIO::POST;
529 if (order == 0) {
530 interpolatorS_ = new NearestInterpolator1D<Double, Float>();
531 interpolatorT_ = new NearestInterpolator1D<Double, Float>();
532 }
533 else {
534 interpolatorS_ = new PolynomialInterpolator1D<Double, Float>();
535 interpolatorT_ = new PolynomialInterpolator1D<Double, Float>();
536 interpolatorS_->setOrder(order);
537 interpolatorT_->setOrder(order);
538 }
539 break;
540 }
541 default:
542 {
543 os_ << "use BufferedLinearInterpolator in time axis" << LogIO::POST;
544 interpolatorS_ = new BufferedLinearInterpolator1D<Double, Float>();
545 interpolatorT_ = new BufferedLinearInterpolator1D<Double, Float>();
546 break;
547 }
548 }
549
550 switch (iFreq_) {
551 case STCalEnum::NearestInterpolation:
552 {
553 os_ << "use NearestInterpolator in frequency axis" << LogIO::POST;
554 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
555 break;
556 }
557 case STCalEnum::LinearInterpolation:
558 {
559 os_ << "use BufferedLinearInterpolator in frequency axis" << LogIO::POST;
560 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
561 break;
562 }
563 case STCalEnum::CubicSplineInterpolation:
564 {
565 os_ << "use CubicSplineInterpolator in frequency axis" << LogIO::POST;
566 interpolatorF_ = new CubicSplineInterpolator1D<Double, Float>();
567 break;
568 }
569 case STCalEnum::PolynomialInterpolation:
570 {
571 os_ << "use PolynomialInterpolator in frequency axis" << LogIO::POST;
572 if (order == 0) {
573 interpolatorF_ = new NearestInterpolator1D<Double, Float>();
574 }
575 else {
576 interpolatorF_ = new PolynomialInterpolator1D<Double, Float>();
577 interpolatorF_->setOrder(order);
578 }
579 break;
580 }
581 default:
582 {
583 os_ << "use LinearInterpolator in frequency axis" << LogIO::POST;
584 interpolatorF_ = new BufferedLinearInterpolator1D<Double, Float>();
585 break;
586 }
587 }
588}
589}
Note: See TracBrowser for help on using the repository browser.