source: trunk/src/STSideBandSep.cpp@ 2804

Last change on this file since 2804 was 2794, checked in by Kana Sugimoto, 12 years ago

New Development: Yes

JIRA Issue: Yes (CAS-4141/TRAC-290)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sbseparator module

Description: Pushed python codes down to cpp.
File size: 43.4 KB
Line 
1// C++ Interface: STSideBandSep
2//
3// Description:
4// A class to invoke sideband separation of Scantable
5//
6// Author: Kana Sugimoto <kana.sugi@nao.ac.jp>, (C) 2012
7//
8// Copyright: See COPYING file that comes with this distribution
9//
10//
11
12// STL
13#include <ctype.h>
14
15// cascore
16#include <casa/OS/File.h>
17#include <casa/Logging/LogIO.h>
18#include <casa/Quanta/QuantumHolder.h>
19
20#include <measures/Measures/MFrequency.h>
21#include <measures/Measures/MCFrequency.h>
22
23#include <tables/Tables/TableRow.h>
24#include <tables/Tables/TableRecord.h>
25#include <tables/Tables/TableVector.h>
26
27// asap
28#include "STGrid.h"
29#include "STMath.h"
30#include "MathUtils.h"
31#include "STSideBandSep.h"
32
33using namespace std ;
34using namespace casa ;
35using namespace asap ;
36
37// #ifndef KS_DEBUG
38// #define KS_DEBUG
39// #endif
40
41namespace asap {
42
43// constructors
44STSideBandSep::STSideBandSep(const vector<string> &names)
45{
46 LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
47 os << "Setting scantable names to process." << LogIO::POST ;
48 // Set file names
49 ntable_ = names.size();
50 infileList_.resize(ntable_);
51 for (unsigned int i = 0; i < ntable_; i++){
52 if (!checkFile(names[i], "d"))
53 throw( AipsError("File does not exist") );
54 infileList_[i] = names[i];
55 }
56 intabList_.resize(0);
57
58 init();
59
60 {// Summary
61 os << ntable_ << " files are set: [";
62 for (unsigned int i = 0; i < ntable_; i++) {
63 os << " '" << infileList_[i] << "' ";
64 if (i != ntable_-1) os << ",";
65 }
66 os << "] " << LogIO::POST;
67 }
68};
69
70STSideBandSep::STSideBandSep(const vector<ScantableWrapper> &tables)
71{
72 LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
73 os << "Setting list of scantables to process." << LogIO::POST ;
74 // Set file names
75 ntable_ = tables.size();
76 intabList_.resize(ntable_);
77 for (unsigned int i = 0; i < ntable_; i++){
78 intabList_[i] = tables[i].getCP();
79 }
80 infileList_.resize(0);
81
82 init();
83 tp_ = intabList_[0]->table().tableType();
84
85 os << ntable_ << " tables are set." << LogIO::POST;
86};
87
88STSideBandSep::~STSideBandSep()
89{
90#ifdef KS_DEBUG
91 cout << "Destructor ~STSideBandSep()" << endl;
92#endif
93};
94
95void STSideBandSep::init()
96{
97 // frequency setup
98 sigIfno_= -1;
99 ftol_ = -1;
100 solFrame_ = MFrequency::N_Types;
101 // shifts
102 initshift();
103 // direction tolerance
104 xtol_ = ytol_ = 9.69627e-6; // 2arcsec
105 // solution parameters
106 otherside_ = false;
107 doboth_ = false;
108 rejlimit_ = 0.2;
109 // LO1 values
110 lo1Freq_ = -1;
111 loTime_ = -1;
112 loDir_ = "";
113 // Default LO frame is TOPO
114 loFrame_ = MFrequency::TOPO;
115 // scantable storage
116 tp_ = Table::Memory;
117};
118
119void STSideBandSep::initshift()
120{
121 // shifts
122 nshift_ = 0;
123 nchan_ = 0;
124 sigShift_.resize(0);
125 imgShift_.resize(0);
126 tableList_.resize(0);
127};
128
129void STSideBandSep::setFrequency(const unsigned int ifno,
130 const string freqtol,
131 const string frame)
132{
133 LogIO os(LogOrigin("STSideBandSep","setFrequency()", WHERE));
134
135 initshift();
136
137 // IFNO
138 sigIfno_ = (int) ifno;
139
140 // Frequency tolerance
141 Quantum<Double> qftol;
142 readQuantity(qftol, String(freqtol));
143 if (!qftol.getUnit().empty()){
144 // make sure the quantity is frequency
145 if (qftol.getFullUnit().getValue() != Unit("Hz").getValue())
146 throw( AipsError("Invalid quantity for frequency tolerance.") );
147 qftol.convert("Hz");
148 }
149 ftol_ = qftol;
150
151 // Frequency Frame
152 if (!frame.empty()){
153 MFrequency::Types mft;
154 if (!MFrequency::getType(mft, frame))
155 throw( AipsError("Invalid frame type.") );
156 solFrame_ = mft;
157 } else {
158 solFrame_ = MFrequency::N_Types;
159 }
160
161 {// Summary
162 const String sframe = ( (solFrame_ == MFrequency::N_Types) ?
163 "table frame" :
164 MFrequency::showType(solFrame_) );
165 os << "Frequency setup to search IF group: "
166 << "IFNO of table[0] = " << sigIfno_
167 << " , Freq tolerance = " << ftol_.getValue() << " [ "
168 << (ftol_.getUnit().empty() ? "channel" : ftol_.getUnit() )
169 << " ] (in " << sframe <<")" << LogIO::POST;
170 }
171};
172
173
174void STSideBandSep::setDirTolerance(const vector<string> dirtol)
175{
176 LogIO os(LogOrigin("STSideBandSep","setDirTolerance()", WHERE));
177 Quantum<Double> qcell;
178 if ( (dirtol.size() == 1) && !dirtol[0].empty() ) {
179 readQuantity(qcell, String(dirtol[0]));
180 if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
181 xtol_ = ytol_ = qcell.getValue("rad");
182 else
183 throw( AipsError("Invalid unit for direction tolerance.") );
184 }
185 else if (dirtol.size() > 1) {
186 if ( dirtol[0].empty() && dirtol[1].empty() )
187 throw( AipsError("Direction tolerance is empty.") );
188 if ( !dirtol[0].empty() ) {
189 readQuantity(qcell, String(dirtol[0]));
190 if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
191 xtol_ = qcell.getValue("rad");
192 else
193 throw( AipsError("Invalid unit for direction tolerance.") );
194 }
195 if ( !dirtol[1].empty() ) {
196 readQuantity(qcell, String(dirtol[1]));
197 if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
198 ytol_ = qcell.getValue("rad");
199 else
200 throw( AipsError("Invalid unit for direction tolerance.") );
201 }
202 else {
203 ytol_ = xtol_;
204 }
205 }
206 else throw( AipsError("Invalid direction tolerance.") );
207
208 os << "Direction tolerance: ( "
209 << xtol_ << " , " << ytol_ << " ) [rad]" << LogIO::POST;
210};
211
212void STSideBandSep::setShift(const vector<double> &shift)
213{
214 LogIO os(LogOrigin("STSideBandSep","setShift()", WHERE));
215 imgShift_.resize(shift.size());
216 for (unsigned int i = 0; i < shift.size(); i++)
217 imgShift_[i] = shift[i];
218
219 if (imgShift_.size() == 0) {
220 os << "Channel shifts are cleared." << LogIO::POST;
221 } else {
222 os << "Channel shifts of image sideband are set: ( ";
223 for (unsigned int i = 0; i < imgShift_.size(); i++) {
224 os << imgShift_[i];
225 if (i != imgShift_.size()-1) os << " , ";
226 }
227 os << " ) [channels]" << LogIO::POST;
228 }
229};
230
231void STSideBandSep::setThreshold(const double limit)
232{
233 LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE));
234 if (limit < 0)
235 throw( AipsError("Rejection limit should be a positive number.") );
236
237 rejlimit_ = limit;
238 os << "Rejection limit is set to " << rejlimit_ << LogIO::POST;
239};
240
241void STSideBandSep::separate(string outname)
242{
243#ifdef KS_DEBUG
244 cout << "STSideBandSep::separate" << endl;
245#endif
246 LogIO os(LogOrigin("STSideBandSep","separate()", WHERE));
247 if (outname.empty())
248 outname = "sbseparated.asap";
249
250 // Set up a goup of IFNOs in the list of scantables within
251 // the frequency tolerance and make them a list.
252 nshift_ = setupShift();
253 if (nshift_ < 2)
254 throw( AipsError("At least 2 IFs are necessary for convolution.") );
255 // Grid scantable and generate output tables
256 ScantableWrapper gridst = gridTable();
257 sigTab_p = gridst.getCP();
258 if (doboth_)
259 imgTab_p = gridst.copy().getCP();
260 vector<unsigned int> remRowIds;
261 remRowIds.resize(0);
262 Matrix<float> specMat(nchan_, nshift_);
263 vector<float> sigSpec(nchan_), imgSpec(nchan_);
264 vector<uInt> tabIdvec;
265
266 //Generate FFTServer
267 fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX);
268 fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL);
269
270 /// Loop over sigTab_p and separate sideband
271 for (int irow = 0; irow < sigTab_p->nrow(); irow++){
272 tabIdvec.resize(0);
273 const int polId = sigTab_p->getPol(irow);
274 const int beamId = sigTab_p->getBeam(irow);
275 const vector<double> dir = sigTab_p->getDirectionVector(irow);
276 // Get a set of spectra to solve
277 if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
278 specMat, tabIdvec)){
279 remRowIds.push_back(irow);
280#ifdef KS_DEBUG
281 cout << "no matching row found. skipping row = " << irow << endl;
282#endif
283 continue;
284 }
285 // Solve signal sideband
286 sigSpec = solve(specMat, tabIdvec, true);
287 sigTab_p->setSpectrum(sigSpec, irow);
288
289 // Solve image sideband
290 if (doboth_) {
291 imgSpec = solve(specMat, tabIdvec, false);
292 imgTab_p->setSpectrum(imgSpec, irow);
293 }
294 } // end of row loop
295
296 // Remove or flag rows without relevant data from gridded tables
297 if (remRowIds.size() > 0) {
298 const size_t nrem = remRowIds.size();
299 if ( sigTab_p->table().canRemoveRow() ) {
300 sigTab_p->table().removeRow(remRowIds);
301 os << "Removing " << nrem << " rows from the signal band table"
302 << LogIO::POST;
303 } else {
304 sigTab_p->flagRow(remRowIds, false);
305 os << "Cannot remove rows from the signal band table. Flagging "
306 << nrem << " rows" << LogIO::POST;
307 }
308
309 if (doboth_) {
310 if ( imgTab_p->table().canRemoveRow() ) {
311 imgTab_p->table().removeRow(remRowIds);
312 os << "Removing " << nrem << " rows from the image band table"
313 << LogIO::POST;
314 } else {
315 imgTab_p->flagRow(remRowIds, false);
316 os << "Cannot remove rows from the image band table. Flagging "
317 << nrem << " rows" << LogIO::POST;
318 }
319 }
320 }
321
322 // Finally, save tables on disk
323 if (outname.size() ==0)
324 outname = "sbseparated.asap";
325 const string sigName = outname + ".signalband";
326 os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST;
327 sigTab_p->makePersistent(sigName);
328 if (doboth_) {
329 solveImageFrequency();
330 const string imgName = outname + ".imageband";
331 os << "Saving IMAGE sideband table: " << sigName << LogIO::POST;
332 imgTab_p->makePersistent(imgName);
333 }
334
335};
336
337unsigned int STSideBandSep::setupShift()
338{
339 LogIO os(LogOrigin("STSideBandSep","setupShift()", WHERE));
340 if (infileList_.size() == 0 && intabList_.size() == 0)
341 throw( AipsError("No scantable has been set. Set a list of scantables first.") );
342
343 const bool byname = (intabList_.size() == 0);
344 // Make sure sigIfno_ exists in the first table.
345 CountedPtr<Scantable> stab;
346 vector<string> coordsav;
347 vector<string> coordtmp(3);
348 os << "Checking IFNO in the first table." << LogIO::POST;
349 if (byname) {
350 if (!checkFile(infileList_[0], "d"))
351 os << LogIO::SEVERE << "Could not find scantable '" << infileList_[0]
352 << "'" << LogIO::EXCEPTION;
353 stab = CountedPtr<Scantable>(new Scantable(infileList_[0]));
354 } else {
355 stab = intabList_[0];
356 }
357 if (sigIfno_ < 0) {
358 sigIfno_ = (int) stab->getIF(0);
359 os << "IFNO to process has not been set. Using the first IF = "
360 << sigIfno_ << LogIO::POST;
361 }
362
363 unsigned int basench;
364 double basech0, baseinc, ftolval, inctolval;
365 coordsav = stab->getCoordInfo();
366 const string stfframe = coordsav[1];
367 coordtmp[0] = "Hz";
368 coordtmp[1] = ( (solFrame_ == MFrequency::N_Types) ?
369 stfframe :
370 MFrequency::showType(solFrame_) );
371 coordtmp[2] = coordsav[2];
372 stab->setCoordInfo(coordtmp);
373 if (!getFreqInfo(stab, (unsigned int) sigIfno_, basech0, baseinc, basench)) {
374 os << LogIO::SEVERE << "No data with IFNO=" << sigIfno_
375 << " found in the first table." << LogIO::EXCEPTION;
376 }
377 else {
378 os << "Found IFNO = " << sigIfno_
379 << " in the first table." << LogIO::POST;
380 }
381 if (ftol_.getUnit().empty()) {
382 // tolerance in unit of channels
383 ftolval = ftol_.getValue() * baseinc;
384 }
385 else {
386 ftolval = ftol_.getValue("Hz");
387 }
388 inctolval = abs(baseinc/(double) basench);
389 const string poltype0 = stab->getPolType();
390
391 // Initialize shift values
392 initshift();
393
394 const bool setImg = ( doboth_ && (imgShift_.size() == 0) );
395 // Select IFs
396 for (unsigned int itab = 0; itab < ntable_; itab++ ){
397 os << "Table " << itab << LogIO::POST;
398 if (itab > 0) {
399 if (byname) {
400 if (!checkFile(infileList_[itab], "d"))
401 os << LogIO::SEVERE << "Could not find scantable '"
402 << infileList_[itab] << "'" << LogIO::EXCEPTION;
403 stab = CountedPtr<Scantable>(new Scantable(infileList_[itab]));
404 } else {
405 stab = intabList_[itab];
406 }
407 //POLTYPE should be the same.
408 if (stab->getPolType() != poltype0 ) {
409 os << LogIO::WARN << "POLTYPE differs from the first table."
410 << " Skipping the table" << LogIO::POST;
411 continue;
412 }
413 // Multi beam data may not handled properly
414 if (stab->nbeam() > 1)
415 os << LogIO::WARN << "Table contains multiple beams. "
416 << "It may not be handled properly." << LogIO::POST;
417
418 coordsav = stab->getCoordInfo();
419 coordtmp[2] = coordsav[2];
420 stab->setCoordInfo(coordtmp);
421 }
422 bool selected = false;
423 vector<uint> ifnos = stab->getIFNos();
424 vector<uint>::iterator iter;
425 const STSelector& basesel = stab->getSelection();
426 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
427 unsigned int nch;
428 double freq0, incr;
429 if ( getFreqInfo(stab, *iter, freq0, incr, nch) ){
430 if ( (nch == basench) && (abs(freq0-basech0) < ftolval)
431 && (abs(incr-baseinc) < inctolval) ){
432 //Found
433 STSelector sel(basesel);
434 sel.setIFs(vector<int>(1,(int) *iter));
435 stab->setSelection(sel);
436 CountedPtr<Scantable> seltab = ( new Scantable(*stab, false) );
437 stab->setSelection(basesel);
438 seltab->setCoordInfo(coordsav);
439 const double chShift = (freq0 - basech0) / baseinc;
440 tableList_.push_back(seltab);
441 sigShift_.push_back(chShift);
442 if (setImg)
443 imgShift_.push_back(-chShift);
444
445 selected = true;
446 os << "- IF" << *iter << " selected: sideband shift = "
447 << chShift << " channels" << LogIO::POST;
448 }
449 }
450 } // ifno loop
451 stab->setCoordInfo(coordsav);
452 if (!selected)
453 os << LogIO::WARN << "No data selected in Table "
454 << itab << LogIO::POST;
455 } // table loop
456 nchan_ = basench;
457
458 os << "Total number of IFs selected = " << tableList_.size()
459 << LogIO::POST;
460 if ( setImg && (sigShift_.size() != imgShift_.size()) ){
461 os << LogIO::SEVERE
462 << "User defined channel shift of image sideband has "
463 << imgShift_.size()
464 << " elements, while selected IFNOs are " << sigShift_.size()
465 << "\nThe frequency tolerance (freqtol) may be too small."
466 << LogIO::EXCEPTION;
467 }
468
469 return tableList_.size();
470};
471
472bool STSideBandSep::getFreqInfo(const CountedPtr<Scantable> &stab,
473 const unsigned int &ifno,
474 double &freq0, double &incr,
475 unsigned int &nchan)
476{
477 vector<uint> ifnos = stab->getIFNos();
478 bool found = false;
479 vector<uint>::iterator iter;
480 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
481 if (*iter == ifno) {
482 found = true;
483 break;
484 }
485 }
486 if (!found)
487 return false;
488
489 const STSelector& basesel = stab->getSelection();
490 STSelector sel(basesel);
491 sel.setIFs(vector<int>(1,(int) ifno));
492 stab->setSelection(sel);
493 vector<double> freqs;
494 freqs = stab->getAbcissa(0);
495 freq0 = freqs[0];
496 incr = freqs[1] - freqs[0];
497 nchan = freqs.size();
498 stab->setSelection(basesel);
499 return true;
500};
501
502ScantableWrapper STSideBandSep::gridTable()
503{
504 LogIO os(LogOrigin("STSideBandSep","gridTable()", WHERE));
505 if (tableList_.size() == 0)
506 throw( AipsError("Internal error. No scantable has been set to grid.") );
507 Double xmin, xmax, ymin, ymax;
508 mapExtent(tableList_, xmin, xmax, ymin, ymax);
509 const Double centx = 0.5 * (xmin + xmax);
510 const Double centy = 0.5 * (ymin + ymax);
511 const int nx = max(1, (int) ceil( (xmax - xmin) * cos(centy) /xtol_ ) );
512 const int ny = max(1, (int) ceil( (ymax - ymin) / ytol_ ) );
513
514 string scellx, scelly;
515 {
516 ostringstream oss;
517 oss << xtol_ << "rad" ;
518 scellx = oss.str();
519 }
520 {
521 ostringstream oss;
522 oss << ytol_ << "rad" ;
523 scelly = oss.str();
524 }
525
526 ScantableWrapper stab0;
527 if (intabList_.size() > 0)
528 stab0 = ScantableWrapper(intabList_[0]);
529 else
530 stab0 = ScantableWrapper(infileList_[0]);
531
532 string scenter;
533 {
534 ostringstream oss;
535 oss << stab0.getCP()->getDirectionRefString() << " "
536 << centx << "rad" << " " << centy << "rad";
537 scenter = oss.str();
538 }
539
540 STGrid2 gridder = STGrid2(stab0);
541 gridder.setIF(sigIfno_);
542 gridder.defineImage(nx, ny, scellx, scelly, scenter);
543 gridder.setFunc("box", 1);
544 gridder.setWeight("uniform");
545#ifdef KS_DEBUG
546 cout << "Grid parameter summary: " << endl;
547 cout << "- IF = " << sigIfno_ << endl;
548 cout << "- center = " << scenter << ")\n"
549 << "- npix = (" << nx << ", " << ny << ")\n"
550 << "- cell = (" << scellx << ", " << scelly << endl;
551#endif
552 gridder.grid();
553 const int itp = (tp_ == Table::Memory ? 0 : 1);
554 ScantableWrapper gtab = gridder.getResultAsScantable(itp);
555 return gtab;
556};
557
558void STSideBandSep::mapExtent(vector< CountedPtr<Scantable> > &tablist,
559 Double &xmin, Double &xmax,
560 Double &ymin, Double &ymax)
561{
562 ROArrayColumn<Double> dirCol_;
563 dirCol_.attach( tablist[0]->table(), "DIRECTION" );
564 Matrix<Double> direction = dirCol_.getColumn();
565 Vector<Double> ra( direction.row(0) );
566 mathutil::rotateRA(ra);
567 minMax( xmin, xmax, ra );
568 minMax( ymin, ymax, direction.row(1) );
569 Double amin, amax, bmin, bmax;
570 const uInt ntab = tablist.size();
571 for ( uInt i = 1 ; i < ntab ; i++ ) {
572 dirCol_.attach( tablist[i]->table(), "DIRECTION" );
573 direction.assign( dirCol_.getColumn() );
574 ra.assign( direction.row(0) );
575 mathutil::rotateRA(ra);
576 minMax( amin, amax, ra );
577 minMax( bmin, bmax, direction.row(1) );
578 xmin = min(xmin, amin);
579 xmax = max(xmax, amax);
580 ymin = min(ymin, bmin);
581 ymax = max(ymax, bmax);
582 }
583};
584
585bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
586 const double dirX, const double dirY,
587 Matrix<float> &specMat, vector<uInt> &tabIdvec)
588{
589 LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
590
591 tabIdvec.resize(0);
592 specMat.resize(nchan_, nshift_);
593 Vector<float> spec;
594 uInt nspec = 0;
595 STMath stm(false); // insitu has no effect for average.
596 for (uInt itab = 0 ; itab < nshift_ ; itab++) {
597 CountedPtr<Scantable> currtab_p = tableList_[itab];
598 // Selection by POLNO and BEAMNO
599 const STSelector& basesel = currtab_p->getSelection();
600 STSelector sel(basesel);
601 sel.setPolarizations(vector<int>(1, polId));
602 sel.setBeams(vector<int>(1, beamId));
603 try {
604 currtab_p->setSelection(sel);
605 } catch (...) {
606#ifdef KS_DEBUG
607 cout << "Table " << itab << " - No spectrum found. skipping the table."
608 << endl;
609#endif
610 continue;
611 }
612 // Selection by direction;
613 vector<int> selrow(0);
614 vector<double> currDir(2, 0.);
615 const int nselrow = currtab_p->nrow();
616 for (int irow = 0 ; irow < nselrow ; irow++) {
617 currDir = currtab_p->getDirectionVector(irow);
618 if ( (abs(currDir[0]-dirX) > xtol_) ||
619 (abs(currDir[1]-dirY) > ytol_) )
620 continue;
621 // within direction tolerance
622 selrow.push_back(irow);
623 } // end of row loop
624
625 if (selrow.size() < 1) {
626 currtab_p->setSelection(basesel);
627
628#ifdef KS_DEBUG
629 cout << "Table " << itab << " - No spectrum found. skipping the table."
630 << endl;
631#endif
632
633 continue;
634 }
635
636 // At least a spectrum is selected in this table
637 CountedPtr<Scantable> seltab_p = ( new Scantable(*currtab_p, false) );
638 currtab_p->setSelection(basesel);
639 STSelector sel2(seltab_p->getSelection());
640 sel2.setRows(selrow);
641 seltab_p->setSelection(sel2);
642 CountedPtr<Scantable> avetab_p;
643 vector<bool> mask;
644 if (seltab_p->nrow() > 1) {
645 avetab_p = stm.average(vector< CountedPtr<Scantable> >(1, seltab_p),
646 vector<bool>(), "TINTSYS", "NONE");
647#ifdef KS_DEBUG
648 cout << "Table " << itab << " - more than a spectrum is selected. averaging rows..."
649 << endl;
650#endif
651 if (avetab_p->nrow() > 1)
652 throw( AipsError("Averaged table has more than a row. Somethigs went wrong.") );
653 } else {
654 avetab_p = seltab_p;
655 }
656 spec.reference(specMat.column(nspec));
657 spec = avetab_p->getSpectrum(0);
658 tabIdvec.push_back((uInt) itab);
659 nspec++;
660 } // end of table loop
661 if (nspec != nshift_){
662 //shiftSpecmat.resize(nchan_, nspec, true);
663#ifdef KS_DEBUG
664 cout << "Could not find corresponding rows in some tables."
665 << endl;
666 cout << "Number of spectra selected = " << nspec << endl;
667#endif
668 }
669 if (nspec < 2) {
670#ifdef KS_DEBUG
671 cout << "At least 2 spectra are necessary for convolution"
672 << endl;
673#endif
674 return false;
675 }
676 return true;
677};
678
679vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
680 const vector<uInt> &tabIdvec,
681 const bool signal)
682{
683 LogIO os(LogOrigin("STSideBandSep","solve()", WHERE));
684 if (tabIdvec.size() == 0)
685 throw(AipsError("Internal error. Table index is not defined."));
686 if (specmat.ncolumn() != tabIdvec.size())
687 throw(AipsError("Internal error. The row number of input matrix is not conformant."));
688 if (specmat.nrow() != nchan_)
689 throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
690
691
692#ifdef KS_DEBUG
693 cout << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband."
694 << endl;
695#endif
696
697 const size_t nspec = tabIdvec.size();
698 vector<double> *thisShift, *otherShift;
699 if (signal == otherside_) {
700 // (solve signal && solveother = T) OR (solve image && solveother = F)
701 thisShift = &imgShift_;
702 otherShift = &sigShift_;
703#ifdef KS_DEBUG
704 cout << "Image sideband will be deconvolved." << endl;
705#endif
706 } else {
707 // (solve signal && solveother = F) OR (solve image && solveother = T)
708 thisShift = &sigShift_;
709 otherShift = &imgShift_;
710#ifdef KS_DEBUG
711 cout << "Signal sideband will be deconvolved." << endl;
712#endif
713 }
714
715 vector<double> spshift(nspec);
716 Matrix<float> shiftSpecmat(nchan_, nspec, 0.);
717 double tempshift;
718 Vector<float> shiftspvec;
719 for (uInt i = 0 ; i < nspec; i++) {
720 spshift[i] = otherShift->at(i) - thisShift->at(i);
721 tempshift = - thisShift->at(i);
722 shiftspvec.reference(shiftSpecmat.column(i));
723 shiftSpectrum(specmat.column(i), tempshift, shiftspvec);
724 }
725
726 Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.);
727 vector<float> thisvec(nchan_, 0.);
728
729 float minval, maxval;
730 minMax(minval, maxval, shiftSpecmat);
731#ifdef KS_DEBUG
732 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
733#endif
734
735#ifdef KS_DEBUG
736 cout << "starting deconvolution" << endl;
737#endif
738 deconvolve(shiftSpecmat, spshift, rejlimit_, convmat);
739#ifdef KS_DEBUG
740 cout << "finished deconvolution" << endl;
741#endif
742
743 minMax(minval, maxval, convmat);
744#ifdef KS_DEBUG
745 cout << "Max/Min of output Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
746#endif
747
748 aggregateMat(convmat, thisvec);
749
750 if (!otherside_) return thisvec;
751
752 // subtract from the other side band.
753 vector<float> othervec(nchan_, 0.);
754 subtractFromOther(shiftSpecmat, thisvec, spshift, othervec);
755 return othervec;
756};
757
758
759void STSideBandSep::shiftSpectrum(const Vector<float> &invec,
760 double shift,
761 Vector<float> &outvec)
762{
763 LogIO os(LogOrigin("STSideBandSep","shiftSpectrum()", WHERE));
764 if (invec.size() != nchan_)
765 throw(AipsError("Internal error. The length of input vector differs from nchan_"));
766 if (outvec.size() != nchan_)
767 throw(AipsError("Internal error. The length of output vector differs from nchan_"));
768
769#ifdef KS_DEBUG
770 cout << "Start shifting spectrum for " << shift << "channels" << endl;
771#endif
772
773 // tweak shift to be in 0 ~ nchan_-1
774 if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
775 if (shift < 0.) shift += nchan_;
776 double rweight = fmod(shift, 1.);
777 if (rweight < 0.) rweight += 1.;
778 double lweight = 1. - rweight;
779 uInt lchan, rchan;
780
781 outvec = 0;
782 for (uInt i = 0 ; i < nchan_ ; i++) {
783 lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
784 if (lchan < 0.) lchan += nchan_;
785 rchan = ( (lchan + 1) % nchan_ );
786 outvec(lchan) += invec(i) * lweight;
787 outvec(rchan) += invec(i) * rweight;
788 }
789};
790
791void STSideBandSep::deconvolve(Matrix<float> &specmat,
792 const vector<double> shiftvec,
793 const double threshold,
794 Matrix<float> &outmat)
795{
796 LogIO os(LogOrigin("STSideBandSep","deconvolve()", WHERE));
797 if (specmat.nrow() != nchan_)
798 throw(AipsError("Internal error. The length of input matrix differs from nchan_"));
799 if (specmat.ncolumn() != shiftvec.size())
800 throw(AipsError("Internal error. The number of input shifts and spectrum differs."));
801
802 float minval, maxval;
803#ifdef KS_DEBUG
804 minMax(minval, maxval, specmat);
805 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
806#endif
807
808 uInt ninsp = shiftvec.size();
809 outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.);
810 Matrix<Complex> fftspmat(nchan_/2+1, ninsp, 0.);
811 Vector<float> rvecref(nchan_, 0.);
812 Vector<Complex> cvecref(nchan_/2+1, 0.);
813 uInt icol = 0;
814 unsigned int nreject = 0;
815
816#ifdef KS_DEBUG
817 cout << "Starting initial FFT. The number of input spectra = " << ninsp << endl;
818 cout << "out matrix has ncolumn = " << outmat.ncolumn() << endl;
819#endif
820
821 for (uInt isp = 0 ; isp < ninsp ; isp++) {
822 rvecref.reference( specmat.column(isp) );
823 cvecref.reference( fftspmat.column(isp) );
824
825#ifdef KS_DEBUG
826 minMax(minval, maxval, rvecref);
827 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
828#endif
829
830 fftsf.fft0(cvecref, rvecref, true);
831
832#ifdef KS_DEBUG
833 double maxr=cvecref[0].real(), minr=cvecref[0].real(),
834 maxi=cvecref[0].imag(), mini=cvecref[0].imag();
835 for (uInt i = 1; i<cvecref.size();i++){
836 maxr = max(maxr, cvecref[i].real());
837 maxi = max(maxi, cvecref[i].imag());
838 minr = min(minr, cvecref[i].real());
839 mini = min(mini, cvecref[i].imag());
840 }
841 cout << "Max/Min of inv FFTed Matrix (size=" << cvecref.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
842#endif
843 }
844
845 //Liberate from reference
846 rvecref.unique();
847
848 Vector<Complex> cspec(nchan_/2+1, 0.);
849 const double PI = 6.0 * asin(0.5);
850 const double nchani = 1. / (float) nchan_;
851 const Complex trans(0., 1.);
852#ifdef KS_DEBUG
853 cout << "starting actual deconvolution" << endl;
854#endif
855 for (uInt j = 0 ; j < ninsp ; j++) {
856 for (uInt k = j+1 ; k < ninsp ; k++) {
857 const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani;
858
859#ifdef KS_DEBUG
860 cout << "icol = " << icol << endl;
861#endif
862
863 for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){
864 cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5;
865 double phase = dx*ichan;
866 if ( fabs( sin(phase) ) > threshold){
867 cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5
868 * trans * sin(phase) / ( 1. - cos(phase) );
869 } else {
870 nreject++;
871 }
872 } // end of channel loop
873
874#ifdef KS_DEBUG
875 cout << "done calculation of cspec" << endl;
876#endif
877
878 Vector<Float> rspec;
879 rspec.reference( outmat.column(icol) );
880
881#ifdef KS_DEBUG
882 cout << "Starting inverse FFT. icol = " << icol << endl;
883 //cout << "- size of complex vector = " << cspec.size() << endl;
884 double maxr=cspec[0].real(), minr=cspec[0].real(),
885 maxi=cspec[0].imag(), mini=cspec[0].imag();
886 for (uInt i = 1; i<cspec.size();i++){
887 maxr = max(maxr, cspec[i].real());
888 maxi = max(maxi, cspec[i].imag());
889 minr = min(minr, cspec[i].real());
890 mini = min(mini, cspec[i].imag());
891 }
892 cout << "Max/Min of conv vector (size=" << cspec.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
893#endif
894
895 fftsi.fft0(rspec, cspec, false);
896
897#ifdef KS_DEBUG
898 //cout << "- size of inversed real vector = " << rspec.size() << endl;
899 minMax(minval, maxval, rspec);
900 cout << "Max/Min of inv FFTed Vector (size=" << rspec.size() << ") = (max: " << maxval << ", min: " << minval << ")" << endl;
901 //cout << "Done inverse FFT. icol = " << icol << endl;
902#endif
903
904 icol++;
905 }
906 }
907
908#ifdef KS_DEBUG
909 minMax(minval, maxval, outmat);
910 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
911#endif
912
913 os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl;
914};
915
916////////////////////////////////////////////////////////////////////
917// void STSideBandSep::cpprfft(std::vector<float> invec)
918// {
919// cout << "c++ method cpprfft" << endl;
920// const unsigned int len = invec.size();
921// Vector<Complex> carr(len/2+1, 0.);
922// Vector<float> inarr = Vector<float>(invec);
923// Vector <float> outarr(len, 0.);
924// FFTServer<Float, Complex> fftsf, fftsi;
925// fftsf.resize(IPosition(1, len), FFTEnums::REALTOCOMPLEX);
926// fftsi.resize(IPosition(1, invec.size()), FFTEnums::COMPLEXTOREAL);
927// cout << "Input float array (length = " << len << "):" << endl;
928// for (uInt i = 0 ; i < len ; i++){
929// cout << (i == 0 ? "( " : " ") << inarr[i] << (i == len-1 ? ")" : ",");
930// }
931// cout << endl;
932// cout << "R->C transform" << endl;
933// fftsf.fft0(carr, inarr, true);
934// cout << "FFTed complex array (length = " << carr.size() << "):" << endl;
935// for (uInt i = 0 ; i < carr.size() ; i++){
936// cout << (i == 0 ? "( " : " ") << carr[i] << ( (i == carr.size()-1) ? ")" : "," );
937// }
938// cout << endl;
939// cout << "C->R transform" << endl;
940// fftsi.fft0(outarr, carr, false);
941// cout << "invFFTed float array (length = " << outarr.size() << "):" << endl;
942// for (uInt i = 0 ; i < outarr.size() ; i++){
943// cout << (i == 0 ? "( " : " ") << outarr[i] << ( (i == outarr.size()-1) ? ")" : "," );
944// }
945// cout << endl;
946// };
947////////////////////////////////////////////////////////////////////
948
949
950void STSideBandSep::aggregateMat(Matrix<float> &inmat,
951 vector<float> &outvec)
952{
953 LogIO os(LogOrigin("STSideBandSep","aggregateMat()", WHERE));
954 if (inmat.nrow() != nchan_)
955 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
956// if (outvec.size() != nchan_)
957// throw(AipsError("Internal error. The size of output vector should be equal to nchan_"));
958
959 os << "Averaging " << inmat.ncolumn() << " spectra in the input matrix."
960 << LogIO::POST;
961
962 const uInt nspec = inmat.ncolumn();
963 const double scale = 1./( (double) nspec );
964 // initialize values with 0
965 outvec.assign(nchan_, 0);
966 for (uInt isp = 0 ; isp < nspec ; isp++) {
967 for (uInt ich = 0 ; ich < nchan_ ; ich++) {
968 outvec[ich] += inmat(ich, isp);
969 }
970 }
971
972 vector<float>::iterator iter;
973 for (iter = outvec.begin(); iter != outvec.end(); iter++){
974 *iter *= scale;
975 }
976};
977
978void STSideBandSep::subtractFromOther(const Matrix<float> &shiftmat,
979 const vector<float> &invec,
980 const vector<double> &shift,
981 vector<float> &outvec)
982{
983 LogIO os(LogOrigin("STSideBandSep","subtractFromOther()", WHERE));
984 if (shiftmat.nrow() != nchan_)
985 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
986 if (invec.size() != nchan_)
987 throw(AipsError("Internal error. The length of input vector should be nchan_"));
988 if (shift.size() != shiftmat.ncolumn())
989 throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift"));
990
991 const uInt nspec = shiftmat.ncolumn();
992 Vector<float> subsp(nchan_, 0.), shiftsub;
993 Matrix<float> submat(nchan_, nspec, 0.);
994 vector<float>::iterator iter;
995 for (uInt isp = 0 ; isp < nspec ; isp++) {
996 for (uInt ich = 0; ich < nchan_ ; ich++) {
997 subsp(ich) = shiftmat(ich, isp) - invec[ich];
998 }
999 shiftsub.reference(submat.column(isp));
1000 shiftSpectrum(subsp, shift[isp], shiftsub);
1001 }
1002
1003 aggregateMat(submat, outvec);
1004};
1005
1006
1007void STSideBandSep::setLO1(const string lo1, const string frame,
1008 const double reftime, const string refdir)
1009{
1010 Quantum<Double> qfreq;
1011 readQuantity(qfreq, String(lo1));
1012 lo1Freq_ = qfreq.getValue("Hz");
1013 MFrequency::getType(loFrame_, frame);
1014 loTime_ = reftime;
1015 loDir_ = refdir;
1016
1017#ifdef KS_DEBUG
1018 cout << "STSideBandSep::setLO1" << endl;
1019 if (lo1Freq_ > 0.)
1020 cout << "lo1 = " << lo1Freq_ << " [Hz] (" << frame << ")" << endl;
1021 if (loTime_ > 0.)
1022 cout << "ref time = " << loTime_ << " [day]" << endl;
1023 if (!loDir_.empty())
1024 cout << "ref direction = " << loDir_ << " [day]" << endl;
1025#endif
1026};
1027
1028void STSideBandSep::setLO1Root(string name)
1029{
1030 LogIO os(LogOrigin("STSideBandSep","setLO1Root()", WHERE));
1031 os << "Searching for '" << name << "'..." << LogIO::POST;
1032 // Check for existance of the file
1033 if (!checkFile(name)) {
1034 throw(AipsError("File does not exist"));
1035 }
1036 if (name[(name.size()-1)] == '/')
1037 name = name.substr(0,(name.size()-2));
1038
1039 if (checkFile(name+"/Receiver.xml", "file") &&
1040 checkFile(name+"/SpectralWindow.xml", "file")){
1041 os << "Found '" << name << "/Receiver.xml' ... got an ASDM name." << LogIO::POST;
1042 asdmName_ = name;
1043 } else if (checkFile(name+"/ASDM_RECEIVER") &&
1044 checkFile(name+"/ASDM_SPECTRALWINDOW")){
1045 os << "Found '" << name << "/ASDM_RECEIVER' ... got a Table name." << LogIO::POST;
1046 asisName_ = name;
1047 } else {
1048 throw(AipsError("Invalid file name. Set an MS or ASDM name."));
1049 }
1050
1051#ifdef KS_DEBUG
1052 cout << "STSideBandSep::setLO1Root" << endl;
1053 if (!asdmName_.empty())
1054 cout << "asdm name = " << asdmName_ << endl;
1055 if (!asisName_.empty())
1056 cout << "MS name = " << asisName_ << endl;
1057#endif
1058};
1059
1060
1061void STSideBandSep::solveImageFrequency()
1062{
1063 LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
1064 os << "Start calculating frequencies of image side band" << LogIO::POST;
1065
1066 if (imgTab_p.null())
1067 throw AipsError("STSideBandSep::solveImageFreqency - an image side band scantable should be set first");
1068
1069 // Convert frequency REFVAL to the value in frame of LO.
1070 // The code assumes that imgTab_p has only an IF and only a FREQ_ID
1071 // is associated to an IFNO
1072 // TODO: More complete Procedure would be
1073 // 1. Get freq IDs associated to sigIfno_
1074 // 2. Get freq information of the freq IDs
1075 // 3. For each freqIDs, get freq infromation in TOPO and an LO1
1076 // frequency and calculate image band frequencies.
1077 STFrequencies freqTab_ = imgTab_p->frequencies();
1078 // get the base frame of table
1079 const MFrequency::Types tabframe = freqTab_.getFrame(true);
1080 TableVector<uInt> freqIdVec( imgTab_p->table(), "FREQ_ID" );
1081 // assuming single freqID per IFNO
1082 uInt freqid = freqIdVec(0);
1083 int nChan = imgTab_p->nchan(imgTab_p->getIF(0));
1084 double refpix, refval, increment ;
1085 freqTab_.getEntry(refpix, refval, increment, freqid);
1086 //MFrequency sigrefval = MFrequency(MVFrequency(refval),tabframe);
1087 // get freq infromation of sigIfno_ in loFrame_
1088 const MPosition mp = imgTab_p->getAntennaPosition();
1089 MEpoch me;
1090 MDirection md;
1091 if (loTime_ < 0.)
1092 me = imgTab_p->getEpoch(-1);
1093 else
1094 me = MEpoch(MVEpoch(loTime_));
1095 if (loDir_.empty()) {
1096 ArrayColumn<Double> srcdirCol_;
1097 srcdirCol_.attach(imgTab_p->table(), "SRCDIRECTION");
1098 // Assuming J2000 and SRCDIRECTION in unit of rad
1099 Quantum<Double> srcra = Quantum<Double>(srcdirCol_(0)(IPosition(1,0)), "rad");
1100 Quantum<Double> srcdec = Quantum<Double>(srcdirCol_(0)(IPosition(1,1)), "rad");
1101 md = MDirection(srcra, srcdec, MDirection::J2000);
1102 //imgTab_p->getDirection(0);
1103 } else {
1104 // parse direction string
1105 string::size_type pos0 = loDir_.find(" ");
1106
1107 if (pos0 == string::npos) {
1108 throw AipsError("bad string format in LO1 direction");
1109 }
1110 string::size_type pos1 = loDir_.find(" ", pos0+1);
1111 String sepoch, sra, sdec;
1112 if (pos1 != string::npos) {
1113 sepoch = loDir_.substr(0, pos0);
1114 sra = loDir_.substr(pos0+1, pos1-pos0);
1115 sdec = loDir_.substr(pos1+1);
1116 }
1117 MDirection::Types epoch;
1118 MDirection::getType(epoch, sepoch);
1119 QuantumHolder qh ;
1120 String err ;
1121 qh.fromString( err, sra);
1122 Quantum<Double> ra = qh.asQuantumDouble() ;
1123 qh.fromString( err, sdec ) ;
1124 Quantum<Double> dec = qh.asQuantumDouble() ;
1125 //md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch);
1126 md = MDirection(ra, dec, epoch);
1127 }
1128 MeasFrame mframe( me, mp, md );
1129 MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe));
1130 MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe));
1131 // Convert refval to loFrame_
1132 double sigrefval;
1133 if (tabframe == loFrame_)
1134 sigrefval = refval;
1135 else
1136 sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
1137
1138 // Check for the availability of LO1
1139 if (lo1Freq_ > 0.) {
1140 os << "Using user defined LO1 frequency." << LogIO::POST;
1141 } else if (!asisName_.empty()) {
1142 // MS name is set.
1143 os << "Using user defined MS (asis): " << asisName_ << LogIO::POST;
1144 if (!getLo1FromAsisTab(asisName_, sigrefval, refpix, increment, nChan)) {
1145 throw AipsError("Failed to get LO1 frequency from MS");
1146 }
1147 } else if (!asdmName_.empty()) {
1148 // ASDM name is set.
1149 os << "Using user defined ASDM: " << asdmName_ << LogIO::POST;
1150 if (!getLo1FromAsdm(asdmName_, sigrefval, refpix, increment, nChan)) {
1151 throw AipsError("Failed to get LO1 frequency from ASDM");
1152 }
1153 } else {
1154 // Try getting ASDM name from scantable header
1155 os << "Try getting information from scantable header" << LogIO::POST;
1156 if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) {
1157 //throw AipsError("Failed to get LO1 frequency from asis table");
1158 os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
1159 os << LogIO::WARN << "Could not fill frequency information of IMAGE sideband properly." << LogIO::POST;
1160 os << LogIO::WARN << "Storing values of SIGNAL sideband in FREQUENCIES table" << LogIO::POST;
1161 return;
1162 }
1163 }
1164
1165 // LO1 should now be ready.
1166 if (lo1Freq_ < 0.)
1167 throw(AipsError("Got negative LO1 Frequency"));
1168
1169 // Print summary
1170 {
1171 // LO1
1172 Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
1173 os << "[LO1 settings]" << LogIO::POST;
1174 os << "- Frequency: " << lo1Freq_ << " [Hz] ("
1175 << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
1176 os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
1177 << " [day]" << LogIO::POST;
1178 os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
1179 << "] (" << md.getRefString() << ") " << LogIO::POST;
1180
1181 // signal sideband
1182 os << "[Signal side band]" << LogIO::POST;
1183 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
1184 << LogIO::POST;
1185 os << "- Reference value: " << refval << " [Hz] ("
1186 << MFrequency::showType(tabframe) << ") = "
1187 << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_)
1188 << ")" << LogIO::POST;
1189 os << "- Reference pixel: " << refpix << LogIO::POST;
1190 os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
1191 }
1192
1193 // Calculate image band incr and refval in loFrame_
1194 Double imgincr = -increment;
1195 Double imgrefval = 2 * lo1Freq_ - sigrefval;
1196 Double imgrefval_tab = imgrefval;
1197 // Convert imgrefval back to table base frame
1198 if (tabframe != loFrame_)
1199 imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue();
1200 // Set new frequencies table
1201 uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr);
1202 // Update FREQ_ID in table.
1203 freqIdVec = fIDnew;
1204
1205 // Print summary (Image sideband)
1206 {
1207 os << "[Image side band]" << LogIO::POST;
1208 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
1209 << ")" << LogIO::POST;
1210 os << "- Reference value: " << imgrefval << " [Hz] ("
1211 << MFrequency::showType(tabframe) << ") = "
1212 << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_)
1213 << ")" << LogIO::POST;
1214 os << "- Reference pixel: " << refpix << LogIO::POST;
1215 os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
1216 }
1217};
1218
1219Bool STSideBandSep::checkFile(const string name, string type)
1220{
1221 File file(name);
1222 if (!file.exists()){
1223 return false;
1224 } else if (type.empty()) {
1225 return true;
1226 } else {
1227 // Check for file type
1228 switch (tolower(type[0])) {
1229 case 'f':
1230 return file.isRegular(True);
1231 case 'd':
1232 return file.isDirectory(True);
1233 case 's':
1234 return file.isSymLink();
1235 default:
1236 throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
1237 }
1238 }
1239};
1240
1241bool STSideBandSep::getLo1FromAsdm(const string asdmname,
1242 const double refval,
1243 const double refpix,
1244 const double increment,
1245 const int nChan)
1246{
1247 // Check for relevant tables.
1248 string spwname = asdmname + "/SpectralWindow.xml";
1249 string recname = asdmname + "/Receiver.xml";
1250 if (!checkFile(spwname) || !checkFile(recname)) {
1251 throw(AipsError("Could not find subtables in ASDM"));
1252 }
1253
1254 return false;
1255
1256};
1257
1258
1259bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab,
1260 const double refval,
1261 const double refpix,
1262 const double increment,
1263 const int nChan)
1264{
1265 LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE));
1266 // Check for relevant tables.
1267 const TableRecord &rec = scantab->table().keywordSet() ;
1268 String spwname, recname;
1269 if (rec.isDefined("ASDM_SPECTRALWINDOW") && rec.isDefined("ASDM_RECEIVER")){
1270 spwname = rec.asString("ASDM_SPECTRALWINDOW");
1271 recname = rec.asString("ASDM_RECEIVER");
1272 }
1273 else {
1274 // keywords are not there
1275 os << LogIO::WARN
1276 << "Could not find necessary table names in scantable header."
1277 << LogIO::POST;
1278 return false;
1279 }
1280 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1281 throw(AipsError("Could not find relevant subtables in MS"));
1282 }
1283 // Get root MS name
1284 string msname;
1285 const String recsuff = "/ASDM_RECEIVER";
1286 String::size_type pos;
1287 pos = recname.size()-recsuff.size();
1288 if (recname.substr(pos) == recsuff)
1289 msname = recname.substr(0, pos);
1290 else
1291 throw(AipsError("Internal error in parsing table name from a scantable keyword."));
1292
1293 if (!checkFile(msname))
1294 throw(AipsError("Internal error in parsing MS name from a scantable keyword."));
1295
1296 return getLo1FromAsisTab(msname, refval, refpix, increment, nChan);
1297
1298};
1299
1300bool STSideBandSep::getLo1FromAsisTab(const string msname,
1301 const double refval,
1302 const double refpix,
1303 const double increment,
1304 const int nChan)
1305{
1306 LogIO os(LogOrigin("STSideBandSep","getLo1FromAsisTab()", WHERE));
1307 os << "Searching an LO1 frequency in '" << msname << "'" << LogIO::POST;
1308 // Check for relevant tables.
1309 const string spwname = msname + "/ASDM_SPECTRALWINDOW";
1310 const string recname = msname + "/ASDM_RECEIVER";
1311 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1312 throw(AipsError("Could not find relevant tables in MS"));
1313 }
1314
1315 Table spwtab_ = Table(spwname);
1316 String asdmSpw;
1317 ROTableRow spwrow(spwtab_);
1318 const Double rtol = 0.01;
1319 for (uInt idx = 0; idx < spwtab_.nrow(); idx++){
1320 const TableRecord& rec = spwrow.get(idx);
1321 // Compare nchan
1322 if (rec.asInt("numChan") != (Int) nChan)
1323 continue;
1324 // Compare increment
1325 Double asdminc;
1326 Array<Double> incarr = rec.asArrayDouble("chanWidthArray");
1327 if (incarr.size() > 0)
1328 asdminc = incarr(IPosition(1, (uInt) refpix));
1329 else
1330 asdminc = rec.asDouble("chanWidth");
1331 if (abs(asdminc - abs(increment)) > rtol * abs(increment))
1332 continue;
1333 // Compare refval
1334 Double asdmrefv;
1335 Array<Double> refvarr = rec.asArrayDouble("chanFreqArray");
1336 if (refvarr.size() > 0){
1337 const uInt iref = (uInt) refpix;
1338 const Double ratio = refpix - (Double) iref;
1339 asdmrefv = refvarr(IPosition(1, iref))*(1.-ratio)
1340 + refvarr(IPosition(1,iref+1))*ratio;
1341 }
1342 else {
1343 const Double ch0 = rec.asDouble("chanFreqStart");
1344 const Double chstep = rec.asDouble("chanFreqStep");
1345 asdmrefv = ch0 + chstep * refpix;
1346 }
1347 if (abs(asdmrefv - refval) < 0.5*abs(asdminc)){
1348 asdmSpw = rec.asString("spectralWindowId");
1349 break;
1350 }
1351 }
1352
1353 if (asdmSpw.empty()){
1354 os << LogIO::WARN << "Could not find relevant SPW ID in " << spwname << LogIO::POST;
1355 return false;
1356 }
1357 else {
1358 os << asdmSpw << " in " << spwname
1359 << " matches the freqeuncies of signal side band." << LogIO::POST;
1360 }
1361
1362 Table rectab_ = Table(recname);
1363 ROTableRow recrow(rectab_);
1364 for (uInt idx = 0; idx < rectab_.nrow(); idx++){
1365 const TableRecord& rec = recrow.get(idx);
1366 if (rec.asString("spectralWindowId") == asdmSpw){
1367 const Array<Double> loarr = rec.asArrayDouble("freqLO");
1368 lo1Freq_ = loarr(IPosition(1,0));
1369 os << "Found LO1 Frequency in " << recname << ": "
1370 << lo1Freq_ << " [Hz]" << LogIO::POST;
1371 return true;
1372 }
1373 }
1374 os << LogIO::WARN << "Could not find " << asdmSpw << " in " << recname
1375 << LogIO::POST;
1376 return false;
1377};
1378
1379} //namespace asap
Note: See TracBrowser for help on using the repository browser.