source: trunk/src/Scantable.cpp@ 2016

Last change on this file since 2016 was 2012, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2373, CAS-2620)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: For Scantable::polyBaseline(), parameters and return value have been changed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline, sd.linefinder

Description: (1) CAS-2373-related:

(1.1) Modified Scantable::polyBaseline() to have the row-based loop inside.

Now it fits and subtracts baseline for all rows and also output info
about the fitting result to logger and text file, while in the
previous version this method just did baseline fitting/subtraction
for one row only and had to be put inside a row-based loop at the
python side ("poly_baseline()" in scantable.py) and result output had
also to be controlled at the python side. Using a test data from NRO
45m telescope (348,000 rows, 512 channels), the processing time of
scantable.poly_baseline() has reduced from 130 minutes to 5-6 minutes.

(1.2) For accelerating the another polynomial fitting function, namely

scantable.auto_poly_baseline(), added a method
Scantable::autoPolyBaseline(). This basically does the same thing
with Scantable::polyBaseline(), but uses linefinfer also to
automatically flag the line regions.

(1.3) To make linefinder usable in Scantable class, added a method

linefinder.setdata(). This makes it possible to apply linefinder
for a float-list data given as a parameter, without setting scantable,
while it was indispensable to set scantable to use linefinger previously.

(2) CAS-2620-related:

Added Scantable::cubicSplineBaseline() and autoCubicSplineBaseline() for
fit baseline using the cubic spline function. Parameters include npiece
(number of polynomial pieces), clipthresh (clipping threshold), and
clipniter (maximum iteration number).



  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.6 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13#include <fstream>
14
15#include <casa/aips.h>
16#include <casa/iostream.h>
17#include <casa/iomanip.h>
18#include <casa/OS/Path.h>
19#include <casa/OS/File.h>
20#include <casa/Arrays/Array.h>
21#include <casa/Arrays/ArrayMath.h>
22#include <casa/Arrays/MaskArrMath.h>
23#include <casa/Arrays/ArrayLogical.h>
24#include <casa/Arrays/ArrayAccessor.h>
25#include <casa/Arrays/Vector.h>
26#include <casa/Arrays/VectorSTLIterator.h>
27#include <casa/Arrays/Slice.h>
28#include <casa/BasicMath/Math.h>
29#include <casa/BasicSL/Constants.h>
30#include <casa/Quanta/MVAngle.h>
31#include <casa/Containers/RecordField.h>
32#include <casa/Utilities/GenSort.h>
33#include <casa/Logging/LogIO.h>
34
35#include <tables/Tables/TableParse.h>
36#include <tables/Tables/TableDesc.h>
37#include <tables/Tables/TableCopy.h>
38#include <tables/Tables/SetupNewTab.h>
39#include <tables/Tables/ScaColDesc.h>
40#include <tables/Tables/ArrColDesc.h>
41#include <tables/Tables/TableRow.h>
42#include <tables/Tables/TableVector.h>
43#include <tables/Tables/TableIter.h>
44
45#include <tables/Tables/ExprNode.h>
46#include <tables/Tables/TableRecord.h>
47#include <casa/Quanta/MVTime.h>
48#include <casa/Quanta/MVAngle.h>
49#include <measures/Measures/MeasRef.h>
50#include <measures/Measures/MeasTable.h>
51// needed to avoid error in .tcc
52#include <measures/Measures/MCDirection.h>
53//
54#include <measures/Measures/MDirection.h>
55#include <measures/Measures/MFrequency.h>
56#include <measures/Measures/MEpoch.h>
57#include <measures/TableMeasures/TableMeasRefDesc.h>
58#include <measures/TableMeasures/TableMeasValueDesc.h>
59#include <measures/TableMeasures/TableMeasDesc.h>
60#include <measures/TableMeasures/ScalarMeasColumn.h>
61#include <coordinates/Coordinates/CoordinateUtil.h>
62
63#include <atnf/PKSIO/SrcType.h>
64#include "Scantable.h"
65#include "STPolLinear.h"
66#include "STPolCircular.h"
67#include "STPolStokes.h"
68#include "STAttr.h"
69#include "STLineFinder.h"
70#include "MathUtils.h"
71
72using namespace casa;
73
74namespace asap {
75
76std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
77
78void Scantable::initFactories() {
79 if ( factories_.empty() ) {
80 Scantable::factories_["linear"] = &STPolLinear::myFactory;
81 Scantable::factories_["circular"] = &STPolCircular::myFactory;
82 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
83 }
84}
85
86Scantable::Scantable(Table::TableType ttype) :
87 type_(ttype)
88{
89 initFactories();
90 setupMainTable();
91 freqTable_ = STFrequencies(*this);
92 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
93 weatherTable_ = STWeather(*this);
94 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
95 focusTable_ = STFocus(*this);
96 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
97 tcalTable_ = STTcal(*this);
98 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
99 moleculeTable_ = STMolecules(*this);
100 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
101 historyTable_ = STHistory(*this);
102 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
103 fitTable_ = STFit(*this);
104 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
105 table_.tableInfo().setType( "Scantable" ) ;
106 originalTable_ = table_;
107 attach();
108}
109
110Scantable::Scantable(const std::string& name, Table::TableType ttype) :
111 type_(ttype)
112{
113 initFactories();
114
115 Table tab(name, Table::Update);
116 uInt version = tab.keywordSet().asuInt("VERSION");
117 if (version != version_) {
118 throw(AipsError("Unsupported version of ASAP file."));
119 }
120 if ( type_ == Table::Memory ) {
121 table_ = tab.copyToMemoryTable(generateName());
122 } else {
123 table_ = tab;
124 }
125 table_.tableInfo().setType( "Scantable" ) ;
126
127 attachSubtables();
128 originalTable_ = table_;
129 attach();
130}
131/*
132Scantable::Scantable(const std::string& name, Table::TableType ttype) :
133 type_(ttype)
134{
135 initFactories();
136 Table tab(name, Table::Update);
137 uInt version = tab.keywordSet().asuInt("VERSION");
138 if (version != version_) {
139 throw(AipsError("Unsupported version of ASAP file."));
140 }
141 if ( type_ == Table::Memory ) {
142 table_ = tab.copyToMemoryTable(generateName());
143 } else {
144 table_ = tab;
145 }
146
147 attachSubtables();
148 originalTable_ = table_;
149 attach();
150}
151*/
152
153Scantable::Scantable( const Scantable& other, bool clear )
154{
155 // with or without data
156 String newname = String(generateName());
157 type_ = other.table_.tableType();
158 if ( other.table_.tableType() == Table::Memory ) {
159 if ( clear ) {
160 table_ = TableCopy::makeEmptyMemoryTable(newname,
161 other.table_, True);
162 } else
163 table_ = other.table_.copyToMemoryTable(newname);
164 } else {
165 other.table_.deepCopy(newname, Table::New, False,
166 other.table_.endianFormat(),
167 Bool(clear));
168 table_ = Table(newname, Table::Update);
169 table_.markForDelete();
170 }
171 table_.tableInfo().setType( "Scantable" ) ;
172 /// @todo reindex SCANNO, recompute nbeam, nif, npol
173 if ( clear ) copySubtables(other);
174 attachSubtables();
175 originalTable_ = table_;
176 attach();
177}
178
179void Scantable::copySubtables(const Scantable& other) {
180 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
181 TableCopy::copyRows(t, other.freqTable_.table());
182 t = table_.rwKeywordSet().asTable("FOCUS");
183 TableCopy::copyRows(t, other.focusTable_.table());
184 t = table_.rwKeywordSet().asTable("WEATHER");
185 TableCopy::copyRows(t, other.weatherTable_.table());
186 t = table_.rwKeywordSet().asTable("TCAL");
187 TableCopy::copyRows(t, other.tcalTable_.table());
188 t = table_.rwKeywordSet().asTable("MOLECULES");
189 TableCopy::copyRows(t, other.moleculeTable_.table());
190 t = table_.rwKeywordSet().asTable("HISTORY");
191 TableCopy::copyRows(t, other.historyTable_.table());
192 t = table_.rwKeywordSet().asTable("FIT");
193 TableCopy::copyRows(t, other.fitTable_.table());
194}
195
196void Scantable::attachSubtables()
197{
198 freqTable_ = STFrequencies(table_);
199 focusTable_ = STFocus(table_);
200 weatherTable_ = STWeather(table_);
201 tcalTable_ = STTcal(table_);
202 moleculeTable_ = STMolecules(table_);
203 historyTable_ = STHistory(table_);
204 fitTable_ = STFit(table_);
205}
206
207Scantable::~Scantable()
208{
209 //cout << "~Scantable() " << this << endl;
210}
211
212void Scantable::setupMainTable()
213{
214 TableDesc td("", "1", TableDesc::Scratch);
215 td.comment() = "An ASAP Scantable";
216 td.rwKeywordSet().define("VERSION", uInt(version_));
217
218 // n Cycles
219 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
220 // new index every nBeam x nIF x nPol
221 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
222
223 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
224 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
225 // linear, circular, stokes
226 td.rwKeywordSet().define("POLTYPE", String("linear"));
227 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
228
229 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
230 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
231
232 ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
233 refbeamnoColumn.setDefault(Int(-1));
234 td.addColumn(refbeamnoColumn);
235
236 ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
237 flagrowColumn.setDefault(uInt(0));
238 td.addColumn(flagrowColumn);
239
240 td.addColumn(ScalarColumnDesc<Double>("TIME"));
241 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
242 TableMeasValueDesc measVal(td, "TIME");
243 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
244 mepochCol.write(td);
245
246 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
247
248 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
249 // Type of source (on=0, off=1, other=-1)
250 ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
251 stypeColumn.setDefault(Int(-1));
252 td.addColumn(stypeColumn);
253 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
254
255 //The actual Data Vectors
256 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
257 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
258 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
259
260 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
261 IPosition(1,2),
262 ColumnDesc::Direct));
263 TableMeasRefDesc mdirRef(MDirection::J2000); // default
264 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
265 // the TableMeasDesc gives the column a type
266 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
267 // a uder set table type e.g. GALCTIC, B1950 ...
268 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
269 // writing create the measure column
270 mdirCol.write(td);
271 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
272 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
273 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
274
275 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
276 ScalarColumnDesc<Int> fitColumn("FIT_ID");
277 fitColumn.setDefault(Int(-1));
278 td.addColumn(fitColumn);
279
280 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
281 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
282
283 // columns which just get dragged along, as they aren't used in asap
284 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
285 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
286 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
287 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
288
289 td.rwKeywordSet().define("OBSMODE", String(""));
290
291 // Now create Table SetUp from the description.
292 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
293 table_ = Table(aNewTab, type_, 0);
294 originalTable_ = table_;
295}
296
297void Scantable::attach()
298{
299 timeCol_.attach(table_, "TIME");
300 srcnCol_.attach(table_, "SRCNAME");
301 srctCol_.attach(table_, "SRCTYPE");
302 specCol_.attach(table_, "SPECTRA");
303 flagsCol_.attach(table_, "FLAGTRA");
304 tsysCol_.attach(table_, "TSYS");
305 cycleCol_.attach(table_,"CYCLENO");
306 scanCol_.attach(table_, "SCANNO");
307 beamCol_.attach(table_, "BEAMNO");
308 ifCol_.attach(table_, "IFNO");
309 polCol_.attach(table_, "POLNO");
310 integrCol_.attach(table_, "INTERVAL");
311 azCol_.attach(table_, "AZIMUTH");
312 elCol_.attach(table_, "ELEVATION");
313 dirCol_.attach(table_, "DIRECTION");
314 fldnCol_.attach(table_, "FIELDNAME");
315 rbeamCol_.attach(table_, "REFBEAMNO");
316
317 mweatheridCol_.attach(table_,"WEATHER_ID");
318 mfitidCol_.attach(table_,"FIT_ID");
319 mfreqidCol_.attach(table_, "FREQ_ID");
320 mtcalidCol_.attach(table_, "TCAL_ID");
321 mfocusidCol_.attach(table_, "FOCUS_ID");
322 mmolidCol_.attach(table_, "MOLECULE_ID");
323
324 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
325 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
326
327}
328
329template<class T, class T2>
330void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
331 const String& colName,
332 const T2& defValue)
333{
334 try {
335 col.attach(table_, colName);
336 } catch (TableError& err) {
337 String errMesg = err.getMesg();
338 if (errMesg == "Table column " + colName + " is unknown") {
339 table_.addColumn(ScalarColumnDesc<T>(colName));
340 col.attach(table_, colName);
341 col.fillColumn(static_cast<T>(defValue));
342 } else {
343 throw;
344 }
345 } catch (...) {
346 throw;
347 }
348}
349
350template<class T, class T2>
351void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
352 const String& colName,
353 const Array<T2>& defValue)
354{
355 try {
356 col.attach(table_, colName);
357 } catch (TableError& err) {
358 String errMesg = err.getMesg();
359 if (errMesg == "Table column " + colName + " is unknown") {
360 table_.addColumn(ArrayColumnDesc<T>(colName));
361 col.attach(table_, colName);
362
363 int size = 0;
364 ArrayIterator<T2>& it = defValue.begin();
365 while (it != defValue.end()) {
366 ++size;
367 ++it;
368 }
369 IPosition ip(1, size);
370 Array<T>& arr(ip);
371 for (int i = 0; i < size; ++i)
372 arr[i] = static_cast<T>(defValue[i]);
373
374 col.fillColumn(arr);
375 } else {
376 throw;
377 }
378 } catch (...) {
379 throw;
380 }
381}
382
383void Scantable::setHeader(const STHeader& sdh)
384{
385 table_.rwKeywordSet().define("nIF", sdh.nif);
386 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
387 table_.rwKeywordSet().define("nPol", sdh.npol);
388 table_.rwKeywordSet().define("nChan", sdh.nchan);
389 table_.rwKeywordSet().define("Observer", sdh.observer);
390 table_.rwKeywordSet().define("Project", sdh.project);
391 table_.rwKeywordSet().define("Obstype", sdh.obstype);
392 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
393 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
394 table_.rwKeywordSet().define("Equinox", sdh.equinox);
395 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
396 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
397 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
398 table_.rwKeywordSet().define("UTC", sdh.utc);
399 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
400 table_.rwKeywordSet().define("Epoch", sdh.epoch);
401 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
402}
403
404STHeader Scantable::getHeader() const
405{
406 STHeader sdh;
407 table_.keywordSet().get("nBeam",sdh.nbeam);
408 table_.keywordSet().get("nIF",sdh.nif);
409 table_.keywordSet().get("nPol",sdh.npol);
410 table_.keywordSet().get("nChan",sdh.nchan);
411 table_.keywordSet().get("Observer", sdh.observer);
412 table_.keywordSet().get("Project", sdh.project);
413 table_.keywordSet().get("Obstype", sdh.obstype);
414 table_.keywordSet().get("AntennaName", sdh.antennaname);
415 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
416 table_.keywordSet().get("Equinox", sdh.equinox);
417 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
418 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
419 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
420 table_.keywordSet().get("UTC", sdh.utc);
421 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
422 table_.keywordSet().get("Epoch", sdh.epoch);
423 table_.keywordSet().get("POLTYPE", sdh.poltype);
424 return sdh;
425}
426
427void Scantable::setSourceType( int stype )
428{
429 if ( stype < 0 || stype > 1 )
430 throw(AipsError("Illegal sourcetype."));
431 TableVector<Int> tabvec(table_, "SRCTYPE");
432 tabvec = Int(stype);
433}
434
435bool Scantable::conformant( const Scantable& other )
436{
437 return this->getHeader().conformant(other.getHeader());
438}
439
440
441
442std::string Scantable::formatSec(Double x) const
443{
444 Double xcop = x;
445 MVTime mvt(xcop/24./3600.); // make days
446
447 if (x < 59.95)
448 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
449 else if (x < 3599.95)
450 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
451 else {
452 ostringstream oss;
453 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
454 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
455 return String(oss);
456 }
457};
458
459std::string Scantable::formatDirection(const MDirection& md) const
460{
461 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
462 Int prec = 7;
463
464 MVAngle mvLon(t[0]);
465 String sLon = mvLon.string(MVAngle::TIME,prec);
466 uInt tp = md.getRef().getType();
467 if (tp == MDirection::GALACTIC ||
468 tp == MDirection::SUPERGAL ) {
469 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
470 }
471 MVAngle mvLat(t[1]);
472 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
473 return sLon + String(" ") + sLat;
474}
475
476
477std::string Scantable::getFluxUnit() const
478{
479 return table_.keywordSet().asString("FluxUnit");
480}
481
482void Scantable::setFluxUnit(const std::string& unit)
483{
484 String tmp(unit);
485 Unit tU(tmp);
486 if (tU==Unit("K") || tU==Unit("Jy")) {
487 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
488 } else {
489 throw AipsError("Illegal unit - must be compatible with Jy or K");
490 }
491}
492
493void Scantable::setInstrument(const std::string& name)
494{
495 bool throwIt = true;
496 // create an Instrument to see if this is valid
497 STAttr::convertInstrument(name, throwIt);
498 String nameU(name);
499 nameU.upcase();
500 table_.rwKeywordSet().define(String("AntennaName"), nameU);
501}
502
503void Scantable::setFeedType(const std::string& feedtype)
504{
505 if ( Scantable::factories_.find(feedtype) == Scantable::factories_.end() ) {
506 std::string msg = "Illegal feed type "+ feedtype;
507 throw(casa::AipsError(msg));
508 }
509 table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
510}
511
512MPosition Scantable::getAntennaPosition() const
513{
514 Vector<Double> antpos;
515 table_.keywordSet().get("AntennaPosition", antpos);
516 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
517 return MPosition(mvpos);
518}
519
520void Scantable::makePersistent(const std::string& filename)
521{
522 String inname(filename);
523 Path path(inname);
524 /// @todo reindex SCANNO, recompute nbeam, nif, npol
525 inname = path.expandedName();
526 // WORKAROUND !!! for Table bug
527 // Remove when fixed in casacore
528 if ( table_.tableType() == Table::Memory && !selector_.empty() ) {
529 Table tab = table_.copyToMemoryTable(generateName());
530 tab.deepCopy(inname, Table::New);
531 tab.markForDelete();
532
533 } else {
534 table_.deepCopy(inname, Table::New);
535 }
536}
537
538int Scantable::nbeam( int scanno ) const
539{
540 if ( scanno < 0 ) {
541 Int n;
542 table_.keywordSet().get("nBeam",n);
543 return int(n);
544 } else {
545 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
546 Table t = table_(table_.col("SCANNO") == scanno);
547 ROTableRow row(t);
548 const TableRecord& rec = row.get(0);
549 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
550 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
551 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
552 ROTableVector<uInt> v(subt, "BEAMNO");
553 return int(v.nelements());
554 }
555 return 0;
556}
557
558int Scantable::nif( int scanno ) const
559{
560 if ( scanno < 0 ) {
561 Int n;
562 table_.keywordSet().get("nIF",n);
563 return int(n);
564 } else {
565 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
566 Table t = table_(table_.col("SCANNO") == scanno);
567 ROTableRow row(t);
568 const TableRecord& rec = row.get(0);
569 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
570 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
571 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
572 if ( subt.nrow() == 0 ) return 0;
573 ROTableVector<uInt> v(subt, "IFNO");
574 return int(v.nelements());
575 }
576 return 0;
577}
578
579int Scantable::npol( int scanno ) const
580{
581 if ( scanno < 0 ) {
582 Int n;
583 table_.keywordSet().get("nPol",n);
584 return n;
585 } else {
586 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
587 Table t = table_(table_.col("SCANNO") == scanno);
588 ROTableRow row(t);
589 const TableRecord& rec = row.get(0);
590 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
591 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
592 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
593 if ( subt.nrow() == 0 ) return 0;
594 ROTableVector<uInt> v(subt, "POLNO");
595 return int(v.nelements());
596 }
597 return 0;
598}
599
600int Scantable::ncycle( int scanno ) const
601{
602 if ( scanno < 0 ) {
603 Block<String> cols(2);
604 cols[0] = "SCANNO";
605 cols[1] = "CYCLENO";
606 TableIterator it(table_, cols);
607 int n = 0;
608 while ( !it.pastEnd() ) {
609 ++n;
610 ++it;
611 }
612 return n;
613 } else {
614 Table t = table_(table_.col("SCANNO") == scanno);
615 ROTableRow row(t);
616 const TableRecord& rec = row.get(0);
617 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
618 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
619 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
620 if ( subt.nrow() == 0 ) return 0;
621 return int(subt.nrow());
622 }
623 return 0;
624}
625
626
627int Scantable::nrow( int scanno ) const
628{
629 return int(table_.nrow());
630}
631
632int Scantable::nchan( int ifno ) const
633{
634 if ( ifno < 0 ) {
635 Int n;
636 table_.keywordSet().get("nChan",n);
637 return int(n);
638 } else {
639 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
640 // vary with these
641 Table t = table_(table_.col("IFNO") == ifno);
642 if ( t.nrow() == 0 ) return 0;
643 ROArrayColumn<Float> v(t, "SPECTRA");
644 return v.shape(0)(0);
645 }
646 return 0;
647}
648
649int Scantable::nscan() const {
650 Vector<uInt> scannos(scanCol_.getColumn());
651 uInt nout = genSort( scannos, Sort::Ascending,
652 Sort::QuickSort|Sort::NoDuplicates );
653 return int(nout);
654}
655
656int Scantable::getChannels(int whichrow) const
657{
658 return specCol_.shape(whichrow)(0);
659}
660
661int Scantable::getBeam(int whichrow) const
662{
663 return beamCol_(whichrow);
664}
665
666std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
667{
668 Vector<uInt> nos(col.getColumn());
669 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
670 nos.resize(n, True);
671 std::vector<uint> stlout;
672 nos.tovector(stlout);
673 return stlout;
674}
675
676int Scantable::getIF(int whichrow) const
677{
678 return ifCol_(whichrow);
679}
680
681int Scantable::getPol(int whichrow) const
682{
683 return polCol_(whichrow);
684}
685
686std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
687{
688 return formatTime(me, showdate, 0);
689}
690
691std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
692{
693 MVTime mvt(me.getValue());
694 if (showdate)
695 //mvt.setFormat(MVTime::YMD);
696 mvt.setFormat(MVTime::YMD, prec);
697 else
698 //mvt.setFormat(MVTime::TIME);
699 mvt.setFormat(MVTime::TIME, prec);
700 ostringstream oss;
701 oss << mvt;
702 return String(oss);
703}
704
705void Scantable::calculateAZEL()
706{
707 MPosition mp = getAntennaPosition();
708 MEpoch::ROScalarColumn timeCol(table_, "TIME");
709 ostringstream oss;
710 oss << "Computed azimuth/elevation using " << endl
711 << mp << endl;
712 for (Int i=0; i<nrow(); ++i) {
713 MEpoch me = timeCol(i);
714 MDirection md = getDirection(i);
715 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
716 << endl << " => ";
717 MeasFrame frame(mp, me);
718 Vector<Double> azel =
719 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
720 frame)
721 )().getAngle("rad").getValue();
722 azCol_.put(i,Float(azel[0]));
723 elCol_.put(i,Float(azel[1]));
724 oss << "azel: " << azel[0]/C::pi*180.0 << " "
725 << azel[1]/C::pi*180.0 << " (deg)" << endl;
726 }
727 pushLog(String(oss));
728}
729
730void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
731{
732 for (uInt i=0; i<table_.nrow(); ++i) {
733 Vector<uChar> flgs = flagsCol_(i);
734 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
735 flagsCol_.put(i, flgs);
736 }
737}
738
739std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
740{
741 Vector<uChar> flags;
742 flagsCol_.get(uInt(whichrow), flags);
743 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
744 Vector<Bool> bflag(flags.shape());
745 convertArray(bflag, flags);
746 //bflag = !bflag;
747
748 std::vector<bool> mask;
749 bflag.tovector(mask);
750 return mask;
751}
752
753void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
754 Vector<uChar> flgs)
755{
756 Vector<Float> spcs = specCol_(whichrow);
757 uInt nchannel = nchan();
758 if (spcs.nelements() != nchannel) {
759 throw(AipsError("Data has incorrect number of channels"));
760 }
761 uChar userflag = 1 << 7;
762 if (unflag) {
763 userflag = 0 << 7;
764 }
765 if (clipoutside) {
766 for (uInt j = 0; j < nchannel; ++j) {
767 Float spc = spcs(j);
768 if ((spc >= uthres) || (spc <= dthres)) {
769 flgs(j) = userflag;
770 }
771 }
772 } else {
773 for (uInt j = 0; j < nchannel; ++j) {
774 Float spc = spcs(j);
775 if ((spc < uthres) && (spc > dthres)) {
776 flgs(j) = userflag;
777 }
778 }
779 }
780}
781
782
783void Scantable::flag( int whichrow, const std::vector<bool>& msk, bool unflag ) {
784 std::vector<bool>::const_iterator it;
785 uInt ntrue = 0;
786 if (whichrow >= int(table_.nrow()) ) {
787 throw(AipsError("Invalid row number"));
788 }
789 for (it = msk.begin(); it != msk.end(); ++it) {
790 if ( *it ) {
791 ntrue++;
792 }
793 }
794 //if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
795 if ( whichrow == -1 && !unflag && selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
796 throw(AipsError("Trying to flag whole scantable."));
797 uChar userflag = 1 << 7;
798 if ( unflag ) {
799 userflag = 0 << 7;
800 }
801 if (whichrow > -1 ) {
802 applyChanFlag(uInt(whichrow), msk, userflag);
803 } else {
804 for ( uInt i=0; i<table_.nrow(); ++i) {
805 applyChanFlag(i, msk, userflag);
806 }
807 }
808}
809
810void Scantable::applyChanFlag( uInt whichrow, const std::vector<bool>& msk, uChar flagval )
811{
812 if (whichrow >= table_.nrow() ) {
813 throw( casa::indexError<int>( whichrow, "asap::Scantable::applyChanFlag: Invalid row number" ) );
814 }
815 Vector<uChar> flgs = flagsCol_(whichrow);
816 if ( msk.size() == 0 ) {
817 flgs = flagval;
818 flagsCol_.put(whichrow, flgs);
819 return;
820 }
821 if ( int(msk.size()) != nchan() ) {
822 throw(AipsError("Mask has incorrect number of channels."));
823 }
824 if ( flgs.nelements() != msk.size() ) {
825 throw(AipsError("Mask has incorrect number of channels."
826 " Probably varying with IF. Please flag per IF"));
827 }
828 std::vector<bool>::const_iterator it;
829 uInt j = 0;
830 for (it = msk.begin(); it != msk.end(); ++it) {
831 if ( *it ) {
832 flgs(j) = flagval;
833 }
834 ++j;
835 }
836 flagsCol_.put(whichrow, flgs);
837}
838
839void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
840{
841 if ( selector_.empty() && (rows.size() == table_.nrow()) )
842 throw(AipsError("Trying to flag whole scantable."));
843
844 uInt rowflag = (unflag ? 0 : 1);
845 std::vector<uInt>::const_iterator it;
846 for (it = rows.begin(); it != rows.end(); ++it)
847 flagrowCol_.put(*it, rowflag);
848}
849
850std::vector<bool> Scantable::getMask(int whichrow) const
851{
852 Vector<uChar> flags;
853 flagsCol_.get(uInt(whichrow), flags);
854 Vector<Bool> bflag(flags.shape());
855 convertArray(bflag, flags);
856 bflag = !bflag;
857 std::vector<bool> mask;
858 bflag.tovector(mask);
859 return mask;
860}
861
862std::vector<float> Scantable::getSpectrum( int whichrow,
863 const std::string& poltype ) const
864{
865 String ptype = poltype;
866 if (poltype == "" ) ptype = getPolType();
867 if ( whichrow < 0 || whichrow >= nrow() )
868 throw(AipsError("Illegal row number."));
869 std::vector<float> out;
870 Vector<Float> arr;
871 uInt requestedpol = polCol_(whichrow);
872 String basetype = getPolType();
873 if ( ptype == basetype ) {
874 specCol_.get(whichrow, arr);
875 } else {
876 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
877 basetype));
878 uInt row = uInt(whichrow);
879 stpol->setSpectra(getPolMatrix(row));
880 Float fang,fhand,parang;
881 fang = focusTable_.getTotalAngle(mfocusidCol_(row));
882 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
883 stpol->setPhaseCorrections(fang, fhand);
884 arr = stpol->getSpectrum(requestedpol, ptype);
885 }
886 if ( arr.nelements() == 0 )
887 pushLog("Not enough polarisations present to do the conversion.");
888 arr.tovector(out);
889 return out;
890}
891
892void Scantable::setSpectrum( const std::vector<float>& spec,
893 int whichrow )
894{
895 Vector<Float> spectrum(spec);
896 Vector<Float> arr;
897 specCol_.get(whichrow, arr);
898 if ( spectrum.nelements() != arr.nelements() )
899 throw AipsError("The spectrum has incorrect number of channels.");
900 specCol_.put(whichrow, spectrum);
901}
902
903
904String Scantable::generateName()
905{
906 return (File::newUniqueName("./","temp")).baseName();
907}
908
909const casa::Table& Scantable::table( ) const
910{
911 return table_;
912}
913
914casa::Table& Scantable::table( )
915{
916 return table_;
917}
918
919std::string Scantable::getPolType() const
920{
921 return table_.keywordSet().asString("POLTYPE");
922}
923
924void Scantable::unsetSelection()
925{
926 table_ = originalTable_;
927 attach();
928 selector_.reset();
929}
930
931void Scantable::setSelection( const STSelector& selection )
932{
933 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
934 if ( tab.nrow() == 0 ) {
935 throw(AipsError("Selection contains no data. Not applying it."));
936 }
937 table_ = tab;
938 attach();
939 selector_ = selection;
940}
941
942std::string Scantable::summary( bool verbose )
943{
944 // Format header info
945 ostringstream oss;
946 oss << endl;
947 oss << asap::SEPERATOR << endl;
948 oss << " Scan Table Summary" << endl;
949 oss << asap::SEPERATOR << endl;
950 oss.flags(std::ios_base::left);
951 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
952 << setw(15) << "IFs:" << setw(4) << nif() << endl
953 << setw(15) << "Polarisations:" << setw(4) << npol()
954 << "(" << getPolType() << ")" << endl
955 << setw(15) << "Channels:" << nchan() << endl;
956 String tmp;
957 oss << setw(15) << "Observer:"
958 << table_.keywordSet().asString("Observer") << endl;
959 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
960 table_.keywordSet().get("Project", tmp);
961 oss << setw(15) << "Project:" << tmp << endl;
962 table_.keywordSet().get("Obstype", tmp);
963 oss << setw(15) << "Obs. Type:" << tmp << endl;
964 table_.keywordSet().get("AntennaName", tmp);
965 oss << setw(15) << "Antenna Name:" << tmp << endl;
966 table_.keywordSet().get("FluxUnit", tmp);
967 oss << setw(15) << "Flux Unit:" << tmp << endl;
968 //Vector<Double> vec(moleculeTable_.getRestFrequencies());
969 int nid = moleculeTable_.nrow();
970 Bool firstline = True;
971 oss << setw(15) << "Rest Freqs:";
972 for (int i=0; i<nid; i++) {
973 Table t = table_(table_.col("MOLECULE_ID") == i);
974 if (t.nrow() > 0) {
975 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
976 if (vec.nelements() > 0) {
977 if (firstline) {
978 oss << setprecision(10) << vec << " [Hz]" << endl;
979 firstline=False;
980 }
981 else{
982 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
983 }
984 } else {
985 oss << "none" << endl;
986 }
987 }
988 }
989
990 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
991 oss << selector_.print() << endl;
992 oss << endl;
993 // main table
994 String dirtype = "Position ("
995 + getDirectionRefString()
996 + ")";
997 oss << setw(5) << "Scan" << setw(15) << "Source"
998 << setw(10) << "Time" << setw(18) << "Integration"
999 << setw(15) << "Source Type" << endl;
1000 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
1001 oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
1002 << setw(8) << "Frame" << setw(16)
1003 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
1004 << setw(7) << "Channels"
1005 << endl;
1006 oss << asap::SEPERATOR << endl;
1007 TableIterator iter(table_, "SCANNO");
1008 while (!iter.pastEnd()) {
1009 Table subt = iter.table();
1010 ROTableRow row(subt);
1011 MEpoch::ROScalarColumn timeCol(subt,"TIME");
1012 const TableRecord& rec = row.get(0);
1013 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1014 << std::left << setw(1) << ""
1015 << setw(15) << rec.asString("SRCNAME")
1016 << setw(10) << formatTime(timeCol(0), false);
1017 // count the cycles in the scan
1018 TableIterator cyciter(subt, "CYCLENO");
1019 int nint = 0;
1020 while (!cyciter.pastEnd()) {
1021 ++nint;
1022 ++cyciter;
1023 }
1024 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
1025 << setw(11) << formatSec(rec.asFloat("INTERVAL")) << setw(1) << ""
1026 << setw(15) << SrcType::getName(rec.asInt("SRCTYPE")) << endl;
1027
1028 TableIterator biter(subt, "BEAMNO");
1029 while (!biter.pastEnd()) {
1030 Table bsubt = biter.table();
1031 ROTableRow brow(bsubt);
1032 const TableRecord& brec = brow.get(0);
1033 uInt row0 = bsubt.rowNumbers(table_)[0];
1034 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
1035 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
1036 TableIterator iiter(bsubt, "IFNO");
1037 while (!iiter.pastEnd()) {
1038 Table isubt = iiter.table();
1039 ROTableRow irow(isubt);
1040 const TableRecord& irec = irow.get(0);
1041 oss << setw(9) << "";
1042 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
1043 << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1044 << setw(3) << "" << nchan(irec.asuInt("IFNO"))
1045 << endl;
1046
1047 ++iiter;
1048 }
1049 ++biter;
1050 }
1051 ++iter;
1052 }
1053 /// @todo implement verbose mode
1054 return String(oss);
1055}
1056
1057// std::string Scantable::getTime(int whichrow, bool showdate) const
1058// {
1059// MEpoch::ROScalarColumn timeCol(table_, "TIME");
1060// MEpoch me;
1061// if (whichrow > -1) {
1062// me = timeCol(uInt(whichrow));
1063// } else {
1064// Double tm;
1065// table_.keywordSet().get("UTC",tm);
1066// me = MEpoch(MVEpoch(tm));
1067// }
1068// return formatTime(me, showdate);
1069// }
1070
1071std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
1072{
1073 MEpoch me;
1074 me = getEpoch(whichrow);
1075 return formatTime(me, showdate, prec);
1076}
1077
1078MEpoch Scantable::getEpoch(int whichrow) const
1079{
1080 if (whichrow > -1) {
1081 return timeCol_(uInt(whichrow));
1082 } else {
1083 Double tm;
1084 table_.keywordSet().get("UTC",tm);
1085 return MEpoch(MVEpoch(tm));
1086 }
1087}
1088
1089std::string Scantable::getDirectionString(int whichrow) const
1090{
1091 return formatDirection(getDirection(uInt(whichrow)));
1092}
1093
1094
1095SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1096 const MPosition& mp = getAntennaPosition();
1097 const MDirection& md = getDirection(whichrow);
1098 const MEpoch& me = timeCol_(whichrow);
1099 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1100 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1101 return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1102 mfreqidCol_(whichrow));
1103}
1104
1105std::vector< double > Scantable::getAbcissa( int whichrow ) const
1106{
1107 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
1108 std::vector<double> stlout;
1109 int nchan = specCol_(whichrow).nelements();
1110 String us = freqTable_.getUnitString();
1111 if ( us == "" || us == "pixel" || us == "channel" ) {
1112 for (int i=0; i<nchan; ++i) {
1113 stlout.push_back(double(i));
1114 }
1115 return stlout;
1116 }
1117 SpectralCoordinate spc = getSpectralCoordinate(whichrow);
1118 Vector<Double> pixel(nchan);
1119 Vector<Double> world;
1120 indgen(pixel);
1121 if ( Unit(us) == Unit("Hz") ) {
1122 for ( int i=0; i < nchan; ++i) {
1123 Double world;
1124 spc.toWorld(world, pixel[i]);
1125 stlout.push_back(double(world));
1126 }
1127 } else if ( Unit(us) == Unit("km/s") ) {
1128 Vector<Double> world;
1129 spc.pixelToVelocity(world, pixel);
1130 world.tovector(stlout);
1131 }
1132 return stlout;
1133}
1134void Scantable::setDirectionRefString( const std::string & refstr )
1135{
1136 MDirection::Types mdt;
1137 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1138 throw(AipsError("Illegal Direction frame."));
1139 }
1140 if ( refstr == "" ) {
1141 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1142 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1143 } else {
1144 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1145 }
1146}
1147
1148std::string Scantable::getDirectionRefString( ) const
1149{
1150 return table_.keywordSet().asString("DIRECTIONREF");
1151}
1152
1153MDirection Scantable::getDirection(int whichrow ) const
1154{
1155 String usertype = table_.keywordSet().asString("DIRECTIONREF");
1156 String type = MDirection::showType(dirCol_.getMeasRef().getType());
1157 if ( usertype != type ) {
1158 MDirection::Types mdt;
1159 if (!MDirection::getType(mdt, usertype)) {
1160 throw(AipsError("Illegal Direction frame."));
1161 }
1162 return dirCol_.convert(uInt(whichrow), mdt);
1163 } else {
1164 return dirCol_(uInt(whichrow));
1165 }
1166}
1167
1168std::string Scantable::getAbcissaLabel( int whichrow ) const
1169{
1170 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
1171 const MPosition& mp = getAntennaPosition();
1172 const MDirection& md = getDirection(whichrow);
1173 const MEpoch& me = timeCol_(whichrow);
1174 //const Double& rf = mmolidCol_(whichrow);
1175 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1176 SpectralCoordinate spc =
1177 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1178
1179 String s = "Channel";
1180 Unit u = Unit(freqTable_.getUnitString());
1181 if (u == Unit("km/s")) {
1182 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
1183 } else if (u == Unit("Hz")) {
1184 Vector<String> wau(1);wau = u.getName();
1185 spc.setWorldAxisUnits(wau);
1186 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
1187 }
1188 return s;
1189
1190}
1191
1192/**
1193void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1194 const std::string& unit )
1195**/
1196void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1197 const std::string& unit )
1198
1199{
1200 ///@todo lookup in line table to fill in name and formattedname
1201 Unit u(unit);
1202 //Quantum<Double> urf(rf, u);
1203 Quantum<Vector<Double> >urf(rf, u);
1204 Vector<String> formattedname(0);
1205 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1206
1207 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1208 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1209 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1210 tabvec = id;
1211}
1212
1213/**
1214void asap::Scantable::setRestFrequencies( const std::string& name )
1215{
1216 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1217 ///@todo implement
1218}
1219**/
1220
1221void Scantable::setRestFrequencies( const vector<std::string>& name )
1222{
1223 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1224 ///@todo implement
1225}
1226
1227std::vector< unsigned int > Scantable::rownumbers( ) const
1228{
1229 std::vector<unsigned int> stlout;
1230 Vector<uInt> vec = table_.rowNumbers();
1231 vec.tovector(stlout);
1232 return stlout;
1233}
1234
1235
1236Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1237{
1238 ROTableRow row(table_);
1239 const TableRecord& rec = row.get(whichrow);
1240 Table t =
1241 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1242 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1243 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1244 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1245 ROArrayColumn<Float> speccol(t, "SPECTRA");
1246 return speccol.getColumn();
1247}
1248
1249std::vector< std::string > Scantable::columnNames( ) const
1250{
1251 Vector<String> vec = table_.tableDesc().columnNames();
1252 return mathutil::tovectorstring(vec);
1253}
1254
1255MEpoch::Types Scantable::getTimeReference( ) const
1256{
1257 return MEpoch::castType(timeCol_.getMeasRef().getType());
1258}
1259
1260void Scantable::addFit( const STFitEntry& fit, int row )
1261{
1262 //cout << mfitidCol_(uInt(row)) << endl;
1263 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1264 os << mfitidCol_(uInt(row)) << LogIO::POST ;
1265 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1266 mfitidCol_.put(uInt(row), id);
1267}
1268
1269void Scantable::shift(int npix)
1270{
1271 Vector<uInt> fids(mfreqidCol_.getColumn());
1272 genSort( fids, Sort::Ascending,
1273 Sort::QuickSort|Sort::NoDuplicates );
1274 for (uInt i=0; i<fids.nelements(); ++i) {
1275 frequencies().shiftRefPix(npix, fids[i]);
1276 }
1277}
1278
1279String Scantable::getAntennaName() const
1280{
1281 String out;
1282 table_.keywordSet().get("AntennaName", out);
1283 String::size_type pos1 = out.find("@") ;
1284 String::size_type pos2 = out.find("//") ;
1285 if ( pos2 != String::npos )
1286 out = out.substr(pos2+2,pos1) ;
1287 else if ( pos1 != String::npos )
1288 out = out.substr(0,pos1) ;
1289 return out;
1290}
1291
1292int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1293{
1294 String tbpath;
1295 int ret = 0;
1296 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1297 table_.keywordSet().get("GBT_GO", tbpath);
1298 Table t(tbpath,Table::Old);
1299 // check each scan if other scan of the pair exist
1300 int nscan = scanlist.size();
1301 for (int i = 0; i < nscan; i++) {
1302 Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1303 if (subt.nrow()==0) {
1304 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1305 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1306 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1307 ret = 1;
1308 break;
1309 }
1310 ROTableRow row(subt);
1311 const TableRecord& rec = row.get(0);
1312 int scan1seqn = rec.asuInt("PROCSEQN");
1313 int laston1 = rec.asuInt("LASTON");
1314 if ( rec.asuInt("PROCSIZE")==2 ) {
1315 if ( i < nscan-1 ) {
1316 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1317 if ( subt2.nrow() == 0) {
1318 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1319
1320 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1321 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1322 ret = 1;
1323 break;
1324 }
1325 ROTableRow row2(subt2);
1326 const TableRecord& rec2 = row2.get(0);
1327 int scan2seqn = rec2.asuInt("PROCSEQN");
1328 int laston2 = rec2.asuInt("LASTON");
1329 if (scan1seqn == 1 && scan2seqn == 2) {
1330 if (laston1 == laston2) {
1331 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1332 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1333 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1334 i +=1;
1335 }
1336 else {
1337 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1338 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1339 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1340 }
1341 }
1342 else if (scan1seqn==2 && scan2seqn == 1) {
1343 if (laston1 == laston2) {
1344 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1345 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1346 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1347 ret = 1;
1348 break;
1349 }
1350 }
1351 else {
1352 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1353 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1354 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1355 ret = 1;
1356 break;
1357 }
1358 }
1359 }
1360 else {
1361 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1362 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1363 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1364 }
1365 //if ( i >= nscan ) break;
1366 }
1367 }
1368 else {
1369 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1370 //cerr<<"No reference to GBT_GO table."<<endl;
1371 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1372 ret = 1;
1373 }
1374 return ret;
1375}
1376
1377std::vector<double> Scantable::getDirectionVector(int whichrow) const
1378{
1379 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1380 std::vector<double> dir;
1381 Dir.tovector(dir);
1382 return dir;
1383}
1384
1385void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1386 throw( casa::AipsError )
1387{
1388 // assumed that all rows have same nChan
1389 Vector<Float> arr = specCol_( 0 ) ;
1390 int nChan = arr.nelements() ;
1391
1392 // if nmin < 0 or nmax < 0, nothing to do
1393 if ( nmin < 0 ) {
1394 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1395 }
1396 if ( nmax < 0 ) {
1397 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1398 }
1399
1400 // if nmin > nmax, exchange values
1401 if ( nmin > nmax ) {
1402 int tmp = nmax ;
1403 nmax = nmin ;
1404 nmin = tmp ;
1405 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1406 os << "Swap values. Applied range is ["
1407 << nmin << ", " << nmax << "]" << LogIO::POST ;
1408 }
1409
1410 // if nmin exceeds nChan, nothing to do
1411 if ( nmin >= nChan ) {
1412 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1413 }
1414
1415 // if nmax exceeds nChan, reset nmax to nChan
1416 if ( nmax >= nChan ) {
1417 if ( nmin == 0 ) {
1418 // nothing to do
1419 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1420 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1421 return ;
1422 }
1423 else {
1424 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1425 os << "Specified maximum exceeds nChan. Applied range is ["
1426 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1427 nmax = nChan - 1 ;
1428 }
1429 }
1430
1431 // reshape specCol_ and flagCol_
1432 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1433 reshapeSpectrum( nmin, nmax, irow ) ;
1434 }
1435
1436 // update FREQUENCIES subtable
1437 Double refpix ;
1438 Double refval ;
1439 Double increment ;
1440 int freqnrow = freqTable_.table().nrow() ;
1441 Vector<uInt> oldId( freqnrow ) ;
1442 Vector<uInt> newId( freqnrow ) ;
1443 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1444 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1445 /***
1446 * need to shift refpix to nmin
1447 * note that channel nmin in old index will be channel 0 in new one
1448 ***/
1449 refval = refval - ( refpix - nmin ) * increment ;
1450 refpix = 0 ;
1451 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1452 }
1453
1454 // update nchan
1455 int newsize = nmax - nmin + 1 ;
1456 table_.rwKeywordSet().define( "nChan", newsize ) ;
1457
1458 // update bandwidth
1459 // assumed all spectra in the scantable have same bandwidth
1460 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1461
1462 return ;
1463}
1464
1465void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1466{
1467 // reshape specCol_ and flagCol_
1468 Vector<Float> oldspec = specCol_( irow ) ;
1469 Vector<uChar> oldflag = flagsCol_( irow ) ;
1470 uInt newsize = nmax - nmin + 1 ;
1471 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1472 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1473
1474 return ;
1475}
1476
1477void asap::Scantable::regridChannel( int nChan, double dnu )
1478{
1479 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1480 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1481 // assumed that all rows have same nChan
1482 Vector<Float> arr = specCol_( 0 ) ;
1483 int oldsize = arr.nelements() ;
1484
1485 // if oldsize == nChan, nothing to do
1486 if ( oldsize == nChan ) {
1487 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1488 return ;
1489 }
1490
1491 // if oldChan < nChan, unphysical operation
1492 if ( oldsize < nChan ) {
1493 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1494 return ;
1495 }
1496
1497 // change channel number for specCol_ and flagCol_
1498 Vector<Float> newspec( nChan, 0 ) ;
1499 Vector<uChar> newflag( nChan, false ) ;
1500 vector<string> coordinfo = getCoordInfo() ;
1501 string oldinfo = coordinfo[0] ;
1502 coordinfo[0] = "Hz" ;
1503 setCoordInfo( coordinfo ) ;
1504 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1505 regridChannel( nChan, dnu, irow ) ;
1506 }
1507 coordinfo[0] = oldinfo ;
1508 setCoordInfo( coordinfo ) ;
1509
1510
1511 // NOTE: this method does not update metadata such as
1512 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1513
1514 return ;
1515}
1516
1517void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1518{
1519 // logging
1520 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1521 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1522
1523 Vector<Float> oldspec = specCol_( irow ) ;
1524 Vector<uChar> oldflag = flagsCol_( irow ) ;
1525 Vector<Float> newspec( nChan, 0 ) ;
1526 Vector<uChar> newflag( nChan, false ) ;
1527
1528 // regrid
1529 vector<double> abcissa = getAbcissa( irow ) ;
1530 int oldsize = abcissa.size() ;
1531 double olddnu = abcissa[1] - abcissa[0] ;
1532 //int refChan = 0 ;
1533 //double frac = 0.0 ;
1534 //double wedge = 0.0 ;
1535 //double pile = 0.0 ;
1536 int ichan = 0 ;
1537 double wsum = 0.0 ;
1538 Vector<Float> z( nChan ) ;
1539 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
1540 for ( int ii = 1 ; ii < nChan ; ii++ )
1541 z[ii] = z[ii-1] + dnu ;
1542 Vector<Float> zi( nChan+1 ) ;
1543 Vector<Float> yi( oldsize + 1 ) ;
1544 zi[0] = z[0] - 0.5 * dnu ;
1545 zi[1] = z[0] + 0.5 * dnu ;
1546 for ( int ii = 2 ; ii < nChan ; ii++ )
1547 zi[ii] = zi[ii-1] + dnu ;
1548 zi[nChan] = z[nChan-1] + 0.5 * dnu ;
1549 yi[0] = abcissa[0] - 0.5 * olddnu ;
1550 yi[1] = abcissa[1] + 0.5 * olddnu ;
1551 for ( int ii = 2 ; ii < oldsize ; ii++ )
1552 yi[ii] = abcissa[ii-1] + olddnu ;
1553 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
1554 if ( dnu > 0.0 ) {
1555 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1556 double zl = zi[ii] ;
1557 double zr = zi[ii+1] ;
1558 for ( int j = ichan ; j < oldsize ; j++ ) {
1559 double yl = yi[j] ;
1560 double yr = yi[j+1] ;
1561 if ( yl <= zl ) {
1562 if ( yr <= zl ) {
1563 continue ;
1564 }
1565 else if ( yr <= zr ) {
1566 newspec[ii] += oldspec[j] * ( yr - zl ) ;
1567 newflag[ii] = newflag[ii] || oldflag[j] ;
1568 wsum += ( yr - zl ) ;
1569 }
1570 else {
1571 newspec[ii] += oldspec[j] * dnu ;
1572 newflag[ii] = newflag[ii] || oldflag[j] ;
1573 wsum += dnu ;
1574 ichan = j ;
1575 break ;
1576 }
1577 }
1578 else if ( yl < zr ) {
1579 if ( yr <= zr ) {
1580 newspec[ii] += oldspec[j] * ( yr - yl ) ;
1581 newflag[ii] = newflag[ii] || oldflag[j] ;
1582 wsum += ( yr - yl ) ;
1583 }
1584 else {
1585 newspec[ii] += oldspec[j] * ( zr - yl ) ;
1586 newflag[ii] = newflag[ii] || oldflag[j] ;
1587 wsum += ( zr - yl ) ;
1588 ichan = j ;
1589 break ;
1590 }
1591 }
1592 else {
1593 ichan = j - 1 ;
1594 break ;
1595 }
1596 }
1597 newspec[ii] /= wsum ;
1598 wsum = 0.0 ;
1599 }
1600 }
1601 else if ( dnu < 0.0 ) {
1602 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1603 double zl = zi[ii] ;
1604 double zr = zi[ii+1] ;
1605 for ( int j = ichan ; j < oldsize ; j++ ) {
1606 double yl = yi[j] ;
1607 double yr = yi[j+1] ;
1608 if ( yl >= zl ) {
1609 if ( yr >= zl ) {
1610 continue ;
1611 }
1612 else if ( yr >= zr ) {
1613 newspec[ii] += oldspec[j] * abs( yr - zl ) ;
1614 newflag[ii] = newflag[ii] || oldflag[j] ;
1615 wsum += abs( yr - zl ) ;
1616 }
1617 else {
1618 newspec[ii] += oldspec[j] * abs( dnu ) ;
1619 newflag[ii] = newflag[ii] || oldflag[j] ;
1620 wsum += abs( dnu ) ;
1621 ichan = j ;
1622 break ;
1623 }
1624 }
1625 else if ( yl > zr ) {
1626 if ( yr >= zr ) {
1627 newspec[ii] += oldspec[j] * abs( yr - yl ) ;
1628 newflag[ii] = newflag[ii] || oldflag[j] ;
1629 wsum += abs( yr - yl ) ;
1630 }
1631 else {
1632 newspec[ii] += oldspec[j] * abs( zr - yl ) ;
1633 newflag[ii] = newflag[ii] || oldflag[j] ;
1634 wsum += abs( zr - yl ) ;
1635 ichan = j ;
1636 break ;
1637 }
1638 }
1639 else {
1640 ichan = j - 1 ;
1641 break ;
1642 }
1643 }
1644 newspec[ii] /= wsum ;
1645 wsum = 0.0 ;
1646 }
1647 }
1648// * ichan = 0
1649// ***/
1650// //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
1651// pile += dnu ;
1652// wedge = olddnu * ( refChan + 1 ) ;
1653// while ( wedge < pile ) {
1654// newspec[0] += olddnu * oldspec[refChan] ;
1655// newflag[0] = newflag[0] || oldflag[refChan] ;
1656// //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
1657// refChan++ ;
1658// wedge += olddnu ;
1659// wsum += olddnu ;
1660// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1661// }
1662// frac = ( wedge - pile ) / olddnu ;
1663// wsum += ( 1.0 - frac ) * olddnu ;
1664// newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1665// newflag[0] = newflag[0] || oldflag[refChan] ;
1666// //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
1667// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1668// newspec[0] /= wsum ;
1669// //ofs << "newspec[0] = " << newspec[0] << endl ;
1670// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1671
1672// /***
1673// * ichan = 1 - nChan-2
1674// ***/
1675// for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
1676// pile += dnu ;
1677// newspec[ichan] += frac * olddnu * oldspec[refChan] ;
1678// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1679// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
1680// refChan++ ;
1681// wedge += olddnu ;
1682// wsum = frac * olddnu ;
1683// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1684// while ( wedge < pile ) {
1685// newspec[ichan] += olddnu * oldspec[refChan] ;
1686// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1687// //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
1688// refChan++ ;
1689// wedge += olddnu ;
1690// wsum += olddnu ;
1691// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1692// }
1693// frac = ( wedge - pile ) / olddnu ;
1694// wsum += ( 1.0 - frac ) * olddnu ;
1695// newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1696// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1697// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
1698// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1699// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1700// newspec[ichan] /= wsum ;
1701// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
1702// }
1703
1704// /***
1705// * ichan = nChan-1
1706// ***/
1707// // NOTE: Assumed that all spectra have the same bandwidth
1708// pile += dnu ;
1709// newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
1710// newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
1711// //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1712// refChan++ ;
1713// wedge += olddnu ;
1714// wsum = frac * olddnu ;
1715// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1716// for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
1717// newspec[nChan-1] += olddnu * oldspec[jchan] ;
1718// newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
1719// wsum += olddnu ;
1720// //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1721// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1722// }
1723// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1724// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1725// newspec[nChan-1] /= wsum ;
1726// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
1727
1728// specCol_.put( irow, newspec ) ;
1729// flagsCol_.put( irow, newflag ) ;
1730
1731// // ofs.close() ;
1732
1733
1734 return ;
1735}
1736
1737std::vector<float> Scantable::getWeather(int whichrow) const
1738{
1739 std::vector<float> out(5);
1740 //Float temperature, pressure, humidity, windspeed, windaz;
1741 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
1742 mweatheridCol_(uInt(whichrow)));
1743
1744
1745 return out;
1746}
1747
1748bool Scantable::getFlagtraFast(int whichrow)
1749{
1750 uChar flag;
1751 Vector<uChar> flags;
1752 flagsCol_.get(uInt(whichrow), flags);
1753 flag = flags[0];
1754 for (int i = 1; i < flags.size(); ++i) {
1755 flag &= flags[i];
1756 }
1757 return ((flag >> 7) == 1);
1758}
1759
1760void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool outLogger, const std::string& blfile) {
1761 ofstream ofs;
1762 String coordInfo;
1763 bool hasSameNchan;
1764 int firstIF;
1765 bool outTextFile = false;
1766
1767 if (blfile != "") {
1768 ofs.open(blfile.c_str(), ios::out | ios::app);
1769 if (ofs) outTextFile = true;
1770 }
1771
1772 if (outLogger || outTextFile) {
1773 coordInfo = getCoordInfo()[0];
1774 if (coordInfo == "") coordInfo = "channel";
1775 hasSameNchan = hasSameNchanOverIFs();
1776 firstIF = getIF(0);
1777 }
1778
1779 //Fitter fitter = Fitter();
1780 //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
1781
1782 int nRow = nrow();
1783 std::vector<bool> chanMask;
1784
1785 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
1786
1787 chanMask = getCompositeChanMask(whichrow, mask);
1788 //fitBaseline(chanMask, whichrow, fitter);
1789 //setSpectrum(fitter.getResidual(), whichrow);
1790 std::vector<int> pieceRanges;
1791 std::vector<float> params;
1792 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
1793 setSpectrum(res, whichrow);
1794
1795 if (outLogger || outTextFile) {
1796 std::vector<bool> fixed;
1797 fixed.clear();
1798 float rms = getRms(chanMask, whichrow);
1799 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
1800
1801 if (outLogger) {
1802 LogIO ols(LogOrigin("Scantable", "cubicSplineBaseline()", WHERE));
1803 ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
1804 }
1805 if (outTextFile) {
1806 ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
1807 }
1808 }
1809
1810 }
1811
1812 if (outTextFile) ofs.close();
1813}
1814
1815
1816void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
1817{
1818 ofstream ofs;
1819 String coordInfo;
1820 bool hasSameNchan;
1821 int firstIF;
1822 bool outTextFile = false;
1823
1824 if (blfile != "") {
1825 ofs.open(blfile.c_str(), ios::out | ios::app);
1826 if (ofs) outTextFile = true;
1827 }
1828
1829 if (outLogger || outTextFile) {
1830 coordInfo = getCoordInfo()[0];
1831 if (coordInfo == "") coordInfo = "channel";
1832 hasSameNchan = hasSameNchanOverIFs();
1833 firstIF = getIF(0);
1834 }
1835
1836 //Fitter fitter = Fitter();
1837 //fitter.setExpression("cspline", nPiece, thresClip, nIterClip);
1838
1839 int nRow = nrow();
1840 std::vector<bool> chanMask;
1841 int minEdgeSize = getIFNos().size()*2;
1842 STLineFinder lineFinder = STLineFinder();
1843 lineFinder.setOptions(threshold, 3, chanAvgLimit);
1844
1845 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
1846
1847 //-------------------------------------------------------
1848 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
1849 //-------------------------------------------------------
1850 int edgeSize = edge.size();
1851 std::vector<int> currentEdge;
1852 if (edgeSize >= 2) {
1853 int idx = 0;
1854 if (edgeSize > 2) {
1855 if (edgeSize < minEdgeSize) {
1856 throw(AipsError("Length of edge element info is less than that of IFs"));
1857 }
1858 idx = 2 * getIF(whichrow);
1859 }
1860 currentEdge.push_back(edge[idx]);
1861 currentEdge.push_back(edge[idx+1]);
1862 } else {
1863 throw(AipsError("Wrong length of edge element"));
1864 }
1865 lineFinder.setData(getSpectrum(whichrow));
1866 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
1867 chanMask = lineFinder.getMask();
1868 //-------------------------------------------------------
1869
1870
1871 //fitBaseline(chanMask, whichrow, fitter);
1872 //setSpectrum(fitter.getResidual(), whichrow);
1873 std::vector<int> pieceRanges;
1874 std::vector<float> params;
1875 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, pieceRanges, params, nPiece, thresClip, nIterClip);
1876 setSpectrum(res, whichrow);
1877
1878 if (outLogger || outTextFile) {
1879 std::vector<bool> fixed;
1880 fixed.clear();
1881 float rms = getRms(chanMask, whichrow);
1882 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
1883
1884 if (outLogger) {
1885 LogIO ols(LogOrigin("Scantable", "autoCubicSplineBaseline()", WHERE));
1886 ols << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
1887 }
1888 if (outTextFile) {
1889 ofs << formatPiecewiseBaselineParams(pieceRanges, params, fixed, rms, masklist, whichrow, true) << flush;
1890 }
1891 }
1892
1893 }
1894
1895 if (outTextFile) ofs.close();
1896}
1897
1898
1899std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, std::vector<int>& sectionRanges, std::vector<float>& params, int nPiece, float thresClip, int nIterClip) {
1900 if (nPiece < 1) {
1901 throw(AipsError("wrong number of the sections for fitting"));
1902 }
1903 if (data.size() != mask.size()) {
1904 throw(AipsError("data and mask have different size"));
1905 }
1906
1907 int nChan = data.size();
1908 std::vector<int> maskArray;
1909 std::vector<int> x;
1910 for (int i = 0; i < nChan; ++i) {
1911 maskArray.push_back(mask[i] ? 1 : 0);
1912 if (mask[i]) {
1913 x.push_back(i);
1914 }
1915 }
1916
1917 int nData = x.size();
1918 int nElement = (int)(floor(floor((double)(nData/nPiece))+0.5));
1919
1920 std::vector<double> sectionList0, sectionList1, sectionListR;
1921 sectionList0.push_back((double)x[0]);
1922 sectionRanges.clear();
1923 sectionRanges.push_back(x[0]);
1924 for (int i = 1; i < nPiece; ++i) {
1925 double valX = (double)x[nElement*i];
1926 sectionList0.push_back(valX);
1927 sectionList1.push_back(valX);
1928 sectionListR.push_back(1.0/valX);
1929
1930 sectionRanges.push_back(x[nElement*i]-1);
1931 sectionRanges.push_back(x[nElement*i]);
1932 }
1933 sectionList1.push_back((double)(x[x.size()-1]+1));
1934 sectionRanges.push_back(x[x.size()-1]);
1935
1936 std::vector<double> x1, x2, x3, z1, x1z1, x2z1, x3z1, r1;
1937 for (int i = 0; i < nChan; ++i) {
1938 x1.push_back((double)i);
1939 x2.push_back((double)(i*i));
1940 x3.push_back((double)(i*i*i));
1941 z1.push_back((double)data[i]);
1942 x1z1.push_back(((double)i)*(double)data[i]);
1943 x2z1.push_back(((double)(i*i))*(double)data[i]);
1944 x3z1.push_back(((double)(i*i*i))*(double)data[i]);
1945 r1.push_back(0.0);
1946 }
1947
1948 int currentNData = nData;
1949 int nDOF = nPiece + 3; //number of parameters to solve, namely, 4+(nPiece-1).
1950
1951 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
1952 //xMatrix : horizontal concatenation of the least-sq. matrix (left) and an
1953 //identity matrix (right).
1954 //the right part is used to calculate the inverse matrix of the left part.
1955 double xMatrix[nDOF][2*nDOF];
1956 double zMatrix[nDOF];
1957 for (int i = 0; i < nDOF; ++i) {
1958 for (int j = 0; j < 2*nDOF; ++j) {
1959 xMatrix[i][j] = 0.0;
1960 }
1961 xMatrix[i][nDOF+i] = 1.0;
1962 zMatrix[i] = 0.0;
1963 }
1964
1965 for (int n = 0; n < nPiece; ++n) {
1966 for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
1967 if (maskArray[i] == 0) continue;
1968 xMatrix[0][0] += 1.0;
1969 xMatrix[0][1] += x1[i];
1970 xMatrix[0][2] += x2[i];
1971 xMatrix[0][3] += x3[i];
1972 xMatrix[1][1] += x2[i];
1973 xMatrix[1][2] += x3[i];
1974 xMatrix[1][3] += x2[i]*x2[i];
1975 xMatrix[2][2] += x2[i]*x2[i];
1976 xMatrix[2][3] += x3[i]*x2[i];
1977 xMatrix[3][3] += x3[i]*x3[i];
1978 zMatrix[0] += z1[i];
1979 zMatrix[1] += x1z1[i];
1980 zMatrix[2] += x2z1[i];
1981 zMatrix[3] += x3z1[i];
1982 for (int j = 0; j < n; ++j) {
1983 double q = 1.0 - x1[i]*sectionListR[j];
1984 q = q*q*q;
1985 xMatrix[0][j+4] += q;
1986 xMatrix[1][j+4] += q*x1[i];
1987 xMatrix[2][j+4] += q*x2[i];
1988 xMatrix[3][j+4] += q*x3[i];
1989 for (int k = 0; k < j; ++k) {
1990 double r = 1.0 - x1[i]*sectionListR[k];
1991 r = r*r*r;
1992 xMatrix[k+4][j+4] += r*q;
1993 }
1994 xMatrix[j+4][j+4] += q*q;
1995 zMatrix[j+4] += q*z1[i];
1996 }
1997 }
1998 }
1999
2000 for (int i = 0; i < nDOF; ++i) {
2001 for (int j = 0; j < i; ++j) {
2002 xMatrix[i][j] = xMatrix[j][i];
2003 }
2004 }
2005
2006 std::vector<double> invDiag;
2007 for (int i = 0; i < nDOF; ++i) {
2008 invDiag.push_back(1.0/xMatrix[i][i]);
2009 for (int j = 0; j < nDOF; ++j) {
2010 xMatrix[i][j] *= invDiag[i];
2011 }
2012 }
2013
2014 for (int k = 0; k < nDOF; ++k) {
2015 for (int i = 0; i < nDOF; ++i) {
2016 if (i != k) {
2017 double factor1 = xMatrix[k][k];
2018 double factor2 = xMatrix[i][k];
2019 for (int j = k; j < 2*nDOF; ++j) {
2020 xMatrix[i][j] *= factor1;
2021 xMatrix[i][j] -= xMatrix[k][j]*factor2;
2022 xMatrix[i][j] /= factor1;
2023 }
2024 }
2025 }
2026 double xDiag = xMatrix[k][k];
2027 for (int j = k; j < 2*nDOF; ++j) {
2028 xMatrix[k][j] /= xDiag;
2029 }
2030 }
2031
2032 for (int i = 0; i < nDOF; ++i) {
2033 for (int j = 0; j < nDOF; ++j) {
2034 xMatrix[i][nDOF+j] *= invDiag[j];
2035 }
2036 }
2037 //compute a vector y which consists of the coefficients of the best-fit spline curves
2038 //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
2039 //cubic terms for the other pieces (in case nPiece>1).
2040 std::vector<double> y;
2041 for (int i = 0; i < nDOF; ++i) {
2042 y.push_back(0.0);
2043 for (int j = 0; j < nDOF; ++j) {
2044 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
2045 }
2046 }
2047
2048 double a0 = y[0];
2049 double a1 = y[1];
2050 double a2 = y[2];
2051 double a3 = y[3];
2052 params.clear();
2053
2054 for (int n = 0; n < nPiece; ++n) {
2055 for (int i = sectionList0[n]; i < sectionList1[n]; ++i) {
2056 r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
2057 }
2058 params.push_back(a0);
2059 params.push_back(a1);
2060 params.push_back(a2);
2061 params.push_back(a3);
2062
2063 if (n == nPiece-1) break;
2064
2065 double d = y[4+n];
2066 a0 += d;
2067 a1 -= 3.0*d*sectionListR[n];
2068 a2 += 3.0*d*sectionListR[n]*sectionListR[n];
2069 a3 -= d*sectionListR[n]*sectionListR[n]*sectionListR[n];
2070 }
2071
2072 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
2073 break;
2074 } else {
2075 std::vector<double> rz;
2076 double stdDev = 0.0;
2077 for (int i = 0; i < nChan; ++i) {
2078 double val = abs(r1[i] - z1[i]);
2079 rz.push_back(val);
2080 stdDev += val*val*(double)maskArray[i];
2081 }
2082 stdDev = sqrt(stdDev/(double)nData);
2083
2084 double thres = stdDev * thresClip;
2085 int newNData = 0;
2086 for (int i = 0; i < nChan; ++i) {
2087 if (rz[i] >= thres) {
2088 maskArray[i] = 0;
2089 }
2090 if (maskArray[i] > 0) {
2091 newNData++;
2092 }
2093 }
2094 if (newNData == currentNData) {
2095 break; //no additional flag. finish iteration.
2096 } else {
2097 currentNData = newNData;
2098 }
2099 }
2100 }
2101
2102 std::vector<float> residual;
2103 for (int i = 0; i < nChan; ++i) {
2104 residual.push_back((float)(z1[i] - r1[i]));
2105 }
2106 return residual;
2107
2108}
2109
2110void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
2111{
2112 std::vector<double> abcsd = getAbcissa(whichrow);
2113 std::vector<float> abcs;
2114 for (int i = 0; i < abcsd.size(); ++i) {
2115 abcs.push_back((float)abcsd[i]);
2116 }
2117 std::vector<float> spec = getSpectrum(whichrow);
2118
2119 fitter.setData(abcs, spec, mask);
2120 fitter.lfit();
2121}
2122
2123std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
2124{
2125 std::vector<bool> chanMask = getMask(whichrow);
2126 int chanMaskSize = chanMask.size();
2127 if (chanMaskSize != inMask.size()) {
2128 throw(AipsError("different mask sizes"));
2129 }
2130 for (int i = 0; i < chanMaskSize; ++i) {
2131 chanMask[i] = chanMask[i] && inMask[i];
2132 }
2133
2134 return chanMask;
2135}
2136
2137/*
2138std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
2139{
2140 int edgeSize = edge.size();
2141 std::vector<int> currentEdge;
2142 if (edgeSize >= 2) {
2143 int idx = 0;
2144 if (edgeSize > 2) {
2145 if (edgeSize < minEdgeSize) {
2146 throw(AipsError("Length of edge element info is less than that of IFs"));
2147 }
2148 idx = 2 * getIF(whichrow);
2149 }
2150 currentEdge.push_back(edge[idx]);
2151 currentEdge.push_back(edge[idx+1]);
2152 } else {
2153 throw(AipsError("Wrong length of edge element"));
2154 }
2155
2156 lineFinder.setData(getSpectrum(whichrow));
2157 lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
2158
2159 return lineFinder.getMask();
2160}
2161*/
2162
2163void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool outLogger, const std::string& blfile)
2164{
2165 ofstream ofs;
2166 String coordInfo;
2167 bool hasSameNchan;
2168 int firstIF;
2169 bool outTextFile = false;
2170
2171 if (blfile != "") {
2172 ofs.open(blfile.c_str(), ios::out | ios::app);
2173 if (ofs) outTextFile = true;
2174 }
2175
2176 if (outLogger || outTextFile) {
2177 coordInfo = getCoordInfo()[0];
2178 if (coordInfo == "") coordInfo = "channel";
2179 hasSameNchan = hasSameNchanOverIFs();
2180 firstIF = getIF(0);
2181 }
2182
2183 Fitter fitter = Fitter();
2184 fitter.setExpression("poly", order);
2185
2186 int nRow = nrow();
2187 std::vector<bool> chanMask;
2188
2189 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2190
2191 chanMask = getCompositeChanMask(whichrow, mask);
2192 fitBaseline(chanMask, whichrow, fitter);
2193 setSpectrum(fitter.getResidual(), whichrow);
2194
2195 if (outLogger || outTextFile) {
2196 std::vector<float> params = fitter.getParameters();
2197 std::vector<bool> fixed = fitter.getFixedParameters();
2198 float rms = getRms(chanMask, whichrow);
2199 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
2200
2201 if (outLogger) {
2202 LogIO ols(LogOrigin("Scantable", "polyBaseline()", WHERE));
2203 ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
2204 }
2205 if (outTextFile) {
2206 ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
2207 }
2208 }
2209
2210 }
2211
2212 if (outTextFile) ofs.close();
2213}
2214
2215void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool outLogger, const std::string& blfile)
2216{
2217 ofstream ofs;
2218 String coordInfo;
2219 bool hasSameNchan;
2220 int firstIF;
2221 bool outTextFile = false;
2222
2223 if (blfile != "") {
2224 ofs.open(blfile.c_str(), ios::out | ios::app);
2225 if (ofs) outTextFile = true;
2226 }
2227
2228 if (outLogger || outTextFile) {
2229 coordInfo = getCoordInfo()[0];
2230 if (coordInfo == "") coordInfo = "channel";
2231 hasSameNchan = hasSameNchanOverIFs();
2232 firstIF = getIF(0);
2233 }
2234
2235 Fitter fitter = Fitter();
2236 fitter.setExpression("poly", order);
2237
2238 int nRow = nrow();
2239 std::vector<bool> chanMask;
2240 int minEdgeSize = getIFNos().size()*2;
2241 STLineFinder lineFinder = STLineFinder();
2242 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2243
2244 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2245
2246 //-------------------------------------------------------
2247 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2248 //-------------------------------------------------------
2249 int edgeSize = edge.size();
2250 std::vector<int> currentEdge;
2251 if (edgeSize >= 2) {
2252 int idx = 0;
2253 if (edgeSize > 2) {
2254 if (edgeSize < minEdgeSize) {
2255 throw(AipsError("Length of edge element info is less than that of IFs"));
2256 }
2257 idx = 2 * getIF(whichrow);
2258 }
2259 currentEdge.push_back(edge[idx]);
2260 currentEdge.push_back(edge[idx+1]);
2261 } else {
2262 throw(AipsError("Wrong length of edge element"));
2263 }
2264 lineFinder.setData(getSpectrum(whichrow));
2265 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2266 chanMask = lineFinder.getMask();
2267 //-------------------------------------------------------
2268
2269
2270 fitBaseline(chanMask, whichrow, fitter);
2271 setSpectrum(fitter.getResidual(), whichrow);
2272
2273 if (outLogger || outTextFile) {
2274 std::vector<float> params = fitter.getParameters();
2275 std::vector<bool> fixed = fitter.getFixedParameters();
2276 float rms = getRms(chanMask, whichrow);
2277 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan, firstIF, true);
2278
2279 if (outLogger) {
2280 LogIO ols(LogOrigin("Scantable", "autoPolyBaseline()", WHERE));
2281 ols << formatBaselineParams(params, fixed, rms, masklist, whichrow, false) << LogIO::POST ;
2282 }
2283 if (outTextFile) {
2284 ofs << formatBaselineParams(params, fixed, rms, masklist, whichrow, true) << flush;
2285 }
2286 }
2287
2288 }
2289
2290 if (outTextFile) ofs.close();
2291}
2292
2293
2294float Scantable::getRms(const std::vector<bool>& mask, int whichrow) {
2295 Vector<Float> spec;
2296 specCol_.get(whichrow, spec);
2297
2298 float mean = 0.0;
2299 float smean = 0.0;
2300 int n = 0;
2301 for (int i = 0; i < spec.nelements(); ++i) {
2302 if (mask[i]) {
2303 mean += spec[i];
2304 smean += spec[i]*spec[i];
2305 n++;
2306 }
2307 }
2308
2309 mean /= (float)n;
2310 smean /= (float)n;
2311
2312 return sqrt(smean - mean*mean);
2313}
2314
2315
2316std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
2317{
2318 ostringstream oss;
2319
2320 if (verbose) {
2321 for (int i = 0; i < 60; ++i) {
2322 oss << "-";
2323 }
2324 oss << endl;
2325 oss << " Scan[" << getScan(whichrow) << "]";
2326 oss << " Beam[" << getBeam(whichrow) << "]";
2327 oss << " IF[" << getIF(whichrow) << "]";
2328 oss << " Pol[" << getPol(whichrow) << "]";
2329 oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
2330 oss << "Fitter range = " << masklist << endl;
2331 oss << "Baseline parameters" << endl;
2332 oss << flush;
2333 }
2334
2335 return String(oss);
2336}
2337
2338std::string Scantable::formatBaselineParamsFooter(float rms, bool verbose) const
2339{
2340 ostringstream oss;
2341
2342 if (verbose) {
2343 oss << endl;
2344 oss << "Results of baseline fit" << endl;
2345 oss << " rms = " << setprecision(6) << rms << endl;
2346 }
2347
2348 return String(oss);
2349}
2350
2351std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
2352{
2353 ostringstream oss;
2354 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
2355
2356 if ((params.size() > 0) && (ranges.size() > 0)) {
2357 if (((ranges.size() % 2) == 0) && ((params.size() % (ranges.size() / 2)) == 0)) {
2358 int nParam = params.size() / (ranges.size() / 2);
2359 for (uInt j = 0; j < ranges.size(); j+=2) {
2360 oss << "[" << ranges[j] << "," << ranges[j+1] << "]";
2361 for (uInt i = 0; i < nParam; ++i) {
2362 if (i > 0) {
2363 oss << ",";
2364 }
2365 String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
2366 oss << " p" << i << fix << "= " << setprecision(6) << params[j/2*nParam+i];
2367 }
2368 }
2369 } else {
2370 oss << " ";
2371 }
2372 } else {
2373 oss << " Not fitted";
2374 }
2375
2376 oss << formatBaselineParamsFooter(rms, verbose);
2377 oss << flush;
2378
2379 return String(oss);
2380}
2381
2382std::string Scantable::formatBaselineParams(const std::vector<float>& params, const std::vector<bool>& fixed, float rms, const std::string& masklist, int whichrow, bool verbose) const
2383{
2384 ostringstream oss;
2385 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
2386
2387 if (params.size() > 0) {
2388 for (uInt i = 0; i < params.size(); ++i) {
2389 if (i > 0) {
2390 oss << ",";
2391 }
2392 String fix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
2393 oss << " p" << i << fix << "= " << setprecision(6) << params[i];
2394 }
2395 } else {
2396 oss << " Not fitted";
2397 }
2398
2399 oss << formatBaselineParamsFooter(rms, verbose);
2400 oss << flush;
2401
2402 return String(oss);
2403}
2404
2405std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, int firstIF, bool silent)
2406{
2407 if (mask.size() < 2) {
2408 throw(AipsError("The mask elements should be > 1"));
2409 }
2410 if (mask.size() != nchan()) {
2411 throw(AipsError("Number of channels in scantable != number of mask elements"));
2412 }
2413
2414 std::vector<int> startIdx = getMaskEdgeIndices(mask, true);
2415 std::vector<int> endIdx = getMaskEdgeIndices(mask, false);
2416
2417 if (startIdx.size() != endIdx.size()) {
2418 throw(AipsError("Numbers of mask start and mask end are not identical"));
2419 }
2420 for (int i = 0; i < startIdx.size(); ++i) {
2421 if (startIdx[i] > endIdx[i]) {
2422 throw(AipsError("Mask start index > mask end index"));
2423 }
2424 }
2425
2426 if (!silent) {
2427 LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
2428 logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
2429 if (!hasSameNchan) {
2430 logOs << endl << "This mask is only valid for IF=" << firstIF;
2431 }
2432 logOs << LogIO::POST;
2433 }
2434
2435 std::vector<double> abcissa = getAbcissa(whichrow);
2436 ostringstream oss;
2437 oss.setf(ios::fixed);
2438 oss << setprecision(1) << "[";
2439 for (int i = 0; i < startIdx.size(); ++i) {
2440 std::vector<float> aMaskRange;
2441 aMaskRange.push_back((float)abcissa[startIdx[i]]);
2442 aMaskRange.push_back((float)abcissa[endIdx[i]]);
2443 sort(aMaskRange.begin(), aMaskRange.end());
2444 if (i > 0) oss << ",";
2445 oss << "[" << aMaskRange[0] << "," << aMaskRange[1] << "]";
2446 }
2447 oss << "]" << flush;
2448
2449 return String(oss);
2450}
2451
2452bool Scantable::hasSameNchanOverIFs()
2453{
2454 int nIF = nif(-1);
2455 int nCh;
2456 int totalPositiveNChan = 0;
2457 int nPositiveNChan = 0;
2458
2459 for (int i = 0; i < nIF; ++i) {
2460 nCh = nchan(i);
2461 if (nCh > 0) {
2462 totalPositiveNChan += nCh;
2463 nPositiveNChan++;
2464 }
2465 }
2466
2467 return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
2468}
2469
2470std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask, bool getStartIndices)
2471{
2472 if (mask.size() < 2) {
2473 throw(AipsError("The mask elements should be > 1"));
2474 }
2475 if (mask.size() != nchan()) {
2476 throw(AipsError("Number of channels in scantable != number of mask elements"));
2477 }
2478
2479 std::vector<int> out;
2480 int endIdx = mask.size() - 1;
2481
2482 if (getStartIndices) {
2483 if (mask[0]) {
2484 out.push_back(0);
2485 }
2486 for (int i = 0; i < endIdx; ++i) {
2487 if ((!mask[i]) && mask[i+1]) {
2488 out.push_back(i+1);
2489 }
2490 }
2491 } else {
2492 for (int i = 0; i < endIdx; ++i) {
2493 if (mask[i] && (!mask[i+1])) {
2494 out.push_back(i);
2495 }
2496 }
2497 if (mask[endIdx]) {
2498 out.push_back(endIdx);
2499 }
2500 }
2501
2502 return out;
2503}
2504
2505
2506/*
2507STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
2508{
2509 Fitter fitter = Fitter();
2510 fitter.setExpression("poly", order);
2511
2512 std::vector<bool> fmask = getMask(rowno);
2513 if (fmask.size() != mask.size()) {
2514 throw(AipsError("different mask sizes"));
2515 }
2516 for (int i = 0; i < fmask.size(); ++i) {
2517 fmask[i] = fmask[i] && mask[i];
2518 }
2519
2520 fitBaseline(fmask, rowno, fitter);
2521 setSpectrum(fitter.getResidual(), rowno);
2522 return fitter.getFitEntry();
2523}
2524*/
2525
2526}
2527//namespace asap
Note: See TracBrowser for help on using the repository browser.