source: trunk/src/STSideBandSep.cpp@ 2857

Last change on this file since 2857 was 2857, checked in by Kana Sugimoto, 11 years ago

New Development: No

JIRA Issue: Yes

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description: a bug fix to flag of image sideband.
File size: 46.8 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 Matrix<bool> flagMat(nchan_, nshift_);
264 vector<float> sigSpec(nchan_), imgSpec(nchan_);
265 Vector<bool> flagVec(nchan_);
266 vector<uInt> tabIdvec;
267
268 //Generate FFTServer
269 fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX);
270 fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL);
271
272 /// Loop over sigTab_p and separate sideband
273 for (int irow = 0; irow < sigTab_p->nrow(); irow++){
274 tabIdvec.resize(0);
275 const int polId = sigTab_p->getPol(irow);
276 const int beamId = sigTab_p->getBeam(irow);
277 const vector<double> dir = sigTab_p->getDirectionVector(irow);
278 // Get a set of spectra to solve
279 if (!getSpectraToSolve(polId, beamId, dir[0], dir[1],
280 specMat, flagMat, tabIdvec)){
281 remRowIds.push_back(irow);
282#ifdef KS_DEBUG
283 cout << "no matching row found. skipping row = " << irow << endl;
284#endif
285 continue;
286 }
287 // Solve signal sideband
288 sigSpec = solve(specMat, tabIdvec, true);
289 sigTab_p->setSpectrum(sigSpec, irow);
290 if (sigTab_p->isAllChannelsFlagged(irow)){
291 // unflag the spectrum since there should be some valid data
292 sigTab_p->flagRow(vector<uInt>(irow), true);
293 // need to unflag whole channels anyway
294 sigTab_p->flag(irow, vector<bool>(), true);
295 }
296 // apply channel flag
297 flagVec = collapseFlag(flagMat, tabIdvec, true);
298 //boolVec = !boolVec; // flag
299 vector<bool> tmpflag;
300 flagVec.tovector(tmpflag);
301 sigTab_p->flag(irow, tmpflag, false);
302
303 // Solve image sideband
304 if (doboth_) {
305 imgSpec = solve(specMat, tabIdvec, false);
306 imgTab_p->setSpectrum(imgSpec, irow);
307 if (imgTab_p->isAllChannelsFlagged(irow)){
308 // unflag the spectrum since there should be some valid data
309 imgTab_p->flagRow(vector<uInt>(irow), true);
310 // need to unflag whole channels anyway
311 imgTab_p->flag(irow, vector<bool>(), true);
312 }
313 // apply channel flag
314 flagVec = collapseFlag(flagMat, tabIdvec, false);
315 //boolVec = !boolVec; // flag
316 flagVec.tovector(tmpflag);
317 imgTab_p->flag(irow, tmpflag, false);
318 }
319 } // end of row loop
320
321 // Remove or flag rows without relevant data from gridded tables
322 if (remRowIds.size() > 0) {
323 const size_t nrem = remRowIds.size();
324 if ( sigTab_p->table().canRemoveRow() ) {
325 sigTab_p->table().removeRow(remRowIds);
326 os << "Removing " << nrem << " rows from the signal band table"
327 << LogIO::POST;
328 } else {
329 sigTab_p->flagRow(remRowIds, false);
330 os << "Cannot remove rows from the signal band table. Flagging "
331 << nrem << " rows" << LogIO::POST;
332 }
333
334 if (doboth_) {
335 if ( imgTab_p->table().canRemoveRow() ) {
336 imgTab_p->table().removeRow(remRowIds);
337 os << "Removing " << nrem << " rows from the image band table"
338 << LogIO::POST;
339 } else {
340 imgTab_p->flagRow(remRowIds, false);
341 os << "Cannot remove rows from the image band table. Flagging "
342 << nrem << " rows" << LogIO::POST;
343 }
344 }
345 }
346
347 // Finally, save tables on disk
348 if (outname.size() ==0)
349 outname = "sbseparated.asap";
350 const string sigName = outname + ".signalband";
351 os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST;
352 sigTab_p->makePersistent(sigName);
353 if (doboth_) {
354 solveImageFrequency();
355 const string imgName = outname + ".imageband";
356 os << "Saving IMAGE sideband table: " << sigName << LogIO::POST;
357 imgTab_p->makePersistent(imgName);
358 }
359
360};
361
362unsigned int STSideBandSep::setupShift()
363{
364 LogIO os(LogOrigin("STSideBandSep","setupShift()", WHERE));
365 if (infileList_.size() == 0 && intabList_.size() == 0)
366 throw( AipsError("No scantable has been set. Set a list of scantables first.") );
367
368 const bool byname = (intabList_.size() == 0);
369 // Make sure sigIfno_ exists in the first table.
370 CountedPtr<Scantable> stab;
371 vector<string> coordsav;
372 vector<string> coordtmp(3);
373 os << "Checking IFNO in the first table." << LogIO::POST;
374 if (byname) {
375 if (!checkFile(infileList_[0], "d"))
376 os << LogIO::SEVERE << "Could not find scantable '" << infileList_[0]
377 << "'" << LogIO::EXCEPTION;
378 stab = CountedPtr<Scantable>(new Scantable(infileList_[0]));
379 } else {
380 stab = intabList_[0];
381 }
382 if (sigIfno_ < 0) {
383 sigIfno_ = (int) stab->getIF(0);
384 os << "IFNO to process has not been set. Using the first IF = "
385 << sigIfno_ << LogIO::POST;
386 }
387
388 unsigned int basench;
389 double basech0, baseinc, ftolval, inctolval;
390 coordsav = stab->getCoordInfo();
391 const string stfframe = coordsav[1];
392 coordtmp[0] = "Hz";
393 coordtmp[1] = ( (solFrame_ == MFrequency::N_Types) ?
394 stfframe :
395 MFrequency::showType(solFrame_) );
396 coordtmp[2] = coordsav[2];
397 stab->setCoordInfo(coordtmp);
398 if (!getFreqInfo(stab, (unsigned int) sigIfno_, basech0, baseinc, basench)) {
399 os << LogIO::SEVERE << "No data with IFNO=" << sigIfno_
400 << " found in the first table." << LogIO::EXCEPTION;
401 }
402 else {
403 os << "Found IFNO = " << sigIfno_
404 << " in the first table." << LogIO::POST;
405 }
406 if (ftol_.getUnit().empty()) {
407 // tolerance in unit of channels
408 ftolval = ftol_.getValue() * baseinc;
409 }
410 else {
411 ftolval = ftol_.getValue("Hz");
412 }
413 inctolval = abs(baseinc/(double) basench);
414 const string poltype0 = stab->getPolType();
415
416 // Initialize shift values
417 initshift();
418
419 const bool setImg = ( doboth_ && (imgShift_.size() == 0) );
420 // Select IFs
421 for (unsigned int itab = 0; itab < ntable_; itab++ ){
422 os << "Table " << itab << LogIO::POST;
423 if (itab > 0) {
424 if (byname) {
425 if (!checkFile(infileList_[itab], "d"))
426 os << LogIO::SEVERE << "Could not find scantable '"
427 << infileList_[itab] << "'" << LogIO::EXCEPTION;
428 stab = CountedPtr<Scantable>(new Scantable(infileList_[itab]));
429 } else {
430 stab = intabList_[itab];
431 }
432 //POLTYPE should be the same.
433 if (stab->getPolType() != poltype0 ) {
434 os << LogIO::WARN << "POLTYPE differs from the first table."
435 << " Skipping the table" << LogIO::POST;
436 continue;
437 }
438 // Multi beam data may not handled properly
439 if (stab->nbeam() > 1)
440 os << LogIO::WARN << "Table contains multiple beams. "
441 << "It may not be handled properly." << LogIO::POST;
442
443 coordsav = stab->getCoordInfo();
444 coordtmp[2] = coordsav[2];
445 stab->setCoordInfo(coordtmp);
446 }
447 bool selected = false;
448 vector<uint> ifnos = stab->getIFNos();
449 vector<uint>::iterator iter;
450 const STSelector& basesel = stab->getSelection();
451 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
452 unsigned int nch;
453 double freq0, incr;
454 if ( getFreqInfo(stab, *iter, freq0, incr, nch) ){
455 if ( (nch == basench) && (abs(freq0-basech0) < ftolval)
456 && (abs(incr-baseinc) < inctolval) ){
457 //Found
458 STSelector sel(basesel);
459 sel.setIFs(vector<int>(1,(int) *iter));
460 stab->setSelection(sel);
461 CountedPtr<Scantable> seltab = ( new Scantable(*stab, false) );
462 stab->setSelection(basesel);
463 seltab->setCoordInfo(coordsav);
464 const double chShift = (freq0 - basech0) / baseinc;
465 tableList_.push_back(seltab);
466 sigShift_.push_back(-chShift);
467 if (setImg)
468 imgShift_.push_back(chShift);
469
470 selected = true;
471 os << "- IF" << *iter << " selected: sideband shift = "
472 << chShift << " channels" << LogIO::POST;
473 }
474 }
475 } // ifno loop
476 stab->setCoordInfo(coordsav);
477 if (!selected)
478 os << LogIO::WARN << "No data selected in Table "
479 << itab << LogIO::POST;
480 } // table loop
481 nchan_ = basench;
482
483 os << "Total number of IFs selected = " << tableList_.size()
484 << LogIO::POST;
485 if ( setImg && (sigShift_.size() != imgShift_.size()) ){
486 os << LogIO::SEVERE
487 << "User defined channel shift of image sideband has "
488 << imgShift_.size()
489 << " elements, while selected IFNOs are " << sigShift_.size()
490 << "\nThe frequency tolerance (freqtol) may be too small."
491 << LogIO::EXCEPTION;
492 }
493
494 return tableList_.size();
495};
496
497bool STSideBandSep::getFreqInfo(const CountedPtr<Scantable> &stab,
498 const unsigned int &ifno,
499 double &freq0, double &incr,
500 unsigned int &nchan)
501{
502 vector<uint> ifnos = stab->getIFNos();
503 bool found = false;
504 vector<uint>::iterator iter;
505 for (iter = ifnos.begin(); iter != ifnos.end(); iter++){
506 if (*iter == ifno) {
507 found = true;
508 break;
509 }
510 }
511 if (!found)
512 return false;
513
514 const STSelector& basesel = stab->getSelection();
515 STSelector sel(basesel);
516 sel.setIFs(vector<int>(1,(int) ifno));
517 stab->setSelection(sel);
518 vector<double> freqs;
519 freqs = stab->getAbcissa(0);
520 freq0 = freqs[0];
521 incr = freqs[1] - freqs[0];
522 nchan = freqs.size();
523 stab->setSelection(basesel);
524 return true;
525};
526
527ScantableWrapper STSideBandSep::gridTable()
528{
529 LogIO os(LogOrigin("STSideBandSep","gridTable()", WHERE));
530 if (tableList_.size() == 0)
531 throw( AipsError("Internal error. No scantable has been set to grid.") );
532 Double xmin, xmax, ymin, ymax;
533 mapExtent(tableList_, xmin, xmax, ymin, ymax);
534 const Double centx = 0.5 * (xmin + xmax);
535 const Double centy = 0.5 * (ymin + ymax);
536 const int nx = max(1, (int) ceil( (xmax - xmin) * cos(centy) /xtol_ ) );
537 const int ny = max(1, (int) ceil( (ymax - ymin) / ytol_ ) );
538
539 string scellx, scelly;
540 {
541 ostringstream oss;
542 oss << xtol_ << "rad" ;
543 scellx = oss.str();
544 }
545 {
546 ostringstream oss;
547 oss << ytol_ << "rad" ;
548 scelly = oss.str();
549 }
550
551 ScantableWrapper stab0;
552 if (intabList_.size() > 0)
553 stab0 = ScantableWrapper(intabList_[0]);
554 else
555 stab0 = ScantableWrapper(infileList_[0]);
556
557 string scenter;
558 {
559 ostringstream oss;
560 oss << stab0.getCP()->getDirectionRefString() << " "
561 << centx << "rad" << " " << centy << "rad";
562 scenter = oss.str();
563 }
564
565 STGrid2 gridder = STGrid2(stab0);
566 gridder.setIF(sigIfno_);
567 gridder.defineImage(nx, ny, scellx, scelly, scenter);
568 // gridder.setFunc("box", 1); // convsupport=1 fails
569 gridder.setFunc("box");
570 gridder.setWeight("uniform");
571#ifdef KS_DEBUG
572 cout << "Grid parameter summary: " << endl;
573 cout << "- IF = " << sigIfno_ << endl;
574 cout << "- center = " << scenter << "\n"
575 << "- npix = (" << nx << ", " << ny << ")\n"
576 << "- cell = (" << scellx << ", " << scelly << ")" << endl;
577#endif
578 gridder.grid();
579 const int itp = (tp_ == Table::Memory ? 0 : 1);
580 ScantableWrapper gtab = gridder.getResultAsScantable(itp);
581 return gtab;
582};
583
584void STSideBandSep::mapExtent(vector< CountedPtr<Scantable> > &tablist,
585 Double &xmin, Double &xmax,
586 Double &ymin, Double &ymax)
587{
588 ROArrayColumn<Double> dirCol_;
589 dirCol_.attach( tablist[0]->table(), "DIRECTION" );
590 Matrix<Double> direction = dirCol_.getColumn();
591 Vector<Double> ra( direction.row(0) );
592 mathutil::rotateRA(ra);
593 minMax( xmin, xmax, ra );
594 minMax( ymin, ymax, direction.row(1) );
595 Double amin, amax, bmin, bmax;
596 const uInt ntab = tablist.size();
597 for ( uInt i = 1 ; i < ntab ; i++ ) {
598 dirCol_.attach( tablist[i]->table(), "DIRECTION" );
599 direction.assign( dirCol_.getColumn() );
600 ra.assign( direction.row(0) );
601 mathutil::rotateRA(ra);
602 minMax( amin, amax, ra );
603 minMax( bmin, bmax, direction.row(1) );
604 xmin = min(xmin, amin);
605 xmax = max(xmax, amax);
606 ymin = min(ymin, bmin);
607 ymax = max(ymax, bmax);
608 }
609};
610
611// bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
612// const double dirX, const double dirY,
613// Matrix<float> &specMat, vector<uInt> &tabIdvec)
614bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId,
615 const double dirX, const double dirY,
616 Matrix<float> &specMat,
617 Matrix<bool> &flagMat,
618 vector<uInt> &tabIdvec)
619{
620 LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE));
621
622 tabIdvec.resize(0);
623 specMat.resize(nchan_, nshift_);
624 Vector<float> spec;
625 Vector<bool> boolVec;
626 uInt nspec = 0;
627 STMath stm(false); // insitu has no effect for average.
628 for (uInt itab = 0 ; itab < nshift_ ; itab++) {
629 CountedPtr<Scantable> currtab_p = tableList_[itab];
630 // Selection by POLNO and BEAMNO
631 const STSelector& basesel = currtab_p->getSelection();
632 STSelector sel(basesel);
633 sel.setPolarizations(vector<int>(1, polId));
634 sel.setBeams(vector<int>(1, beamId));
635 try {
636 currtab_p->setSelection(sel);
637 } catch (...) {
638#ifdef KS_DEBUG
639 cout << "Table " << itab << " - No spectrum found. skipping the table."
640 << endl;
641#endif
642 continue;
643 }
644 // Selection by direction;
645 vector<int> selrow(0);
646 vector<double> currDir(2, 0.);
647 const int nselrow = currtab_p->nrow();
648 for (int irow = 0 ; irow < nselrow ; irow++) {
649 currDir = currtab_p->getDirectionVector(irow);
650 if ( (abs(currDir[0]-dirX) > xtol_) ||
651 (abs(currDir[1]-dirY) > ytol_) )
652 continue;
653 // within direction tolerance
654 selrow.push_back(irow);
655 } // end of row loop
656
657 if (selrow.size() < 1) {
658 currtab_p->setSelection(basesel);
659
660#ifdef KS_DEBUG
661 cout << "Table " << itab << " - No spectrum found. skipping the table."
662 << endl;
663#endif
664
665 continue;
666 }
667
668 // At least a spectrum is selected in this table
669 CountedPtr<Scantable> seltab_p = ( new Scantable(*currtab_p, false) );
670 currtab_p->setSelection(basesel);
671 STSelector sel2(seltab_p->getSelection());
672 sel2.setRows(selrow);
673 seltab_p->setSelection(sel2);
674 CountedPtr<Scantable> avetab_p;
675 if (seltab_p->nrow() > 1) {
676 // STMath::average also merges FLAGTRA and FLAGROW
677 avetab_p = stm.average(vector< CountedPtr<Scantable> >(1, seltab_p),
678 vector<bool>(), "TINTSYS", "NONE");
679#ifdef KS_DEBUG
680 cout << "Table " << itab << " - more than a spectrum is selected. averaging rows..."
681 << endl;
682#endif
683 if (avetab_p->nrow() > 1)
684 throw( AipsError("Averaged table has more than a row. Somethigs went wrong.") );
685 } else {
686 avetab_p = seltab_p;
687 }
688 // Check FLAGTRA and FLAGROW if there's any valid channel in the spectrum
689 if (avetab_p->getFlagRow(0) || avetab_p->isAllChannelsFlagged(0)) {
690#ifdef KS_DEBUG
691 cout << "Table " << itab << " - All data are flagged. skipping the table."
692 << endl;
693#endif
694 continue;
695 }
696 // Interpolate flagged channels of the spectrum.
697 Vector<Float> tmpSpec = avetab_p->getSpectrum(0);
698 // Mask is true if the data is valid (!flag)
699 vector<bool> mask = avetab_p->getMask(0);
700 mathutil::doZeroOrderInterpolation(tmpSpec, mask);
701 spec.reference(specMat.column(nspec));
702 spec = tmpSpec;
703 boolVec.reference(flagMat.column(nspec));
704 boolVec = mask; // cast std::vector to casa::Vector
705 boolVec = !boolVec;
706 tabIdvec.push_back((uInt) itab);
707 nspec++;
708 //Liberate from reference
709 spec.unique();
710 boolVec.unique();
711 } // end of table loop
712 // Check the number of selected spectra and resize matrix.
713 if (nspec != nshift_){
714 //shiftSpecmat.resize(nchan_, nspec, true);
715 specMat.resize(nchan_, nspec, true);
716 flagMat.resize(nchan_, nspec, true);
717#ifdef KS_DEBUG
718 cout << "Could not find corresponding rows in some tables."
719 << endl;
720 cout << "Number of spectra selected = " << nspec << endl;
721#endif
722 }
723 if (nspec < 2) {
724#ifdef KS_DEBUG
725 cout << "At least 2 spectra are necessary for convolution"
726 << endl;
727#endif
728 return false;
729 }
730 return true;
731};
732
733
734Vector<bool> STSideBandSep::collapseFlag(const Matrix<bool> &flagMat,
735 const vector<uInt> &tabIdvec,
736 const bool signal)
737{
738 LogIO os(LogOrigin("STSideBandSep","collapseFlag()", WHERE));
739 if (tabIdvec.size() == 0)
740 throw(AipsError("Internal error. Table index is not defined."));
741 if (flagMat.ncolumn() != tabIdvec.size())
742 throw(AipsError("Internal error. The row number of input matrix is not conformant."));
743 if (flagMat.nrow() != nchan_)
744 throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
745
746 const size_t nspec = tabIdvec.size();
747 vector<double> *thisShift;
748 if (signal == otherside_) {
749 // (solve signal && solveother = T) OR (solve image && solveother = F)
750 thisShift = &imgShift_;
751 } else {
752 // (solve signal && solveother = F) OR (solve image && solveother = T)
753 thisShift = &sigShift_;
754 }
755
756 Vector<bool> outflag(nchan_, false);
757 double tempshift;
758 Vector<bool> shiftvec(nchan_, false);
759 Vector<bool> accflag(nchan_, false);
760 uInt shiftId;
761 for (uInt i = 0 ; i < nspec; ++i) {
762 shiftId = tabIdvec[i];
763 tempshift = - thisShift->at(shiftId);
764 shiftFlag(flagMat.column(i), tempshift, shiftvec);
765 // Now accumulate Flag
766 for (uInt j = 0 ; j < nchan_ ; ++j)
767 accflag[j] |= shiftvec[j];
768 }
769 // Shift back Flag
770 shiftFlag(accflag, thisShift->at(0), outflag);
771
772 return outflag;
773}
774
775
776vector<float> STSideBandSep::solve(const Matrix<float> &specmat,
777 const vector<uInt> &tabIdvec,
778 const bool signal)
779{
780 LogIO os(LogOrigin("STSideBandSep","solve()", WHERE));
781 if (tabIdvec.size() == 0)
782 throw(AipsError("Internal error. Table index is not defined."));
783 if (specmat.ncolumn() != tabIdvec.size())
784 throw(AipsError("Internal error. The row number of input matrix is not conformant."));
785 if (specmat.nrow() != nchan_)
786 throw(AipsError("Internal error. The channel size of input matrix is not conformant."));
787
788
789#ifdef KS_DEBUG
790 cout << "Solving " << (signal ? "SIGNAL" : "IMAGE") << " sideband."
791 << endl;
792#endif
793
794 const size_t nspec = tabIdvec.size();
795 vector<double> *thisShift, *otherShift;
796 if (signal == otherside_) {
797 // (solve signal && solveother = T) OR (solve image && solveother = F)
798 thisShift = &imgShift_;
799 otherShift = &sigShift_;
800#ifdef KS_DEBUG
801 cout << "Image sideband will be deconvolved." << endl;
802#endif
803 } else {
804 // (solve signal && solveother = F) OR (solve image && solveother = T)
805 thisShift = &sigShift_;
806 otherShift = &imgShift_;
807#ifdef KS_DEBUG
808 cout << "Signal sideband will be deconvolved." << endl;
809#endif
810 }
811
812 vector<double> spshift(nspec);
813 Matrix<float> shiftSpecmat(nchan_, nspec, 0.);
814 double tempshift;
815 Vector<float> shiftspvec;
816 uInt shiftId;
817 for (uInt i = 0 ; i < nspec; i++) {
818 shiftId = tabIdvec[i];
819 spshift[i] = otherShift->at(shiftId) - thisShift->at(shiftId);
820 tempshift = - thisShift->at(shiftId);
821 shiftspvec.reference(shiftSpecmat.column(i));
822 shiftSpectrum(specmat.column(i), tempshift, shiftspvec);
823 }
824
825 Matrix<float> convmat(nchan_, nspec*(nspec-1)/2, 0.);
826 vector<float> thisvec(nchan_, 0.);
827
828 float minval, maxval;
829 minMax(minval, maxval, shiftSpecmat);
830#ifdef KS_DEBUG
831 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
832#endif
833
834#ifdef KS_DEBUG
835 cout << "starting deconvolution" << endl;
836#endif
837 deconvolve(shiftSpecmat, spshift, rejlimit_, convmat);
838#ifdef KS_DEBUG
839 cout << "finished deconvolution" << endl;
840#endif
841
842 minMax(minval, maxval, convmat);
843#ifdef KS_DEBUG
844 cout << "Max/Min of output Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
845#endif
846
847 aggregateMat(convmat, thisvec);
848
849 if (!otherside_) return thisvec;
850
851 // subtract from the other side band.
852 vector<float> othervec(nchan_, 0.);
853 subtractFromOther(shiftSpecmat, thisvec, spshift, othervec);
854 return othervec;
855};
856
857
858void STSideBandSep::shiftSpectrum(const Vector<float> &invec,
859 double shift,
860 Vector<float> &outvec)
861{
862 LogIO os(LogOrigin("STSideBandSep","shiftSpectrum()", WHERE));
863 if (invec.size() != nchan_)
864 throw(AipsError("Internal error. The length of input vector differs from nchan_"));
865 if (outvec.size() != nchan_)
866 throw(AipsError("Internal error. The length of output vector differs from nchan_"));
867
868#ifdef KS_DEBUG
869 cout << "Start shifting spectrum for " << shift << "channels" << endl;
870#endif
871
872 // tweak shift to be in 0 ~ nchan_-1
873 if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
874 if (shift < 0.) shift += nchan_;
875 double rweight = fmod(shift, 1.);
876 if (rweight < 0.) rweight += 1.;
877 double lweight = 1. - rweight;
878 uInt lchan, rchan;
879
880 outvec = 0;
881 for (uInt i = 0 ; i < nchan_ ; i++) {
882 lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
883 if (lchan < 0.) lchan += nchan_;
884 rchan = ( (lchan + 1) % nchan_ );
885 outvec(lchan) += invec(i) * lweight;
886 outvec(rchan) += invec(i) * rweight;
887 }
888};
889
890
891void STSideBandSep::shiftFlag(const Vector<bool> &invec,
892 double shift,
893 Vector<bool> &outvec)
894{
895 LogIO os(LogOrigin("STSideBandSep","shiftFlag()", WHERE));
896 if (invec.size() != nchan_)
897 throw(AipsError("Internal error. The length of input vector differs from nchan_"));
898 if (outvec.size() != nchan_)
899 throw(AipsError("Internal error. The length of output vector differs from nchan_"));
900
901#ifdef KS_DEBUG
902 cout << "Start shifting flag for " << shift << "channels" << endl;
903#endif
904
905 // shift is almost integer think it as int.
906 // tolerance should be in 0 - 1
907 double tolerance = 0.01;
908 // tweak shift to be in 0 ~ nchan_-1
909 if ( fabs(shift) > nchan_ ) shift = fmod(shift, nchan_);
910 if (shift < 0.) shift += nchan_;
911 double rweight = fmod(shift, 1.);
912 bool ruse(true), luse(true);
913 if (rweight < 0.) rweight += 1.;
914 if (rweight < tolerance){
915 ruse = true;
916 luse = false;
917 }
918 if (rweight > 1-tolerance){
919 ruse = false;
920 luse = true;
921 }
922 uInt lchan, rchan;
923
924 outvec = false;
925 for (uInt i = 0 ; i < nchan_ ; i++) {
926 lchan = uInt( floor( fmod( (i + shift), nchan_ ) ) );
927 if (lchan < 0.) lchan += nchan_;
928 rchan = ( (lchan + 1) % nchan_ );
929 outvec(lchan) |= (invec(i) && luse);
930 outvec(rchan) |= (invec(i) && ruse);
931 }
932};
933
934
935void STSideBandSep::deconvolve(Matrix<float> &specmat,
936 const vector<double> shiftvec,
937 const double threshold,
938 Matrix<float> &outmat)
939{
940 LogIO os(LogOrigin("STSideBandSep","deconvolve()", WHERE));
941 if (specmat.nrow() != nchan_)
942 throw(AipsError("Internal error. The length of input matrix differs from nchan_"));
943 if (specmat.ncolumn() != shiftvec.size())
944 throw(AipsError("Internal error. The number of input shifts and spectrum differs."));
945
946#ifdef KS_DEBUG
947 float minval, maxval;
948#endif
949#ifdef KS_DEBUG
950 minMax(minval, maxval, specmat);
951 cout << "Max/Min of input Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
952#endif
953
954 uInt ninsp = shiftvec.size();
955 outmat.resize(nchan_, ninsp*(ninsp-1)/2, 0.);
956 Matrix<Complex> fftspmat(nchan_/2+1, ninsp, 0.);
957 Vector<float> rvecref(nchan_, 0.);
958 Vector<Complex> cvecref(nchan_/2+1, 0.);
959 uInt icol = 0;
960 unsigned int nreject = 0;
961
962#ifdef KS_DEBUG
963 cout << "Starting initial FFT. The number of input spectra = " << ninsp << endl;
964 cout << "out matrix has ncolumn = " << outmat.ncolumn() << endl;
965#endif
966
967 for (uInt isp = 0 ; isp < ninsp ; isp++) {
968 rvecref.reference( specmat.column(isp) );
969 cvecref.reference( fftspmat.column(isp) );
970
971#ifdef KS_DEBUG
972 minMax(minval, maxval, rvecref);
973 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
974#endif
975
976 fftsf.fft0(cvecref, rvecref, true);
977
978#ifdef KS_DEBUG
979 double maxr=cvecref[0].real(), minr=cvecref[0].real(),
980 maxi=cvecref[0].imag(), mini=cvecref[0].imag();
981 for (uInt i = 1; i<cvecref.size();i++){
982 maxr = max(maxr, cvecref[i].real());
983 maxi = max(maxi, cvecref[i].imag());
984 minr = min(minr, cvecref[i].real());
985 mini = min(mini, cvecref[i].imag());
986 }
987 cout << "Max/Min of inv FFTed Matrix (size=" << cvecref.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
988#endif
989 }
990
991 //Liberate from reference
992 rvecref.unique();
993
994 Vector<Complex> cspec(nchan_/2+1, 0.);
995 const double PI = 6.0 * asin(0.5);
996 const double nchani = 1. / (float) nchan_;
997 const Complex trans(0., 1.);
998#ifdef KS_DEBUG
999 cout << "starting actual deconvolution" << endl;
1000#endif
1001 for (uInt j = 0 ; j < ninsp ; j++) {
1002 for (uInt k = j+1 ; k < ninsp ; k++) {
1003 const double dx = (shiftvec[k] - shiftvec[j]) * 2. * PI * nchani;
1004
1005#ifdef KS_DEBUG
1006 cout << "icol = " << icol << endl;
1007#endif
1008
1009 for (uInt ichan = 0 ; ichan < cspec.size() ; ichan++){
1010 cspec[ichan] = ( fftspmat(ichan, j) + fftspmat(ichan, k) )*0.5;
1011 double phase = dx*ichan;
1012 if ( fabs( sin(phase) ) > threshold){
1013 cspec[ichan] += ( fftspmat(ichan, j) - fftspmat(ichan, k) ) * 0.5
1014 * trans * sin(phase) / ( 1. - cos(phase) );
1015 } else {
1016 nreject++;
1017 }
1018 } // end of channel loop
1019
1020#ifdef KS_DEBUG
1021 cout << "done calculation of cspec" << endl;
1022#endif
1023
1024 Vector<Float> rspec;
1025 rspec.reference( outmat.column(icol) );
1026
1027#ifdef KS_DEBUG
1028 cout << "Starting inverse FFT. icol = " << icol << endl;
1029 //cout << "- size of complex vector = " << cspec.size() << endl;
1030 double maxr=cspec[0].real(), minr=cspec[0].real(),
1031 maxi=cspec[0].imag(), mini=cspec[0].imag();
1032 for (uInt i = 1; i<cspec.size();i++){
1033 maxr = max(maxr, cspec[i].real());
1034 maxi = max(maxi, cspec[i].imag());
1035 minr = min(minr, cspec[i].real());
1036 mini = min(mini, cspec[i].imag());
1037 }
1038 cout << "Max/Min of conv vector (size=" << cspec.size() << ") = (max: " << maxr << " + " << maxi << "j , min: " << minr << " + " << mini << "j)" << endl;
1039#endif
1040
1041 fftsi.fft0(rspec, cspec, false);
1042
1043#ifdef KS_DEBUG
1044 //cout << "- size of inversed real vector = " << rspec.size() << endl;
1045 minMax(minval, maxval, rspec);
1046 cout << "Max/Min of inv FFTed Vector (size=" << rspec.size() << ") = (max: " << maxval << ", min: " << minval << ")" << endl;
1047 //cout << "Done inverse FFT. icol = " << icol << endl;
1048#endif
1049
1050 icol++;
1051 }
1052 }
1053
1054#ifdef KS_DEBUG
1055 minMax(minval, maxval, outmat);
1056 cout << "Max/Min of inv FFTed Matrix = (max: " << maxval << ", min: " << minval << ")" << endl;
1057#endif
1058
1059 os << "Threshold = " << threshold << ", Rejected channels = " << nreject << endl;
1060};
1061
1062
1063void STSideBandSep::aggregateMat(Matrix<float> &inmat,
1064 vector<float> &outvec)
1065{
1066 LogIO os(LogOrigin("STSideBandSep","aggregateMat()", WHERE));
1067 if (inmat.nrow() != nchan_)
1068 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
1069// if (outvec.size() != nchan_)
1070// throw(AipsError("Internal error. The size of output vector should be equal to nchan_"));
1071
1072 os << "Averaging " << inmat.ncolumn() << " spectra in the input matrix."
1073 << LogIO::POST;
1074
1075 const uInt nspec = inmat.ncolumn();
1076 const double scale = 1./( (double) nspec );
1077 // initialize values with 0
1078 outvec.assign(nchan_, 0);
1079 for (uInt isp = 0 ; isp < nspec ; isp++) {
1080 for (uInt ich = 0 ; ich < nchan_ ; ich++) {
1081 outvec[ich] += inmat(ich, isp);
1082 }
1083 }
1084
1085 vector<float>::iterator iter;
1086 for (iter = outvec.begin(); iter != outvec.end(); iter++){
1087 *iter *= scale;
1088 }
1089};
1090
1091void STSideBandSep::subtractFromOther(const Matrix<float> &shiftmat,
1092 const vector<float> &invec,
1093 const vector<double> &shift,
1094 vector<float> &outvec)
1095{
1096 LogIO os(LogOrigin("STSideBandSep","subtractFromOther()", WHERE));
1097 if (shiftmat.nrow() != nchan_)
1098 throw(AipsError("Internal error. The row numbers of input matrix differs from nchan_"));
1099 if (invec.size() != nchan_)
1100 throw(AipsError("Internal error. The length of input vector should be nchan_"));
1101 if (shift.size() != shiftmat.ncolumn())
1102 throw(AipsError("Internal error. The column numbers of input matrix != the number of elements in shift"));
1103
1104 const uInt nspec = shiftmat.ncolumn();
1105 Vector<float> subsp(nchan_, 0.), shiftsub;
1106 Matrix<float> submat(nchan_, nspec, 0.);
1107 vector<float>::iterator iter;
1108 for (uInt isp = 0 ; isp < nspec ; isp++) {
1109 for (uInt ich = 0; ich < nchan_ ; ich++) {
1110 subsp(ich) = shiftmat(ich, isp) - invec[ich];
1111 }
1112 shiftsub.reference(submat.column(isp));
1113 shiftSpectrum(subsp, shift[isp], shiftsub);
1114 }
1115
1116 aggregateMat(submat, outvec);
1117};
1118
1119
1120void STSideBandSep::setLO1(const string lo1, const string frame,
1121 const double reftime, const string refdir)
1122{
1123 Quantum<Double> qfreq;
1124 readQuantity(qfreq, String(lo1));
1125 lo1Freq_ = qfreq.getValue("Hz");
1126 MFrequency::getType(loFrame_, frame);
1127 loTime_ = reftime;
1128 loDir_ = refdir;
1129
1130#ifdef KS_DEBUG
1131 cout << "STSideBandSep::setLO1" << endl;
1132 if (lo1Freq_ > 0.)
1133 cout << "lo1 = " << lo1Freq_ << " [Hz] (" << frame << ")" << endl;
1134 if (loTime_ > 0.)
1135 cout << "ref time = " << loTime_ << " [day]" << endl;
1136 if (!loDir_.empty())
1137 cout << "ref direction = " << loDir_ << " [day]" << endl;
1138#endif
1139};
1140
1141void STSideBandSep::setLO1Root(string name)
1142{
1143 LogIO os(LogOrigin("STSideBandSep","setLO1Root()", WHERE));
1144 os << "Searching for '" << name << "'..." << LogIO::POST;
1145 // Check for existance of the file
1146 if (!checkFile(name)) {
1147 throw(AipsError("File does not exist"));
1148 }
1149 if (name[(name.size()-1)] == '/')
1150 name = name.substr(0,(name.size()-2));
1151
1152 if (checkFile(name+"/Receiver.xml", "file") &&
1153 checkFile(name+"/SpectralWindow.xml", "file")){
1154 os << "Found '" << name << "/Receiver.xml' ... got an ASDM name." << LogIO::POST;
1155 asdmName_ = name;
1156 } else if (checkFile(name+"/ASDM_RECEIVER") &&
1157 checkFile(name+"/ASDM_SPECTRALWINDOW")){
1158 os << "Found '" << name << "/ASDM_RECEIVER' ... got a Table name." << LogIO::POST;
1159 asisName_ = name;
1160 } else {
1161 throw(AipsError("Invalid file name. Set an MS or ASDM name."));
1162 }
1163
1164#ifdef KS_DEBUG
1165 cout << "STSideBandSep::setLO1Root" << endl;
1166 if (!asdmName_.empty())
1167 cout << "asdm name = " << asdmName_ << endl;
1168 if (!asisName_.empty())
1169 cout << "MS name = " << asisName_ << endl;
1170#endif
1171};
1172
1173
1174void STSideBandSep::solveImageFrequency()
1175{
1176 LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
1177 os << "Start calculating frequencies of image side band" << LogIO::POST;
1178
1179 if (imgTab_p.null())
1180 throw AipsError("STSideBandSep::solveImageFreqency - an image side band scantable should be set first");
1181
1182 // Convert frequency REFVAL to the value in frame of LO.
1183 // The code assumes that imgTab_p has only an IF and only a FREQ_ID
1184 // is associated to an IFNO
1185 // TODO: More complete Procedure would be
1186 // 1. Get freq IDs associated to sigIfno_
1187 // 2. Get freq information of the freq IDs
1188 // 3. For each freqIDs, get freq infromation in TOPO and an LO1
1189 // frequency and calculate image band frequencies.
1190 STFrequencies freqTab_ = imgTab_p->frequencies();
1191 // get the base frame of table
1192 const MFrequency::Types tabframe = freqTab_.getFrame(true);
1193 TableVector<uInt> freqIdVec( imgTab_p->table(), "FREQ_ID" );
1194 // assuming single freqID per IFNO
1195 uInt freqid = freqIdVec(0);
1196 int nChan = imgTab_p->nchan(imgTab_p->getIF(0));
1197 double refpix, refval, increment ;
1198 freqTab_.getEntry(refpix, refval, increment, freqid);
1199 //MFrequency sigrefval = MFrequency(MVFrequency(refval),tabframe);
1200 // get freq infromation of sigIfno_ in loFrame_
1201 const MPosition mp = imgTab_p->getAntennaPosition();
1202 MEpoch me;
1203 MDirection md;
1204 if (loTime_ < 0.)
1205 me = imgTab_p->getEpoch(-1);
1206 else
1207 me = MEpoch(MVEpoch(loTime_));
1208 if (loDir_.empty()) {
1209 ArrayColumn<Double> srcdirCol_;
1210 srcdirCol_.attach(imgTab_p->table(), "SRCDIRECTION");
1211 // Assuming J2000 and SRCDIRECTION in unit of rad
1212 Quantum<Double> srcra = Quantum<Double>(srcdirCol_(0)(IPosition(1,0)), "rad");
1213 Quantum<Double> srcdec = Quantum<Double>(srcdirCol_(0)(IPosition(1,1)), "rad");
1214 md = MDirection(srcra, srcdec, MDirection::J2000);
1215 //imgTab_p->getDirection(0);
1216 } else {
1217 // parse direction string
1218 string::size_type pos0 = loDir_.find(" ");
1219
1220 if (pos0 == string::npos) {
1221 throw AipsError("bad string format in LO1 direction");
1222 }
1223 string::size_type pos1 = loDir_.find(" ", pos0+1);
1224 String sepoch, sra, sdec;
1225 if (pos1 != string::npos) {
1226 sepoch = loDir_.substr(0, pos0);
1227 sra = loDir_.substr(pos0+1, pos1-pos0);
1228 sdec = loDir_.substr(pos1+1);
1229 }
1230 MDirection::Types epoch;
1231 MDirection::getType(epoch, sepoch);
1232 QuantumHolder qh ;
1233 String err ;
1234 qh.fromString( err, sra);
1235 Quantum<Double> ra = qh.asQuantumDouble() ;
1236 qh.fromString( err, sdec ) ;
1237 Quantum<Double> dec = qh.asQuantumDouble() ;
1238 //md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch);
1239 md = MDirection(ra, dec, epoch);
1240 }
1241 MeasFrame mframe( me, mp, md );
1242 MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe));
1243 MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe));
1244 // Convert refval to loFrame_
1245 double sigrefval;
1246 if (tabframe == loFrame_)
1247 sigrefval = refval;
1248 else
1249 sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
1250
1251 // Check for the availability of LO1
1252 if (lo1Freq_ > 0.) {
1253 os << "Using user defined LO1 frequency." << LogIO::POST;
1254 } else if (!asisName_.empty()) {
1255 // MS name is set.
1256 os << "Using user defined MS (asis): " << asisName_ << LogIO::POST;
1257 if (!getLo1FromAsisTab(asisName_, sigrefval, refpix, increment, nChan)) {
1258 throw AipsError("Failed to get LO1 frequency from MS");
1259 }
1260 } else if (!asdmName_.empty()) {
1261 // ASDM name is set.
1262 os << "Using user defined ASDM: " << asdmName_ << LogIO::POST;
1263 if (!getLo1FromAsdm(asdmName_, sigrefval, refpix, increment, nChan)) {
1264 throw AipsError("Failed to get LO1 frequency from ASDM");
1265 }
1266 } else {
1267 // Try getting ASDM name from scantable header
1268 os << "Try getting information from scantable header" << LogIO::POST;
1269 if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) {
1270 //throw AipsError("Failed to get LO1 frequency from asis table");
1271 os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
1272 os << LogIO::WARN << "Could not fill frequency information of IMAGE sideband properly." << LogIO::POST;
1273 os << LogIO::WARN << "Storing values of SIGNAL sideband in FREQUENCIES table" << LogIO::POST;
1274 return;
1275 }
1276 }
1277
1278 // LO1 should now be ready.
1279 if (lo1Freq_ < 0.)
1280 throw(AipsError("Got negative LO1 Frequency"));
1281
1282 // Print summary
1283 {
1284 // LO1
1285 Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
1286 os << "[LO1 settings]" << LogIO::POST;
1287 os << "- Frequency: " << lo1Freq_ << " [Hz] ("
1288 << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
1289 os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
1290 << " [day]" << LogIO::POST;
1291 os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
1292 << "] (" << md.getRefString() << ") " << LogIO::POST;
1293
1294 // signal sideband
1295 os << "[Signal side band]" << LogIO::POST;
1296 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
1297 << LogIO::POST;
1298 os << "- Reference value: " << refval << " [Hz] ("
1299 << MFrequency::showType(tabframe) << ") = "
1300 << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_)
1301 << ")" << LogIO::POST;
1302 os << "- Reference pixel: " << refpix << LogIO::POST;
1303 os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
1304 }
1305
1306 // Calculate image band incr and refval in loFrame_
1307 Double imgincr = -increment;
1308 Double imgrefval = 2 * lo1Freq_ - sigrefval;
1309 Double imgrefval_tab = imgrefval;
1310 // Convert imgrefval back to table base frame
1311 if (tabframe != loFrame_)
1312 imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue();
1313 // Set new frequencies table
1314 uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr);
1315 // Update FREQ_ID in table.
1316 freqIdVec = fIDnew;
1317
1318 // Print summary (Image sideband)
1319 {
1320 os << "[Image side band]" << LogIO::POST;
1321 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
1322 << ")" << LogIO::POST;
1323 os << "- Reference value: " << imgrefval << " [Hz] ("
1324 << MFrequency::showType(tabframe) << ") = "
1325 << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_)
1326 << ")" << LogIO::POST;
1327 os << "- Reference pixel: " << refpix << LogIO::POST;
1328 os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
1329 }
1330};
1331
1332Bool STSideBandSep::checkFile(const string name, string type)
1333{
1334 File file(name);
1335 if (!file.exists()){
1336 return false;
1337 } else if (type.empty()) {
1338 return true;
1339 } else {
1340 // Check for file type
1341 switch (tolower(type[0])) {
1342 case 'f':
1343 return file.isRegular(True);
1344 case 'd':
1345 return file.isDirectory(True);
1346 case 's':
1347 return file.isSymLink();
1348 default:
1349 throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
1350 }
1351 }
1352};
1353
1354bool STSideBandSep::getLo1FromAsdm(const string asdmname,
1355 const double /*refval*/,
1356 const double /*refpix*/,
1357 const double /*increment*/,
1358 const int /*nChan*/)
1359{
1360 // Check for relevant tables.
1361 string spwname = asdmname + "/SpectralWindow.xml";
1362 string recname = asdmname + "/Receiver.xml";
1363 if (!checkFile(spwname) || !checkFile(recname)) {
1364 throw(AipsError("Could not find subtables in ASDM"));
1365 }
1366
1367 return false;
1368
1369};
1370
1371
1372bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab,
1373 const double refval,
1374 const double refpix,
1375 const double increment,
1376 const int nChan)
1377{
1378 LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE));
1379 // Check for relevant tables.
1380 const TableRecord &rec = scantab->table().keywordSet() ;
1381 String spwname, recname;
1382 if (rec.isDefined("ASDM_SPECTRALWINDOW") && rec.isDefined("ASDM_RECEIVER")){
1383 spwname = rec.asString("ASDM_SPECTRALWINDOW");
1384 recname = rec.asString("ASDM_RECEIVER");
1385 }
1386 else {
1387 // keywords are not there
1388 os << LogIO::WARN
1389 << "Could not find necessary table names in scantable header."
1390 << LogIO::POST;
1391 return false;
1392 }
1393 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1394 throw(AipsError("Could not find relevant subtables in MS"));
1395 }
1396 // Get root MS name
1397 string msname;
1398 const String recsuff = "/ASDM_RECEIVER";
1399 String::size_type pos;
1400 pos = recname.size()-recsuff.size();
1401 if (recname.substr(pos) == recsuff)
1402 msname = recname.substr(0, pos);
1403 else
1404 throw(AipsError("Internal error in parsing table name from a scantable keyword."));
1405
1406 if (!checkFile(msname))
1407 throw(AipsError("Internal error in parsing MS name from a scantable keyword."));
1408
1409 return getLo1FromAsisTab(msname, refval, refpix, increment, nChan);
1410
1411};
1412
1413bool STSideBandSep::getLo1FromAsisTab(const string msname,
1414 const double refval,
1415 const double refpix,
1416 const double increment,
1417 const int nChan)
1418{
1419 LogIO os(LogOrigin("STSideBandSep","getLo1FromAsisTab()", WHERE));
1420 os << "Searching an LO1 frequency in '" << msname << "'" << LogIO::POST;
1421 // Check for relevant tables.
1422 const string spwname = msname + "/ASDM_SPECTRALWINDOW";
1423 const string recname = msname + "/ASDM_RECEIVER";
1424 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
1425 throw(AipsError("Could not find relevant tables in MS"));
1426 }
1427
1428 Table spwtab_ = Table(spwname);
1429 String asdmSpw;
1430 ROTableRow spwrow(spwtab_);
1431 const Double rtol = 0.01;
1432 for (uInt idx = 0; idx < spwtab_.nrow(); idx++){
1433 const TableRecord& rec = spwrow.get(idx);
1434 // Compare nchan
1435 if (rec.asInt("numChan") != (Int) nChan)
1436 continue;
1437 // Compare increment
1438 Double asdminc;
1439 Array<Double> incarr = rec.asArrayDouble("chanWidthArray");
1440 if (incarr.size() > 0)
1441 asdminc = incarr(IPosition(1, (uInt) refpix));
1442 else
1443 asdminc = rec.asDouble("chanWidth");
1444 if (abs(asdminc - abs(increment)) > rtol * abs(increment))
1445 continue;
1446 // Compare refval
1447 Double asdmrefv;
1448 Array<Double> refvarr = rec.asArrayDouble("chanFreqArray");
1449 if (refvarr.size() > 0){
1450 const uInt iref = (uInt) refpix;
1451 const Double ratio = refpix - (Double) iref;
1452 asdmrefv = refvarr(IPosition(1, iref))*(1.-ratio)
1453 + refvarr(IPosition(1,iref+1))*ratio;
1454 }
1455 else {
1456 const Double ch0 = rec.asDouble("chanFreqStart");
1457 const Double chstep = rec.asDouble("chanFreqStep");
1458 asdmrefv = ch0 + chstep * refpix;
1459 }
1460 if (abs(asdmrefv - refval) < 0.5*abs(asdminc)){
1461 asdmSpw = rec.asString("spectralWindowId");
1462 break;
1463 }
1464 }
1465
1466 if (asdmSpw.empty()){
1467 os << LogIO::WARN << "Could not find relevant SPW ID in " << spwname << LogIO::POST;
1468 return false;
1469 }
1470 else {
1471 os << asdmSpw << " in " << spwname
1472 << " matches the freqeuncies of signal side band." << LogIO::POST;
1473 }
1474
1475 Table rectab_ = Table(recname);
1476 ROTableRow recrow(rectab_);
1477 for (uInt idx = 0; idx < rectab_.nrow(); idx++){
1478 const TableRecord& rec = recrow.get(idx);
1479 if (rec.asString("spectralWindowId") == asdmSpw){
1480 const Array<Double> loarr = rec.asArrayDouble("freqLO");
1481 lo1Freq_ = loarr(IPosition(1,0));
1482 os << "Found LO1 Frequency in " << recname << ": "
1483 << lo1Freq_ << " [Hz]" << LogIO::POST;
1484 return true;
1485 }
1486 }
1487 os << LogIO::WARN << "Could not find " << asdmSpw << " in " << recname
1488 << LogIO::POST;
1489 return false;
1490};
1491
1492} //namespace asap
Note: See TracBrowser for help on using the repository browser.