source: trunk/src/STSideBandSep.cpp@ 2854

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

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