source: trunk/src/STSideBandSep.cpp @ 2856

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

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