source: trunk/src/STMath.cpp@ 3045

Last change on this file since 3045 was 3010, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6599

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: modified so that sd.scantable.stats() returns None for flagged rows and rows with all channels flagged.


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