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