source: trunk/src/STMath.cpp@ 2941

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

New Development: No

JIRA Issue: Yes CAS-5875

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix for freqtol.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 146.3 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <sstream>
14
15#include <casa/iomanip.h>
16#include <casa/Arrays/MaskArrLogi.h>
17#include <casa/Arrays/MaskArrMath.h>
18#include <casa/Arrays/ArrayLogical.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/Slice.h>
21#include <casa/Arrays/Slicer.h>
22#include <casa/BasicSL/String.h>
23#include <casa/Containers/Block.h>
24#include <casa/Containers/RecordField.h>
25#include <casa/Exceptions/Error.h>
26#include <casa/Logging/LogIO.h>
27#include <casa/Quanta/Quantum.h>
28
29#include <coordinates/Coordinates/CoordinateSystem.h>
30#include <coordinates/Coordinates/CoordinateUtil.h>
31#include <coordinates/Coordinates/FrequencyAligner.h>
32#include <coordinates/Coordinates/SpectralCoordinate.h>
33
34#include <lattices/Lattices/LatticeUtilities.h>
35
36#include <scimath/Functionals/Polynomial.h>
37#include <scimath/Mathematics/Convolver.h>
38#include <scimath/Mathematics/VectorKernel.h>
39
40#include <tables/Tables/ExprNode.h>
41#include <tables/Tables/ReadAsciiTable.h>
42#include <tables/Tables/TableCopy.h>
43#include <tables/Tables/TableIter.h>
44#include <tables/Tables/TableParse.h>
45#include <tables/Tables/TableRecord.h>
46#include <tables/Tables/TableRow.h>
47#include <tables/Tables/TableVector.h>
48#include <tables/Tables/TabVecMath.h>
49
50#include <atnf/PKSIO/SrcType.h>
51
52#include "RowAccumulator.h"
53#include "STAttr.h"
54#include "STMath.h"
55#include "STSelector.h"
56#include "Accelerator.h"
57#include "STIdxIter.h"
58
59#include "CalibrationHelper.h"
60
61using namespace casa;
62using namespace asap;
63
64// 2012/02/17 TN
65// Since STGrid is implemented, average doesn't consider direction
66// when accumulating
67// tolerance for direction comparison (rad)
68// #define TOL_OTF 1.0e-15
69// #define TOL_POINT 2.9088821e-4 // 1 arcmin
70
71STMath::STMath(bool insitu) :
72 insitu_(insitu)
73{
74}
75
76
77STMath::~STMath()
78{
79}
80
81CountedPtr<Scantable>
82STMath::average( const std::vector<CountedPtr<Scantable> >& in,
83 const std::vector<bool>& mask,
84 const std::string& weight,
85 const std::string& avmode)
86{
87// double t0, t1 ;
88// t0 = mathutil::gettimeofday_sec() ;
89
90 LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
91 if ( avmode == "SCAN" && in.size() != 1 )
92 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
93 "Use merge first."));
94 WeightType wtype = stringToWeight(weight);
95
96 // 2012/02/17 TN
97 // Since STGrid is implemented, average doesn't consider direction
98 // when accumulating
99 // check if OTF observation
100// String obstype = in[0]->getHeader().obstype ;
101// Double tol = 0.0 ;
102// if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
103// tol = TOL_OTF ;
104// }
105// else {
106// tol = TOL_POINT ;
107// }
108
109 // output
110 // clone as this is non insitu
111 bool insitu = insitu_;
112 setInsitu(false);
113 CountedPtr< Scantable > out = getScantable(in[0], true);
114 setInsitu(insitu);
115 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
116 ++stit;
117 while ( stit != in.end() ) {
118 out->appendToHistoryTable((*stit)->history());
119 ++stit;
120 }
121
122 Table& tout = out->table();
123
124 /// @todo check if all scantables are conformant
125
126 ArrayColumn<Float> specColOut(tout,"SPECTRA");
127 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
128 ArrayColumn<Float> tsysColOut(tout,"TSYS");
129 ScalarColumn<Double> mjdColOut(tout,"TIME");
130 ScalarColumn<Double> intColOut(tout,"INTERVAL");
131 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
132 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
133 ScalarColumn<uInt> flagRowColOut(tout,"FLAGROW");
134
135 // set up the output table rows. These are based on the structure of the
136 // FIRST scantable in the vector
137 const Table& baset = in[0]->table();
138
139 RowAccumulator acc(wtype);
140 Vector<Bool> cmask(mask);
141 acc.setUserMask(cmask);
142// ROTableRow row(tout);
143 ROArrayColumn<Float> specCol, tsysCol;
144 ROArrayColumn<uChar> flagCol;
145 ROScalarColumn<Double> mjdCol, intCol;
146 ROScalarColumn<Int> scanIDCol;
147 ROScalarColumn<uInt> flagRowCol;
148
149 //Vector<uInt> rowstodelete;
150 Block<uInt> rowstodelB( in[0]->nrow() ) ;
151 uInt nrowdel = 0 ;
152
153// Block<String> cols(3);
154 vector<string> cols(3) ;
155 cols[0] = String("BEAMNO");
156 cols[1] = String("IFNO");
157 cols[2] = String("POLNO");
158 if ( avmode == "SOURCE" ) {
159 cols.resize(4);
160 cols[3] = String("SRCNAME");
161 }
162 if ( avmode == "SCAN" && in.size() == 1) {
163 //cols.resize(4);
164 //cols[3] = String("SCANNO");
165 cols.resize(5);
166 cols[3] = String("SRCNAME");
167 cols[4] = String("SCANNO");
168 }
169 uInt outrowCount = 0;
170 // use STIdxIter2 instead of TableIterator
171 STIdxIter2 iter( in[0], cols ) ;
172// double t2 = 0 ;
173// double t3 = 0 ;
174// double t4 = 0 ;
175// double t5 = 0 ;
176// TableIterator iter(baset, cols);
177// int count = 0 ;
178 while (!iter.pastEnd()) {
179 Vector<uInt> rows = iter.getRows( SHARE ) ;
180 if ( rows.nelements() == 0 ) {
181 iter.next() ;
182 continue ;
183 }
184 Record current = iter.currentValue() ;
185 //Table subt = iter.table();
186 // copy the first row of this selection into the new table
187 tout.addRow();
188// t4 = mathutil::gettimeofday_sec() ;
189 // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
190 // overwritten in the following process
191 copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
192// t5 += mathutil::gettimeofday_sec() - t4 ;
193 // re-index to 0
194 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
195 scanColOut.put(outrowCount, uInt(0));
196 }
197
198 // 2012/02/17 TN
199 // Since STGrid is implemented, average doesn't consider direction
200 // when accumulating
201// MDirection::ScalarColumn dircol ;
202// dircol.attach( subt, "DIRECTION" ) ;
203// Int length = subt.nrow() ;
204// vector< Vector<Double> > dirs ;
205// vector<int> indexes ;
206// for ( Int i = 0 ; i < length ; i++ ) {
207// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
208// //os << << count++ << ": " ;
209// //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
210// bool adddir = true ;
211// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
212// //if ( allTrue( t == dirs[j] ) ) {
213// Double dx = t[0] - dirs[j][0] ;
214// Double dy = t[1] - dirs[j][1] ;
215// Double dd = sqrt( dx * dx + dy * dy ) ;
216// //if ( allNearAbs( t, dirs[j], tol ) ) {
217// if ( dd <= tol ) {
218// adddir = false ;
219// break ;
220// }
221// }
222// if ( adddir ) {
223// dirs.push_back( t ) ;
224// indexes.push_back( i ) ;
225// }
226// }
227// uInt rowNum = dirs.size() ;
228// tout.addRow( rowNum ) ;
229// for ( uInt i = 0 ; i < rowNum ; i++ ) {
230// TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
231// // re-index to 0
232// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
233// scanColOut.put(outrowCount+i, uInt(0));
234// }
235// }
236// outrowCount += rowNum ;
237
238 // merge loop
239 uInt i = outrowCount ;
240 // in[0] is already selected by iterator
241 specCol.attach(baset,"SPECTRA");
242 flagCol.attach(baset,"FLAGTRA");
243 tsysCol.attach(baset,"TSYS");
244 intCol.attach(baset,"INTERVAL");
245 mjdCol.attach(baset,"TIME");
246 flagRowCol.attach(baset,"FLAGROW");
247 Vector<Float> spec,tsys;
248 Vector<uChar> flag;
249 Double inter,time;
250 uInt flagRow;
251
252 for (uInt l = 0; l < rows.nelements(); ++l ) {
253 uInt k = rows[l] ;
254 flagCol.get(k, flag);
255 Vector<Bool> bflag(flag.shape());
256 flagRowCol.get(k, flagRow);
257 if (flagRow > 0)
258 bflag = true;
259 else
260 convertArray(bflag, flag);
261 /*
262 if ( allEQ(bflag, True) ) {
263 continue;//don't accumulate
264 }
265 */
266 specCol.get(k, spec);
267 tsysCol.get(k, tsys);
268 intCol.get(k, inter);
269 mjdCol.get(k, time);
270 // spectrum has to be added last to enable weighting by the other values
271// t2 = mathutil::gettimeofday_sec() ;
272 acc.add(spec, !bflag, tsys, inter, time);
273// t3 += mathutil::gettimeofday_sec() - t2 ;
274
275 }
276
277
278 // in[0] is already selected by TableIterator so that index is
279 // started from 1
280 for ( int j=1; j < int(in.size()); ++j ) {
281 const Table& tin = in[j]->table();
282 //const TableRecord& rec = row.get(i);
283 ROScalarColumn<Double> tmp(tin, "TIME");
284 Double td;tmp.get(0,td);
285
286#if 1
287 static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
288 //uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
289 uInt const values1[] = { current.asuInt("IFNO"), current.asuInt("BEAMNO"), current.asuInt("POLNO") };
290 SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
291 CustomTableExprNodeRep myNodeRep(tin, myPred);
292 myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
293 CustomTableExprNode myExpr(myNodeRep);
294 Table basesubt = tin(myExpr);
295#else
296// Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
297// && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
298// && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
299 Table basesubt = tin( tin.col("BEAMNO") == current.asuInt("BEAMNO")
300 && tin.col("IFNO") == current.asuInt("IFNO")
301 && tin.col("POLNO") == current.asuInt("POLNO") );
302#endif
303 Table subt;
304 if ( avmode == "SOURCE") {
305// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
306 subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME") );
307
308 } else if (avmode == "SCAN") {
309// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
310// && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
311 subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME")
312 && basesubt.col("SCANNO") == current.asuInt("SCANNO") );
313 } else {
314 subt = basesubt;
315 }
316
317 // 2012/02/17 TN
318 // Since STGrid is implemented, average doesn't consider direction
319 // when accumulating
320// vector<uInt> removeRows ;
321// uInt nrsubt = subt.nrow() ;
322// for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
323// //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
324// Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
325// Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
326// double dx = x0[0] - x1[0];
327// double dy = x0[1] - x1[1];
328// Double dd = sqrt( dx * dx + dy * dy ) ;
329// //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
330// if ( dd > tol ) {
331// removeRows.push_back( irow ) ;
332// }
333// }
334// if ( removeRows.size() != 0 ) {
335// subt.removeRow( removeRows ) ;
336// }
337
338// if ( nrsubt == removeRows.size() )
339// throw(AipsError("Averaging data is empty.")) ;
340
341 specCol.attach(subt,"SPECTRA");
342 flagCol.attach(subt,"FLAGTRA");
343 tsysCol.attach(subt,"TSYS");
344 intCol.attach(subt,"INTERVAL");
345 mjdCol.attach(subt,"TIME");
346 flagRowCol.attach(subt,"FLAGROW");
347 for (uInt k = 0; k < subt.nrow(); ++k ) {
348 flagCol.get(k, flag);
349 Vector<Bool> bflag(flag.shape());
350 flagRowCol.get(k, flagRow);
351 if (flagRow > 0)
352 bflag = true;
353 else
354 convertArray(bflag, flag);
355 /*
356 if ( allEQ(bflag, True) ) {
357 continue;//don't accumulate
358 }
359 */
360 specCol.get(k, spec);
361 //tsysCol.get(k, tsys);
362 tsys.assign( tsysCol(k) );
363 intCol.get(k, inter);
364 mjdCol.get(k, time);
365 // spectrum has to be added last to enable weighting by the other values
366// t2 = mathutil::gettimeofday_sec() ;
367 acc.add(spec, !bflag, tsys, inter, time);
368// t3 += mathutil::gettimeofday_sec() - t2 ;
369 }
370
371 }
372 const Vector<Bool>& msk = acc.getMask();
373 if ( allEQ(msk, False) ) {
374 rowstodelB[nrowdel] = i ;
375 nrowdel++ ;
376 continue;
377 }
378 //write out
379 if (acc.state()) {
380 // If there exists a channel at which all the input spectra are masked,
381 // spec has 'nan' values for that channel and it may affect the following
382 // processes. To avoid this, replacing 'nan' values in spec with
383 // weighted-mean of all spectra in the following line.
384 // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
385 acc.replaceNaN();
386
387 Vector<uChar> flg(msk.shape());
388 convertArray(flg, !msk);
389 for (uInt k = 0; k < flg.nelements(); ++k) {
390 uChar userFlag = 1 << 7;
391 if (msk[k]==True) userFlag = 0 << 7;
392 flg(k) = userFlag;
393 }
394
395 flagColOut.put(i, flg);
396 specColOut.put(i, acc.getSpectrum());
397 tsysColOut.put(i, acc.getTsys());
398 intColOut.put(i, acc.getInterval());
399 mjdColOut.put(i, acc.getTime());
400 // we should only have one cycle now -> reset it to be 0
401 // frequency switched data has different CYCLENO for different IFNO
402 // which requires resetting this value
403 cycColOut.put(i, uInt(0));
404 // completely flagged rows are removed anyway
405 flagRowColOut.put(i, uInt(0));
406 } else {
407 os << "For output row="<<i<<", all input rows of data are flagged. no averaging" << LogIO::POST;
408 }
409 acc.reset();
410
411 // merge with while loop for preparing out table
412 ++outrowCount;
413// ++iter ;
414 iter.next() ;
415 }
416
417 if ( nrowdel > 0 ) {
418 Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ;
419 os << rowstodelete << LogIO::POST ;
420 tout.removeRow(rowstodelete);
421 if (tout.nrow() == 0) {
422 throw(AipsError("Can't average fully flagged data."));
423 }
424 }
425
426// t1 = mathutil::gettimeofday_sec() ;
427// cout << "elapsed time for average(): " << t1-t0 << " sec" << endl ;
428// cout << " elapsed time for acc.add(): " << t3 << " sec" << endl ;
429// cout << " elapsed time for copyRows(): " << t5 << " sec" << endl ;
430
431 return out;
432}
433
434CountedPtr< Scantable >
435STMath::averageChannel( const CountedPtr < Scantable > & in,
436 const std::string & mode,
437 const std::string& avmode )
438{
439 (void) mode; // currently unused
440 // 2012/02/17 TN
441 // Since STGrid is implemented, average doesn't consider direction
442 // when accumulating
443 // check if OTF observation
444// String obstype = in->getHeader().obstype ;
445// Double tol = 0.0 ;
446// if ( obstype.find( "OTF" ) != String::npos ) {
447// tol = TOL_OTF ;
448// }
449// else {
450// tol = TOL_POINT ;
451// }
452
453 // clone as this is non insitu
454 bool insitu = insitu_;
455 setInsitu(false);
456 CountedPtr< Scantable > out = getScantable(in, true);
457 setInsitu(insitu);
458 Table& tout = out->table();
459 ArrayColumn<Float> specColOut(tout,"SPECTRA");
460 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
461 ArrayColumn<Float> tsysColOut(tout,"TSYS");
462 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
463 ScalarColumn<Double> intColOut(tout, "INTERVAL");
464 Table tmp = in->table().sort("BEAMNO");
465 Block<String> cols(3);
466 cols[0] = String("BEAMNO");
467 cols[1] = String("IFNO");
468 cols[2] = String("POLNO");
469 if ( avmode == "SCAN") {
470 cols.resize(4);
471 cols[3] = String("SCANNO");
472 }
473 uInt outrowCount = 0;
474 uChar userflag = 1 << 7;
475 TableIterator iter(tmp, cols);
476 while (!iter.pastEnd()) {
477 Table subt = iter.table();
478 ROArrayColumn<Float> specCol, tsysCol;
479 ROArrayColumn<uChar> flagCol;
480 ROScalarColumn<Double> intCol(subt, "INTERVAL");
481 specCol.attach(subt,"SPECTRA");
482 flagCol.attach(subt,"FLAGTRA");
483 tsysCol.attach(subt,"TSYS");
484
485 tout.addRow();
486 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
487 if ( avmode != "SCAN") {
488 scanColOut.put(outrowCount, uInt(0));
489 }
490 Vector<Float> tmp;
491 specCol.get(0, tmp);
492 uInt nchan = tmp.nelements();
493 // have to do channel by channel here as MaskedArrMath
494 // doesn't have partialMedians
495 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
496 Vector<Float> outspec(nchan);
497 Vector<uChar> outflag(nchan,0);
498 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
499 for (uInt i=0; i<nchan; ++i) {
500 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
501 MaskedArray<Float> ma = maskedArray(specs,flags);
502 outspec[i] = median(ma);
503 if ( allEQ(ma.getMask(), False) )
504 outflag[i] = userflag;// flag data
505 }
506 outtsys[0] = median(tsysCol.getColumn());
507 specColOut.put(outrowCount, outspec);
508 flagColOut.put(outrowCount, outflag);
509 tsysColOut.put(outrowCount, outtsys);
510 Double intsum = sum(intCol.getColumn());
511 intColOut.put(outrowCount, intsum);
512 ++outrowCount;
513 ++iter;
514
515 // 2012/02/17 TN
516 // Since STGrid is implemented, average doesn't consider direction
517 // when accumulating
518// MDirection::ScalarColumn dircol ;
519// dircol.attach( subt, "DIRECTION" ) ;
520// Int length = subt.nrow() ;
521// vector< Vector<Double> > dirs ;
522// vector<int> indexes ;
523// // Handle MX mode averaging
524// if (in->nbeam() > 1 ) {
525// length = 1;
526// }
527// for ( Int i = 0 ; i < length ; i++ ) {
528// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
529// bool adddir = true ;
530// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
531// //if ( allTrue( t == dirs[j] ) ) {
532// Double dx = t[0] - dirs[j][0] ;
533// Double dy = t[1] - dirs[j][1] ;
534// Double dd = sqrt( dx * dx + dy * dy ) ;
535// //if ( allNearAbs( t, dirs[j], tol ) ) {
536// if ( dd <= tol ) {
537// adddir = false ;
538// break ;
539// }
540// }
541// if ( adddir ) {
542// dirs.push_back( t ) ;
543// indexes.push_back( i ) ;
544// }
545// }
546// uInt rowNum = dirs.size() ;
547// tout.addRow( rowNum );
548// for ( uInt i = 0 ; i < rowNum ; i++ ) {
549// TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
550// // Handle MX mode averaging
551// if ( avmode != "SCAN") {
552// scanColOut.put(outrowCount+i, uInt(0));
553// }
554// }
555// MDirection::ScalarColumn dircolOut ;
556// dircolOut.attach( tout, "DIRECTION" ) ;
557// for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
558// Vector<Double> t = \
559// dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
560// Vector<Float> tmp;
561// specCol.get(0, tmp);
562// uInt nchan = tmp.nelements();
563// // have to do channel by channel here as MaskedArrMath
564// // doesn't have partialMedians
565// Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
566// // mask spectra for different DIRECTION
567// for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
568// Vector<Double> direction = \
569// dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
570// //if ( t[0] != direction[0] || t[1] != direction[1] ) {
571// Double dx = t[0] - direction[0];
572// Double dy = t[1] - direction[1];
573// Double dd = sqrt(dx*dx + dy*dy);
574// //if ( !allNearAbs( t, direction, tol ) ) {
575// if ( dd > tol && in->nbeam() < 2 ) {
576// flags[jrow] = userflag ;
577// }
578// }
579// Vector<Float> outspec(nchan);
580// Vector<uChar> outflag(nchan,0);
581// Vector<Float> outtsys(1);/// @fixme when tsys is channel based
582// for (uInt i=0; i<nchan; ++i) {
583// Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
584// MaskedArray<Float> ma = maskedArray(specs,flags);
585// outspec[i] = median(ma);
586// if ( allEQ(ma.getMask(), False) )
587// outflag[i] = userflag;// flag data
588// }
589// outtsys[0] = median(tsysCol.getColumn());
590// specColOut.put(outrowCount+irow, outspec);
591// flagColOut.put(outrowCount+irow, outflag);
592// tsysColOut.put(outrowCount+irow, outtsys);
593// Vector<Double> integ = intCol.getColumn() ;
594// MaskedArray<Double> mi = maskedArray( integ, flags ) ;
595// Double intsum = sum(mi);
596// intColOut.put(outrowCount+irow, intsum);
597// }
598// outrowCount += rowNum ;
599// ++iter;
600 }
601 return out;
602}
603
604CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
605 bool droprows)
606{
607 if (insitu_) {
608 return in;
609 }
610 else {
611 // clone
612 return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
613 }
614}
615
616CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
617 float val,
618 const std::string& mode,
619 bool tsys )
620{
621 CountedPtr< Scantable > out = getScantable(in, false);
622 Table& tab = out->table();
623 ArrayColumn<Float> specCol(tab,"SPECTRA");
624 ArrayColumn<Float> tsysCol(tab,"TSYS");
625 if (mode=="DIV") val = 1.0/val ;
626 else if (mode=="SUB") val *= -1.0 ;
627 for (uInt i=0; i<tab.nrow(); ++i) {
628 Vector<Float> spec;
629 Vector<Float> ts;
630 specCol.get(i, spec);
631 tsysCol.get(i, ts);
632 if (mode == "MUL" || mode == "DIV") {
633 //if (mode == "DIV") val = 1.0/val;
634 spec *= val;
635 specCol.put(i, spec);
636 if ( tsys ) {
637 ts *= val;
638 tsysCol.put(i, ts);
639 }
640 } else if ( mode == "ADD" || mode == "SUB") {
641 //if (mode == "SUB") val *= -1.0;
642 spec += val;
643 specCol.put(i, spec);
644 if ( tsys ) {
645 ts += val;
646 tsysCol.put(i, ts);
647 }
648 }
649 }
650 return out;
651}
652
653CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
654 const std::vector<float> val,
655 const std::string& mode,
656 const std::string& opmode,
657 bool tsys )
658{
659 CountedPtr< Scantable > out ;
660 if ( opmode == "channel" ) {
661 out = arrayOperateChannel( in, val, mode, tsys ) ;
662 }
663 else if ( opmode == "row" ) {
664 out = arrayOperateRow( in, val, mode, tsys ) ;
665 }
666 else {
667 throw( AipsError( "Unknown array operation mode." ) ) ;
668 }
669 return out ;
670}
671
672CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
673 const std::vector<float> val,
674 const std::string& mode,
675 bool tsys )
676{
677 if ( val.size() == 1 ){
678 return unaryOperate( in, val[0], mode, tsys ) ;
679 }
680
681 // conformity of SPECTRA and TSYS
682 if ( tsys ) {
683 TableIterator titer(in->table(), "IFNO");
684 while ( !titer.pastEnd() ) {
685 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
686 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
687 Array<Float> spec = specCol.getColumn() ;
688 Array<Float> ts = tsysCol.getColumn() ;
689 if ( !spec.conform( ts ) ) {
690 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
691 }
692 titer.next() ;
693 }
694 }
695
696 // check if all spectra in the scantable have the same number of channel
697 vector<uInt> nchans;
698 vector<uInt> ifnos = in->getIFNos() ;
699 for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
700 nchans.push_back( in->nchan( ifnos[i] ) ) ;
701 }
702 Vector<uInt> mchans( nchans ) ;
703 if ( anyNE( mchans, mchans[0] ) ) {
704 throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
705 }
706
707 // check if vector size is equal to nchan
708 Vector<Float> fact( val ) ;
709 if ( fact.nelements() != mchans[0] ) {
710 throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
711 }
712
713 // check divided by zero
714 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
715 throw( AipsError("Divided by zero is not recommended." ) ) ;
716 }
717
718 CountedPtr< Scantable > out = getScantable(in, false);
719 Table& tab = out->table();
720 ArrayColumn<Float> specCol(tab,"SPECTRA");
721 ArrayColumn<Float> tsysCol(tab,"TSYS");
722 if (mode == "DIV") fact = (float)1.0 / fact;
723 else if (mode == "SUB") fact *= (float)-1.0 ;
724 for (uInt i=0; i<tab.nrow(); ++i) {
725 Vector<Float> spec;
726 Vector<Float> ts;
727 specCol.get(i, spec);
728 tsysCol.get(i, ts);
729 if (mode == "MUL" || mode == "DIV") {
730 //if (mode == "DIV") fact = (float)1.0 / fact;
731 spec *= fact;
732 specCol.put(i, spec);
733 if ( tsys ) {
734 ts *= fact;
735 tsysCol.put(i, ts);
736 }
737 } else if ( mode == "ADD" || mode == "SUB") {
738 //if (mode == "SUB") fact *= (float)-1.0 ;
739 spec += fact;
740 specCol.put(i, spec);
741 if ( tsys ) {
742 ts += fact;
743 tsysCol.put(i, ts);
744 }
745 }
746 }
747 return out;
748}
749
750CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
751 const std::vector<float> val,
752 const std::string& mode,
753 bool tsys )
754{
755 if ( val.size() == 1 ) {
756 return unaryOperate( in, val[0], mode, tsys ) ;
757 }
758
759 // conformity of SPECTRA and TSYS
760 if ( tsys ) {
761 TableIterator titer(in->table(), "IFNO");
762 while ( !titer.pastEnd() ) {
763 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
764 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
765 Array<Float> spec = specCol.getColumn() ;
766 Array<Float> ts = tsysCol.getColumn() ;
767 if ( !spec.conform( ts ) ) {
768 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
769 }
770 titer.next() ;
771 }
772 }
773
774 // check if vector size is equal to nrow
775 Vector<Float> fact( val ) ;
776 if (fact.nelements() != uInt(in->nrow())) {
777 throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
778 }
779
780 // check divided by zero
781 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
782 throw( AipsError("Divided by zero is not recommended." ) ) ;
783 }
784
785 CountedPtr< Scantable > out = getScantable(in, false);
786 Table& tab = out->table();
787 ArrayColumn<Float> specCol(tab,"SPECTRA");
788 ArrayColumn<Float> tsysCol(tab,"TSYS");
789 if (mode == "DIV") fact = (float)1.0 / fact;
790 if (mode == "SUB") fact *= (float)-1.0 ;
791 for (uInt i=0; i<tab.nrow(); ++i) {
792 Vector<Float> spec;
793 Vector<Float> ts;
794 specCol.get(i, spec);
795 tsysCol.get(i, ts);
796 if (mode == "MUL" || mode == "DIV") {
797 spec *= fact[i];
798 specCol.put(i, spec);
799 if ( tsys ) {
800 ts *= fact[i];
801 tsysCol.put(i, ts);
802 }
803 } else if ( mode == "ADD" || mode == "SUB") {
804 spec += fact[i];
805 specCol.put(i, spec);
806 if ( tsys ) {
807 ts += fact[i];
808 tsysCol.put(i, ts);
809 }
810 }
811 }
812 return out;
813}
814
815CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
816 const std::vector< std::vector<float> > val,
817 const std::string& mode,
818 bool tsys )
819{
820 // conformity of SPECTRA and TSYS
821 if ( tsys ) {
822 TableIterator titer(in->table(), "IFNO");
823 while ( !titer.pastEnd() ) {
824 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
825 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
826 Array<Float> spec = specCol.getColumn() ;
827 Array<Float> ts = tsysCol.getColumn() ;
828 if ( !spec.conform( ts ) ) {
829 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
830 }
831 titer.next() ;
832 }
833 }
834
835 // some checks
836 vector<uInt> nchans;
837 for (Int i = 0 ; i < in->nrow() ; i++) {
838 nchans.push_back((in->getSpectrum(i)).size());
839 }
840 //Vector<uInt> mchans( nchans ) ;
841 vector< Vector<Float> > facts ;
842 for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
843 Vector<Float> tmp( val[i] ) ;
844 // check divided by zero
845 if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
846 throw( AipsError("Divided by zero is not recommended." ) ) ;
847 }
848 // conformity check
849 if ( tmp.nelements() != nchans[i] ) {
850 stringstream ss ;
851 ss << "Row " << i << ": Vector size must be same as number of channel." ;
852 throw( AipsError( ss.str() ) ) ;
853 }
854 facts.push_back( tmp ) ;
855 }
856
857
858 CountedPtr< Scantable > out = getScantable(in, false);
859 Table& tab = out->table();
860 ArrayColumn<Float> specCol(tab,"SPECTRA");
861 ArrayColumn<Float> tsysCol(tab,"TSYS");
862 for (uInt i=0; i<tab.nrow(); ++i) {
863 Vector<Float> fact = facts[i] ;
864 Vector<Float> spec;
865 Vector<Float> ts;
866 specCol.get(i, spec);
867 tsysCol.get(i, ts);
868 if (mode == "MUL" || mode == "DIV") {
869 if (mode == "DIV") fact = (float)1.0 / fact;
870 spec *= fact;
871 specCol.put(i, spec);
872 if ( tsys ) {
873 ts *= fact;
874 tsysCol.put(i, ts);
875 }
876 } else if ( mode == "ADD" || mode == "SUB") {
877 if (mode == "SUB") fact *= (float)-1.0 ;
878 spec += fact;
879 specCol.put(i, spec);
880 if ( tsys ) {
881 ts += fact;
882 tsysCol.put(i, ts);
883 }
884 }
885 }
886 return out;
887}
888
889CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
890 const CountedPtr<Scantable>& right,
891 const std::string& mode)
892{
893 bool insitu = insitu_;
894 if ( ! left->conformant(*right) ) {
895 throw(AipsError("'left' and 'right' scantables are not conformant."));
896 }
897 setInsitu(false);
898 CountedPtr< Scantable > out = getScantable(left, false);
899 setInsitu(insitu);
900 Table& tout = out->table();
901 Block<String> coln(5);
902 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
903 coln[3] = "IFNO"; coln[4] = "POLNO";
904 Table tmpl = tout.sort(coln);
905 Table tmpr = right->table().sort(coln);
906 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
907 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
908 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
909 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
910
911 for (uInt i=0; i<tout.nrow(); ++i) {
912 Vector<Float> lspecvec, rspecvec;
913 Vector<uChar> lflagvec, rflagvec;
914 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
915 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
916 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
917 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
918 if (mode == "ADD") {
919 mleft += mright;
920 } else if ( mode == "SUB") {
921 mleft -= mright;
922 } else if ( mode == "MUL") {
923 mleft *= mright;
924 } else if ( mode == "DIV") {
925 mleft /= mright;
926 } else {
927 throw(AipsError("Illegal binary operator"));
928 }
929 lspecCol.put(i, mleft.getArray());
930 }
931 return out;
932}
933
934
935
936MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
937 const Vector<uChar>& f)
938{
939 Vector<Bool> mask;
940 mask.resize(f.shape());
941 convertArray(mask, f);
942 return MaskedArray<Float>(s,!mask);
943}
944
945MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
946 const Vector<uChar>& f)
947{
948 Vector<Bool> mask;
949 mask.resize(f.shape());
950 convertArray(mask, f);
951 return MaskedArray<Double>(s,!mask);
952}
953
954Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
955{
956 const Vector<Bool>& m = ma.getMask();
957 Vector<uChar> flags(m.shape());
958 convertArray(flags, !m);
959 return flags;
960}
961
962CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
963 const std::string & mode,
964 bool preserve )
965{
966 /// @todo make other modes available
967 /// modes should be "nearest", "pair"
968 // make this operation non insitu
969 (void) mode; //currently unused
970 const Table& tin = in->table();
971 Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
972 Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
973 if ( offs.nrow() == 0 )
974 throw(AipsError("No 'off' scans present."));
975 // put all "on" scans into output table
976
977 bool insitu = insitu_;
978 setInsitu(false);
979 CountedPtr< Scantable > out = getScantable(in, true);
980 setInsitu(insitu);
981 Table& tout = out->table();
982
983 TableCopy::copyRows(tout, ons);
984 TableRow row(tout);
985 ROScalarColumn<Double> offtimeCol(offs, "TIME");
986 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
987 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
988 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
989 for (uInt i=0; i < tout.nrow(); ++i) {
990 const TableRecord& rec = row.get(i);
991 Double ontime = rec.asDouble("TIME");
992 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
993 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
994 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
995 ROScalarColumn<Double> offtimeCol(presel, "TIME");
996
997 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
998 // Timestamp may vary within a cycle ???!!!
999 // increase this by 0.01 sec in case of rounding errors...
1000 // There might be a better way to do this.
1001 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
1002 mindeltat += 0.01/24./60./60.;
1003 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
1004
1005 if ( sel.nrow() < 1 ) {
1006 throw(AipsError("No closest in time found... This could be a rounding "
1007 "issue. Try quotient instead."));
1008 }
1009 TableRow offrow(sel);
1010 const TableRecord& offrec = offrow.get(0);//should only be one row
1011 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1012 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1013 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1014 /// @fixme this assumes tsys is a scalar not vector
1015 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1016 Vector<Float> specon, tsyson;
1017 outtsysCol.get(i, tsyson);
1018 outspecCol.get(i, specon);
1019 Vector<uChar> flagon;
1020 outflagCol.get(i, flagon);
1021 MaskedArray<Float> mon = maskedArray(specon, flagon);
1022 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
1023 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1024 if (preserve) {
1025 quot -= tsysoffscalar;
1026 } else {
1027 quot -= tsyson[0];
1028 }
1029 outspecCol.put(i, quot.getArray());
1030 outflagCol.put(i, flagsFromMA(quot));
1031 }
1032 // renumber scanno
1033 TableIterator it(tout, "SCANNO");
1034 uInt i = 0;
1035 while ( !it.pastEnd() ) {
1036 Table t = it.table();
1037 TableVector<uInt> vec(t, "SCANNO");
1038 vec = i;
1039 ++i;
1040 ++it;
1041 }
1042 return out;
1043}
1044
1045
1046CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
1047 const CountedPtr< Scantable > & off,
1048 bool preserve )
1049{
1050 bool insitu = insitu_;
1051 if ( ! on->conformant(*off) ) {
1052 throw(AipsError("'on' and 'off' scantables are not conformant."));
1053 }
1054 setInsitu(false);
1055 CountedPtr< Scantable > out = getScantable(on, false);
1056 setInsitu(insitu);
1057 Table& tout = out->table();
1058 const Table& toff = off->table();
1059 TableIterator sit(tout, "SCANNO");
1060 TableIterator s2it(toff, "SCANNO");
1061 while ( !sit.pastEnd() ) {
1062 Table ton = sit.table();
1063 TableRow row(ton);
1064 Table t = s2it.table();
1065 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1066 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
1067 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1068 for (uInt i=0; i < ton.nrow(); ++i) {
1069 const TableRecord& rec = row.get(i);
1070 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1071 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
1072 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
1073 if ( offsel.nrow() == 0 )
1074 throw AipsError("STMath::quotient: no matching off");
1075 TableRow offrow(offsel);
1076 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
1077 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1078 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1079 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1080 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1081 Vector<Float> specon, tsyson;
1082 outtsysCol.get(i, tsyson);
1083 outspecCol.get(i, specon);
1084 Vector<uChar> flagon;
1085 outflagCol.get(i, flagon);
1086 MaskedArray<Float> mon = maskedArray(specon, flagon);
1087 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
1088 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1089 if (preserve) {
1090 quot -= tsysoffscalar;
1091 } else {
1092 quot -= tsyson[0];
1093 }
1094 outspecCol.put(i, quot.getArray());
1095 outflagCol.put(i, flagsFromMA(quot));
1096 }
1097 ++sit;
1098 ++s2it;
1099 // take the first off for each on scan which doesn't have a
1100 // matching off scan
1101 // non <= noff: matching pairs, non > noff matching pairs then first off
1102 if ( s2it.pastEnd() ) s2it.reset();
1103 }
1104 return out;
1105}
1106
1107// dototalpower (migration of GBTIDL procedure dototalpower.pro)
1108// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
1109// do it for each cycles in a specific scan.
1110CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
1111 const CountedPtr< Scantable >& caloff, Float tcal )
1112{
1113 if ( ! calon->conformant(*caloff) ) {
1114 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
1115 }
1116 setInsitu(false);
1117 CountedPtr< Scantable > out = getScantable(caloff, false);
1118 Table& tout = out->table();
1119 const Table& tcon = calon->table();
1120 Vector<Float> tcalout;
1121
1122 std::map<uInt,uInt> tcalIdToRecNoMap;
1123 const Table& calOffTcalTable = caloff->tcal().table();
1124 {
1125 ROScalarColumn<uInt> calOffTcalTable_IDcol(calOffTcalTable, "ID");
1126 const Vector<uInt> tcalIds(calOffTcalTable_IDcol.getColumn());
1127 size_t tcalIdsEnd = tcalIds.nelements();
1128 for (uInt i = 0; i < tcalIdsEnd; i++) {
1129 tcalIdToRecNoMap[tcalIds[i]] = i;
1130 }
1131 }
1132 ROArrayColumn<Float> calOffTcalTable_TCALcol(calOffTcalTable, "TCAL");
1133
1134 if ( tout.nrow() != tcon.nrow() ) {
1135 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
1136 }
1137 // iteration by scanno or cycle no.
1138 TableIterator sit(tout, "SCANNO");
1139 TableIterator s2it(tcon, "SCANNO");
1140 while ( !sit.pastEnd() ) {
1141 Table toff = sit.table();
1142 TableRow row(toff);
1143 Table t = s2it.table();
1144 ScalarColumn<Double> outintCol(toff, "INTERVAL");
1145 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
1146 ArrayColumn<Float> outtsysCol(toff, "TSYS");
1147 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
1148 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
1149 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
1150 ROScalarColumn<Double> onintCol(t, "INTERVAL");
1151 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
1152 ROArrayColumn<Float> ontsysCol(t, "TSYS");
1153 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
1154 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
1155
1156 for (uInt i=0; i < toff.nrow(); ++i) {
1157 //skip these checks -> assumes the data order are the same between the cal on off pairs
1158 //
1159 Vector<Float> specCalon, specCaloff;
1160 // to store scalar (mean) tsys
1161 Vector<Float> tsysout(1);
1162 uInt tcalId, polno;
1163 Double offint, onint;
1164 outpolCol.get(i, polno);
1165 outspecCol.get(i, specCaloff);
1166 onspecCol.get(i, specCalon);
1167 Vector<uChar> flagCaloff, flagCalon;
1168 outflagCol.get(i, flagCaloff);
1169 onflagCol.get(i, flagCalon);
1170 outtcalIdCol.get(i, tcalId);
1171 outintCol.get(i, offint);
1172 onintCol.get(i, onint);
1173 // caluculate mean Tsys
1174 uInt nchan = specCaloff.nelements();
1175 // percentage of edge cut off
1176 uInt pc = 10;
1177 uInt bchan = nchan/pc;
1178 uInt echan = nchan-bchan;
1179
1180 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
1181 Vector<Float> testsubsp = specCaloff(chansl);
1182 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
1183 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
1184 MaskedArray<Float> spdiff = spon-spoff;
1185 uInt noff = spoff.nelementsValid();
1186 //uInt non = spon.nelementsValid();
1187 uInt ndiff = spdiff.nelementsValid();
1188 Float meantsys;
1189
1190/**
1191 Double subspec, subdiff;
1192 uInt usednchan;
1193 subspec = 0;
1194 subdiff = 0;
1195 usednchan = 0;
1196 for(uInt k=(bchan-1); k<echan; k++) {
1197 subspec += specCaloff[k];
1198 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
1199 ++usednchan;
1200 }
1201**/
1202 // get tcal if input tcal <= 0
1203 Float tcalUsed;
1204 tcalUsed = tcal;
1205 if ( tcal <= 0.0 ) {
1206 uInt tcalRecNo = tcalIdToRecNoMap[tcalId];
1207 calOffTcalTable_TCALcol.get(tcalRecNo, tcalout);
1208// if (polno<=3) {
1209// tcalUsed = tcalout[polno];
1210// }
1211// else {
1212// tcalUsed = tcalout[0];
1213// }
1214 if ( tcalout.size() == 1 )
1215 tcalUsed = tcalout[0] ;
1216 else if ( tcalout.size() == nchan )
1217 tcalUsed = mean(tcalout) ;
1218 else {
1219 uInt ipol = polno ;
1220 if ( ipol > 3 ) ipol = 0 ;
1221 tcalUsed = tcalout[ipol] ;
1222 }
1223 }
1224
1225 Float meanoff;
1226 Float meandiff;
1227 if (noff && ndiff) {
1228 //Debug
1229 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
1230 //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ;
1231 //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
1232 meanoff = sum(spoff)/noff;
1233 meandiff = sum(spdiff)/ndiff;
1234 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
1235 }
1236 else {
1237 meantsys=1;
1238 }
1239
1240 tsysout[0] = Float(meantsys);
1241 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
1242 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
1243 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
1244 //uInt ncaloff = mcaloff.nelementsValid();
1245 //uInt ncalon = mcalon.nelementsValid();
1246
1247 outintCol.put(i, offint+onint);
1248 outspecCol.put(i, sig.getArray());
1249 outflagCol.put(i, flagsFromMA(sig));
1250 outtsysCol.put(i, tsysout);
1251 }
1252 ++sit;
1253 ++s2it;
1254 }
1255 return out;
1256}
1257
1258//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
1259// observatiions.
1260// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
1261// no smoothing).
1262// output: resultant scantable [= (sig-ref/ref)*tsys]
1263CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
1264 const CountedPtr < Scantable >& ref,
1265 int smoothref,
1266 casa::Float tsysv,
1267 casa::Float tau )
1268{
1269 LogIO os( casa::LogOrigin( "STMath", "dosigref()"));
1270if ( ! ref->conformant(*sig) ) {
1271 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
1272 }
1273 setInsitu(false);
1274 CountedPtr< Scantable > out = getScantable(sig, false);
1275 CountedPtr< Scantable > smref;
1276 if ( smoothref > 1 ) {
1277 float fsmoothref = static_cast<float>(smoothref);
1278 std::string inkernel = "boxcar";
1279 smref = smooth(ref, inkernel, fsmoothref );
1280 ostringstream oss;
1281 os <<"Applied smoothing of "<<fsmoothref<<" on the reference."
1282 << LogIO::POST;
1283 }
1284 else {
1285 smref = ref;
1286 }
1287 Table& tout = out->table();
1288 const Table& tref = smref->table();
1289 if ( tout.nrow() != tref.nrow() ) {
1290 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
1291 }
1292 // iteration by scanno? or cycle no.
1293 TableIterator sit(tout, "SCANNO");
1294 TableIterator s2it(tref, "SCANNO");
1295 while ( !sit.pastEnd() ) {
1296 Table ton = sit.table();
1297 Table t = s2it.table();
1298 ScalarColumn<Double> outintCol(ton, "INTERVAL");
1299 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1300 ArrayColumn<Float> outtsysCol(ton, "TSYS");
1301 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1302 ArrayColumn<Float> refspecCol(t, "SPECTRA");
1303 ROScalarColumn<Double> refintCol(t, "INTERVAL");
1304 ROArrayColumn<Float> reftsysCol(t, "TSYS");
1305 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
1306 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
1307 for (uInt i=0; i < ton.nrow(); ++i) {
1308
1309 Double onint, refint;
1310 Vector<Float> specon, specref;
1311 // to store scalar (mean) tsys
1312 Vector<Float> tsysref;
1313 outintCol.get(i, onint);
1314 refintCol.get(i, refint);
1315 outspecCol.get(i, specon);
1316 refspecCol.get(i, specref);
1317 Vector<uChar> flagref, flagon;
1318 outflagCol.get(i, flagon);
1319 refflagCol.get(i, flagref);
1320 reftsysCol.get(i, tsysref);
1321
1322 Float tsysrefscalar;
1323 if ( tsysv > 0.0 ) {
1324 ostringstream oss;
1325 Float elev;
1326 refelevCol.get(i, elev);
1327 os << "user specified Tsys = " << tsysv;
1328 // do recalc elevation if EL = 0
1329 if ( elev == 0 ) {
1330 throw(AipsError("EL=0, elevation data is missing."));
1331 } else {
1332 if ( tau <= 0.0 ) {
1333 throw(AipsError("Valid tau is not supplied."));
1334 } else {
1335 tsysrefscalar = tsysv * exp(tau/elev);
1336 }
1337 }
1338 os << ", corrected (for El) tsys= "<<tsysrefscalar;
1339 }
1340 else {
1341 tsysrefscalar = tsysref[0];
1342 }
1343 //get quotient spectrum
1344 MaskedArray<Float> mref = maskedArray(specref, flagref);
1345 MaskedArray<Float> mon = maskedArray(specon, flagon);
1346 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
1347 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
1348
1349 //Debug
1350 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
1351 //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
1352 //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
1353 // fill the result, replay signal tsys by reference tsys
1354 outintCol.put(i, resint);
1355 outspecCol.put(i, specres.getArray());
1356 outflagCol.put(i, flagsFromMA(specres));
1357 outtsysCol.put(i, tsysref);
1358 }
1359 ++sit;
1360 ++s2it;
1361 }
1362
1363 out->setFluxUnit("K");
1364
1365 return out;
1366}
1367
1368CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
1369 const std::vector<int>& scans,
1370 int smoothref,
1371 casa::Float tsysv,
1372 casa::Float tau,
1373 casa::Float tcal )
1374
1375{
1376 setInsitu(false);
1377 LogIO os( casa::LogOrigin( "STMath", "donod()"));
1378 STSelector sel;
1379 std::vector<int> scan1, scan2, beams, types;
1380 std::vector< vector<int> > scanpair;
1381 //std::vector<string> calstate;
1382 std::vector<int> calstate;
1383 String msg;
1384
1385 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
1386 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
1387
1388 std::vector< CountedPtr< Scantable > > sctables;
1389 sctables.push_back(s1b1on);
1390 sctables.push_back(s1b1off);
1391 sctables.push_back(s1b2on);
1392 sctables.push_back(s1b2off);
1393 sctables.push_back(s2b1on);
1394 sctables.push_back(s2b1off);
1395 sctables.push_back(s2b2on);
1396 sctables.push_back(s2b2off);
1397
1398 //check scanlist
1399 int n=s->checkScanInfo(scans);
1400 if (n==1) {
1401 throw(AipsError("Incorrect scan pairs. "));
1402 }
1403
1404 // Assume scans contain only a pair of consecutive scan numbers.
1405 // It is assumed that first beam, b1, is on target.
1406 // There is no check if the first beam is on or not.
1407 if ( scans.size()==1 ) {
1408 scan1.push_back(scans[0]);
1409 scan2.push_back(scans[0]+1);
1410 } else if ( scans.size()==2 ) {
1411 scan1.push_back(scans[0]);
1412 scan2.push_back(scans[1]);
1413 } else {
1414 if ( scans.size()%2 == 0 ) {
1415 for (uInt i=0; i<scans.size(); i++) {
1416 if (i%2 == 0) {
1417 scan1.push_back(scans[i]);
1418 }
1419 else {
1420 scan2.push_back(scans[i]);
1421 }
1422 }
1423 } else {
1424 throw(AipsError("Odd numbers of scans, cannot form pairs."));
1425 }
1426 }
1427 scanpair.push_back(scan1);
1428 scanpair.push_back(scan2);
1429 //calstate.push_back("*calon");
1430 //calstate.push_back("*[^calon]");
1431 calstate.push_back(SrcType::NODCAL);
1432 calstate.push_back(SrcType::NOD);
1433 CountedPtr< Scantable > ws = getScantable(s, false);
1434 uInt l=0;
1435 while ( l < sctables.size() ) {
1436 for (uInt i=0; i < 2; i++) {
1437 for (uInt j=0; j < 2; j++) {
1438 for (uInt k=0; k < 2; k++) {
1439 sel.reset();
1440 sel.setScans(scanpair[i]);
1441 //sel.setName(calstate[k]);
1442 types.clear();
1443 types.push_back(calstate[k]);
1444 sel.setTypes(types);
1445 beams.clear();
1446 beams.push_back(j);
1447 sel.setBeams(beams);
1448 ws->setSelection(sel);
1449 sctables[l]= getScantable(ws, false);
1450 l++;
1451 }
1452 }
1453 }
1454 }
1455
1456 // replace here by splitData or getData functionality
1457 CountedPtr< Scantable > sig1;
1458 CountedPtr< Scantable > ref1;
1459 CountedPtr< Scantable > sig2;
1460 CountedPtr< Scantable > ref2;
1461 CountedPtr< Scantable > calb1;
1462 CountedPtr< Scantable > calb2;
1463
1464 msg=String("Processing dototalpower for subset of the data");
1465 os << msg << LogIO::POST;
1466 // Debug for IRC CS data
1467 //float tcal1=7.0;
1468 //float tcal2=4.0;
1469 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
1470 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
1471 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
1472 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
1473
1474 // correction of user-specified tsys for elevation here
1475
1476 // dosigref calibration
1477 msg=String("Processing dosigref for subset of the data");
1478 os << msg << endl;
1479 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
1480 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
1481
1482 // iteration by scanno or cycle no.
1483 Table& tcalb1 = calb1->table();
1484 Table& tcalb2 = calb2->table();
1485 TableIterator sit(tcalb1, "SCANNO");
1486 TableIterator s2it(tcalb2, "SCANNO");
1487 while ( !sit.pastEnd() ) {
1488 Table t1 = sit.table();
1489 Table t2= s2it.table();
1490 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
1491 ArrayColumn<Float> outtsysCol(t1, "TSYS");
1492 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
1493 ScalarColumn<Double> outintCol(t1, "INTERVAL");
1494 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
1495 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
1496 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
1497 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
1498 for (uInt i=0; i < t1.nrow(); ++i) {
1499 Vector<Float> spec1, spec2;
1500 // to store scalar (mean) tsys
1501 Vector<Float> tsys1, tsys2;
1502 Vector<uChar> flag1, flag2;
1503 Double tint1, tint2;
1504 outspecCol.get(i, spec1);
1505 t2specCol.get(i, spec2);
1506 outflagCol.get(i, flag1);
1507 t2flagCol.get(i, flag2);
1508 outtsysCol.get(i, tsys1);
1509 t2tsysCol.get(i, tsys2);
1510 outintCol.get(i, tint1);
1511 t2intCol.get(i, tint2);
1512 // average
1513 // assume scalar tsys for weights
1514 Float wt1, wt2, tsyssq1, tsyssq2;
1515 tsyssq1 = tsys1[0]*tsys1[0];
1516 tsyssq2 = tsys2[0]*tsys2[0];
1517 wt1 = Float(tint1)/tsyssq1;
1518 wt2 = Float(tint2)/tsyssq2;
1519 Float invsumwt=1/(wt1+wt2);
1520 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
1521 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
1522 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
1523 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
1524 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
1525 // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
1526 // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
1527 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
1528 Array<Float> avtsys = tsys1;
1529
1530 outspecCol.put(i, avspec.getArray());
1531 outflagCol.put(i, flagsFromMA(avspec));
1532 outtsysCol.put(i, avtsys);
1533 }
1534 ++sit;
1535 ++s2it;
1536 }
1537
1538 calb1->setFluxUnit("K");
1539
1540 return calb1;
1541}
1542
1543//GBTIDL version of frequency switched data calibration
1544CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
1545 const std::vector<int>& scans,
1546 int smoothref,
1547 casa::Float tsysv,
1548 casa::Float tau,
1549 casa::Float tcal )
1550{
1551
1552
1553 (void) scans; //currently unused
1554 STSelector sel;
1555 CountedPtr< Scantable > ws = getScantable(s, false);
1556 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
1557 CountedPtr< Scantable > calsig, calref, out, out1, out2;
1558 Bool nofold=False;
1559 vector<int> types ;
1560
1561 //split the data
1562 //sel.setName("*_fs");
1563 types.push_back( SrcType::FSON ) ;
1564 sel.setTypes( types ) ;
1565 ws->setSelection(sel);
1566 sig = getScantable(ws,false);
1567 sel.reset();
1568 types.clear() ;
1569 //sel.setName("*_fs_calon");
1570 types.push_back( SrcType::FONCAL ) ;
1571 sel.setTypes( types ) ;
1572 ws->setSelection(sel);
1573 sigwcal = getScantable(ws,false);
1574 sel.reset();
1575 types.clear() ;
1576 //sel.setName("*_fsr");
1577 types.push_back( SrcType::FSOFF ) ;
1578 sel.setTypes( types ) ;
1579 ws->setSelection(sel);
1580 ref = getScantable(ws,false);
1581 sel.reset();
1582 types.clear() ;
1583 //sel.setName("*_fsr_calon");
1584 types.push_back( SrcType::FOFFCAL ) ;
1585 sel.setTypes( types ) ;
1586 ws->setSelection(sel);
1587 refwcal = getScantable(ws,false);
1588 sel.reset() ;
1589 types.clear() ;
1590
1591 calsig = dototalpower(sigwcal, sig, tcal=tcal);
1592 calref = dototalpower(refwcal, ref, tcal=tcal);
1593
1594 out1=dosigref(calsig,calref,smoothref,tsysv,tau);
1595 out2=dosigref(calref,calsig,smoothref,tsysv,tau);
1596
1597 Table& tabout1=out1->table();
1598 Table& tabout2=out2->table();
1599 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
1600 ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
1601 ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
1602 Vector<Float> spec; specCol.get(0, spec);
1603 uInt nchan = spec.nelements();
1604 uInt freqid1; freqidCol1.get(0,freqid1);
1605 uInt freqid2; freqidCol2.get(0,freqid2);
1606 Double rp1, rp2, rv1, rv2, inc1, inc2;
1607 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1608 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1609 //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1610 //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1611 //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
1612 if (rp1==rp2) {
1613 Double foffset = rv1 - rv2;
1614 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1615 if (choffset >= nchan) {
1616 //cerr<<"out-band frequency switching, no folding"<<endl;
1617 LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1618 os<<"out-band frequency switching, no folding"<<LogIO::POST;
1619 nofold = True;
1620 }
1621 }
1622
1623 if (nofold) {
1624 std::vector< CountedPtr< Scantable > > tabs;
1625 tabs.push_back(out1);
1626 tabs.push_back(out2);
1627 out = merge(tabs);
1628 }
1629 else {
1630 //out = out1;
1631 Double choffset = ( rv1 - rv2 ) / inc2 ;
1632 out = dofold( out1, out2, choffset ) ;
1633 }
1634
1635 out->setFluxUnit("K");
1636
1637 return out;
1638}
1639
1640CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1641 const CountedPtr<Scantable> &ref,
1642 Double choffset,
1643 Double choffset2 )
1644{
1645 LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
1646 os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
1647
1648 // output scantable
1649 CountedPtr<Scantable> out = getScantable( sig, false ) ;
1650
1651 // separate choffset to integer part and decimal part
1652 Int ioffset = (Int)choffset ;
1653 Double doffset = choffset - ioffset ;
1654 Int ioffset2 = (Int)choffset2 ;
1655 Double doffset2 = choffset2 - ioffset2 ;
1656 os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
1657 os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ;
1658
1659 // get column
1660 ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1661 ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1662 ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1663 ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1664 ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1665 ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1666 ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1667 ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1668 ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1669 ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1670
1671 // check
1672 if ( ioffset == 0 ) {
1673 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1674 os << "channel offset is zero, no folding" << LogIO::POST ;
1675 return out ;
1676 }
1677 int nchan = ref->nchan() ;
1678 if ( abs(ioffset) >= nchan ) {
1679 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1680 os << "out-band frequency switching, no folding" << LogIO::POST ;
1681 return out ;
1682 }
1683
1684 // attach column for output scantable
1685 ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1686 ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1687 ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1688 ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1689 ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1690 ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
1691
1692 // for each row
1693 // assume that the data order are same between sig and ref
1694 RowAccumulator acc( asap::W_TINTSYS ) ;
1695 for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1696 // get values
1697 Vector<Float> spsig ;
1698 specCol1.get( i, spsig ) ;
1699 Vector<Float> spref ;
1700 specCol2.get( i, spref ) ;
1701 Vector<Float> tsyssig ;
1702 tsysCol1.get( i, tsyssig ) ;
1703 Vector<Float> tsysref ;
1704 tsysCol2.get( i, tsysref ) ;
1705 Vector<uChar> flagsig ;
1706 flagCol1.get( i, flagsig ) ;
1707 Vector<uChar> flagref ;
1708 flagCol2.get( i, flagref ) ;
1709 Double timesig ;
1710 mjdCol1.get( i, timesig ) ;
1711 Double timeref ;
1712 mjdCol2.get( i, timeref ) ;
1713 Double intsig ;
1714 intervalCol1.get( i, intsig ) ;
1715 Double intref ;
1716 intervalCol2.get( i, intref ) ;
1717
1718 // shift reference spectra
1719 int refchan = spref.nelements() ;
1720 Vector<Float> sspref( spref.nelements() ) ;
1721 Vector<Float> stsysref( tsysref.nelements() ) ;
1722 Vector<uChar> sflagref( flagref.nelements() ) ;
1723 if ( ioffset > 0 ) {
1724 // SPECTRA and FLAGTRA
1725 for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1726 sspref[j] = spref[j+ioffset] ;
1727 sflagref[j] = flagref[j+ioffset] ;
1728 }
1729 for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1730 sspref[j] = spref[j-refchan+ioffset] ;
1731 sflagref[j] = flagref[j-refchan+ioffset] ;
1732 }
1733 spref = sspref.copy() ;
1734 flagref = sflagref.copy() ;
1735 for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1736 sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
1737 sflagref[j] = flagref[j+1] + flagref[j] ;
1738 }
1739 sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
1740 sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
1741
1742 // TSYS
1743 if ( spref.nelements() == tsysref.nelements() ) {
1744 for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1745 stsysref[j] = tsysref[j+ioffset] ;
1746 }
1747 for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1748 stsysref[j] = tsysref[j-refchan+ioffset] ;
1749 }
1750 tsysref = stsysref.copy() ;
1751 for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1752 stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
1753 }
1754 stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
1755 }
1756 }
1757 else {
1758 // SPECTRA and FLAGTRA
1759 for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1760 sspref[j] = spref[refchan+ioffset+j] ;
1761 sflagref[j] = flagref[refchan+ioffset+j] ;
1762 }
1763 for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1764 sspref[j] = spref[j+ioffset] ;
1765 sflagref[j] = flagref[j+ioffset] ;
1766 }
1767 spref = sspref.copy() ;
1768 flagref = sflagref.copy() ;
1769 sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
1770 sflagref[0] = flagref[0] + flagref[refchan-1] ;
1771 for ( int j = 1 ; j < refchan ; j++ ) {
1772 sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
1773 sflagref[j] = flagref[j-1] + flagref[j] ;
1774 }
1775 // TSYS
1776 if ( spref.nelements() == tsysref.nelements() ) {
1777 for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1778 stsysref[j] = tsysref[refchan+ioffset+j] ;
1779 }
1780 for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1781 stsysref[j] = tsysref[j+ioffset] ;
1782 }
1783 tsysref = stsysref.copy() ;
1784 stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
1785 for ( int j = 1 ; j < refchan ; j++ ) {
1786 stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
1787 }
1788 }
1789 }
1790
1791 // shift signal spectra if necessary (only for APEX?)
1792 if ( choffset2 != 0.0 ) {
1793 int sigchan = spsig.nelements() ;
1794 Vector<Float> sspsig( spsig.nelements() ) ;
1795 Vector<Float> stsyssig( tsyssig.nelements() ) ;
1796 Vector<uChar> sflagsig( flagsig.nelements() ) ;
1797 if ( ioffset2 > 0 ) {
1798 // SPECTRA and FLAGTRA
1799 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1800 sspsig[j] = spsig[j+ioffset2] ;
1801 sflagsig[j] = flagsig[j+ioffset2] ;
1802 }
1803 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1804 sspsig[j] = spsig[j-sigchan+ioffset2] ;
1805 sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
1806 }
1807 spsig = sspsig.copy() ;
1808 flagsig = sflagsig.copy() ;
1809 for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1810 sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
1811 sflagsig[j] = flagsig[j+1] || flagsig[j] ;
1812 }
1813 sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
1814 sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ;
1815 // TSTS
1816 if ( spsig.nelements() == tsyssig.nelements() ) {
1817 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1818 stsyssig[j] = tsyssig[j+ioffset2] ;
1819 }
1820 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1821 stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
1822 }
1823 tsyssig = stsyssig.copy() ;
1824 for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1825 stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1826 }
1827 stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
1828 }
1829 }
1830 else {
1831 // SPECTRA and FLAGTRA
1832 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1833 sspsig[j] = spsig[sigchan+ioffset2+j] ;
1834 sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
1835 }
1836 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1837 sspsig[j] = spsig[j+ioffset2] ;
1838 sflagsig[j] = flagsig[j+ioffset2] ;
1839 }
1840 spsig = sspsig.copy() ;
1841 flagsig = sflagsig.copy() ;
1842 sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
1843 sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
1844 for ( int j = 1 ; j < sigchan ; j++ ) {
1845 sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
1846 sflagsig[j] = flagsig[j-1] + flagsig[j] ;
1847 }
1848 // TSYS
1849 if ( spsig.nelements() == tsyssig.nelements() ) {
1850 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1851 stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
1852 }
1853 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1854 stsyssig[j] = tsyssig[j+ioffset2] ;
1855 }
1856 tsyssig = stsyssig.copy() ;
1857 stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
1858 for ( int j = 1 ; j < sigchan ; j++ ) {
1859 stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1860 }
1861 }
1862 }
1863 }
1864
1865 // folding
1866 acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1867 acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
1868
1869 // put result
1870 specColOut.put( i, acc.getSpectrum() ) ;
1871 const Vector<Bool> &msk = acc.getMask() ;
1872 Vector<uChar> flg( msk.shape() ) ;
1873 convertArray( flg, !msk ) ;
1874 flagColOut.put( i, flg ) ;
1875 tsysColOut.put( i, acc.getTsys() ) ;
1876 intervalColOut.put( i, acc.getInterval() ) ;
1877 mjdColOut.put( i, acc.getTime() ) ;
1878 // change FREQ_ID to unshifted IF setting (only for APEX?)
1879 if ( choffset2 != 0.0 ) {
1880 uInt freqid = fidColOut( 0 ) ; // assume single-IF data
1881 double refpix, refval, increment ;
1882 out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
1883 refval -= choffset * increment ;
1884 uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
1885 Vector<uInt> freqids = fidColOut.getColumn() ;
1886 for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
1887 if ( freqids[j] == freqid )
1888 freqids[j] = newfreqid ;
1889 }
1890 fidColOut.putColumn( freqids ) ;
1891 }
1892
1893 acc.reset() ;
1894 }
1895
1896 return out ;
1897}
1898
1899
1900CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1901{
1902 // make copy or reference
1903 CountedPtr< Scantable > out = getScantable(in, false);
1904 Table& tout = out->table();
1905 Block<String> cols(4);
1906 cols[0] = String("SCANNO");
1907 cols[1] = String("CYCLENO");
1908 cols[2] = String("BEAMNO");
1909 cols[3] = String("POLNO");
1910 TableIterator iter(tout, cols);
1911 while (!iter.pastEnd()) {
1912 Table subt = iter.table();
1913 // this should leave us with two rows for the two IFs....if not ignore
1914 if (subt.nrow() != 2 ) {
1915 continue;
1916 }
1917 ArrayColumn<Float> specCol(subt, "SPECTRA");
1918 ArrayColumn<Float> tsysCol(subt, "TSYS");
1919 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1920 Vector<Float> onspec,offspec, ontsys, offtsys;
1921 Vector<uChar> onflag, offflag;
1922 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1923 specCol.get(0, onspec); specCol.get(1, offspec);
1924 flagCol.get(0, onflag); flagCol.get(1, offflag);
1925 MaskedArray<Float> on = maskedArray(onspec, onflag);
1926 MaskedArray<Float> off = maskedArray(offspec, offflag);
1927 MaskedArray<Float> oncopy = on.copy();
1928
1929 on /= off; on -= 1.0f;
1930 on *= ontsys[0];
1931 off /= oncopy; off -= 1.0f;
1932 off *= offtsys[0];
1933 specCol.put(0, on.getArray());
1934 const Vector<Bool>& m0 = on.getMask();
1935 Vector<uChar> flags0(m0.shape());
1936 convertArray(flags0, !m0);
1937 flagCol.put(0, flags0);
1938
1939 specCol.put(1, off.getArray());
1940 const Vector<Bool>& m1 = off.getMask();
1941 Vector<uChar> flags1(m1.shape());
1942 convertArray(flags1, !m1);
1943 flagCol.put(1, flags1);
1944 ++iter;
1945 }
1946
1947 return out;
1948}
1949
1950std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1951 const std::vector< bool > & mask,
1952 const std::string& which )
1953{
1954
1955 Vector<Bool> m(mask);
1956 const Table& tab = in->table();
1957 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1958 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1959 std::vector<float> out;
1960 for (uInt i=0; i < tab.nrow(); ++i ) {
1961 Vector<Float> spec; specCol.get(i, spec);
1962 Vector<uChar> flag; flagCol.get(i, flag);
1963 MaskedArray<Float> ma = maskedArray(spec, flag);
1964 float outstat = 0.0;
1965 if ( spec.nelements() == m.nelements() ) {
1966 outstat = mathutil::statistics(which, ma(m));
1967 } else {
1968 outstat = mathutil::statistics(which, ma);
1969 }
1970 out.push_back(outstat);
1971 }
1972 return out;
1973}
1974
1975std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
1976 const std::vector< bool > & mask,
1977 const std::string& which,
1978 int row )
1979{
1980
1981 Vector<Bool> m(mask);
1982 const Table& tab = in->table();
1983 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1984 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1985 std::vector<float> out;
1986
1987 Vector<Float> spec; specCol.get(row, spec);
1988 Vector<uChar> flag; flagCol.get(row, flag);
1989 MaskedArray<Float> ma = maskedArray(spec, flag);
1990 float outstat = 0.0;
1991 if ( spec.nelements() == m.nelements() ) {
1992 outstat = mathutil::statistics(which, ma(m));
1993 } else {
1994 outstat = mathutil::statistics(which, ma);
1995 }
1996 out.push_back(outstat);
1997
1998 return out;
1999}
2000
2001std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
2002 const std::vector< bool > & mask,
2003 const std::string& which )
2004{
2005
2006 Vector<Bool> m(mask);
2007 const Table& tab = in->table();
2008 ROArrayColumn<Float> specCol(tab, "SPECTRA");
2009 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2010 std::vector<int> out;
2011 for (uInt i=0; i < tab.nrow(); ++i ) {
2012 Vector<Float> spec; specCol.get(i, spec);
2013 Vector<uChar> flag; flagCol.get(i, flag);
2014 MaskedArray<Float> ma = maskedArray(spec, flag);
2015 if (ma.ndim() != 1) {
2016 throw (ArrayError(
2017 "std::vector<int> STMath::minMaxChan("
2018 "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
2019 " std::string &which)"
2020 " - MaskedArray is not 1D"));
2021 }
2022 IPosition outpos(1,0);
2023 if ( spec.nelements() == m.nelements() ) {
2024 outpos = mathutil::minMaxPos(which, ma(m));
2025 } else {
2026 outpos = mathutil::minMaxPos(which, ma);
2027 }
2028 out.push_back(outpos[0]);
2029 }
2030 return out;
2031}
2032
2033CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
2034 int width )
2035{
2036 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
2037 CountedPtr< Scantable > out = getScantable(in, false);
2038 Table& tout = out->table();
2039 out->frequencies().rescale(width, "BIN");
2040 ArrayColumn<Float> specCol(tout, "SPECTRA");
2041 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
2042 ArrayColumn<Float> tsysCol(tout, "TSYS");
2043
2044 for (uInt i=0; i < tout.nrow(); ++i ) {
2045 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
2046 MaskedArray<Float> maout;
2047 LatticeUtilities::bin(maout, main, 0, Int(width));
2048 specCol.put(i, maout.getArray());
2049 flagCol.put(i, flagsFromMA(maout));
2050 if (tsysCol(i).nelements() == specCol(i).nelements()) {
2051 MaskedArray<Float> matsysin = maskedArray(tsysCol(i), flagCol(i));
2052 MaskedArray<Float> matsysout;
2053 LatticeUtilities::bin(matsysout, matsysin, 0, Int(width));
2054 tsysCol.put(i, matsysout.getArray());
2055 }
2056 // take only the first binned spectrum's length for the deprecated
2057 // global header item nChan
2058 if (i==0) tout.rwKeywordSet().define(String("nChan"),
2059 Int(maout.getArray().nelements()));
2060 }
2061 return out;
2062}
2063
2064CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
2065 const std::string& method,
2066 float width )
2067//
2068// Should add the possibility of width being specified in km/s. This means
2069// that for each freqID (SpectralCoordinate) we will need to convert to an
2070// average channel width (say at the reference pixel). Then we would need
2071// to be careful to make sure each spectrum (of different freqID)
2072// is the same length.
2073//
2074{
2075 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
2076 Int interpMethod(stringToIMethod(method));
2077
2078 CountedPtr< Scantable > out = getScantable(in, false);
2079 Table& tout = out->table();
2080
2081// Resample SpectralCoordinates (one per freqID)
2082 out->frequencies().rescale(width, "RESAMPLE");
2083 TableIterator iter(tout, "IFNO");
2084 TableRow row(tout);
2085 while ( !iter.pastEnd() ) {
2086 Table tab = iter.table();
2087 ArrayColumn<Float> specCol(tab, "SPECTRA");
2088 //ArrayColumn<Float> tsysCol(tout, "TSYS");
2089 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2090 Vector<Float> spec;
2091 Vector<uChar> flag;
2092 specCol.get(0,spec); // the number of channels should be constant per IF
2093 uInt nChanIn = spec.nelements();
2094 Vector<Float> xIn(nChanIn); indgen(xIn);
2095 Int fac = Int(nChanIn/width);
2096 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
2097 uInt k = 0;
2098 Float x = 0.0;
2099 while (x < Float(nChanIn) ) {
2100 xOut(k) = x;
2101 k++;
2102 x += width;
2103 }
2104 uInt nChanOut = k;
2105 xOut.resize(nChanOut, True);
2106 // process all rows for this IFNO
2107 Vector<Float> specOut;
2108 Vector<Bool> maskOut;
2109 Vector<uChar> flagOut;
2110 for (uInt i=0; i < tab.nrow(); ++i) {
2111 specCol.get(i, spec);
2112 flagCol.get(i, flag);
2113 Vector<Bool> mask(flag.nelements());
2114 convertArray(mask, flag);
2115
2116 IPosition shapeIn(spec.shape());
2117 //sh.nchan = nChanOut;
2118 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2119 xIn, spec, mask,
2120 interpMethod, True, True);
2121 /// @todo do the same for channel based Tsys
2122 flagOut.resize(maskOut.nelements());
2123 convertArray(flagOut, maskOut);
2124 specCol.put(i, specOut);
2125 flagCol.put(i, flagOut);
2126 }
2127 ++iter;
2128 }
2129
2130 return out;
2131}
2132
2133STMath::imethod STMath::stringToIMethod(const std::string& in)
2134{
2135 static STMath::imap lookup;
2136
2137 // initialize the lookup table if necessary
2138 if ( lookup.empty() ) {
2139 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
2140 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2141 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
2142 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
2143 }
2144
2145 STMath::imap::const_iterator iter = lookup.find(in);
2146
2147 if ( lookup.end() == iter ) {
2148 std::string message = in;
2149 message += " is not a valid interpolation mode";
2150 throw(AipsError(message));
2151 }
2152 return iter->second;
2153}
2154
2155WeightType STMath::stringToWeight(const std::string& in)
2156{
2157 static std::map<std::string, WeightType> lookup;
2158
2159 // initialize the lookup table if necessary
2160 if ( lookup.empty() ) {
2161 lookup["NONE"] = asap::W_NONE;
2162 lookup["TINT"] = asap::W_TINT;
2163 lookup["TINTSYS"] = asap::W_TINTSYS;
2164 lookup["TSYS"] = asap::W_TSYS;
2165 lookup["VAR"] = asap::W_VAR;
2166 }
2167
2168 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2169
2170 if ( lookup.end() == iter ) {
2171 std::string message = in;
2172 message += " is not a valid weighting mode";
2173 throw(AipsError(message));
2174 }
2175 return iter->second;
2176}
2177
2178CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2179 const vector< float > & coeff,
2180 const std::string & filename,
2181 const std::string& method)
2182{
2183 LogIO os( LogOrigin( "STMath", "gainElevation", WHERE ) ) ;
2184 // Get elevation data from Scantable and convert to degrees
2185 CountedPtr< Scantable > out = getScantable(in, false);
2186 Table& tab = out->table();
2187 ROScalarColumn<Float> elev(tab, "ELEVATION");
2188 Vector<Float> x = elev.getColumn();
2189 x *= Float(180 / C::pi); // Degrees
2190
2191 Vector<Float> coeffs(coeff);
2192 const uInt nc = coeffs.nelements();
2193 if ( filename.length() > 0 && nc > 0 ) {
2194 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2195 }
2196
2197 // Correct
2198 if ( nc > 0 || filename.length() == 0 ) {
2199 // Find instrument
2200 Bool throwit = True;
2201 Instrument inst =
2202 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2203 throwit);
2204
2205 // Set polynomial
2206 Polynomial<Float>* ppoly = 0;
2207 Vector<Float> coeff;
2208 String msg;
2209 if ( nc > 0 ) {
2210 ppoly = new Polynomial<Float>(nc-1);
2211 coeff = coeffs;
2212 msg = String("user");
2213 } else {
2214 STAttr sdAttr;
2215 coeff = sdAttr.gainElevationPoly(inst);
2216 ppoly = new Polynomial<Float>(coeff.nelements()-1);
2217 msg = String("built in");
2218 }
2219
2220 if ( coeff.nelements() > 0 ) {
2221 ppoly->setCoefficients(coeff);
2222 } else {
2223 delete ppoly;
2224 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2225 }
2226 os << "Making polynomial correction with " << msg << " coefficients:" << endl;
2227 os << " " << coeff << LogIO::POST;
2228 const uInt nrow = tab.nrow();
2229 Vector<Float> factor(nrow);
2230 for ( uInt i=0; i < nrow; ++i ) {
2231 factor[i] = 1.0 / (*ppoly)(x[i]);
2232 }
2233 delete ppoly;
2234 scaleByVector(tab, factor, true);
2235
2236 } else {
2237 // Read and correct
2238 os << "Making correction from ascii Table" << LogIO::POST;
2239 scaleFromAsciiTable(tab, filename, method, x, true);
2240 }
2241 return out;
2242}
2243
2244void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2245 const std::string& method,
2246 const Vector<Float>& xout, bool dotsys)
2247{
2248
2249// Read gain-elevation ascii file data into a Table.
2250
2251 String formatString;
2252 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2253 scaleFromTable(in, tbl, method, xout, dotsys);
2254}
2255
2256void STMath::scaleFromTable(Table& in,
2257 const Table& table,
2258 const std::string& method,
2259 const Vector<Float>& xout, bool dotsys)
2260{
2261
2262 ROScalarColumn<Float> geElCol(table, "ELEVATION");
2263 ROScalarColumn<Float> geFacCol(table, "FACTOR");
2264 Vector<Float> xin = geElCol.getColumn();
2265 Vector<Float> yin = geFacCol.getColumn();
2266 Vector<Bool> maskin(xin.nelements(),True);
2267
2268 // Interpolate (and extrapolate) with desired method
2269
2270 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2271
2272 Vector<Float> yout;
2273 Vector<Bool> maskout;
2274 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2275 xin, yin, maskin, interp,
2276 True, True);
2277
2278 scaleByVector(in, Float(1.0)/yout, dotsys);
2279}
2280
2281void STMath::scaleByVector( Table& in,
2282 const Vector< Float >& factor,
2283 bool dotsys )
2284{
2285 uInt nrow = in.nrow();
2286 if ( factor.nelements() != nrow ) {
2287 throw(AipsError("factors.nelements() != table.nelements()"));
2288 }
2289 ArrayColumn<Float> specCol(in, "SPECTRA");
2290 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2291 ArrayColumn<Float> tsysCol(in, "TSYS");
2292 for (uInt i=0; i < nrow; ++i) {
2293 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2294 ma *= factor[i];
2295 specCol.put(i, ma.getArray());
2296 flagCol.put(i, flagsFromMA(ma));
2297 if ( dotsys ) {
2298 Vector<Float> tsys = tsysCol(i);
2299 tsys *= factor[i];
2300 tsysCol.put(i,tsys);
2301 }
2302 }
2303}
2304
2305CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2306 float d, float etaap,
2307 float jyperk )
2308{
2309 LogIO os( LogOrigin( "STMath", "convertFlux", WHERE ) ) ;
2310
2311 CountedPtr< Scantable > out = getScantable(in, false);
2312 Table& tab = in->table();
2313 Table& outtab = out->table();
2314 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2315 Unit K(String("K"));
2316 Unit JY(String("Jy"));
2317
2318 bool tokelvin = true;
2319 Double cfac = 1.0;
2320
2321 if ( fluxUnit == JY ) {
2322 os << "Converting to K" << LogIO::POST;
2323 Quantum<Double> t(1.0,fluxUnit);
2324 Quantum<Double> t2 = t.get(JY);
2325 cfac = (t2 / t).getValue(); // value to Jy
2326
2327 tokelvin = true;
2328 out->setFluxUnit("K");
2329 } else if ( fluxUnit == K ) {
2330 os << "Converting to Jy" << LogIO::POST;
2331 Quantum<Double> t(1.0,fluxUnit);
2332 Quantum<Double> t2 = t.get(K);
2333 cfac = (t2 / t).getValue(); // value to K
2334
2335 tokelvin = false;
2336 out->setFluxUnit("Jy");
2337 } else {
2338 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2339 }
2340 // Make sure input values are converted to either Jy or K first...
2341 Float factor = cfac;
2342
2343 // Select method
2344 if (jyperk > 0.0) {
2345 factor *= jyperk;
2346 if ( tokelvin ) factor = 1.0 / jyperk;
2347 os << "Jy/K = " << jyperk << LogIO::POST;
2348 Vector<Float> factors(outtab.nrow(), factor);
2349 scaleByVector(outtab,factors, false);
2350 } else if ( etaap > 0.0) {
2351 if (d < 0) {
2352 Instrument inst =
2353 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2354 True);
2355 STAttr sda;
2356 d = sda.diameter(inst);
2357 }
2358 jyperk = STAttr::findJyPerK(etaap, d);
2359 os << "Jy/K = " << jyperk << LogIO::POST;
2360 factor *= jyperk;
2361 if ( tokelvin ) {
2362 factor = 1.0 / factor;
2363 }
2364 Vector<Float> factors(outtab.nrow(), factor);
2365 scaleByVector(outtab, factors, False);
2366 } else {
2367
2368 // OK now we must deal with automatic look up of values.
2369 // We must also deal with the fact that the factors need
2370 // to be computed per IF and may be different and may
2371 // change per integration.
2372
2373 os <<"Looking up conversion factors" << LogIO::POST;
2374 convertBrightnessUnits(out, tokelvin, cfac);
2375 }
2376
2377 return out;
2378}
2379
2380void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2381 bool tokelvin, float cfac )
2382{
2383 Table& table = in->table();
2384 Instrument inst =
2385 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2386 TableIterator iter(table, "FREQ_ID");
2387 STFrequencies stfreqs = in->frequencies();
2388 STAttr sdAtt;
2389 while (!iter.pastEnd()) {
2390 Table tab = iter.table();
2391 ArrayColumn<Float> specCol(tab, "SPECTRA");
2392 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2393 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2394 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2395
2396 uInt freqid; freqidCol.get(0, freqid);
2397 Vector<Float> tmpspec; specCol.get(0, tmpspec);
2398 // STAttr.JyPerK has a Vector interface... change sometime.
2399 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2400 for ( uInt i=0; i<tab.nrow(); ++i) {
2401 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2402 Float factor = cfac * jyperk;
2403 if ( tokelvin ) factor = Float(1.0) / factor;
2404 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2405 ma *= factor;
2406 specCol.put(i, ma.getArray());
2407 flagCol.put(i, flagsFromMA(ma));
2408 }
2409 ++iter;
2410 }
2411}
2412
2413CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2414 const std::vector<float>& tau )
2415{
2416 CountedPtr< Scantable > out = getScantable(in, false);
2417
2418 Table outtab = out->table();
2419
2420 const Int ntau = uInt(tau.size());
2421 std::vector<float>::const_iterator tauit = tau.begin();
2422 AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
2423 AipsError);
2424 TableIterator iiter(outtab, "IFNO");
2425 while ( !iiter.pastEnd() ) {
2426 Table itab = iiter.table();
2427 TableIterator piter(itab, "POLNO");
2428 while ( !piter.pastEnd() ) {
2429 Table tab = piter.table();
2430 ROScalarColumn<Float> elev(tab, "ELEVATION");
2431 ArrayColumn<Float> specCol(tab, "SPECTRA");
2432 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2433 ArrayColumn<Float> tsysCol(tab, "TSYS");
2434 for ( uInt i=0; i<tab.nrow(); ++i) {
2435 Float zdist = Float(C::pi_2) - elev(i);
2436 Float factor = exp(*tauit/cos(zdist));
2437 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2438 ma *= factor;
2439 specCol.put(i, ma.getArray());
2440 flagCol.put(i, flagsFromMA(ma));
2441 Vector<Float> tsys;
2442 tsysCol.get(i, tsys);
2443 tsys *= factor;
2444 tsysCol.put(i, tsys);
2445 }
2446 if (ntau == in->nif()*in->npol() ) {
2447 tauit++;
2448 }
2449 piter++;
2450 }
2451 if (ntau >= in->nif() ) {
2452 tauit++;
2453 }
2454 iiter++;
2455 }
2456 return out;
2457}
2458
2459CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2460 const std::string& kernel,
2461 float width, int order)
2462{
2463 CountedPtr< Scantable > out = getScantable(in, false);
2464 Table table = out->table();
2465
2466 TableIterator iter(table, "IFNO");
2467 while (!iter.pastEnd()) {
2468 Table tab = iter.table();
2469 ArrayColumn<Float> specCol(tab, "SPECTRA");
2470 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2471 Vector<Float> spec;
2472 Vector<uChar> flag;
2473 for (uInt i = 0; i < tab.nrow(); ++i) {
2474 specCol.get(i, spec);
2475 flagCol.get(i, flag);
2476 Vector<Bool> mask(flag.nelements());
2477 convertArray(mask, flag);
2478 Vector<Float> specout;
2479 Vector<Bool> maskout;
2480 if (kernel == "hanning") {
2481 mathutil::hanning(specout, maskout, spec, !mask);
2482 } else if (kernel == "rmedian") {
2483 mathutil::runningMedian(specout, maskout, spec , mask, width);
2484 } else if (kernel == "poly") {
2485 mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2486 }
2487
2488 for (uInt j = 0; j < flag.nelements(); ++j) {
2489 uChar userFlag = 1 << 7;
2490 if (maskout[j]==True) userFlag = 0 << 7;
2491 flag(j) = userFlag;
2492 }
2493
2494 flagCol.put(i, flag);
2495 specCol.put(i, specout);
2496 }
2497 ++iter;
2498 }
2499 return out;
2500}
2501
2502CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2503 const std::string& kernel, float width,
2504 int order)
2505{
2506 if (kernel == "rmedian" || kernel == "hanning" || kernel == "poly") {
2507 return smoothOther(in, kernel, width, order);
2508 }
2509 CountedPtr< Scantable > out = getScantable(in, false);
2510 Table& table = out->table();
2511 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2512 // same IFNO should have same no of channels
2513 // this saves overhead
2514 TableIterator iter(table, "IFNO");
2515 while (!iter.pastEnd()) {
2516 Table tab = iter.table();
2517 ArrayColumn<Float> specCol(tab, "SPECTRA");
2518 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2519 Vector<Float> spec = specCol( 0 );
2520 uInt nchan = spec.nelements();
2521 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2522 Convolver<Float> conv(kvec, IPosition(1,nchan));
2523 Vector<uChar> flag;
2524 Vector<Bool> mask(nchan);
2525 for ( uInt i=0; i<tab.nrow(); ++i) {
2526 specCol.get(i, spec);
2527 flagCol.get(i, flag);
2528 convertArray(mask, flag);
2529 Vector<Float> specout;
2530 //mathutil::replaceMaskByZero(specout, mask);
2531 mathutil::replaceMaskByZero(spec, !mask);
2532 //std::vector<bool> vmask;
2533 //(!mask).tovector(vmask);
2534 //mathutil::doZeroOrderInterpolation(spec, vmask);
2535 conv.linearConv(specout, spec);
2536 specCol.put(i, specout);
2537 }
2538 ++iter;
2539 }
2540 return out;
2541}
2542
2543CountedPtr< Scantable >
2544STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
2545 const std::string &freqTol )
2546{
2547 Double freqTolInHz = 0.0; // default is 0.0Hz (merge only when exact match)
2548 if (freqTol.size() > 0) {
2549 Quantum<Double> freqTolInQuantity;
2550 if (!Quantum<Double>::read(freqTolInQuantity, freqTol)) {
2551 throw(AipsError("Failed to convert freqTol string to quantity"));
2552 }
2553 if (!freqTolInQuantity.isConform("Hz")) {
2554 throw(AipsError("Invalid freqTol string"));
2555 }
2556 freqTolInHz = freqTolInQuantity.getValue("Hz");
2557 LogIO os(LogOrigin("STMath", "merge", WHERE));
2558 os << "frequency tolerance = " << freqTolInHz << "Hz" << LogIO::POST;
2559 }
2560
2561 if ( in.size() < 2 ) {
2562 throw(AipsError("Need at least two scantables to perform a merge."));
2563 }
2564 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2565 bool insitu = insitu_;
2566 setInsitu(false);
2567 CountedPtr< Scantable > out = getScantable(*it, false);
2568 setInsitu(insitu);
2569 Table& tout = out->table();
2570 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2571 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2572 ScalarColumn<uInt> ifnocol(tout, "IFNO");
2573 // Renumber SCANNO to be 0-based
2574 Vector<uInt> scannos = scannocol.getColumn();
2575 uInt offset = min(scannos);
2576 scannos -= offset;
2577 scannocol.putColumn(scannos);
2578 uInt newscanno = max(scannos)+1;
2579 ++it;
2580
2581 // new IFNO
2582 uInt ifnoCounter = max(ifnocol.getColumn()) + 1;
2583
2584 // Here we assume that each IFNO has unique MOLECULE_ID
2585 // molIdMap:
2586 // KEY: IFNO
2587 // VALUE: MOLECULE_ID
2588 map<uInt, uInt> molIdMap;
2589 {
2590 TableIterator ifit(tout, "IFNO");
2591 while (!ifit.pastEnd()) {
2592 ROTableRow row(ifit.table());
2593 const TableRecord& rec = row.get(0);
2594 molIdMap[rec.asuInt("IFNO")] = rec.asuInt("MOLECULE_ID");
2595 ifit.next();
2596 }
2597 }
2598
2599 while ( it != in.end() ){
2600 // Check FREQUENCIES/BASEFRAME
2601 if ( out->frequencies().getFrame(true) != (*it)->frequencies().getFrame(true) ) {
2602 throw(AipsError("BASEFRAME is not identical"));
2603 }
2604
2605 if ( ! (*it)->conformant(*out) ) {
2606 // non conformant.
2607 LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2608 os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2609 }
2610 out->appendToHistoryTable((*it)->history());
2611 const Table& tab = (*it)->table();
2612
2613 Block<String> cols(3);
2614 cols[0] = String("FREQ_ID");
2615 cols[1] = String("MOLECULE_ID");
2616 cols[2] = String("FOCUS_ID");
2617
2618 TableIterator scanit(tab, "SCANNO");
2619 while (!scanit.pastEnd()) {
2620 TableIterator subit(scanit.table(), cols);
2621 while ( !subit.pastEnd() ) {
2622 uInt nrow = tout.nrow();
2623 Table thetab = subit.table();
2624 ROTableRow row(thetab);
2625 Vector<uInt> thecolvals(thetab.nrow());
2626 // The selected subset of table should have
2627 // the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2628 const TableRecord& rec = row.get(0);
2629
2630 // Set the proper FREQ_ID
2631 Double rv,rp,inc;
2632 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2633 uInt id;
2634
2635 // default value is new unique IFNO
2636 uInt newifno = ifnoCounter;
2637 Bool isIfMerged = False;
2638 uInt nchan = rec.asArrayFloat("SPECTRA").shape()[0];
2639 //id = out->frequencies().addEntry(rp, rv, inc);
2640 if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
2641 // add new entry to FREQUENCIES table
2642 id = out->frequencies().addEntry(rp, rv, inc);
2643
2644 // increment counter for IFNO
2645 ifnoCounter++;
2646 }
2647 else {
2648 // should renumber IFNO to be same as existing rows that have same FREQ_ID
2649 LogIO os(LogOrigin("STMath", "merge", WHERE));
2650 Table outFreqIdSelected = tout(tout.col("FREQ_ID") == id);
2651 TableIterator _iter(outFreqIdSelected, "IFNO");
2652 map<uInt, uInt> nchanMap;
2653 while (!_iter.pastEnd()) {
2654 const Table _table = _iter.table();
2655 ROTableRow _row(_table);
2656 const TableRecord &_rec = _row.get(0);
2657 uInt nchan = _rec.asArrayFloat("SPECTRA").shape()[0];
2658 if (nchanMap.find(nchan) != nchanMap.end()) {
2659 throw(AipsError("There are non-unique IFNOs assigned to spectra that have same FREQ_ID and same nchan. Something wrong."));
2660 }
2661 nchanMap[nchan] = _rec.asuInt("IFNO");
2662 _iter.next();
2663 }
2664
2665 os << LogIO::DEBUGGING << "nchanMap for " << id << ":" << LogIO::POST;
2666 for (map<uInt, uInt>::iterator i = nchanMap.begin(); i != nchanMap.end(); ++i) {
2667 os << LogIO::DEBUGGING << "nchanMap[" << i->first << "] = " << i->second << LogIO::POST;
2668 }
2669
2670 if (nchanMap.find(nchan) == nchanMap.end()) {
2671 // increment counter for IFNO
2672 ifnoCounter++;
2673 }
2674 else {
2675 // renumber IFNO to be same as existing value that corresponds to nchan
2676 newifno = nchanMap[nchan];
2677 isIfMerged = True;
2678 }
2679 os << LogIO::DEBUGGING << "newifno = " << newifno << LogIO::POST;
2680 }
2681 // copy rows to output table
2682 tout.addRow(thetab.nrow());
2683 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2684
2685 Slicer slice(IPosition(1, nrow), IPosition(1, thetab.nrow()),
2686 Slicer::endIsLength);
2687
2688 thecolvals = id;
2689 freqidcol.putColumnRange(slice, thecolvals);
2690
2691 thecolvals = newifno;
2692 ifnocol.putColumnRange(slice, thecolvals);
2693
2694 // Set the proper MOLECULE_ID
2695 Vector<String> name,fname;Vector<Double> rf;
2696 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2697 id = out->molecules().addEntry(rf, name, fname);
2698 if (molIdMap.find(newifno) == molIdMap.end()) {
2699 // add new entry to molIdMap
2700 molIdMap[newifno] = id;
2701 }
2702 if (isIfMerged) {
2703 id = molIdMap[newifno];
2704 }
2705 thecolvals = id;
2706 molidcol.putColumnRange(slice, thecolvals);
2707
2708 // Set the proper FOCUS_ID
2709 Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2710 (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2711 fxy, fxyp, rec.asuInt("FOCUS_ID"));
2712 id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2713 fxy, fxyp);
2714 thecolvals = id;
2715 focusidcol.putColumnRange(slice, thecolvals);
2716
2717 // Set the proper SCANNO
2718 thecolvals = newscanno;
2719 scannocol.putColumnRange(slice, thecolvals);
2720
2721 ++subit;
2722 }
2723 ++newscanno;
2724 ++scanit;
2725 }
2726 ++it;
2727 }
2728 return out;
2729}
2730
2731CountedPtr< Scantable >
2732 STMath::invertPhase( const CountedPtr < Scantable >& in )
2733{
2734 return applyToPol(in, &STPol::invertPhase, Float(0.0));
2735}
2736
2737CountedPtr< Scantable >
2738 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2739{
2740 return applyToPol(in, &STPol::rotatePhase, Float(phase));
2741}
2742
2743CountedPtr< Scantable >
2744 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2745{
2746 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2747}
2748
2749CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2750 STPol::polOperation fptr,
2751 Float phase )
2752{
2753 CountedPtr< Scantable > out = getScantable(in, false);
2754 Table& tout = out->table();
2755 Block<String> cols(4);
2756 cols[0] = String("SCANNO");
2757 cols[1] = String("BEAMNO");
2758 cols[2] = String("IFNO");
2759 cols[3] = String("CYCLENO");
2760 TableIterator iter(tout, cols);
2761 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2762 out->getPolType() );
2763 while (!iter.pastEnd()) {
2764 Table t = iter.table();
2765 ArrayColumn<Float> speccol(t, "SPECTRA");
2766 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2767 Matrix<Float> pols(speccol.getColumn());
2768 try {
2769 stpol->setSpectra(pols);
2770 Float fang,fhand;
2771 fang = in->focusTable_.getTotalAngle(focidcol(0));
2772 fhand = in->focusTable_.getFeedHand(focidcol(0));
2773 stpol->setPhaseCorrections(fang, fhand);
2774 // use a member function pointer in STPol. This only works on
2775 // the STPol pointer itself, not the Counted Pointer so
2776 // derefernce it.
2777 (&(*(stpol))->*fptr)(phase);
2778 speccol.putColumn(stpol->getSpectra());
2779 } catch (AipsError& e) {
2780 //delete stpol;stpol=0;
2781 throw(e);
2782 }
2783 ++iter;
2784 }
2785 //delete stpol;stpol=0;
2786 return out;
2787}
2788
2789CountedPtr< Scantable >
2790 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2791{
2792 CountedPtr< Scantable > out = getScantable(in, false);
2793 Table& tout = out->table();
2794 Table t0 = tout(tout.col("POLNO") == 0);
2795 Table t1 = tout(tout.col("POLNO") == 1);
2796 if ( t0.nrow() != t1.nrow() )
2797 throw(AipsError("Inconsistent number of polarisations"));
2798 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2799 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2800 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2801 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2802 Matrix<Float> s0 = speccol0.getColumn();
2803 Matrix<uChar> f0 = flagcol0.getColumn();
2804 speccol0.putColumn(speccol1.getColumn());
2805 flagcol0.putColumn(flagcol1.getColumn());
2806 speccol1.putColumn(s0);
2807 flagcol1.putColumn(f0);
2808 return out;
2809}
2810
2811CountedPtr< Scantable >
2812 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2813 const std::vector<bool>& mask,
2814 const std::string& weight )
2815{
2816 if (in->npol() < 2 )
2817 throw(AipsError("averagePolarisations can only be applied to two or more"
2818 "polarisations"));
2819 bool insitu = insitu_;
2820 setInsitu(false);
2821 CountedPtr< Scantable > pols = getScantable(in, true);
2822 setInsitu(insitu);
2823 Table& tout = pols->table();
2824 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2825 Table tab = tableCommand(taql, in->table());
2826 if (tab.nrow() == 0 )
2827 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2828 TableCopy::copyRows(tout, tab);
2829 TableVector<uInt> vec(tout, "POLNO");
2830 vec = 0;
2831 pols->table_.rwKeywordSet().define("nPol", Int(1));
2832 pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2833 //pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2834 std::vector<CountedPtr<Scantable> > vpols;
2835 vpols.push_back(pols);
2836 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2837 return out;
2838}
2839
2840CountedPtr< Scantable >
2841 STMath::averageBeams( const CountedPtr< Scantable > & in,
2842 const std::vector<bool>& mask,
2843 const std::string& weight )
2844{
2845 bool insitu = insitu_;
2846 setInsitu(false);
2847 CountedPtr< Scantable > beams = getScantable(in, false);
2848 setInsitu(insitu);
2849 Table& tout = beams->table();
2850 // give all rows the same BEAMNO
2851 TableVector<uInt> vec(tout, "BEAMNO");
2852 vec = 0;
2853 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2854 std::vector<CountedPtr<Scantable> > vbeams;
2855 vbeams.push_back(beams);
2856 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2857 return out;
2858}
2859
2860
2861CountedPtr< Scantable >
2862 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2863 const std::string & refTime,
2864 const std::string & method)
2865{
2866 LogIO os( casa::LogOrigin("STMath", "frequencyAlign()", WHERE));
2867 // clone as this is not working insitu
2868 bool insitu = insitu_;
2869 setInsitu(false);
2870 CountedPtr< Scantable > out = getScantable(in, false);
2871 setInsitu(insitu);
2872 Table& tout = out->table();
2873 // Get reference Epoch to time of first row or given String
2874 Unit DAY(String("d"));
2875 MEpoch::Ref epochRef(in->getTimeReference());
2876 MEpoch refEpoch;
2877 if (refTime.length()>0) {
2878 Quantum<Double> qt;
2879 if (MVTime::read(qt,refTime)) {
2880 MVEpoch mv(qt);
2881 refEpoch = MEpoch(mv, epochRef);
2882 } else {
2883 throw(AipsError("Invalid format for Epoch string"));
2884 }
2885 } else {
2886 refEpoch = in->timeCol_(0);
2887 }
2888 MPosition refPos = in->getAntennaPosition();
2889
2890 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2891 /*
2892 // Comment from MV.
2893 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2894 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2895 // test if user frame is different to base frame
2896 if ( in->frequencies().getFrameString(true)
2897 == in->frequencies().getFrameString(false) ) {
2898 throw(AipsError("Can't convert as no output frame has been set"
2899 " (use set_freqframe) or it is aligned already."));
2900 }
2901 */
2902 MFrequency::Types system = in->frequencies().getFrame();
2903 MVTime mvt(refEpoch.getValue());
2904 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2905 os << "Aligning at reference Epoch " << epochout
2906 << " in frame " << MFrequency::showType(system) << LogIO::POST;
2907 // set up the iterator
2908 Block<String> cols(4);
2909 // select by constant direction
2910 cols[0] = String("SRCNAME");
2911 cols[1] = String("BEAMNO");
2912 // select by IF ( no of channels varies over this )
2913 cols[2] = String("IFNO");
2914 // select by restfrequency
2915 cols[3] = String("MOLECULE_ID");
2916 TableIterator iter(tout, cols);
2917 while ( !iter.pastEnd() ) {
2918 Table t = iter.table();
2919 ROScalarColumn<String> snCol(t, "SRCNAME");
2920 os << "Aligning to position of source '" << snCol(0) << "'" << LogIO::POST;
2921 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2922 TableIterator fiter(t, "FREQ_ID");
2923 // determine nchan from the first row. This should work as
2924 // we are iterating over BEAMNO and IFNO
2925 // we should have constant direction
2926
2927 ROArrayColumn<Float> sCol(t, "SPECTRA");
2928 const MDirection direction = dirCol(0);
2929 const uInt nchan = sCol(0).nelements();
2930
2931 // skip operations if there is nothing to align
2932 if (fiter.pastEnd()) {
2933 continue;
2934 }
2935
2936 Table ftab = fiter.table();
2937 // align all frequency ids with respect to the first encountered id
2938 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2939 // get the SpectralCoordinate for the freqid, which we are iterating over
2940 SpectralCoordinate sC = \
2941 in->frequencies().getSpectralCoordinate(freqidCol(0));
2942 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2943 direction, refPos, system );
2944 // realign the SpectralCoordinate and put into the output Scantable
2945 Vector<String> units(1);
2946 units = String("Hz");
2947 Bool linear=True;
2948 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2949 sc2.setWorldAxisUnits(units);
2950 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2951 sc2.referenceValue()[0],
2952 sc2.increment()[0]);
2953 while ( !fiter.pastEnd() ) {
2954
2955 ftab = fiter.table();
2956 // spectral coordinate for the current FREQ_ID
2957 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2958 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2959 // create the "global" abcissa for alignment with same FREQ_ID
2960 Vector<Double> abc(nchan);
2961 for (uInt i=0; i<nchan; i++) {
2962 Double w;
2963 sC.toWorld(w,Double(i));
2964 abc[i] = w;
2965 }
2966 TableVector<uInt> tvec(ftab, "FREQ_ID");
2967 // assign new frequency id to all rows
2968 tvec = id;
2969 // cache abcissa for same time stamps, so iterate over those
2970 TableIterator timeiter(ftab, "TIME");
2971 while ( !timeiter.pastEnd() ) {
2972 Table tab = timeiter.table();
2973 ArrayColumn<Float> specCol(tab, "SPECTRA");
2974 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2975 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2976 // use align abcissa cache after the first row
2977 // these rows should be just be POLNO
2978 bool first = true;
2979 for (int i=0; i<int(tab.nrow()); ++i) {
2980 // input values
2981 Vector<uChar> flag = flagCol(i);
2982 Vector<Bool> mask(flag.shape());
2983 Vector<Float> specOut, spec;
2984 spec = specCol(i);
2985 Vector<Bool> maskOut;Vector<uChar> flagOut;
2986 convertArray(mask, flag);
2987 // alignment
2988 Bool ok = fa.align(specOut, maskOut, abc, spec,
2989 mask, timeCol(i), !first,
2990 interp, False);
2991 (void) ok; // unused stop compiler nagging
2992 // back into scantable
2993 flagOut.resize(maskOut.nelements());
2994 convertArray(flagOut, maskOut);
2995 flagCol.put(i, flagOut);
2996 specCol.put(i, specOut);
2997 // start abcissa caching
2998 first = false;
2999 }
3000 // next timestamp
3001 ++timeiter;
3002 }
3003 // next FREQ_ID
3004 ++fiter;
3005 }
3006 // next aligner
3007 ++iter;
3008 }
3009 // set this afterwards to ensure we are doing insitu correctly.
3010 out->frequencies().setFrame(system, true);
3011 return out;
3012}
3013
3014CountedPtr<Scantable>
3015 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
3016 const std::string & newtype )
3017{
3018 if (in->npol() != 2 && in->npol() != 4)
3019 throw(AipsError("Can only convert two or four polarisations."));
3020 if ( in->getPolType() == newtype )
3021 throw(AipsError("No need to convert."));
3022 if ( ! in->selector_.empty() )
3023 throw(AipsError("Can only convert whole scantable. Unset the selection."));
3024 bool insitu = insitu_;
3025 setInsitu(false);
3026 CountedPtr< Scantable > out = getScantable(in, true);
3027 setInsitu(insitu);
3028 Table& tout = out->table();
3029 tout.rwKeywordSet().define("POLTYPE", String(newtype));
3030
3031 Block<String> cols(4);
3032 cols[0] = "SCANNO";
3033 cols[1] = "CYCLENO";
3034 cols[2] = "BEAMNO";
3035 cols[3] = "IFNO";
3036 TableIterator it(in->originalTable_, cols);
3037 String basetype = in->getPolType();
3038 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
3039 try {
3040 while ( !it.pastEnd() ) {
3041 Table tab = it.table();
3042 uInt row = tab.rowNumbers()[0];
3043 stpol->setSpectra(in->getPolMatrix(row));
3044 Float fang,fhand;
3045 fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
3046 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
3047 stpol->setPhaseCorrections(fang, fhand);
3048 Int npolout = 0;
3049 for (uInt i=0; i<tab.nrow(); ++i) {
3050 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
3051 if ( outvec.nelements() > 0 ) {
3052 tout.addRow();
3053 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
3054 ArrayColumn<Float> sCol(tout,"SPECTRA");
3055 ScalarColumn<uInt> pCol(tout,"POLNO");
3056 sCol.put(tout.nrow()-1 ,outvec);
3057 pCol.put(tout.nrow()-1 ,uInt(npolout));
3058 npolout++;
3059 }
3060 }
3061 tout.rwKeywordSet().define("nPol", npolout);
3062 ++it;
3063 }
3064 } catch (AipsError& e) {
3065 delete stpol;
3066 throw(e);
3067 }
3068 delete stpol;
3069 return out;
3070}
3071
3072CountedPtr< Scantable >
3073 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
3074 const std::string & scantype )
3075{
3076 bool insitu = insitu_;
3077 setInsitu(false);
3078 CountedPtr< Scantable > out = getScantable(in, true);
3079 setInsitu(insitu);
3080 Table& tout = out->table();
3081 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
3082 if (scantype == "on") {
3083 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
3084 }
3085 Table tab = tableCommand(taql, in->table());
3086 TableCopy::copyRows(tout, tab);
3087 if (scantype == "on") {
3088 // re-index SCANNO to 0
3089 TableVector<uInt> vec(tout, "SCANNO");
3090 vec = 0;
3091 }
3092 return out;
3093}
3094
3095std::vector<float>
3096 asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
3097 const std::vector<int>& whichrow,
3098 bool getRealImag )
3099{
3100 std::vector<float> res;
3101 Table tab = in->table();
3102 std::vector<bool> mask;
3103
3104 if (whichrow.size() < 1) { // for all rows (by default)
3105 int nrow = int(tab.nrow());
3106 for (int i = 0; i < nrow; ++i) {
3107 res = in->execFFT(i, mask, getRealImag);
3108 }
3109 } else { // for specified rows
3110 for (uInt i = 0; i < whichrow.size(); ++i) {
3111 res = in->execFFT(i, mask, getRealImag);
3112 }
3113 }
3114
3115 return res;
3116}
3117
3118
3119CountedPtr<Scantable>
3120 asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
3121 double start, double end,
3122 const std::string& mode )
3123{
3124 CountedPtr<Scantable> out = getScantable(in, false);
3125 Table& tout = out->table();
3126 TableIterator iter(tout, "FREQ_ID");
3127 FFTServer<Float,Complex> ffts;
3128
3129 while ( !iter.pastEnd() ) {
3130 Table tab = iter.table();
3131 Double rp,rv,inc;
3132 ROTableRow row(tab);
3133 const TableRecord& rec = row.get(0);
3134 uInt freqid = rec.asuInt("FREQ_ID");
3135 out->frequencies().getEntry(rp, rv, inc, freqid);
3136 ArrayColumn<Float> specCol(tab, "SPECTRA");
3137 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3138
3139 for (int i=0; i<int(tab.nrow()); ++i) {
3140 Vector<Float> spec = specCol(i);
3141 Vector<uChar> flag = flagCol(i);
3142 std::vector<bool> mask;
3143 for (uInt j = 0; j < flag.nelements(); ++j) {
3144 mask.push_back(!(flag[j]>0));
3145 }
3146 mathutil::doZeroOrderInterpolation(spec, mask);
3147
3148 Vector<Complex> lags;
3149 ffts.fft0(lags, spec);
3150
3151 Int lag0(start+0.5);
3152 Int lag1(end+0.5);
3153 if (mode == "frequency") {
3154 lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3155 lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3156 }
3157 Int lstart = max(0, lag0);
3158 Int lend = min(Int(lags.nelements()-1), lag1);
3159 if (lstart == lend) {
3160 lags[lstart] = Complex(0.0);
3161 } else {
3162 if (lstart > lend) {
3163 Int tmp = lend;
3164 lend = lstart;
3165 lstart = tmp;
3166 }
3167 for (int j=lstart; j <=lend ;++j) {
3168 lags[j] = Complex(0.0);
3169 }
3170 }
3171
3172 ffts.fft0(spec, lags);
3173
3174 specCol.put(i, spec);
3175 }
3176 ++iter;
3177 }
3178 return out;
3179}
3180
3181// Averaging spectra with different channel/resolution
3182CountedPtr<Scantable>
3183STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3184 const bool& compel,
3185 const std::vector<bool>& mask,
3186 const std::string& weight,
3187 const std::string& avmode )
3188 throw ( casa::AipsError )
3189{
3190 LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3191 if ( avmode == "SCAN" && in.size() != 1 )
3192 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3193 "Use merge first."));
3194
3195 CountedPtr<Scantable> out ; // processed result
3196 if ( compel ) {
3197 std::vector< CountedPtr<Scantable> > newin ; // input for average process
3198 uInt insize = in.size() ; // number of input scantables
3199
3200 // setup newin
3201 bool oldInsitu = insitu_ ;
3202 setInsitu( false ) ;
3203 newin.resize( insize ) ;
3204 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3205 newin[itable] = getScantable( in[itable], false ) ;
3206 }
3207 setInsitu( oldInsitu ) ;
3208
3209 // warning
3210 os << "Average spectra with different spectral resolution" << LogIO::POST ;
3211
3212 // temporarily set coordinfo
3213 vector<string> oldinfo( insize ) ;
3214 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3215 vector<string> coordinfo = in[itable]->getCoordInfo() ;
3216 oldinfo[itable] = coordinfo[0] ;
3217 coordinfo[0] = "Hz" ;
3218 newin[itable]->setCoordInfo( coordinfo ) ;
3219 }
3220
3221 ostringstream oss ;
3222
3223 // check IF frequency coverage
3224 // freqid: list of FREQ_ID, which is used, in each table
3225 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
3226 // each table
3227 // freqid[insize][numIF]
3228 // freqid: [[id00, id01, ...],
3229 // [id10, id11, ...],
3230 // ...
3231 // [idn0, idn1, ...]]
3232 // iffreq[insize][numIF*2]
3233 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3234 // [min_id10, max_id10, min_id11, max_id11, ...],
3235 // ...
3236 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3237 //os << "Check IF settings in each table" << LogIO::POST ;
3238 vector< vector<uInt> > freqid( insize );
3239 vector< vector<double> > iffreq( insize ) ;
3240 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3241 Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
3242 vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
3243 for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
3244 i != uniqueFreqId.end() ; i++ ) {
3245 //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3246 uInt target = 0 ;
3247 while ( freqIds[target] != *i )
3248 target++ ;
3249 vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
3250 freqid[itable].push_back( *i ) ;
3251 double incr = abs( abcissa[1] - abcissa[0] ) ;
3252 iffreq[itable].push_back( (*min_element(abcissa.begin(),abcissa.end()))-0.5*incr ) ;
3253 iffreq[itable].push_back( (*max_element(abcissa.begin(),abcissa.end()))+0.5*incr ) ;
3254 }
3255 }
3256
3257 // debug
3258// os << "IF settings summary:" << endl ;
3259// for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3260// os << " Table" << i << endl ;
3261// for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3262// os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
3263// }
3264// }
3265// os << endl ;
3266// os.post() ;
3267
3268 // IF grouping based on their frequency coverage
3269 // ifgrp: number of member in each IF group
3270 // ifgrp[numgrp]
3271 // ifgrp: [n0, n1, ...]
3272 //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3273
3274 // parameter for IF grouping
3275 // groupmode = OR retrieve all region
3276 // AND only retrieve overlaped region
3277 //string groupmode = "AND" ;
3278 string groupmode = "OR" ;
3279 uInt sizecr = 0 ;
3280 if ( groupmode == "AND" )
3281 sizecr = 1 ;
3282 else if ( groupmode == "OR" )
3283 sizecr = 0 ;
3284
3285 vector<double> sortedfreq ;
3286 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3287 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3288 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3289 sortedfreq.push_back( iffreq[i][j] ) ;
3290 }
3291 }
3292 sort( sortedfreq.begin(), sortedfreq.end() ) ;
3293 vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
3294 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3295 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3296 double range0 = iffreq[itable][2*iif] ;
3297 double range1 = iffreq[itable][2*iif+1] ;
3298 for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
3299 double fmin = max( range0, sortedfreq[j] ) ;
3300 double fmax = min( range1, sortedfreq[j+1] ) ;
3301 if ( fmin < fmax ) {
3302 ifgrp[j]++ ;
3303 }
3304 }
3305 }
3306 }
3307
3308 // Grouping continuous IF groups (without frequency gap)
3309 // freqgrp: list of IF group indexes in each frequency group
3310 // freqgrp[numgrp][nummember]
3311 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3312 // [ifgrp10, ifgrp11, ifgrp12, ...],
3313 // ...
3314 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3315 // grprange[2*numgrp]
3316 // grprange: [fmin0,fmax0,fmin1,fmax1,...]
3317 vector< vector<uInt> > freqgrp ;
3318 vector<double> grprange ;
3319 vector<uInt> grpedge ;
3320 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3321 if ( ifgrp[igrp] <= sizecr ) {
3322 grpedge.push_back( igrp ) ;
3323 }
3324 }
3325 grpedge.push_back( ifgrp.size() ) ;
3326 uInt itmp = 0 ;
3327 for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
3328 int n = grpedge[i] - itmp ;
3329 if ( n > 0 ) {
3330 vector<uInt> members( n ) ;
3331 for ( int j = 0 ; j < n ; j++ ) {
3332 members[j] = itmp+j ;
3333 }
3334 freqgrp.push_back( members ) ;
3335 grprange.push_back( sortedfreq[itmp] ) ;
3336 grprange.push_back( sortedfreq[grpedge[i]] ) ;
3337 }
3338 itmp += n + 1 ;
3339 }
3340
3341 // print frequency group
3342 oss.str("") ;
3343 oss << "Frequency Group summary: " << endl ;
3344 oss << " GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
3345 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3346 oss << " GROUP " << setw( 2 ) << i << ": [" << grprange[2*i] << "," << grprange[2*i+1] << "]" ;
3347 oss << endl ;
3348 }
3349 oss << endl ;
3350 os << oss.str() << LogIO::POST ;
3351
3352 // groups: list of frequency group index whose frequency range overlaps
3353 // with that of each table and IF
3354 // groups[numtable][numIF]
3355 // groups: [[grpx, grpy,...],
3356 // [grpa, grpb,...],
3357 // ...
3358 // [grpk, grpm,...]]
3359 vector< vector<uInt> > groups( insize ) ;
3360 for ( uInt i = 0 ; i < insize ; i++ ) {
3361 groups[i].resize( freqid[i].size() ) ;
3362 }
3363 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3364 for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
3365 double minf = iffreq[itable][2*ifreq] ;
3366 uInt groupid ;
3367 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3368 vector<uInt> memberlist = freqgrp[igrp] ;
3369 if ( (minf >= grprange[2*igrp]) && (minf <= grprange[2*igrp+1]) ) {
3370 groupid = igrp ;
3371 break ;
3372 }
3373 }
3374 groups[itable][ifreq] = groupid ;
3375 }
3376 }
3377
3378
3379 // print membership
3380 oss.str("") ;
3381 for ( uInt i = 0 ; i < insize ; i++ ) {
3382 oss << "Table " << i << endl ;
3383 for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3384 oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3385 oss << setw( 2 ) << groups[i][j] ;
3386 oss << endl ;
3387 }
3388 }
3389 os << oss.str() << LogIO::POST ;
3390
3391 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3392 //os << "All IF number is set to IF group index" << LogIO::POST ;
3393 // reset SCANNO only when avmode != "SCAN"
3394 if ( avmode != "SCAN" ) {
3395 os << "All scan number is set to 0" << LogIO::POST ;
3396 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3397 uInt nrow = newin[itable]->nrow() ;
3398 Vector<uInt> resetScan( nrow, 0 ) ;
3399 newin[itable]->scanCol_.putColumn( resetScan ) ;
3400 }
3401 }
3402
3403 // reset spectra and flagtra: align spectral resolution
3404 //os << "Align spectral resolution" << LogIO::POST ;
3405 // gmaxdnu: the coarsest frequency resolution in the frequency group
3406 // gminfreq: lower frequency edge of the frequency group
3407 // gnchan: number of channels for the frequency group
3408 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3409 vector<double> gminfreq( freqgrp.size() ) ;
3410 vector<double> gnchan( freqgrp.size() ) ;
3411 for ( uInt i = 0 ; i < insize ; i++ ) {
3412 vector<uInt> members = groups[i] ;
3413 for ( uInt j = 0 ; j < members.size() ; j++ ) {
3414 uInt groupid = members[j] ;
3415 Double rp,rv,ic ;
3416 newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
3417 if ( abs(ic) > abs(gmaxdnu[groupid]) )
3418 gmaxdnu[groupid] = ic ;
3419 }
3420 }
3421 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3422 gminfreq[igrp] = grprange[2*igrp] ;
3423 double maxfreq = grprange[2*igrp+1] ;
3424 gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
3425 }
3426
3427 // regrid spectral data and update frequency info
3428 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3429 Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
3430 Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
3431
3432 // update MAIN
3433 for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
3434 uInt groupid = groups[itable][oldFreqId[irow]] ;
3435 newFreqId[irow] = groupid ;
3436 newin[itable]->regridChannel( gnchan[groupid],
3437 gmaxdnu[groupid],
3438 gminfreq[groupid],
3439 irow ) ;
3440 }
3441 newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
3442 newin[itable]->ifCol_.putColumn( newFreqId ) ;
3443
3444 // update FREQUENCIES
3445 Table tab = newin[itable]->frequencies().table() ;
3446 ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
3447 ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
3448 ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
3449 ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
3450 if ( freqgrp.size() > tab.nrow() ) {
3451 tab.addRow( freqgrp.size()-tab.nrow() ) ;
3452 }
3453 for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
3454 Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
3455 Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
3456 Double increment = gmaxdnu[irow] ;
3457 fIdCol.put( irow, irow ) ;
3458 fRefPixCol.put( irow, refpix ) ;
3459 fRefValCol.put( irow, refval ) ;
3460 fIncrCol.put( irow, increment ) ;
3461 }
3462 }
3463
3464 // set back coordinfo
3465 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3466 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3467 coordinfo[0] = oldinfo[itable] ;
3468 newin[itable]->setCoordInfo( coordinfo ) ;
3469 }
3470
3471 // average
3472 out = average( newin, mask, weight, avmode ) ;
3473 }
3474 else {
3475 // simple average
3476 out = average( in, mask, weight, avmode ) ;
3477 }
3478
3479 return out;
3480}
3481
3482CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3483 const String calmode,
3484 const String antname )
3485{
3486 // frequency switch
3487 if ( calmode == "fs" ) {
3488 return cwcalfs( s, antname ) ;
3489 }
3490 else {
3491 vector<bool> masks = s->getMask( 0 ) ;
3492 vector<int> types ;
3493
3494 // save original table selection
3495 Table torg = s->table_ ;
3496
3497 // sky scan
3498 bool insitu = insitu_ ;
3499 insitu_ = false ;
3500 // share calibration scans before average with out
3501 CountedPtr<Scantable> out = getScantable( s, true ) ;
3502 insitu_ = insitu ;
3503 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3504 out->attach() ;
3505 CountedPtr<Scantable> asky = averageWithinSession( out,
3506 masks,
3507 "TINT" ) ;
3508 // hot scan
3509 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3510 out->attach() ;
3511 CountedPtr<Scantable> ahot = averageWithinSession( out,
3512 masks,
3513 "TINT" ) ;
3514 // cold scan
3515 CountedPtr<Scantable> acold ;
3516// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3517// out->attach() ;
3518// CountedPtr<Scantable> acold = averageWithinSession( out,
3519// masks,
3520// "TINT" ) ;
3521
3522 // off scan
3523 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3524 out->attach() ;
3525 CountedPtr<Scantable> aoff = averageWithinSession( out,
3526 masks,
3527 "TINT" ) ;
3528
3529 // on scan
3530 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3531 s->attach() ;
3532 out->table_ = out->originalTable_ ;
3533 out->attach() ;
3534 out->table().addRow( s->nrow() ) ;
3535 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
3536
3537 // iterate throgh STIdxIter2
3538 ChopperWheelCalibrator cal(out, s, aoff, asky, ahot, acold);
3539 STIdxIter2::Iterate<ChopperWheelCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3540
3541 s->table_ = torg ;
3542 s->attach() ;
3543
3544 // flux unit
3545 out->setFluxUnit( "K" ) ;
3546
3547 return out ;
3548 }
3549}
3550
3551CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3552 const String calmode )
3553{
3554 // frequency switch
3555 if ( calmode == "fs" ) {
3556 return almacalfs( s ) ;
3557 }
3558 else {
3559// double t0, t1 ;
3560// t0 = mathutil::gettimeofday_sec() ;
3561 vector<bool> masks = s->getMask( 0 ) ;
3562
3563 // save original table selection
3564 Table torg = s->table_ ;
3565
3566 // off scan
3567 // TODO 2010/01/08 TN
3568 // Grouping by time should be needed before averaging.
3569 // Each group must have own unique SCANNO (should be renumbered).
3570 // See PIPELINE/SDCalibration.py
3571 bool insitu = insitu_ ;
3572 insitu_ = false ;
3573 // share off scan before average with out
3574 CountedPtr<Scantable> out = getScantable( s, true ) ;
3575 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3576 out->attach() ;
3577 insitu_ = insitu ;
3578 CountedPtr<Scantable> aoff = averageWithinSession( out,
3579 masks,
3580 "TINT" ) ;
3581
3582 // on scan
3583// t0 = mathutil::gettimeofday_sec() ;
3584 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3585 s->attach() ;
3586 out->table_ = out->originalTable_ ;
3587 out->attach() ;
3588 out->table().addRow( s->nrow() ) ;
3589 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
3590// t1 = mathutil::gettimeofday_sec() ;
3591// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
3592
3593 // process each on scan
3594// t0 = mathutil::gettimeofday_sec() ;
3595
3596 // iterate throgh STIdxIter2
3597 AlmaCalibrator cal(out, s, aoff);
3598 STIdxIter2::Iterate<AlmaCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3599
3600 // finalize
3601 s->table_ = torg ;
3602 s->attach() ;
3603
3604// t1 = mathutil::gettimeofday_sec() ;
3605// cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
3606
3607 // flux unit
3608 out->setFluxUnit( "K" ) ;
3609
3610 return out ;
3611 }
3612}
3613
3614CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3615 const String antname )
3616{
3617 vector<int> types ;
3618
3619 // APEX calibration mode
3620 int apexcalmode = 1 ;
3621
3622 if ( antname.find( "APEX" ) != string::npos ) {
3623 // check if off scan exists or not
3624 STSelector sel = STSelector() ;
3625 //sel.setName( offstr1 ) ;
3626 types.push_back( SrcType::FLOOFF ) ;
3627 sel.setTypes( types ) ;
3628 try {
3629 s->setSelection( sel ) ;
3630 }
3631 catch ( AipsError &e ) {
3632 apexcalmode = 0 ;
3633 }
3634 sel.reset() ;
3635 }
3636 s->unsetSelection() ;
3637 types.clear() ;
3638
3639 vector<bool> masks = s->getMask( 0 ) ;
3640 CountedPtr<Scantable> ssig, sref ;
3641 //CountedPtr<Scantable> out ;
3642 bool insitu = insitu_ ;
3643 insitu_ = False ;
3644 CountedPtr<Scantable> out = getScantable( s, true ) ;
3645 insitu_ = insitu ;
3646
3647 if ( antname.find( "APEX" ) != string::npos ) {
3648 // APEX calibration
3649 // sky scan
3650 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
3651 out->attach() ;
3652 CountedPtr<Scantable> askylo = averageWithinSession( out,
3653 masks,
3654 "TINT" ) ;
3655 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
3656 out->attach() ;
3657 CountedPtr<Scantable> askyhi = averageWithinSession( out,
3658 masks,
3659 "TINT" ) ;
3660
3661 // hot scan
3662 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
3663 out->attach() ;
3664 CountedPtr<Scantable> ahotlo = averageWithinSession( out,
3665 masks,
3666 "TINT" ) ;
3667 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
3668 out->attach() ;
3669 CountedPtr<Scantable> ahothi = averageWithinSession( out,
3670 masks,
3671 "TINT" ) ;
3672
3673 // cold scan
3674 CountedPtr<Scantable> acoldlo, acoldhi ;
3675// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
3676// out->attach() ;
3677// CountedPtr<Scantable> acoldlo = averageWithinSession( out,
3678// masks,
3679// "TINT" ) ;
3680// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
3681// out->attach() ;
3682// CountedPtr<Scantable> acoldhi = averageWithinSession( out,
3683// masks,
3684// "TINT" ) ;
3685
3686 // ref scan
3687 insitu_ = false ;
3688 sref = getScantable( s, true ) ;
3689 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3690 insitu_ = insitu ;
3691 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
3692 rref->attach() ;
3693 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3694
3695 // sig scan
3696 insitu_ = false ;
3697 ssig = getScantable( s, true ) ;
3698 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3699 insitu_ = insitu ;
3700 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
3701 rsig->attach() ;
3702 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3703
3704 if ( apexcalmode == 0 ) {
3705 // using STIdxIter2
3706 vector<string> cols( 3 ) ;
3707 cols[0] = "BEAMNO" ;
3708 cols[1] = "POLNO" ;
3709 cols[2] = "IFNO" ;
3710 STIdxIter2 iter(ssig, cols) ;
3711 STSelector sel ;
3712 vector< CountedPtr<Scantable> > on( 2 ) ;
3713 on[0] = rsig ;
3714 on[1] = rref ;
3715 vector< CountedPtr<Scantable> > sky( 2 ) ;
3716 sky[0] = askylo ;
3717 sky[1] = askyhi ;
3718 vector< CountedPtr<Scantable> > hot( 2 ) ;
3719 hot[0] = ahotlo ;
3720 hot[1] = ahothi ;
3721 vector< CountedPtr<Scantable> > cold( 2 ) ;
3722 while ( !iter.pastEnd() ) {
3723 Record ids = iter.currentValue() ;
3724 stringstream ss ;
3725 ss << "SELECT FROM $1 WHERE "
3726 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3727 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3728 << "IFNO==" << ids.asuInt("IFNO") ;
3729 //cout << "TaQL string: " << ss.str() << endl ;
3730 sel.setTaQL( ss.str() ) ;
3731 sky[0]->setSelection( sel ) ;
3732 sky[1]->setSelection( sel ) ;
3733 hot[0]->setSelection( sel ) ;
3734 hot[1]->setSelection( sel ) ;
3735 Vector<uInt> rows = iter.getRows( SHARE ) ;
3736 calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
3737 sky[0]->unsetSelection() ;
3738 sky[1]->unsetSelection() ;
3739 hot[0]->unsetSelection() ;
3740 hot[1]->unsetSelection() ;
3741 sel.reset() ;
3742 iter.next() ;
3743 }
3744
3745 }
3746 else if ( apexcalmode == 1 ) {
3747 // APEX fs data with off scan
3748 // off scan
3749 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
3750 out->attach() ;
3751 CountedPtr<Scantable> aofflo = averageWithinSession( out,
3752 masks,
3753 "TINT" ) ;
3754 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
3755 out->attach() ;
3756 CountedPtr<Scantable> aoffhi = averageWithinSession( out,
3757 masks,
3758 "TINT" ) ;
3759
3760 // process each sig and ref scan
3761 // iterate throgh STIdxIter2
3762 ChopperWheelCalibrator cal_sig(ssig, rsig, aofflo, askylo, ahotlo, acoldlo);
3763 STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_sig, "BEAMNO,POLNO,IFNO");
3764 ChopperWheelCalibrator cal_ref(sref, rref, aoffhi, askyhi, ahothi, acoldhi);
3765 STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_ref, "BEAMNO,POLNO,IFNO");
3766 }
3767 }
3768 else {
3769 // non-APEX fs data
3770 // sky scan
3771 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3772 out->attach() ;
3773 CountedPtr<Scantable> asky = averageWithinSession( out,
3774 masks,
3775 "TINT" ) ;
3776 STSelector sel = STSelector() ;
3777
3778 // hot scan
3779 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3780 out->attach() ;
3781 CountedPtr<Scantable> ahot = averageWithinSession( out,
3782 masks,
3783 "TINT" ) ;
3784
3785 // cold scan
3786 CountedPtr<Scantable> acold ;
3787// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3788// out->attach() ;
3789// CountedPtr<Scantable> acold = averageWithinSession( out,
3790// masks,
3791// "TINT" ) ;
3792
3793 // ref scan
3794 bool insitu = insitu_ ;
3795 insitu_ = false ;
3796 sref = getScantable( s, true ) ;
3797 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3798 insitu_ = insitu ;
3799 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3800 rref->attach() ;
3801 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3802
3803 // sig scan
3804 insitu_ = false ;
3805 ssig = getScantable( s, true ) ;
3806 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3807 insitu_ = insitu ;
3808 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3809 rsig->attach() ;
3810 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3811
3812 // process each sig and ref scan
3813 vector<string> cols( 3 ) ;
3814 cols[0] = "BEAMNO" ;
3815 cols[1] = "POLNO" ;
3816 cols[2] = "IFNO" ;
3817 STIdxIter2 iter(ssig, cols);
3818 while ( !iter.pastEnd() ) {
3819 Record ids = iter.currentValue() ;
3820 stringstream ss ;
3821 ss << "SELECT FROM $1 WHERE "
3822 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3823 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3824 << "IFNO==" << ids.asuInt("IFNO") ;
3825 //cout << "TaQL string: " << ss.str() << endl ;
3826 sel.setTaQL( ss.str() ) ;
3827 ahot->setSelection( sel ) ;
3828 asky->setSelection( sel ) ;
3829 Vector<uInt> rows = iter.getRows( SHARE ) ;
3830 // out should be an exact copy of s except that SPECTRA column is empty
3831 calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
3832 ahot->unsetSelection() ;
3833 asky->unsetSelection() ;
3834 sel.reset() ;
3835 iter.next() ;
3836 }
3837 }
3838
3839 // do folding if necessary
3840 Table sigtab = ssig->table() ;
3841 Table reftab = sref->table() ;
3842 ScalarColumn<uInt> reffidCol ;
3843 Int nchan = (Int)ssig->nchan() ;
3844 reffidCol.attach( reftab, "FREQ_ID" ) ;
3845 Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
3846 Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
3847 vector<uInt> sfids_unique ;
3848 vector<uInt> rfids_unique ;
3849 vector<uInt> sifno_unique ;
3850 vector<uInt> rifno_unique ;
3851 for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
3852 if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
3853 sfids_unique.push_back( sfids[i] ) ;
3854 sifno_unique.push_back( ssig->getIF( i ) ) ;
3855 }
3856 if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
3857 rfids_unique.push_back( rfids[i] ) ;
3858 rifno_unique.push_back( sref->getIF( i ) ) ;
3859 }
3860 }
3861 double refpix_sig, refval_sig, increment_sig ;
3862 double refpix_ref, refval_ref, increment_ref ;
3863 vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
3864 for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
3865 ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
3866 sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
3867 if ( refpix_sig == refpix_ref ) {
3868 double foffset = refval_ref - refval_sig ;
3869 int choffset = static_cast<int>(foffset/increment_sig) ;
3870 double doffset = foffset / increment_sig ;
3871 if ( abs(choffset) >= nchan ) {
3872 LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
3873 os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
3874 os << "Just return signal data" << LogIO::POST ;
3875 //std::vector< CountedPtr<Scantable> > tabs ;
3876 //tabs.push_back( ssig ) ;
3877 //tabs.push_back( sref ) ;
3878 //out = merge( tabs ) ;
3879 tmp[i] = ssig ;
3880 }
3881 else {
3882 STSelector sel = STSelector() ;
3883 vector<int> v( 1, sifno_unique[i] ) ;
3884 sel.setIFs( v ) ;
3885 ssig->setSelection( sel ) ;
3886 sel.reset() ;
3887 v[0] = rifno_unique[i] ;
3888 sel.setIFs( v ) ;
3889 sref->setSelection( sel ) ;
3890 sel.reset() ;
3891 if ( antname.find( "APEX" ) != string::npos ) {
3892 tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
3893 //tmp[i] = dofold( ssig, sref, doffset ) ;
3894 }
3895 else {
3896 tmp[i] = dofold( ssig, sref, doffset ) ;
3897 }
3898 ssig->unsetSelection() ;
3899 sref->unsetSelection() ;
3900 }
3901 }
3902 }
3903
3904 if ( tmp.size() > 1 ) {
3905 out = merge( tmp ) ;
3906 }
3907 else {
3908 out = tmp[0] ;
3909 }
3910
3911 // flux unit
3912 out->setFluxUnit( "K" ) ;
3913
3914 return out ;
3915}
3916
3917CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
3918{
3919 (void) s; //currently unused
3920 CountedPtr<Scantable> out ;
3921
3922 return out ;
3923}
3924
3925void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
3926 CountedPtr<Scantable> &ref,
3927 const vector< CountedPtr<Scantable> >& on,
3928 const vector< CountedPtr<Scantable> >& sky,
3929 const vector< CountedPtr<Scantable> >& hot,
3930 const vector< CountedPtr<Scantable> >& cold,
3931 const Vector<uInt> &rows )
3932{
3933 // if rows is empty, just return
3934 if ( rows.nelements() == 0 )
3935 return ;
3936 ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
3937 Vector<Double> timeSkyS = timeCol.getColumn() ;
3938 timeCol.attach( sky[1]->table(), "TIME" ) ;
3939 Vector<Double> timeSkyR = timeCol.getColumn() ;
3940 timeCol.attach( hot[0]->table(), "TIME" ) ;
3941 Vector<Double> timeHotS = timeCol.getColumn() ;
3942 timeCol.attach( hot[1]->table(), "TIME" ) ;
3943 Vector<Double> timeHotR = timeCol.getColumn() ;
3944 timeCol.attach( sig->table(), "TIME" ) ;
3945 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
3946 ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
3947 SpectralData skyspectraS(arrayFloatCol.getColumn());
3948 arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
3949 SpectralData skyspectraR(arrayFloatCol.getColumn());
3950 arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
3951 SpectralData hotspectraS(arrayFloatCol.getColumn());
3952 arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
3953 SpectralData hotspectraR(arrayFloatCol.getColumn());
3954 TcalData tcaldataS(sky[0]);
3955 TsysData tsysdataS(sky[0]);
3956 TcalData tcaldataR(sky[1]);
3957 TsysData tsysdataR(sky[1]);
3958 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
3959 Vector<Float> spec( spsize ) ;
3960 // I know that the data is contiguous
3961 const uInt *p = rows.data() ;
3962 vector<int> ids( 2 ) ;
3963 Block<uInt> flagchan( spsize ) ;
3964 uInt nflag = 0 ;
3965 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
3966 double reftime = timeCol.asdouble(*p) ;
3967 ids = getRowIdFromTime( reftime, timeSkyS ) ;
3968 Vector<Float> spskyS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyS, ids, skyspectraS, "linear");
3969 Vector<Float> tcalS = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyS, ids, tcaldataS, "linear");
3970 Vector<Float> tsysS = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyS, ids, tsysdataS, "linear");
3971 ids = getRowIdFromTime( reftime, timeHotS ) ;
3972 Vector<Float> sphotS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotS, ids, hotspectraS, "linear");
3973 reftime = timeCol2.asdouble(*p) ;
3974 ids = getRowIdFromTime( reftime, timeSkyR ) ;
3975 Vector<Float> spskyR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyR, ids, skyspectraR, "linear");
3976 Vector<Float> tcalR = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyR, ids, tcaldataR, "linear");
3977 Vector<Float> tsysR = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyR, ids, tsysdataR, "linear");
3978 ids = getRowIdFromTime( reftime, timeHotR ) ;
3979 Vector<Float> sphotR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotR, ids, hotspectraR, "linear");
3980 Vector<Float> spsig = on[0]->specCol_( *p ) ;
3981 Vector<Float> spref = on[1]->specCol_( *p ) ;
3982 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
3983 if ( (sphotS[j]-spskyS[j]) == 0.0 || (sphotR[j]-spskyR[j]) == 0.0 ) {
3984 spec[j] = 0.0 ;
3985 flagchan[nflag++] = j ;
3986 }
3987 else {
3988 spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
3989 - tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
3990 }
3991 }
3992 sig->specCol_.put( *p, spec ) ;
3993 sig->tsysCol_.put( *p, tsysS ) ;
3994 spec *= (Float)-1.0 ;
3995 ref->specCol_.put( *p, spec ) ;
3996 ref->tsysCol_.put( *p, tsysR ) ;
3997 if ( nflag > 0 ) {
3998 Vector<uChar> flsig = sig->flagsCol_( *p ) ;
3999 Vector<uChar> flref = ref->flagsCol_( *p ) ;
4000 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4001 flsig[flagchan[j]] = (uChar)True ;
4002 flref[flagchan[j]] = (uChar)True ;
4003 }
4004 sig->flagsCol_.put( *p, flsig ) ;
4005 ref->flagsCol_.put( *p, flref ) ;
4006 }
4007 nflag = 0 ;
4008 p++ ;
4009 }
4010}
4011
4012void STMath::calibrateFS( CountedPtr<Scantable> &sig,
4013 CountedPtr<Scantable> &ref,
4014 const CountedPtr<Scantable>& rsig,
4015 const CountedPtr<Scantable>& rref,
4016 const CountedPtr<Scantable>& sky,
4017 const CountedPtr<Scantable>& hot,
4018 const CountedPtr<Scantable>& cold,
4019 const Vector<uInt> &rows )
4020{
4021 // if rows is empty, just return
4022 if ( rows.nelements() == 0 )
4023 return ;
4024 ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
4025 Vector<Double> timeSky = timeCol.getColumn() ;
4026 timeCol.attach( hot->table(), "TIME" ) ;
4027 Vector<Double> timeHot = timeCol.getColumn() ;
4028 timeCol.attach( sig->table(), "TIME" ) ;
4029 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4030 ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
4031 SpectralData skyspectra(arrayFloatCol.getColumn());
4032 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4033 SpectralData hotspectra(arrayFloatCol.getColumn());
4034 TcalData tcaldata(sky);
4035 TsysData tsysdata(sky);
4036 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4037 Vector<Float> spec( spsize ) ;
4038 // I know that the data is contiguous
4039 const uInt *p = rows.data() ;
4040 vector<int> ids( 2 ) ;
4041 Block<uInt> flagchan( spsize ) ;
4042 uInt nflag = 0 ;
4043 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4044 double reftime = timeCol.asdouble(*p) ;
4045 ids = getRowIdFromTime( reftime, timeSky ) ;
4046 Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4047 Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4048 Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4049 ids = getRowIdFromTime( reftime, timeHot ) ;
4050 Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4051 Vector<Float> spsig = rsig->specCol_( *p ) ;
4052 Vector<Float> spref = rref->specCol_( *p ) ;
4053 // using gain array
4054 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4055 if ( spref[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4056 spec[j] = 0.0 ;
4057 flagchan[nflag++] = j ;
4058 }
4059 else {
4060 spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
4061 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4062 }
4063 }
4064 sig->specCol_.put( *p, spec ) ;
4065 sig->tsysCol_.put( *p, tsys ) ;
4066 if ( nflag > 0 ) {
4067 Vector<uChar> fl = sig->flagsCol_( *p ) ;
4068 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4069 fl[flagchan[j]] = (uChar)True ;
4070 }
4071 sig->flagsCol_.put( *p, fl ) ;
4072 }
4073 nflag = 0 ;
4074
4075 reftime = timeCol2.asdouble(*p) ;
4076 ids = getRowIdFromTime( reftime, timeSky ) ;
4077 spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4078 tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4079 tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4080 ids = getRowIdFromTime( reftime, timeHot ) ;
4081 sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4082 // using gain array
4083 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4084 if ( spsig[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4085 spec[j] = 0.0 ;
4086 flagchan[nflag++] = j ;
4087 }
4088 else {
4089 spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
4090 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4091 }
4092 }
4093 ref->specCol_.put( *p, spec ) ;
4094 ref->tsysCol_.put( *p, tsys ) ;
4095 if ( nflag > 0 ) {
4096 Vector<uChar> fl = ref->flagsCol_( *p ) ;
4097 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4098 fl[flagchan[j]] = (uChar)True ;
4099 }
4100 ref->flagsCol_.put( *p, fl ) ;
4101 }
4102 nflag = 0 ;
4103 p++ ;
4104 }
4105}
4106
4107void STMath::copyRows( Table &out,
4108 const Table &in,
4109 uInt startout,
4110 uInt startin,
4111 uInt nrow,
4112 Bool copySpectra,
4113 Bool copyFlagtra,
4114 Bool copyTsys )
4115{
4116 uInt nexclude = 0 ;
4117 Block<String> excludeColsBlock( 3 ) ;
4118 if ( !copySpectra ) {
4119 excludeColsBlock[nexclude] = "SPECTRA" ;
4120 nexclude++ ;
4121 }
4122 if ( !copyFlagtra ) {
4123 excludeColsBlock[nexclude] = "FLAGTRA" ;
4124 nexclude++ ;
4125 }
4126 if ( !copyTsys ) {
4127 excludeColsBlock[nexclude] = "TSYS" ;
4128 nexclude++ ;
4129 }
4130 // if ( nexclude < 3 ) {
4131 // excludeCols.resize( nexclude, True ) ;
4132 // }
4133 Vector<String> excludeCols( IPosition(1,nexclude),
4134 excludeColsBlock.storage(),
4135 SHARE ) ;
4136// cout << "excludeCols=" << excludeCols << endl ;
4137 TableRow rowout( out, excludeCols, True ) ;
4138 ROTableRow rowin( in, excludeCols, True ) ;
4139 uInt rin = startin ;
4140 uInt rout = startout ;
4141 for ( uInt i = 0 ; i < nrow ; i++ ) {
4142 rowin.get( rin ) ;
4143 rowout.putMatchingFields( rout, rowin.record() ) ;
4144 rin++ ;
4145 rout++ ;
4146 }
4147}
4148
4149CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
4150 vector<bool> &mask,
4151 string weight )
4152{
4153 // prepare output table
4154 bool insitu = insitu_ ;
4155 insitu_ = false ;
4156 CountedPtr<Scantable> a = getScantable( s, true ) ;
4157 insitu_ = insitu ;
4158 Table &atab = a->table() ;
4159 ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
4160
4161 if ( s->nrow() == 0 )
4162 return a ;
4163
4164 // setup RowAccumulator
4165 WeightType wtype = stringToWeight( weight ) ;
4166 RowAccumulator acc( wtype ) ;
4167 Vector<Bool> cmask( mask ) ;
4168 acc.setUserMask( cmask ) ;
4169
4170 vector<string> cols( 3 ) ;
4171 cols[0] = "IFNO" ;
4172 cols[1] = "POLNO" ;
4173 cols[2] = "BEAMNO" ;
4174 STIdxIter2 iter( s, cols ) ;
4175
4176 Table ttab = s->table() ;
4177 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
4178 Vector<Double> timeVec = timeCol->getColumn() ;
4179 delete timeCol ;
4180 Vector<Double> interval = s->integrCol_.getColumn() ;
4181 uInt nrow = timeVec.nelements() ;
4182 uInt outrow = 0 ;
4183
4184 while( !iter.pastEnd() ) {
4185
4186 Vector<uInt> rows = iter.getRows( SHARE ) ;
4187 uInt len = rows.nelements() ;
4188
4189 if ( len == 0 ) {
4190 iter.next() ;
4191 continue ;
4192 }
4193
4194 uInt nchan = s->nchan(s->getIF(rows[0])) ;
4195 Vector<uChar> flag( nchan ) ;
4196 Vector<Bool> bflag( nchan ) ;
4197 Vector<Float> spec( nchan ) ;
4198 Vector<Float> tsys( nchan ) ;
4199
4200 Vector<Double> timeSep( len-1 ) ;
4201 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4202 timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
4203 }
4204
4205 uInt irow ;
4206 uInt jrow ;
4207 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4208 irow = rows[i] ;
4209 jrow = rows[i+1] ;
4210 // accumulate data
4211 s->flagsCol_.get( irow, flag ) ;
4212 convertArray( bflag, flag ) ;
4213 s->specCol_.get( irow, spec ) ;
4214 tsys.assign( s->tsysCol_( irow ) ) ;
4215 if ( !allEQ(bflag,True) )
4216 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4217 double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
4218 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
4219 if ( gap > 1.1 ) {
4220 //cout << "detected gap between " << i << " and " << i+1 << endl ;
4221 // put data to output table
4222 // reset RowAccumulator
4223 if ( acc.state() ) {
4224 atab.addRow() ;
4225 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4226 acc.replaceNaN() ;
4227 const Vector<Bool> &msk = acc.getMask() ;
4228 convertArray( flag, !msk ) ;
4229 for (uInt k = 0; k < nchan; ++k) {
4230 uChar userFlag = 1 << 7;
4231 if (msk[k]==True) userFlag = 0 << 7;
4232 flag(k) = userFlag;
4233 }
4234 a->flagsCol_.put( outrow, flag ) ;
4235 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4236 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4237 a->integrCol_.put( outrow, acc.getInterval() ) ;
4238 timeColOut.put( outrow, acc.getTime() ) ;
4239 a->cycleCol_.put( outrow, 0 ) ;
4240 }
4241 acc.reset() ;
4242 outrow++ ;
4243 }
4244 }
4245
4246 // accumulate and add last data
4247 irow = rows[len-1] ;
4248 s->flagsCol_.get( irow, flag ) ;
4249 convertArray( bflag, flag ) ;
4250 s->specCol_.get( irow, spec ) ;
4251 tsys.assign( s->tsysCol_( irow ) ) ;
4252 if (!allEQ(bflag,True) )
4253 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4254 if ( acc.state() ) {
4255 atab.addRow() ;
4256 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4257 acc.replaceNaN() ;
4258 const Vector<Bool> &msk = acc.getMask() ;
4259 convertArray( flag, !msk ) ;
4260 for (uInt k = 0; k < nchan; ++k) {
4261 uChar userFlag = 1 << 7;
4262 if (msk[k]==True) userFlag = 0 << 7;
4263 flag(k) = userFlag;
4264 }
4265 a->flagsCol_.put( outrow, flag ) ;
4266 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4267 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4268 a->integrCol_.put( outrow, acc.getInterval() ) ;
4269 timeColOut.put( outrow, acc.getTime() ) ;
4270 a->cycleCol_.put( outrow, 0 ) ;
4271 }
4272 acc.reset() ;
4273 outrow++ ;
4274
4275 iter.next() ;
4276 }
4277
4278 return a ;
4279}
Note: See TracBrowser for help on using the repository browser.