source: trunk/src/STSideBandSep.cpp@ 2871

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

New Development: No

JIRA Issue: Yes (CAS-4141)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asap.sbseparator

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