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 "STSideBandSep.h"
|
---|
29 |
|
---|
30 | using namespace std ;
|
---|
31 | using namespace casa ;
|
---|
32 | using namespace asap ;
|
---|
33 |
|
---|
34 | #ifndef KS_DEBUG
|
---|
35 | #define KS_DEBUG
|
---|
36 | #endif
|
---|
37 |
|
---|
38 | namespace asap {
|
---|
39 |
|
---|
40 | // constructors
|
---|
41 | STSideBandSep::STSideBandSep(const vector<string> &names)
|
---|
42 | {
|
---|
43 | LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
|
---|
44 | os << "Setting scantable names to process." << LogIO::POST ;
|
---|
45 | // Set file names
|
---|
46 | ntable_ = names.size();
|
---|
47 | infileList_.resize(ntable_);
|
---|
48 | for (unsigned int i = 0; i < ntable_; i++){
|
---|
49 | if (!checkFile(names[i], "d"))
|
---|
50 | throw( AipsError("File does not exist") );
|
---|
51 | infileList_[i] = names[i];
|
---|
52 | }
|
---|
53 | intabList_.resize(0);
|
---|
54 |
|
---|
55 | init();
|
---|
56 |
|
---|
57 | {// Summary
|
---|
58 | os << ntable_ << " files are set: [";
|
---|
59 | for (unsigned int i = 0; i < ntable_; i++) {
|
---|
60 | os << " '" << infileList_[i] << "' ";
|
---|
61 | if (i != ntable_-1) os << ",";
|
---|
62 | }
|
---|
63 | os << "] " << LogIO::POST;
|
---|
64 | }
|
---|
65 | };
|
---|
66 |
|
---|
67 | STSideBandSep::STSideBandSep(const vector<ScantableWrapper> &tables)
|
---|
68 | {
|
---|
69 | LogIO os(LogOrigin("STSideBandSep","STSideBandSep()", WHERE));
|
---|
70 | os << "Setting list of scantables to process." << LogIO::POST ;
|
---|
71 | // Set file names
|
---|
72 | ntable_ = tables.size();
|
---|
73 | intabList_.resize(ntable_);
|
---|
74 | for (unsigned int i = 0; i < ntable_; i++){
|
---|
75 | intabList_[i] = tables[i].getCP();
|
---|
76 | }
|
---|
77 | infileList_.resize(0);
|
---|
78 |
|
---|
79 | init();
|
---|
80 |
|
---|
81 | os << ntable_ << " tables are set." << LogIO::POST;
|
---|
82 | };
|
---|
83 |
|
---|
84 | STSideBandSep::~STSideBandSep()
|
---|
85 | {
|
---|
86 | #ifdef KS_DEBUG
|
---|
87 | cout << "Destructor ~STSideBandSep()" << endl;
|
---|
88 | #endif
|
---|
89 | };
|
---|
90 |
|
---|
91 | void STSideBandSep::init()
|
---|
92 | {
|
---|
93 | // frequency setup
|
---|
94 | sigIfno_= 0;
|
---|
95 | ftol_ = -1;
|
---|
96 | solFrame_ = MFrequency::N_Types;
|
---|
97 | // shifts
|
---|
98 | initshift();
|
---|
99 | // direction tolerance
|
---|
100 | xtol_ = ytol_ = 9.69627e-6; // 2arcsec
|
---|
101 | // solution parameters
|
---|
102 | otherside_ = false;
|
---|
103 | doboth_ = false;
|
---|
104 | rejlimit_ = 0.2;
|
---|
105 | // LO1 values
|
---|
106 | lo1Freq_ = -1;
|
---|
107 | loTime_ = -1;
|
---|
108 | loDir_ = "";
|
---|
109 | // Default LO frame is TOPO
|
---|
110 | loFrame_ = MFrequency::TOPO;
|
---|
111 | };
|
---|
112 |
|
---|
113 | void STSideBandSep::initshift()
|
---|
114 | {
|
---|
115 | // shifts
|
---|
116 | nshift_ = 0;
|
---|
117 | nchan_ = 0;
|
---|
118 | sigShift_.resize(0);
|
---|
119 | imgShift_.resize(0);
|
---|
120 | tableList_.resize(0);
|
---|
121 | };
|
---|
122 |
|
---|
123 | void STSideBandSep::setFrequency(const unsigned int ifno,
|
---|
124 | const string freqtol,
|
---|
125 | const string frame)
|
---|
126 | {
|
---|
127 | LogIO os(LogOrigin("STSideBandSep","setFrequency()", WHERE));
|
---|
128 |
|
---|
129 | initshift();
|
---|
130 |
|
---|
131 | // IFNO
|
---|
132 | sigIfno_ = ifno;
|
---|
133 |
|
---|
134 | // Frequency tolerance
|
---|
135 | Quantum<Double> qftol;
|
---|
136 | readQuantity(qftol, String(freqtol));
|
---|
137 | if (!qftol.getUnit().empty()){
|
---|
138 | // make sure the quantity is frequency
|
---|
139 | if (qftol.getFullUnit().getValue() != Unit("Hz").getValue())
|
---|
140 | throw( AipsError("Invalid quantity for frequency tolerance.") );
|
---|
141 | qftol.convert("Hz");
|
---|
142 | }
|
---|
143 | ftol_ = qftol;
|
---|
144 |
|
---|
145 | // Frequency Frame
|
---|
146 | if (!frame.empty()){
|
---|
147 | MFrequency::Types mft;
|
---|
148 | if (!MFrequency::getType(mft, frame))
|
---|
149 | throw( AipsError("Invalid frame type.") );
|
---|
150 | solFrame_ = mft;
|
---|
151 | } else {
|
---|
152 | solFrame_ = MFrequency::N_Types;
|
---|
153 | }
|
---|
154 |
|
---|
155 | {// Summary
|
---|
156 | const String sframe = ( (solFrame_ == MFrequency::N_Types) ?
|
---|
157 | "table frame" :
|
---|
158 | MFrequency::showType(solFrame_) );
|
---|
159 | os << "Frequency setup to search IF group: "
|
---|
160 | << "IFNO of table[0] = " << sigIfno_
|
---|
161 | << " , Freq tolerance = " << ftol_.getValue() << " [ "
|
---|
162 | << (ftol_.getUnit().empty() ? "channel" : ftol_.getUnit() )
|
---|
163 | << " ] (in " << sframe <<")" << LogIO::POST;
|
---|
164 | }
|
---|
165 | };
|
---|
166 |
|
---|
167 |
|
---|
168 | void STSideBandSep::setDirTolerance(const vector<string> dirtol)
|
---|
169 | {
|
---|
170 | LogIO os(LogOrigin("STSideBandSep","setDirTolerance()", WHERE));
|
---|
171 | Quantum<Double> qcell;
|
---|
172 | if ( (dirtol.size() == 1) && !dirtol[0].empty() ) {
|
---|
173 | readQuantity(qcell, String(dirtol[0]));
|
---|
174 | if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
|
---|
175 | xtol_ = ytol_ = qcell.getValue("rad");
|
---|
176 | else
|
---|
177 | throw( AipsError("Invalid unit for direction tolerance.") );
|
---|
178 | }
|
---|
179 | else if (dirtol.size() > 1) {
|
---|
180 | if ( dirtol[0].empty() && dirtol[1].empty() )
|
---|
181 | throw( AipsError("Direction tolerance is empty.") );
|
---|
182 | if ( !dirtol[0].empty() ) {
|
---|
183 | readQuantity(qcell, String(dirtol[0]));
|
---|
184 | if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
|
---|
185 | xtol_ = qcell.getValue("rad");
|
---|
186 | else
|
---|
187 | throw( AipsError("Invalid unit for direction tolerance.") );
|
---|
188 | }
|
---|
189 | if ( !dirtol[1].empty() ) {
|
---|
190 | readQuantity(qcell, String(dirtol[1]));
|
---|
191 | if (qcell.getFullUnit().getValue() == Unit("rad").getValue())
|
---|
192 | ytol_ = qcell.getValue("rad");
|
---|
193 | else
|
---|
194 | throw( AipsError("Invalid unit for direction tolerance.") );
|
---|
195 | }
|
---|
196 | else {
|
---|
197 | ytol_ = xtol_;
|
---|
198 | }
|
---|
199 | }
|
---|
200 | else throw( AipsError("Invalid direction tolerance.") );
|
---|
201 |
|
---|
202 | os << "Direction tolerance: ( "
|
---|
203 | << xtol_ << " , " << ytol_ << " ) [rad]" << LogIO::POST;
|
---|
204 | };
|
---|
205 |
|
---|
206 | void STSideBandSep::setShift(const vector<double> &shift)
|
---|
207 | {
|
---|
208 | LogIO os(LogOrigin("STSideBandSep","setShift()", WHERE));
|
---|
209 | imgShift_.resize(shift.size());
|
---|
210 | for (unsigned int i = 0; i < shift.size(); i++)
|
---|
211 | imgShift_[i] = shift[i];
|
---|
212 |
|
---|
213 | if (imgShift_.size() == 0) {
|
---|
214 | os << "Channel shifts are cleared." << LogIO::POST;
|
---|
215 | } else {
|
---|
216 | os << "Channel shifts of image sideband are set: ( ";
|
---|
217 | for (unsigned int i = 0; i < imgShift_.size(); i++) {
|
---|
218 | os << imgShift_[i];
|
---|
219 | if (i != imgShift_.size()-1) os << " , ";
|
---|
220 | }
|
---|
221 | os << " ) [channel]" << LogIO::POST;
|
---|
222 | }
|
---|
223 | };
|
---|
224 |
|
---|
225 | void STSideBandSep::setThreshold(const double limit)
|
---|
226 | {
|
---|
227 | LogIO os(LogOrigin("STSideBandSep","setThreshold()", WHERE));
|
---|
228 | if (limit < 0)
|
---|
229 | throw( AipsError("Rejection limit should be positive number.") );
|
---|
230 |
|
---|
231 | rejlimit_ = limit;
|
---|
232 | os << "Rejection limit is set to " << rejlimit_;
|
---|
233 | };
|
---|
234 |
|
---|
235 | // Temporal function to set scantable of image sideband
|
---|
236 | void STSideBandSep::setImageTable(const ScantableWrapper &s)
|
---|
237 | {
|
---|
238 | #ifdef KS_DEBUG
|
---|
239 | cout << "STSideBandSep::setImageTable" << endl;
|
---|
240 | cout << "got image table nrow = " << s.nrow() << endl;
|
---|
241 | #endif
|
---|
242 | imgTab_p = s.getCP();
|
---|
243 | AlwaysAssert(!imgTab_p.null(),AipsError);
|
---|
244 |
|
---|
245 | };
|
---|
246 |
|
---|
247 | //
|
---|
248 | void STSideBandSep::setLO1(const double lo1, const string frame,
|
---|
249 | const double reftime, const string refdir)
|
---|
250 | {
|
---|
251 | lo1Freq_ = lo1;
|
---|
252 | MFrequency::getType(loFrame_, frame);
|
---|
253 | loTime_ = reftime;
|
---|
254 | loDir_ = refdir;
|
---|
255 |
|
---|
256 | #ifdef KS_DEBUG
|
---|
257 | cout << "STSideBandSep::setLO1" << endl;
|
---|
258 | if (lo1Freq_ > 0.)
|
---|
259 | cout << "lo1 = " << lo1Freq_ << " [Hz] (" << frame << ")" << endl;
|
---|
260 | if (loTime_ > 0.)
|
---|
261 | cout << "ref time = " << loTime_ << " [day]" << endl;
|
---|
262 | if (!loDir_.empty())
|
---|
263 | cout << "ref direction = " << loDir_ << " [day]" << endl;
|
---|
264 | #endif
|
---|
265 | };
|
---|
266 |
|
---|
267 | void STSideBandSep::setLO1Root(string name)
|
---|
268 | {
|
---|
269 | LogIO os(LogOrigin("STSideBandSep","setLO1Root()", WHERE));
|
---|
270 | os << "Searching for '" << name << "'..." << LogIO::POST;
|
---|
271 | // Check for existance of the file
|
---|
272 | if (!checkFile(name)) {
|
---|
273 | throw(AipsError("File does not exist"));
|
---|
274 | }
|
---|
275 | if (name[(name.size()-1)] == '/')
|
---|
276 | name = name.substr(0,(name.size()-2));
|
---|
277 |
|
---|
278 | if (checkFile(name+"/Receiver.xml", "file") &&
|
---|
279 | checkFile(name+"/SpectralWindow.xml", "file")){
|
---|
280 | os << "Found '" << name << "/Receiver.xml' ... got an ASDM name." << LogIO::POST;
|
---|
281 | asdmName_ = name;
|
---|
282 | } else if (checkFile(name+"/ASDM_RECEIVER") &&
|
---|
283 | checkFile(name+"/ASDM_SPECTRALWINDOW")){
|
---|
284 | os << "Found '" << name << "/ASDM_RECEIVER' ... got a Table name." << LogIO::POST;
|
---|
285 | asisName_ = name;
|
---|
286 | } else {
|
---|
287 | throw(AipsError("Invalid file name. Set an MS or ASDM name."));
|
---|
288 | }
|
---|
289 |
|
---|
290 | #ifdef KS_DEBUG
|
---|
291 | cout << "STSideBandSep::setLO1Root" << endl;
|
---|
292 | if (!asdmName_.empty())
|
---|
293 | cout << "asdm name = " << asdmName_ << endl;
|
---|
294 | if (!asisName_.empty())
|
---|
295 | cout << "MS name = " << asisName_ << endl;
|
---|
296 | #endif
|
---|
297 | };
|
---|
298 |
|
---|
299 | ///// TEMPORAL FUNCTION!!! /////
|
---|
300 | void STSideBandSep::setScanTb0(const ScantableWrapper &s){
|
---|
301 | st0_ = s.getCP();
|
---|
302 | };
|
---|
303 | ////////////////////////////////
|
---|
304 |
|
---|
305 | void STSideBandSep::solveImageFreqency()
|
---|
306 | {
|
---|
307 | #ifdef KS_DEBUG
|
---|
308 | cout << "STSideBandSep::solveImageFrequency" << endl;
|
---|
309 | #endif
|
---|
310 | LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
|
---|
311 | os << "Start calculating frequencies of image side band" << LogIO::POST;
|
---|
312 |
|
---|
313 | if (imgTab_p.null())
|
---|
314 | throw AipsError("STSideBandSep::solveImageFreqency - an image side band scantable should be set first");
|
---|
315 |
|
---|
316 | // Convert frequency REFVAL to the value in frame of LO.
|
---|
317 | // The code assumes that imgTab_p has only an IF and only a FREQ_ID
|
---|
318 | // is associated to an IFNO
|
---|
319 | // TODO: More complete Procedure would be
|
---|
320 | // 1. Get freq IDs associated to sigIfno_
|
---|
321 | // 2. Get freq information of the freq IDs
|
---|
322 | // 3. For each freqIDs, get freq infromation in TOPO and an LO1
|
---|
323 | // frequency and calculate image band frequencies.
|
---|
324 | STFrequencies freqTab_ = imgTab_p->frequencies();
|
---|
325 | // get the base frame of table
|
---|
326 | const MFrequency::Types tabframe = freqTab_.getFrame(true);
|
---|
327 | TableVector<uInt> freqIdVec( imgTab_p->table(), "FREQ_ID" );
|
---|
328 | // assuming single freqID per IFNO
|
---|
329 | uInt freqid = freqIdVec(0);
|
---|
330 | int nChan = imgTab_p->nchan(imgTab_p->getIF(0));
|
---|
331 | double refpix, refval, increment ;
|
---|
332 | freqTab_.getEntry(refpix, refval, increment, freqid);
|
---|
333 | //MFrequency sigrefval = MFrequency(MVFrequency(refval),tabframe);
|
---|
334 | // get freq infromation of sigIfno_ in loFrame_
|
---|
335 | const MPosition mp = imgTab_p->getAntennaPosition();
|
---|
336 | MEpoch me;
|
---|
337 | MDirection md;
|
---|
338 | if (loTime_ < 0.)
|
---|
339 | me = imgTab_p->getEpoch(-1);
|
---|
340 | else
|
---|
341 | me = MEpoch(MVEpoch(loTime_));
|
---|
342 | if (loDir_.empty()) {
|
---|
343 | ArrayColumn<Double> srcdirCol_;
|
---|
344 | srcdirCol_.attach(imgTab_p->table(), "SRCDIRECTION");
|
---|
345 | // Assuming J2000 and SRCDIRECTION in unit of rad
|
---|
346 | Quantum<Double> srcra = Quantum<Double>(srcdirCol_(0)(IPosition(1,0)), "rad");
|
---|
347 | Quantum<Double> srcdec = Quantum<Double>(srcdirCol_(0)(IPosition(1,1)), "rad");
|
---|
348 | md = MDirection(srcra, srcdec, MDirection::J2000);
|
---|
349 | //imgTab_p->getDirection(0);
|
---|
350 | } else {
|
---|
351 | // parse direction string
|
---|
352 | string::size_type pos0 = loDir_.find(" ");
|
---|
353 |
|
---|
354 | if (pos0 == string::npos) {
|
---|
355 | throw AipsError("bad string format in LO1 direction");
|
---|
356 | }
|
---|
357 | string::size_type pos1 = loDir_.find(" ", pos0+1);
|
---|
358 | String sepoch, sra, sdec;
|
---|
359 | if (pos1 != string::npos) {
|
---|
360 | sepoch = loDir_.substr(0, pos0);
|
---|
361 | sra = loDir_.substr(pos0+1, pos1-pos0);
|
---|
362 | sdec = loDir_.substr(pos1+1);
|
---|
363 | }
|
---|
364 | MDirection::Types epoch;
|
---|
365 | MDirection::getType(epoch, sepoch);
|
---|
366 | QuantumHolder qh ;
|
---|
367 | String err ;
|
---|
368 | qh.fromString( err, sra);
|
---|
369 | Quantum<Double> ra = qh.asQuantumDouble() ;
|
---|
370 | qh.fromString( err, sdec ) ;
|
---|
371 | Quantum<Double> dec = qh.asQuantumDouble() ;
|
---|
372 | //md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch);
|
---|
373 | md = MDirection(ra, dec, epoch);
|
---|
374 | }
|
---|
375 | MeasFrame mframe( me, mp, md );
|
---|
376 | MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe));
|
---|
377 | MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe));
|
---|
378 | // Convert refval to loFrame_
|
---|
379 | double sigrefval;
|
---|
380 | if (tabframe == loFrame_)
|
---|
381 | sigrefval = refval;
|
---|
382 | else
|
---|
383 | sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
|
---|
384 |
|
---|
385 |
|
---|
386 | // Check for the availability of LO1
|
---|
387 | if (lo1Freq_ > 0.) {
|
---|
388 | os << "Using user defined LO1 frequency." << LogIO::POST;
|
---|
389 | } else if (!asisName_.empty()) {
|
---|
390 | // MS name is set.
|
---|
391 | os << "Using user defined MS (asis): " << asisName_ << LogIO::POST;
|
---|
392 | if (!getLo1FromAsisTab(asisName_, sigrefval, refpix, increment, nChan)) {
|
---|
393 | throw AipsError("Failed to get LO1 frequency from MS");
|
---|
394 | }
|
---|
395 | } else if (!asdmName_.empty()) {
|
---|
396 | // ASDM name is set.
|
---|
397 | os << "Using user defined ASDM: " << asdmName_ << LogIO::POST;
|
---|
398 | if (!getLo1FromAsdm(asdmName_, sigrefval, refpix, increment, nChan)) {
|
---|
399 | throw AipsError("Failed to get LO1 frequency from ASDM");
|
---|
400 | }
|
---|
401 | } else {
|
---|
402 | // Try getting ASDM name from scantable header
|
---|
403 | os << "Try getting information from scantable header" << LogIO::POST;
|
---|
404 | if (!getLo1FromScanTab(st0_, sigrefval, refpix, increment, nChan)) {
|
---|
405 | //throw AipsError("Failed to get LO1 frequency from asis table");
|
---|
406 | os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
|
---|
407 | os << LogIO::WARN << "Could not fill frequency information of IMAGE sideband properly." << LogIO::POST;
|
---|
408 | os << LogIO::WARN << "Storing values of SIGNAL sideband in FREQUENCIES table" << LogIO::POST;
|
---|
409 | return;
|
---|
410 | }
|
---|
411 | }
|
---|
412 | // LO1 should now be ready.
|
---|
413 | if (lo1Freq_ < 0.)
|
---|
414 | throw(AipsError("Got negative LO1 Frequency"));
|
---|
415 |
|
---|
416 | // Print summary (LO1)
|
---|
417 | Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
|
---|
418 | os << "[LO1 settings]" << LogIO::POST;
|
---|
419 | os << "- Frequency: " << lo1Freq_ << " [Hz] ("
|
---|
420 | << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
|
---|
421 | os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
|
---|
422 | << " [day]" << LogIO::POST;
|
---|
423 | os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
|
---|
424 | << "] (" << md.getRefString() << ") " << LogIO::POST;
|
---|
425 |
|
---|
426 | //Print summary (signal)
|
---|
427 | os << "[Signal side band]" << LogIO::POST;
|
---|
428 | os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
|
---|
429 | << LogIO::POST;
|
---|
430 | os << "- Reference value: " << refval << " [Hz] ("
|
---|
431 | << MFrequency::showType(tabframe) << ") = "
|
---|
432 | << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_)
|
---|
433 | << ")" << LogIO::POST;
|
---|
434 | os << "- Reference pixel: " << refpix << LogIO::POST;
|
---|
435 | os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
|
---|
436 |
|
---|
437 | // Calculate image band incr and refval in loFrame_
|
---|
438 | Double imgincr = -increment;
|
---|
439 | Double imgrefval = 2 * lo1Freq_ - sigrefval;
|
---|
440 | Double imgrefval_tab = imgrefval;
|
---|
441 | // Convert imgrefval back to table base frame
|
---|
442 | if (tabframe != loFrame_)
|
---|
443 | imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue();
|
---|
444 | // Set new frequencies table
|
---|
445 | uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr);
|
---|
446 | // Update FREQ_ID in table.
|
---|
447 | freqIdVec = fIDnew;
|
---|
448 |
|
---|
449 | // Print summary (Image side band)
|
---|
450 | os << "[Image side band]" << LogIO::POST;
|
---|
451 | os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
|
---|
452 | << ")" << LogIO::POST;
|
---|
453 | os << "- Reference value: " << imgrefval << " [Hz] ("
|
---|
454 | << MFrequency::showType(tabframe) << ") = "
|
---|
455 | << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_)
|
---|
456 | << ")" << LogIO::POST;
|
---|
457 | os << "- Reference pixel: " << refpix << LogIO::POST;
|
---|
458 | os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
|
---|
459 | };
|
---|
460 |
|
---|
461 | Bool STSideBandSep::checkFile(const string name, string type)
|
---|
462 | {
|
---|
463 | File file(name);
|
---|
464 | if (!file.exists()){
|
---|
465 | return false;
|
---|
466 | } else if (type.empty()) {
|
---|
467 | return true;
|
---|
468 | } else {
|
---|
469 | // Check for file type
|
---|
470 | switch (tolower(type[0])) {
|
---|
471 | case 'f':
|
---|
472 | return file.isRegular(True);
|
---|
473 | case 'd':
|
---|
474 | return file.isDirectory(True);
|
---|
475 | case 's':
|
---|
476 | return file.isSymLink();
|
---|
477 | default:
|
---|
478 | throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
|
---|
479 | }
|
---|
480 | }
|
---|
481 | };
|
---|
482 |
|
---|
483 | bool STSideBandSep::getLo1FromAsdm(const string asdmname,
|
---|
484 | const double refval,
|
---|
485 | const double refpix,
|
---|
486 | const double increment,
|
---|
487 | const int nChan)
|
---|
488 | {
|
---|
489 | // Check for relevant tables.
|
---|
490 | string spwname = asdmname + "/SpectralWindow.xml";
|
---|
491 | string recname = asdmname + "/Receiver.xml";
|
---|
492 | if (!checkFile(spwname) || !checkFile(recname)) {
|
---|
493 | throw(AipsError("Could not find subtables in ASDM"));
|
---|
494 | }
|
---|
495 |
|
---|
496 | return false;
|
---|
497 |
|
---|
498 | };
|
---|
499 |
|
---|
500 |
|
---|
501 | bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab,
|
---|
502 | const double refval,
|
---|
503 | const double refpix,
|
---|
504 | const double increment,
|
---|
505 | const int nChan)
|
---|
506 | {
|
---|
507 | LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE));
|
---|
508 | // Check for relevant tables.
|
---|
509 | const TableRecord &rec = scantab->table().keywordSet() ;
|
---|
510 | String spwname, recname;
|
---|
511 | if (rec.isDefined("ASDM_SPECTRALWINDOW") && rec.isDefined("ASDM_RECEIVER")){
|
---|
512 | spwname = rec.asString("ASDM_SPECTRALWINDOW");
|
---|
513 | recname = rec.asString("ASDM_RECEIVER");
|
---|
514 | }
|
---|
515 | else {
|
---|
516 | // keywords are not there
|
---|
517 | os << LogIO::WARN
|
---|
518 | << "Could not find necessary table names in scantable header."
|
---|
519 | << LogIO::POST;
|
---|
520 | return false;
|
---|
521 | }
|
---|
522 | if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
|
---|
523 | throw(AipsError("Could not find relevant subtables in MS"));
|
---|
524 | }
|
---|
525 | // Get root MS name
|
---|
526 | string msname;
|
---|
527 | const String recsuff = "/ASDM_RECEIVER";
|
---|
528 | String::size_type pos;
|
---|
529 | pos = recname.size()-recsuff.size();
|
---|
530 | if (recname.substr(pos) == recsuff)
|
---|
531 | msname = recname.substr(0, pos);
|
---|
532 | else
|
---|
533 | throw(AipsError("Internal error in parsing table name from a scantable keyword."));
|
---|
534 |
|
---|
535 | if (!checkFile(msname))
|
---|
536 | throw(AipsError("Internal error in parsing MS name from a scantable keyword."));
|
---|
537 |
|
---|
538 | return getLo1FromAsisTab(msname, refval, refpix, increment, nChan);
|
---|
539 |
|
---|
540 | };
|
---|
541 |
|
---|
542 | bool STSideBandSep::getLo1FromAsisTab(const string msname,
|
---|
543 | const double refval,
|
---|
544 | const double refpix,
|
---|
545 | const double increment,
|
---|
546 | const int nChan)
|
---|
547 | {
|
---|
548 | LogIO os(LogOrigin("STSideBandSep","getLo1FromAsisTab()", WHERE));
|
---|
549 | os << "Searching an LO1 frequency in '" << msname << "'" << LogIO::POST;
|
---|
550 | // Check for relevant tables.
|
---|
551 | const string spwname = msname + "/ASDM_SPECTRALWINDOW";
|
---|
552 | const string recname = msname + "/ASDM_RECEIVER";
|
---|
553 | if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
|
---|
554 | throw(AipsError("Could not find relevant tables in MS"));
|
---|
555 | }
|
---|
556 |
|
---|
557 | Table spwtab_ = Table(spwname);
|
---|
558 | String asdmSpw;
|
---|
559 | ROTableRow spwrow(spwtab_);
|
---|
560 | const Double rtol = 0.01;
|
---|
561 | for (uInt idx = 0; idx < spwtab_.nrow(); idx++){
|
---|
562 | const TableRecord& rec = spwrow.get(idx);
|
---|
563 | // Compare nchan
|
---|
564 | if (rec.asInt("numChan") != (Int) nChan)
|
---|
565 | continue;
|
---|
566 | // Compare increment
|
---|
567 | Double asdminc;
|
---|
568 | Array<Double> incarr = rec.asArrayDouble("chanWidthArray");
|
---|
569 | if (incarr.size() > 0)
|
---|
570 | asdminc = incarr(IPosition(1, (uInt) refpix));
|
---|
571 | else
|
---|
572 | asdminc = rec.asDouble("chanWidth");
|
---|
573 | if (abs(asdminc - abs(increment)) > rtol * abs(increment))
|
---|
574 | continue;
|
---|
575 | // Compare refval
|
---|
576 | Double asdmrefv;
|
---|
577 | Array<Double> refvarr = rec.asArrayDouble("chanFreqArray");
|
---|
578 | if (refvarr.size() > 0){
|
---|
579 | const uInt iref = (uInt) refpix;
|
---|
580 | const Double ratio = refpix - (Double) iref;
|
---|
581 | asdmrefv = refvarr(IPosition(1, iref))*(1.-ratio)
|
---|
582 | + refvarr(IPosition(1,iref+1))*ratio;
|
---|
583 | }
|
---|
584 | else {
|
---|
585 | const Double ch0 = rec.asDouble("chanFreqStart");
|
---|
586 | const Double chstep = rec.asDouble("chanFreqStep");
|
---|
587 | asdmrefv = ch0 + chstep * refpix;
|
---|
588 | }
|
---|
589 | if (abs(asdmrefv - refval) < 0.5*abs(asdminc)){
|
---|
590 | asdmSpw = rec.asString("spectralWindowId");
|
---|
591 | break;
|
---|
592 | }
|
---|
593 | }
|
---|
594 |
|
---|
595 | if (asdmSpw.empty()){
|
---|
596 | os << LogIO::WARN << "Could not find relevant SPW ID in " << spwname << LogIO::POST;
|
---|
597 | return false;
|
---|
598 | }
|
---|
599 | else {
|
---|
600 | os << asdmSpw << " in " << spwname
|
---|
601 | << " matches the freqeuncies of signal side band." << LogIO::POST;
|
---|
602 | }
|
---|
603 |
|
---|
604 | Table rectab_ = Table(recname);
|
---|
605 | ROTableRow recrow(rectab_);
|
---|
606 | for (uInt idx = 0; idx < rectab_.nrow(); idx++){
|
---|
607 | const TableRecord& rec = recrow.get(idx);
|
---|
608 | if (rec.asString("spectralWindowId") == asdmSpw){
|
---|
609 | const Array<Double> loarr = rec.asArrayDouble("freqLO");
|
---|
610 | lo1Freq_ = loarr(IPosition(1,0));
|
---|
611 | os << "Found LO1 Frequency in " << recname << ": "
|
---|
612 | << lo1Freq_ << " [Hz]" << LogIO::POST;
|
---|
613 | return true;
|
---|
614 | }
|
---|
615 | }
|
---|
616 | os << LogIO::WARN << "Could not find " << asdmSpw << " in " << recname
|
---|
617 | << LogIO::POST;
|
---|
618 | return false;
|
---|
619 | };
|
---|
620 |
|
---|
621 | // String STSideBandSep::()
|
---|
622 | // {
|
---|
623 |
|
---|
624 | // };
|
---|
625 |
|
---|
626 | } //namespace asap
|
---|