source: trunk/src/STMath.cpp@ 2910

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

New Development: No

JIRA Issue: Yes CAS-5875

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdcoadd

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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