1 | // C++ Interface: STSideBandSep |
---|
2 | // |
---|
3 | // Description: |
---|
4 | // A class to invoke sideband separation of Scantable |
---|
5 | // |
---|
6 | // Author: Kana Sugimoto <kana.sugi@nao.ac.jp>, (C) 2012 |
---|
7 | // |
---|
8 | // Copyright: See COPYING file that comes with this distribution |
---|
9 | // |
---|
10 | // |
---|
11 | |
---|
12 | // STL |
---|
13 | #include <ctype.h> |
---|
14 | |
---|
15 | // cascore |
---|
16 | #include <casa/OS/File.h> |
---|
17 | #include <casa/Logging/LogIO.h> |
---|
18 | #include <casa/Quanta/QuantumHolder.h> |
---|
19 | |
---|
20 | #include <measures/Measures/MFrequency.h> |
---|
21 | #include <measures/Measures/MCFrequency.h> |
---|
22 | |
---|
23 | #include <tables/Tables/TableRow.h> |
---|
24 | #include <tables/Tables/TableRecord.h> |
---|
25 | #include <tables/Tables/TableVector.h> |
---|
26 | |
---|
27 | // asap |
---|
28 | #include "STGrid.h" |
---|
29 | #include "STMath.h" |
---|
30 | #include "MathUtils.h" |
---|
31 | #include "STSideBandSep.h" |
---|
32 | |
---|
33 | using namespace std ; |
---|
34 | using namespace casa ; |
---|
35 | using namespace asap ; |
---|
36 | |
---|
37 | // #ifndef KS_DEBUG |
---|
38 | // #define KS_DEBUG |
---|
39 | // #endif |
---|
40 | |
---|
41 | namespace asap { |
---|
42 | |
---|
43 | // constructors |
---|
44 | STSideBandSep::STSideBandSep(const vector<string> &names) |
---|
45 | { |
---|
46 | LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE)); |
---|
47 | os << "Setting scantable names to process." << LogIO::POST ; |
---|
48 | // Set file names |
---|
49 | ntable_ = names.size(); |
---|
50 | infileList_.resize(ntable_); |
---|
51 | for (unsigned int i = 0; i < ntable_; i++){ |
---|
52 | if (!checkFile(names[i], "d")) |
---|
53 | throw( AipsError("File does not exist") ); |
---|
54 | infileList_[i] = names[i]; |
---|
55 | } |
---|
56 | intabList_.resize(0); |
---|
57 | |
---|
58 | init(); |
---|
59 | |
---|
60 | {// Summary |
---|
61 | os << ntable_ << " files are set: ["; |
---|
62 | for (unsigned int i = 0; i < ntable_; i++) { |
---|
63 | os << " '" << infileList_[i] << "' "; |
---|
64 | if (i != ntable_-1) os << ","; |
---|
65 | } |
---|
66 | os << "] " << LogIO::POST; |
---|
67 | } |
---|
68 | }; |
---|
69 | |
---|
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(); |
---|
83 | tp_ = intabList_[0]->table().tableType(); |
---|
84 | |
---|
85 | os << ntable_ << " tables are set." << LogIO::POST; |
---|
86 | }; |
---|
87 | |
---|
88 | STSideBandSep::~STSideBandSep() |
---|
89 | { |
---|
90 | }; |
---|
91 | |
---|
92 | void STSideBandSep::init() |
---|
93 | { |
---|
94 | // frequency setup |
---|
95 | sigIfno_= -1; |
---|
96 | ftol_ = -1; |
---|
97 | solFrame_ = MFrequency::N_Types; |
---|
98 | // shifts |
---|
99 | initshift(); |
---|
100 | // direction tolerance |
---|
101 | xtol_ = ytol_ = 9.69627e-6; // 2arcsec |
---|
102 | // solution parameters |
---|
103 | otherside_ = false; |
---|
104 | doboth_ = false; |
---|
105 | rejlimit_ = 0.2; |
---|
106 | // LO1 values |
---|
107 | lo1Freq_ = -1; |
---|
108 | loTime_ = -1; |
---|
109 | loDir_ = ""; |
---|
110 | // Default LO frame is TOPO |
---|
111 | loFrame_ = MFrequency::TOPO; |
---|
112 | // scantable storage |
---|
113 | tp_ = Table::Memory; |
---|
114 | }; |
---|
115 | |
---|
116 | void STSideBandSep::initshift() |
---|
117 | { |
---|
118 | // shifts |
---|
119 | nshift_ = 0; |
---|
120 | nchan_ = 0; |
---|
121 | sigShift_.resize(0); |
---|
122 | imgShift_.resize(0); |
---|
123 | tableList_.resize(0); |
---|
124 | }; |
---|
125 | |
---|
126 | void STSideBandSep::setFrequency(const int ifno, |
---|
127 | const string freqtol, |
---|
128 | const string frame) |
---|
129 | { |
---|
130 | LogIO os(LogOrigin("STSideBandSep","setFrequency()", WHERE)); |
---|
131 | |
---|
132 | initshift(); |
---|
133 | |
---|
134 | // IFNO |
---|
135 | sigIfno_ = ifno; |
---|
136 | |
---|
137 | // Frequency tolerance |
---|
138 | Quantum<Double> qftol; |
---|
139 | readQuantity(qftol, String(freqtol)); |
---|
140 | if (!qftol.getUnit().empty()){ |
---|
141 | // make sure the quantity is frequency |
---|
142 | if (qftol.getFullUnit().getValue() != Unit("Hz").getValue()) |
---|
143 | throw( AipsError("Invalid quantity for frequency tolerance.") ); |
---|
144 | qftol.convert("Hz"); |
---|
145 | } |
---|
146 | ftol_ = qftol; |
---|
147 | |
---|
148 | // Frequency Frame |
---|
149 | if (!frame.empty()){ |
---|
150 | MFrequency::Types mft; |
---|
151 | if (!MFrequency::getType(mft, frame)) |
---|
152 | throw( AipsError("Invalid frame type.") ); |
---|
153 | solFrame_ = mft; |
---|
154 | } else { |
---|
155 | solFrame_ = MFrequency::N_Types; |
---|
156 | } |
---|
157 | |
---|
158 | {// Summary |
---|
159 | const String sframe = ( (solFrame_ == MFrequency::N_Types) ? |
---|
160 | "table frame" : |
---|
161 | MFrequency::showType(solFrame_) ); |
---|
162 | os << "Frequency setup to search IF group: " |
---|
163 | << "IFNO of table[0] = " << sigIfno_ |
---|
164 | << " , Freq tolerance = " << ftol_.getValue() << " [ " |
---|
165 | << (ftol_.getUnit().empty() ? "channel" : ftol_.getUnit() ) |
---|
166 | << " ] (in " << sframe <<")" << LogIO::POST; |
---|
167 | } |
---|
168 | }; |
---|
169 | |
---|
170 | |
---|
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 | } |
---|
224 | os << " ) [channels]" << LogIO::POST; |
---|
225 | } |
---|
226 | }; |
---|
227 | |
---|
228 | void STSideBandSep::setThreshold(const double limit) |
---|
229 | { |
---|
230 | LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE)); |
---|
231 | if (limit < 0) |
---|
232 | throw( AipsError("Rejection limit should be a positive number.") ); |
---|
233 | |
---|
234 | rejlimit_ = limit; |
---|
235 | os << "Rejection limit is set to " << rejlimit_ << LogIO::POST; |
---|
236 | }; |
---|
237 | |
---|
238 | void STSideBandSep::separate(string outname) |
---|
239 | { |
---|
240 | LogIO os(LogOrigin("STSideBandSep","separate()", WHERE)); |
---|
241 | if (outname.empty()) |
---|
242 | outname = "sbseparated.asap"; |
---|
243 | |
---|
244 | // Set up a goup of IFNOs in the list of scantables within |
---|
245 | // the frequency tolerance and make them a list. |
---|
246 | nshift_ = setupShift(); |
---|
247 | if (nshift_ < 2) |
---|
248 | throw( AipsError("At least 2 IFs are necessary for convolution.") ); |
---|
249 | // Grid scantable and generate output tables |
---|
250 | ScantableWrapper gridst = gridTable(); |
---|
251 | sigTab_p = gridst.getCP(); |
---|
252 | if (doboth_) |
---|
253 | imgTab_p = gridst.copy().getCP(); |
---|
254 | vector<unsigned int> remRowIds; |
---|
255 | remRowIds.resize(0); |
---|
256 | Matrix<float> specMat(nchan_, nshift_); |
---|
257 | Matrix<bool> flagMat(nchan_, nshift_); |
---|
258 | vector<float> sigSpec(nchan_), imgSpec(nchan_); |
---|
259 | Vector<bool> flagVec(nchan_); |
---|
260 | vector<uInt> tabIdvec; |
---|
261 | |
---|
262 | //Generate FFTServer |
---|
263 | fftsf.resize(IPosition(1, nchan_), FFTEnums::REALTOCOMPLEX); |
---|
264 | fftsi.resize(IPosition(1, nchan_), FFTEnums::COMPLEXTOREAL); |
---|
265 | |
---|
266 | /// Loop over sigTab_p and separate sideband |
---|
267 | for (int irow = 0; irow < sigTab_p->nrow(); irow++){ |
---|
268 | tabIdvec.resize(0); |
---|
269 | const int polId = sigTab_p->getPol(irow); |
---|
270 | const int beamId = sigTab_p->getBeam(irow); |
---|
271 | const vector<double> dir = sigTab_p->getDirectionVector(irow); |
---|
272 | // Get a set of spectra to solve |
---|
273 | if (!getSpectraToSolve(polId, beamId, dir[0], dir[1], |
---|
274 | specMat, flagMat, tabIdvec)){ |
---|
275 | remRowIds.push_back(irow); |
---|
276 | #ifdef KS_DEBUG |
---|
277 | cout << "no matching row found. skipping row = " << irow << endl; |
---|
278 | #endif |
---|
279 | continue; |
---|
280 | } |
---|
281 | // Solve signal sideband |
---|
282 | sigSpec = solve(specMat, tabIdvec, true); |
---|
283 | sigTab_p->setSpectrum(sigSpec, irow); |
---|
284 | if (sigTab_p->isAllChannelsFlagged(irow)){ |
---|
285 | // unflag the spectrum since there should be some valid data |
---|
286 | sigTab_p->flagRow(vector<uInt>(irow), true); |
---|
287 | // need to unflag whole channels anyway |
---|
288 | sigTab_p->flag(irow, vector<bool>(), true); |
---|
289 | } |
---|
290 | // apply channel flag |
---|
291 | flagVec = collapseFlag(flagMat, tabIdvec, true); |
---|
292 | //boolVec = !boolVec; // flag |
---|
293 | vector<bool> tmpflag; |
---|
294 | flagVec.tovector(tmpflag); |
---|
295 | sigTab_p->flag(irow, tmpflag, false); |
---|
296 | |
---|
297 | // Solve image sideband |
---|
298 | if (doboth_) { |
---|
299 | imgSpec = solve(specMat, tabIdvec, false); |
---|
300 | imgTab_p->setSpectrum(imgSpec, irow); |
---|
301 | if (imgTab_p->isAllChannelsFlagged(irow)){ |
---|
302 | // unflag the spectrum since there should be some valid data |
---|
303 | imgTab_p->flagRow(vector<uInt>(irow), true); |
---|
304 | // need to unflag whole channels anyway |
---|
305 | imgTab_p->flag(irow, vector<bool>(), true); |
---|
306 | } |
---|
307 | // apply channel flag |
---|
308 | flagVec = collapseFlag(flagMat, tabIdvec, false); |
---|
309 | //boolVec = !boolVec; // flag |
---|
310 | flagVec.tovector(tmpflag); |
---|
311 | imgTab_p->flag(irow, tmpflag, false); |
---|
312 | } |
---|
313 | } // end of row loop |
---|
314 | |
---|
315 | // Remove or flag rows without relevant data from gridded tables |
---|
316 | if (remRowIds.size() > 0) { |
---|
317 | const size_t nrem = remRowIds.size(); |
---|
318 | if ( sigTab_p->table().canRemoveRow() ) { |
---|
319 | sigTab_p->table().removeRow(remRowIds); |
---|
320 | os << "Removing " << nrem << " rows from the signal band table" |
---|
321 | << LogIO::POST; |
---|
322 | } else { |
---|
323 | sigTab_p->flagRow(remRowIds, false); |
---|
324 | os << "Cannot remove rows from the signal band table. Flagging " |
---|
325 | << nrem << " rows" << LogIO::POST; |
---|
326 | } |
---|
327 | |
---|
328 | if (doboth_) { |
---|
329 | if ( imgTab_p->table().canRemoveRow() ) { |
---|
330 | imgTab_p->table().removeRow(remRowIds); |
---|
331 | os << "Removing " << nrem << " rows from the image band table" |
---|
332 | << LogIO::POST; |
---|
333 | } else { |
---|
334 | imgTab_p->flagRow(remRowIds, false); |
---|
335 | os << "Cannot remove rows from the image band table. Flagging " |
---|
336 | << nrem << " rows" << LogIO::POST; |
---|
337 | } |
---|
338 | } |
---|
339 | } |
---|
340 | |
---|
341 | // Finally, save tables on disk |
---|
342 | if (outname.size() ==0) |
---|
343 | outname = "sbseparated.asap"; |
---|
344 | const string sigName = outname + ".signalband"; |
---|
345 | os << "Saving SIGNAL sideband table: " << sigName << LogIO::POST; |
---|
346 | sigTab_p->makePersistent(sigName); |
---|
347 | if (doboth_) { |
---|
348 | solveImageFrequency(); |
---|
349 | const string imgName = outname + ".imageband"; |
---|
350 | os << "Saving IMAGE sideband table: " << sigName << LogIO::POST; |
---|
351 | imgTab_p->makePersistent(imgName); |
---|
352 | } |
---|
353 | |
---|
354 | }; |
---|
355 | |
---|
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; |
---|
399 | #ifdef KS_DEBUG |
---|
400 | cout << "Frequencies in the first table: freq0 = " << basech0 << ", incr = " << baseinc << ", nch = " << basench << endl; |
---|
401 | #endif |
---|
402 | } |
---|
403 | if (ftol_.getUnit().empty()) { |
---|
404 | // tolerance in unit of channels |
---|
405 | ftolval = abs(ftol_.getValue() * baseinc); |
---|
406 | } |
---|
407 | else { |
---|
408 | ftolval = abs(ftol_.getValue("Hz")); |
---|
409 | } |
---|
410 | inctolval = abs(baseinc/(double) basench); |
---|
411 | const string poltype0 = stab->getPolType(); |
---|
412 | #ifdef KS_DEBUG |
---|
413 | cout << "Looking for pairs that satisfies freqtol = " << ftolval << ", inctol = " << inctolval << endl; |
---|
414 | #endif |
---|
415 | |
---|
416 | // Initialize shift values |
---|
417 | initshift(); |
---|
418 | |
---|
419 | // const bool setImg = ( doboth_ && (imgShift_.size() == 0) ); |
---|
420 | const bool setImg = (imgShift_.size() == 0); |
---|
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) ){ |
---|
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 |
---|
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); |
---|
473 | sigShift_.push_back(-chShift); |
---|
474 | if (setImg) |
---|
475 | imgShift_.push_back(chShift); |
---|
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); |
---|
575 | // gridder.setFunc("box", 1); // convsupport=1 fails |
---|
576 | gridder.setFunc("box"); |
---|
577 | gridder.setWeight("uniform"); |
---|
578 | #ifdef KS_DEBUG |
---|
579 | cout << "Grid parameter summary: " << endl; |
---|
580 | cout << "- IF = " << sigIfno_ << endl; |
---|
581 | cout << "- center = " << scenter << "\n" |
---|
582 | << "- npix = (" << nx << ", " << ny << ")\n" |
---|
583 | << "- cell = (" << scellx << ", " << scelly << ")" << endl; |
---|
584 | #endif |
---|
585 | gridder.grid(); |
---|
586 | const int itp = (tp_ == Table::Memory ? 0 : 1); |
---|
587 | ScantableWrapper gtab = gridder.getResultAsScantable(itp); |
---|
588 | // WORKAROUND : Shift TIME for proper pointing resolution in future imaging. |
---|
589 | shiftTimeInGriddedST(gtab.getCP()); |
---|
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 | |
---|
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(); |
---|
650 | ScalarColumn<Double> mjdCol( tab, "TIME"); |
---|
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 | |
---|
684 | // bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId, |
---|
685 | // const double dirX, const double dirY, |
---|
686 | // Matrix<float> &specMat, vector<uInt> &tabIdvec) |
---|
687 | bool STSideBandSep::getSpectraToSolve(const int polId, const int beamId, |
---|
688 | const double dirX, const double dirY, |
---|
689 | Matrix<float> &specMat, |
---|
690 | Matrix<bool> &flagMat, |
---|
691 | vector<uInt> &tabIdvec) |
---|
692 | { |
---|
693 | LogIO os(LogOrigin("STSideBandSep","getSpectraToSolve()", WHERE)); |
---|
694 | |
---|
695 | tabIdvec.resize(0); |
---|
696 | specMat.resize(nchan_, nshift_); |
---|
697 | Vector<float> spec; |
---|
698 | Vector<bool> boolVec; |
---|
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) { |
---|
749 | // STMath::average also merges FLAGTRA and FLAGROW |
---|
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 | } |
---|
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); |
---|
771 | // Mask is true if the data is valid (!flag) |
---|
772 | vector<bool> mask = avetab_p->getMask(0); |
---|
773 | mathutil::doZeroOrderInterpolation(tmpSpec, mask); |
---|
774 | spec.reference(specMat.column(nspec)); |
---|
775 | spec = tmpSpec; |
---|
776 | boolVec.reference(flagMat.column(nspec)); |
---|
777 | boolVec = mask; // cast std::vector to casa::Vector |
---|
778 | boolVec = !boolVec; |
---|
779 | tabIdvec.push_back((uInt) itab); |
---|
780 | nspec++; |
---|
781 | //Liberate from reference |
---|
782 | spec.unique(); |
---|
783 | boolVec.unique(); |
---|
784 | } // end of table loop |
---|
785 | // Check the number of selected spectra and resize matrix. |
---|
786 | if (nspec != nshift_){ |
---|
787 | //shiftSpecmat.resize(nchan_, nspec, true); |
---|
788 | specMat.resize(nchan_, nspec, true); |
---|
789 | flagMat.resize(nchan_, nspec, true); |
---|
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 | |
---|
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 | } |
---|
842 | outflag = accflag; |
---|
843 | // Shift back Flag |
---|
844 | //cout << "Shifting FLAG back to " << thisShift->at(0) << " channels" << endl; |
---|
845 | //shiftFlag(accflag, thisShift->at(0), outflag); |
---|
846 | |
---|
847 | return outflag; |
---|
848 | } |
---|
849 | |
---|
850 | |
---|
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_; |
---|
881 | otherShift = &imgShift_; |
---|
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; |
---|
891 | uInt shiftId; |
---|
892 | for (uInt i = 0 ; i < nspec; i++) { |
---|
893 | shiftId = tabIdvec[i]; |
---|
894 | spshift[i] = otherShift->at(shiftId) - thisShift->at(shiftId); |
---|
895 | tempshift = - thisShift->at(shiftId); |
---|
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 |
---|
944 | cout << "Start shifting spectrum for " << shift << " channels" << endl; |
---|
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; |
---|
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 |
---|
971 | } |
---|
972 | }; |
---|
973 | |
---|
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){ |
---|
999 | // the shift is almost lchan |
---|
1000 | ruse = false; |
---|
1001 | luse = true; |
---|
1002 | } |
---|
1003 | if (rweight > 1-tolerance){ |
---|
1004 | // the shift is almost rchan |
---|
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); |
---|
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 |
---|
1026 | } |
---|
1027 | }; |
---|
1028 | |
---|
1029 | |
---|
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 | |
---|
1041 | #ifdef KS_DEBUG |
---|
1042 | float minval, maxval; |
---|
1043 | #endif |
---|
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, |
---|
1216 | const double reftime, const string refdir) |
---|
1217 | { |
---|
1218 | Quantum<Double> qfreq; |
---|
1219 | readQuantity(qfreq, String(lo1)); |
---|
1220 | lo1Freq_ = qfreq.getValue("Hz"); |
---|
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 | |
---|
1236 | void STSideBandSep::setLO1Root(string name) |
---|
1237 | { |
---|
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)) { |
---|
1242 | throw(AipsError("File does not exist")); |
---|
1243 | } |
---|
1244 | if (name[(name.size()-1)] == '/') |
---|
1245 | name = name.substr(0,(name.size()-2)); |
---|
1246 | |
---|
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 | } |
---|
1258 | |
---|
1259 | #ifdef KS_DEBUG |
---|
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; |
---|
1265 | #endif |
---|
1266 | }; |
---|
1267 | |
---|
1268 | |
---|
1269 | void STSideBandSep::solveImageFrequency() |
---|
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 | |
---|
1277 | // Convert frequency REFVAL to the value in frame of LO. |
---|
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); |
---|
1291 | int nChan = imgTab_p->nchan(imgTab_p->getIF(0)); |
---|
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() ; |
---|
1333 | //md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch); |
---|
1334 | md = MDirection(ra, dec, epoch); |
---|
1335 | } |
---|
1336 | MeasFrame mframe( me, mp, md ); |
---|
1337 | MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe)); |
---|
1338 | MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe)); |
---|
1339 | // Convert refval to loFrame_ |
---|
1340 | double sigrefval; |
---|
1341 | if (tabframe == loFrame_) |
---|
1342 | sigrefval = refval; |
---|
1343 | else |
---|
1344 | sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue(); |
---|
1345 | |
---|
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; |
---|
1364 | if (!getLo1FromScanTab(tableList_[0], sigrefval, refpix, increment, nChan)) { |
---|
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 | } |
---|
1372 | |
---|
1373 | // LO1 should now be ready. |
---|
1374 | if (lo1Freq_ < 0.) |
---|
1375 | throw(AipsError("Got negative LO1 Frequency")); |
---|
1376 | |
---|
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; |
---|
1388 | |
---|
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 | } |
---|
1400 | |
---|
1401 | // Calculate image band incr and refval in loFrame_ |
---|
1402 | Double imgincr = -increment; |
---|
1403 | Double imgrefval = 2 * lo1Freq_ - sigrefval; |
---|
1404 | Double imgrefval_tab = imgrefval; |
---|
1405 | // Convert imgrefval back to table base frame |
---|
1406 | if (tabframe != loFrame_) |
---|
1407 | imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue(); |
---|
1408 | // Set new frequencies table |
---|
1409 | uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr); |
---|
1410 | // Update FREQ_ID in table. |
---|
1411 | freqIdVec = fIDnew; |
---|
1412 | |
---|
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 | } |
---|
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 |
---|
1436 | switch (tolower(type[0])) { |
---|
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 | |
---|
1449 | bool STSideBandSep::getLo1FromAsdm(const string asdmname, |
---|
1450 | const double /*refval*/, |
---|
1451 | const double /*refpix*/, |
---|
1452 | const double /*increment*/, |
---|
1453 | const int /*nChan*/) |
---|
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 | |
---|
1467 | bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab, |
---|
1468 | const double refval, |
---|
1469 | const double refpix, |
---|
1470 | const double increment, |
---|
1471 | const int nChan) |
---|
1472 | { |
---|
1473 | LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE)); |
---|
1474 | // Check for relevant tables. |
---|
1475 | const TableRecord &rec = scantab->table().keywordSet() ; |
---|
1476 | String spwname, recname; |
---|
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 { |
---|
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 | } |
---|
1488 | if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) { |
---|
1489 | throw(AipsError("Could not find relevant subtables in MS")); |
---|
1490 | } |
---|
1491 | // Get root MS name |
---|
1492 | string msname; |
---|
1493 | const String recsuff = "/ASDM_RECEIVER"; |
---|
1494 | String::size_type pos; |
---|
1495 | pos = recname.size()-recsuff.size(); |
---|
1496 | if (recname.substr(pos) == recsuff) |
---|
1497 | msname = recname.substr(0, pos); |
---|
1498 | else |
---|
1499 | throw(AipsError("Internal error in parsing table name from a scantable keyword.")); |
---|
1500 | |
---|
1501 | if (!checkFile(msname)) |
---|
1502 | throw(AipsError("Internal error in parsing MS name from a scantable keyword.")); |
---|
1503 | |
---|
1504 | return getLo1FromAsisTab(msname, refval, refpix, increment, nChan); |
---|
1505 | |
---|
1506 | }; |
---|
1507 | |
---|
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 | } |
---|
1522 | |
---|
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 | |
---|
1587 | } //namespace asap |
---|