source: trunk/src/STApplyCal.cpp@ 2736

Last change on this file since 2736 was 2735, 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...

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