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.