source: trunk/src/STMath.cpp@ 2917

Last change on this file since 2917 was 2917, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: test_sdcal, test_sdcal2

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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