source: trunk/src/CalibrationHelper.h@ 3117

Last change on this file since 3117 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

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


Check-in asap modifications from Jim regarding casacore namespace conversion.
File size: 16.7 KB
Line 
1#include <iostream>
2#include <vector>
3
4#include <casa/Arrays/Vector.h>
5#include <casa/BasicSL/String.h>
6#include <casa/Utilities/CountedPtr.h>
7
8#include "Scantable.h"
9#include "STTcal.h"
10#include "STIdxIter.h"
11#include "STSelector.h"
12
13using namespace casacore;
14using namespace asap;
15
16namespace {
17// Interpolation Helper
18class TcalData
19{
20public:
21 TcalData(CountedPtr<Scantable> s)
22 : table_(s)
23 {}
24 ~TcalData() {}
25 const String method_name() const {return "getTcalFromTime";}
26 uInt nrow() const {return table_->nrow();}
27 Vector<Float> GetEntry(int idx) const
28 {
29 String time;
30 uInt tcalid = table_->getTcalId(idx);
31 Vector<Float> return_value;
32 table_->tcal().getEntry(time, return_value, tcalid);
33 return return_value;
34 }
35private:
36 CountedPtr<Scantable> table_;
37};
38
39class TsysData
40{
41public:
42 TsysData(CountedPtr<Scantable> s)
43 : tsyscolumn_(s->table(), "TSYS")
44 {}
45 ~TsysData() {}
46 const String method_name() const {return "getTsysFromTime";}
47 uInt nrow() const {return tsyscolumn_.nrow();}
48 Vector<Float> GetEntry(int idx) const {return tsyscolumn_(idx);}
49private:
50 ROArrayColumn<Float> tsyscolumn_;
51};
52
53class SpectralData
54{
55public:
56 SpectralData(Matrix<Float> s)
57 : data_(s)
58 {}
59 ~SpectralData() {}
60 const String method_name() const {return "getSpectraFromTime";}
61 uInt nrow() const {return data_.ncolumn();}
62 Vector<Float> GetEntry(int idx) const {return data_.column(idx);}
63private:
64 Matrix<Float> data_;
65};
66
67
68vector<int> getRowIdFromTime(double reftime, const Vector<Double> &t)
69{
70 // double reft = reftime ;
71 double dtmin = 1.0e100 ;
72 double dtmax = -1.0e100 ;
73 // vector<double> dt ;
74 int just_before = -1 ;
75 int just_after = -1 ;
76 Vector<Double> dt = t - reftime ;
77 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
78 if ( dt[i] > 0.0 ) {
79 // after reftime
80 if ( dt[i] < dtmin ) {
81 just_after = i ;
82 dtmin = dt[i] ;
83 }
84 }
85 else if ( dt[i] < 0.0 ) {
86 // before reftime
87 if ( dt[i] > dtmax ) {
88 just_before = i ;
89 dtmax = dt[i] ;
90 }
91 }
92 else {
93 // just a reftime
94 just_before = i ;
95 just_after = i ;
96 dtmax = 0 ;
97 dtmin = 0 ;
98 break ;
99 }
100 }
101
102 vector<int> v(2) ;
103 v[0] = just_before ;
104 v[1] = just_after ;
105
106 return v ;
107}
108
109vector<int> getRowIdFromTime2(double reftime,
110 const Vector<Double> &t,
111 const Vector<uInt> &flagrow,
112 const Matrix<uChar> &flagtra)
113{
114 unsigned int nchan = flagtra[0].nelements();
115 vector<int> v(2*nchan);
116
117 for (unsigned int j = 0; j < nchan; ++j) {
118 // double reft = reftime ;
119 double dtmin = 1.0e100 ;
120 double dtmax = -1.0e100 ;
121 // vector<double> dt ;
122 int just_before = -1 ;
123 int just_after = -1 ;
124 Vector<Double> dt = t - reftime ;
125 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
126 if ( flagrow[i] > 0 ) continue;
127 if ( flagtra.column(i)[j] == 1 << 7) continue;
128
129 if ( dt[i] > 0.0 ) {
130 // after reftime
131 if ( dt[i] < dtmin ) {
132 just_after = i ;
133 dtmin = dt[i] ;
134 }
135 }
136 else if ( dt[i] < 0.0 ) {
137 // before reftime
138 if ( dt[i] > dtmax ) {
139 just_before = i ;
140 dtmax = dt[i] ;
141 }
142 }
143 else {
144 // just a reftime
145 just_before = i ;
146 just_after = i ;
147 dtmax = 0 ;
148 dtmin = 0 ;
149 break ;
150 }
151 }
152
153 v[j*2] = just_before ;
154 v[j*2+1] = just_after ;
155 }
156
157 return v ;
158}
159
160template<class T>
161class SimpleInterpolationHelper
162{
163 public:
164 static Vector<Float> GetFromTime(double reftime,
165 const Vector<Double> &timeVec,
166 const vector<int> &idx,
167 const T &data,
168 const string mode)
169 {
170 Vector<Float> return_value(idx.size()/2);
171 LogIO os_;
172 LogIO os( LogOrigin( "STMath", data.method_name(), WHERE ) ) ;
173 if ( data.nrow() == 0 ) {
174 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
175 }
176 else if ( data.nrow() == 1 ) {
177 return_value = data.GetEntry(0);
178 }
179 else {
180 for (unsigned int i = 0; i < idx.size()/2; ++i) {
181 unsigned int idx0 = 2*i;
182 unsigned int idx1 = 2*i + 1;
183 //no off data available. calibration impossible.
184 if ( ( idx[idx0] == -1 ) && ( idx[idx1] == -1 ) ) continue;
185
186 if ( mode == "before" ) {
187 int id = -1 ;
188 if ( idx[idx0] != -1 ) {
189 id = idx[idx0] ;
190 }
191 else if ( idx[idx1] != -1 ) {
192 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
193 id = idx[idx1] ;
194 }
195
196 return_value[i] = data.GetEntry(id)[i];
197 }
198 else if ( mode == "after" ) {
199 int id = -1 ;
200 if ( idx[idx1] != -1 ) {
201 id = idx[idx1] ;
202 }
203 else if ( idx[idx0] != -1 ) {
204 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
205 id = idx[idx1] ;
206 }
207
208 return_value[i] = data.GetEntry(id)[i];
209 }
210 else if ( mode == "nearest" ) {
211 int id = -1 ;
212 if ( idx[idx0] == -1 ) {
213 id = idx[idx1] ;
214 }
215 else if ( idx[idx1] == -1 ) {
216 id = idx[idx0] ;
217 }
218 else if ( idx[idx0] == idx[idx1] ) {
219 id = idx[idx0] ;
220 }
221 else {
222 double t0 = timeVec[idx[idx0]] ;
223 double t1 = timeVec[idx[idx1]] ;
224 if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
225 id = idx[idx1] ;
226 }
227 else {
228 id = idx[idx0] ;
229 }
230 }
231 return_value[i] = data.GetEntry(id)[i];
232 }
233 else if ( mode == "linear" ) {
234 if ( idx[idx0] == -1 ) {
235 // use after
236 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
237 int id = idx[idx1] ;
238 return_value[i] = data.GetEntry(id)[i];
239 }
240 else if ( idx[idx1] == -1 ) {
241 // use before
242 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
243 int id = idx[idx0] ;
244 return_value[i] = data.GetEntry(id)[i];
245 }
246 else if ( idx[idx0] == idx[idx1] ) {
247 // use before
248 //os << "No need to interporate." << LogIO::POST ;
249 int id = idx[idx0] ;
250 return_value[i] = data.GetEntry(id)[i];
251 }
252 else {
253 // do interpolation
254 double t0 = timeVec[idx[idx0]] ;
255 double t1 = timeVec[idx[idx1]] ;
256 Vector<Float> value0 = data.GetEntry(idx[idx0]);
257 Vector<Float> value1 = data.GetEntry(idx[idx1]);
258 double tfactor = (reftime - t0) / (t1 - t0) ;
259 return_value[i] = ( value1[i] - value0[i] ) * tfactor + value0[i] ;
260 }
261 }
262 else {
263 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
264 }
265 }
266 }
267 return return_value ;
268 }
269};
270
271// Calibration Helper
272class CalibrationHelper
273{
274public:
275 static void CalibrateALMA( CountedPtr<Scantable>& out,
276 const CountedPtr<Scantable>& on,
277 const CountedPtr<Scantable>& off,
278 const Vector<uInt>& rows )
279 {
280 // 2012/05/22 TN
281 // Assume that out has empty SPECTRA column
282
283 // if rows is empty, just return
284 if ( rows.nelements() == 0 )
285 return ;
286
287 ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
288 ROArrayColumn<Float> in_tsys_column(on->table(), "TSYS");
289 ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
290 ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
291 ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
292
293 Vector<Double> timeVec = GetScalarColumn<Double>(off->table(), "TIME");
294 Vector<Double> refTimeVec = GetScalarColumn<Double>(on->table(), "TIME");
295 Vector<uInt> flagrowVec = GetScalarColumn<uInt>(off->table(), "FLAGROW");
296 Vector<uInt> refFlagrowVec = GetScalarColumn<uInt>(on->table(), "FLAGROW");
297 Matrix<uChar> flagtraMtx = GetArrayColumn<uChar>(off->table(), "FLAGTRA");
298 SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
299 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
300 vector<int> ids( 2 * spsize ) ;
301
302 for ( unsigned int irow = 0 ; irow < rows.nelements() ; irow++ ) {
303 uInt row = rows[irow];
304 double reftime = refTimeVec[row];
305 ids = getRowIdFromTime2( reftime, timeVec, flagrowVec, flagtraMtx ) ;
306 Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
307 Vector<Float> spec = in_spectra_column(row);
308 Vector<Float> tsys = in_tsys_column(row);
309 Vector<uChar> flag = in_flagtra_column(row);
310
311 // ALMA Calibration
312 //
313 // Ta* = Tsys * ( ON - OFF ) / OFF
314 //
315 // 2010/01/07 Takeshi Nakazato
316 unsigned int tsyssize = tsys.nelements() ;
317 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
318 //if there is no off data available for a channel, just flag the channel.(2014/7/18 WK)
319 if ((ids[2*j] == -1)&&(ids[2*j+1] == -1)) {
320 flag[j] = 1 << 7;
321 continue;
322 }
323
324 if (refFlagrowVec[row] == 0) {
325 if ( spoff[j] == 0.0 ) {
326 spec[j] = 0.0 ;
327 flag[j] = (uChar)True;
328 }
329 else {
330 spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
331 }
332 spec[j] *= (tsyssize == spsize) ? tsys[j] : tsys[0];
333 }
334 }
335 out_spectra_column.put(row, spec);
336 out_flagtra_column.put(row, flag);
337 }
338 }
339 static void CalibrateChopperWheel( CountedPtr<Scantable> &out,
340 const CountedPtr<Scantable>& on,
341 const CountedPtr<Scantable>& off,
342 const CountedPtr<Scantable>& sky,
343 const CountedPtr<Scantable>& hot,
344 const CountedPtr<Scantable>& /*cold*/,
345 const Vector<uInt> &rows )
346 {
347 // 2012/05/22 TN
348 // Assume that out has empty SPECTRA column
349
350 // if rows is empty, just return
351 if ( rows.nelements() == 0 )
352 return ;
353
354 string antenna_name = out->getAntennaName();
355 ROArrayColumn<Float> in_spectra_column(on->table(), "SPECTRA");
356 ROArrayColumn<uChar> in_flagtra_column(on->table(), "FLAGTRA");
357 ArrayColumn<Float> out_spectra_column(out->table(), "SPECTRA");
358 ArrayColumn<uChar> out_flagtra_column(out->table(), "FLAGTRA");
359 ArrayColumn<Float> out_tsys_column(out->table(), "TSYS");
360
361 Vector<Double> timeOff = GetScalarColumn<Double>(off->table(), "TIME");
362 Vector<Double> timeSky = GetScalarColumn<Double>(sky->table(), "TIME");
363 Vector<Double> timeHot = GetScalarColumn<Double>(hot->table(), "TIME");
364 Vector<Double> timeOn = GetScalarColumn<Double>(on->table(), "TIME");
365 Vector<uInt> flagrowOff = GetScalarColumn<uInt>(off->table(), "FLAGROW");
366 Vector<uInt> flagrowSky = GetScalarColumn<uInt>(sky->table(), "FLAGROW");
367 Vector<uInt> flagrowHot = GetScalarColumn<uInt>(hot->table(), "FLAGROW");
368 Vector<uInt> flagrowOn = GetScalarColumn<uInt>(on->table(), "FLAGROW");
369 Matrix<uChar> flagtraOff = GetArrayColumn<uChar>(off->table(), "FLAGTRA");
370 Matrix<uChar> flagtraSky = GetArrayColumn<uChar>(sky->table(), "FLAGTRA");
371 Matrix<uChar> flagtraHot = GetArrayColumn<uChar>(hot->table(), "FLAGTRA");
372 SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
373 SpectralData skyspectra(Matrix<Float>(GetArrayColumn<Float>(sky->table(), "SPECTRA")));
374 SpectralData hotspectra(Matrix<Float>(GetArrayColumn<Float>(hot->table(), "SPECTRA")));
375 TcalData tcaldata(sky);
376 TsysData tsysdata(sky);
377 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
378 vector<int> idsOff( 2 * spsize ) ;
379 vector<int> idsSky( 2 * spsize ) ;
380 vector<int> idsHot( 2 * spsize ) ;
381 for ( unsigned int irow = 0 ; irow < rows.nelements() ; irow++ ) {
382 uInt row = rows[irow];
383 double reftime = timeOn[row];
384 idsOff = getRowIdFromTime2( reftime, timeOff, flagrowOff, flagtraOff ) ;
385 Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, idsOff, offspectra, "linear");
386 idsSky = getRowIdFromTime2( reftime, timeSky, flagrowSky, flagtraSky ) ;
387 Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, idsSky, skyspectra, "linear");
388 Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, idsSky, tcaldata, "linear");
389 Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, idsSky, tsysdata, "linear");
390 idsHot = getRowIdFromTime2( reftime, timeHot, flagrowHot, flagtraHot ) ;
391 Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, idsHot, hotspectra, "linear");
392 Vector<Float> spec = in_spectra_column(row);
393 Vector<uChar> flag = in_flagtra_column(row);
394 if ( antenna_name.find( "APEX" ) != String::npos ) {
395 // using gain array
396 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
397 //if at least one of off/sky/hot data unavailable, just flag the channel.
398 if (((idsOff[2*j] == -1)&&(idsOff[2*j+1] == -1))||
399 ((idsSky[2*j] == -1)&&(idsSky[2*j+1] == -1))||
400 ((idsHot[2*j] == -1)&&(idsHot[2*j+1] == -1))) {
401 flag[j] = (uChar)True;
402 continue;
403 }
404 if (flagrowOn[row] == 0) {
405 if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
406 spec[j] = 0.0 ;
407 flag[j] = (uChar)True;
408 }
409 else {
410 spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
411 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
412 }
413 }
414 }
415 }
416 else {
417 // Chopper-Wheel calibration (Ulich & Haas 1976)
418 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
419 //if at least one of off/sky/hot data unavailable, just flag the channel.
420 if (((idsOff[2*j] == -1)&&(idsOff[2*j+1] == -1))||
421 ((idsSky[2*j] == -1)&&(idsSky[2*j+1] == -1))||
422 ((idsHot[2*j] == -1)&&(idsHot[2*j+1] == -1))) {
423 flag[j] = (uChar)True;
424 continue;
425 }
426 if (flagrowOn[row] == 0) {
427 if ( (sphot[j]-spsky[j]) == 0.0 ) {
428 spec[j] = 0.0 ;
429 flag[j] = (uChar)True;
430 }
431 else {
432 spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
433 }
434 }
435 }
436 }
437 out_spectra_column.put(row, spec);
438 out_flagtra_column.put(row, flag);
439 out_tsys_column.put(row, tsys);
440 }
441 }
442 static void GetSelector(STSelector &sel, const vector<string> &names, const Record &values)
443 {
444 stringstream ss ;
445 ss << "SELECT FROM $1 WHERE ";
446 string separator = "";
447 for (vector<string>::const_iterator i = names.begin(); i != names.end(); ++i) {
448 ss << separator << *i << "==";
449 switch (values.dataType(*i)) {
450 case TpUInt:
451 ss << values.asuInt(*i);
452 break;
453 case TpInt:
454 ss << values.asInt(*i);
455 break;
456 case TpFloat:
457 ss << values.asFloat(*i);
458 break;
459 case TpDouble:
460 ss << values.asDouble(*i);
461 break;
462 case TpComplex:
463 ss << values.asComplex(*i);
464 break;
465 case TpString:
466 ss << values.asString(*i);
467 break;
468 default:
469 break;
470 }
471 separator = "&&";
472 }
473 sel.setTaQL(ss.str());
474 }
475private:
476 template<class T>
477 static Vector<T> GetScalarColumn(const Table &table, const String &name)
478 {
479 ROScalarColumn<T> column(table, name);
480 return column.getColumn();
481 }
482 template<class T>
483 static Array<T> GetArrayColumn(const Table &table, const String &name)
484 {
485 ROArrayColumn<T> column(table, name);
486 return column.getColumn();
487 }
488};
489
490class AlmaCalibrator
491{
492public:
493 AlmaCalibrator(CountedPtr<Scantable> &out,
494 const CountedPtr<Scantable> &on,
495 const CountedPtr<Scantable> &off)
496 : target_(out),
497 selector_(),
498 on_(on),
499 off_(off)
500 {}
501 ~AlmaCalibrator() {}
502 CountedPtr<Scantable> target() const {return target_;}
503 void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
504 CalibrationHelper::GetSelector(selector_, cols, values);
505 off_->setSelection(selector_);
506 CalibrationHelper::CalibrateALMA(target_, on_, off_, rows);
507 off_->unsetSelection();
508 }
509private:
510 CountedPtr<Scantable> target_;
511 STSelector selector_;
512 const CountedPtr<Scantable> on_;
513 const CountedPtr<Scantable> off_;
514};
515
516class ChopperWheelCalibrator
517{
518public:
519 ChopperWheelCalibrator(CountedPtr<Scantable> &out,
520 const CountedPtr<Scantable> &on,
521 const CountedPtr<Scantable> &sky,
522 const CountedPtr<Scantable> &off,
523 const CountedPtr<Scantable> &hot,
524 const CountedPtr<Scantable> &cold)
525 : target_(out),
526 selector_(),
527 on_(on),
528 off_(off),
529 sky_(sky),
530 hot_(hot),
531 cold_(cold)
532 {}
533 ~ChopperWheelCalibrator() {}
534 CountedPtr<Scantable> target() const {return target_;}
535 void Process(const vector<string> &cols, const Record &values, const Vector<uInt> &rows) {
536 CalibrationHelper::GetSelector(selector_, cols, values);
537 off_->setSelection(selector_);
538 sky_->setSelection(selector_);
539 hot_->setSelection(selector_);
540 CalibrationHelper::CalibrateChopperWheel(target_, on_, off_, sky_, hot_, cold_, rows);
541 off_->unsetSelection();
542 sky_->unsetSelection();
543 hot_->unsetSelection();
544 }
545private:
546 CountedPtr<Scantable> target_;
547 STSelector selector_;
548 const CountedPtr<Scantable> on_;
549 const CountedPtr<Scantable> off_;
550 const CountedPtr<Scantable> sky_;
551 const CountedPtr<Scantable> hot_;
552 const CountedPtr<Scantable> cold_;
553};
554
555} // anonymous namespace
Note: See TracBrowser for help on using the repository browser.