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