source: trunk/src/STSideBandSep.cpp@ 2715

Last change on this file since 2715 was 2712, checked in by Kana Sugimoto, 12 years ago

New Development: No

JIRA Issue: Yes (CAS-4141)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: added temporal methods to set scantable to STSideBandSep for a workaround.

Test Programs:

Put in Release Notes: No

Module(s): STSideBandSep and sbseparator

Description:

Enabled getting LO1 frequency based on information in scantable header
(if necessary information exists).


File size: 14.7 KB
Line 
1// C++ Interface: STSideBandSep
2//
3// Description:
4// A class to invoke sideband separation of Scantable
5//
6// Author: Kana Sugimoto <kana.sugi@nao.ac.jp>, (C) 2012
7//
8// Copyright: See COPYING file that comes with this distribution
9//
10//
11
12// STL
13#include <ctype.h>
14
15// cascore
16#include <casa/OS/File.h>
17#include <casa/Logging/LogIO.h>
18#include <casa/Quanta/QuantumHolder.h>
19
20#include <measures/Measures/MFrequency.h>
21#include <measures/Measures/MCFrequency.h>
22
23#include <tables/Tables/TableRow.h>
24#include <tables/Tables/TableRecord.h>
25#include <tables/Tables/TableVector.h>
26
27// asap
28#include "STSideBandSep.h"
29
30using namespace std ;
31using namespace casa ;
32using namespace asap ;
33
34//#ifndef KS_DEBUG
35//#define KS_DEBUG
36//#endif
37
38namespace asap {
39
40// constructors
41STSideBandSep::STSideBandSep()
42 : sigIfno_(0), ftol_(-1),
43 loTime_(-1), lo1Freq_(-1), loDir_("")
44{
45#ifdef KS_DEBUG
46 cout << "Default constructor STSideBandSep()" << endl;
47#endif
48 // Default LO frame is TOPO
49 loFrame_ = MFrequency::TOPO;
50};
51
52STSideBandSep::~STSideBandSep()
53{
54#ifdef KS_DEBUG
55 cout << "Destructor ~STSideBandSep()" << endl;
56#endif
57};
58
59
60void STSideBandSep::setFrequency(const unsigned int ifno, const double freqtol, string frame)
61{
62#ifdef KS_DEBUG
63 cout << "STSideBandSep::setFrequency" << endl;
64 cout << "IFNO = " << ifno << endl;
65 cout << "freq tol = " << freqtol << " [Hz]" << endl;
66 cout << "frame = " << frame << endl;
67#endif
68};
69
70// Temporal function to set scantable of image sideband
71void STSideBandSep::setImageTable(const ScantableWrapper &s)
72{
73#ifdef KS_DEBUG
74 cout << "STSideBandSep::setImageTable" << endl;
75 cout << "got image table nrow = " << s.nrow() << endl;
76#endif
77 imgTab_p = s.getCP();
78 AlwaysAssert(!imgTab_p.null(),AipsError);
79
80};
81
82//
83void STSideBandSep::setLO1(const double lo1, string frame, double reftime, string refdir)
84{
85 lo1Freq_ = lo1;
86 MFrequency::getType(loFrame_, frame);
87 loTime_ = reftime;
88 loDir_ = refdir;
89
90#ifdef KS_DEBUG
91 cout << "STSideBandSep::setLO1" << endl;
92 if (lo1Freq_ > 0.)
93 cout << "lo1 = " << lo1Freq_ << " [Hz] (" << frame << ")" << endl;
94 if (loTime_ > 0.)
95 cout << "ref time = " << loTime_ << " [day]" << endl;
96 if (!loDir_.empty())
97 cout << "ref direction = " << loDir_ << " [day]" << endl;
98#endif
99};
100
101void STSideBandSep::setLO1Root(string name)
102{
103 LogIO os(LogOrigin("STSideBandSep","setLO1Root()", WHERE));
104 os << "Searching for '" << name << "'..." << LogIO::POST;
105 // Check for existance of the file
106 if (!checkFile(name)) {
107 throw(AipsError("File does not exist"));
108 }
109 if (name[(name.size()-1)] == '/')
110 name = name.substr(0,(name.size()-2));
111
112 if (checkFile(name+"/Receiver.xml", "file") &&
113 checkFile(name+"/SpectralWindow.xml", "file")){
114 os << "Found '" << name << "/Receiver.xml' ... got an ASDM name." << LogIO::POST;
115 asdmName_ = name;
116 } else if (checkFile(name+"/ASDM_RECEIVER") &&
117 checkFile(name+"/ASDM_SPECTRALWINDOW")){
118 os << "Found '" << name << "/ASDM_RECEIVER' ... got a Table name." << LogIO::POST;
119 asisName_ = name;
120 } else {
121 throw(AipsError("Invalid file name. Set an MS or ASDM name."));
122 }
123
124#ifdef KS_DEBUG
125 cout << "STSideBandSep::setLO1Root" << endl;
126 if (!asdmName_.empty())
127 cout << "asdm name = " << asdmName_ << endl;
128 if (!asisName_.empty())
129 cout << "MS name = " << asisName_ << endl;
130#endif
131};
132
133///// TEMPORAL FUNCTION!!! /////
134void STSideBandSep::setScanTb0(const ScantableWrapper &s){
135 st0_ = s.getCP();
136};
137////////////////////////////////
138
139void STSideBandSep::solveImageFreqency()
140{
141#ifdef KS_DEBUG
142 cout << "STSideBandSep::solveImageFrequency" << endl;
143#endif
144 LogIO os(LogOrigin("STSideBandSep","solveImageFreqency()", WHERE));
145 os << "Start calculating frequencies of image side band" << LogIO::POST;
146
147 if (imgTab_p.null())
148 throw AipsError("STSideBandSep::solveImageFreqency - an image side band scantable should be set first");
149
150 // Convert frequency REFVAL to the value in frame of LO.
151 // The code assumes that imgTab_p has only an IF and only a FREQ_ID
152 // is associated to an IFNO
153 // TODO: More complete Procedure would be
154 // 1. Get freq IDs associated to sigIfno_
155 // 2. Get freq information of the freq IDs
156 // 3. For each freqIDs, get freq infromation in TOPO and an LO1
157 // frequency and calculate image band frequencies.
158 STFrequencies freqTab_ = imgTab_p->frequencies();
159 // get the base frame of table
160 const MFrequency::Types tabframe = freqTab_.getFrame(true);
161 TableVector<uInt> freqIdVec( imgTab_p->table(), "FREQ_ID" );
162 // assuming single freqID per IFNO
163 uInt freqid = freqIdVec(0);
164 int nChan = imgTab_p->nchan(imgTab_p->getIF(0));
165 double refpix, refval, increment ;
166 freqTab_.getEntry(refpix, refval, increment, freqid);
167 //MFrequency sigrefval = MFrequency(MVFrequency(refval),tabframe);
168 // get freq infromation of sigIfno_ in loFrame_
169 const MPosition mp = imgTab_p->getAntennaPosition();
170 MEpoch me;
171 MDirection md;
172 if (loTime_ < 0.)
173 me = imgTab_p->getEpoch(-1);
174 else
175 me = MEpoch(MVEpoch(loTime_));
176 if (loDir_.empty()) {
177 ArrayColumn<Double> srcdirCol_;
178 srcdirCol_.attach(imgTab_p->table(), "SRCDIRECTION");
179 // Assuming J2000 and SRCDIRECTION in unit of rad
180 Quantum<Double> srcra = Quantum<Double>(srcdirCol_(0)(IPosition(1,0)), "rad");
181 Quantum<Double> srcdec = Quantum<Double>(srcdirCol_(0)(IPosition(1,1)), "rad");
182 md = MDirection(srcra, srcdec, MDirection::J2000);
183 //imgTab_p->getDirection(0);
184 } else {
185 // parse direction string
186 string::size_type pos0 = loDir_.find(" ");
187
188 if (pos0 == string::npos) {
189 throw AipsError("bad string format in LO1 direction");
190 }
191 string::size_type pos1 = loDir_.find(" ", pos0+1);
192 String sepoch, sra, sdec;
193 if (pos1 != string::npos) {
194 sepoch = loDir_.substr(0, pos0);
195 sra = loDir_.substr(pos0+1, pos1-pos0);
196 sdec = loDir_.substr(pos1+1);
197 }
198 MDirection::Types epoch;
199 MDirection::getType(epoch, sepoch);
200 QuantumHolder qh ;
201 String err ;
202 qh.fromString( err, sra);
203 Quantum<Double> ra = qh.asQuantumDouble() ;
204 qh.fromString( err, sdec ) ;
205 Quantum<Double> dec = qh.asQuantumDouble() ;
206 md = MDirection(ra.getValue("rad"), dec.getValue("rad"),epoch);
207 }
208 MeasFrame mframe( me, mp, md );
209 MFrequency::Convert tobframe(loFrame_, MFrequency::Ref(tabframe, mframe));
210 MFrequency::Convert toloframe(tabframe, MFrequency::Ref(loFrame_, mframe));
211 // Convert refval to loFrame_
212 double sigrefval;
213 if (tabframe == loFrame_)
214 sigrefval = refval;
215 else
216 sigrefval = toloframe(Quantum<Double>(refval, "Hz")).get("Hz").getValue();
217
218
219 // Check for the availability of LO1
220 if (lo1Freq_ > 0.) {
221 os << "Using user defined LO1 frequency." << LogIO::POST;
222 } else if (!asisName_.empty()) {
223 // MS name is set.
224 os << "Using user defined MS (asis): " << asisName_ << LogIO::POST;
225 if (!getLo1FromAsisTab(asisName_, sigrefval, refpix, increment, nChan)) {
226 throw AipsError("Failed to get LO1 frequency from MS");
227 }
228 } else if (!asdmName_.empty()) {
229 // ASDM name is set.
230 os << "Using user defined ASDM: " << asdmName_ << LogIO::POST;
231 if (!getLo1FromAsdm(asdmName_, sigrefval, refpix, increment, nChan)) {
232 throw AipsError("Failed to get LO1 frequency from ASDM");
233 }
234 } else {
235 // Try getting ASDM name from scantable header
236 os << "Try getting information from scantable header" << LogIO::POST;
237 if (!getLo1FromScanTab(st0_, sigrefval, refpix, increment, nChan)) {
238 //throw AipsError("Failed to get LO1 frequency from asis table");
239 os << LogIO::WARN << "Failed to get LO1 frequency using information in scantable." << LogIO::POST;
240 os << LogIO::WARN << "Could not fill frequency information of IMAGE sideband properly." << LogIO::POST;
241 os << LogIO::WARN << "Storing values of SIGNAL sideband in FREQUENCIES table" << LogIO::POST;
242 return;
243 }
244 }
245 // LO1 should now be ready.
246 if (lo1Freq_ < 0.)
247 throw(AipsError("Got negative LO1 Frequency"));
248
249 // Print summary (LO1)
250 Vector<Double> dirvec = md.getAngle(Unit(String("rad"))).getValue();
251 os << "[LO1 settings]" << LogIO::POST;
252 os << "- Frequency: " << lo1Freq_ << " [Hz] ("
253 << MFrequency::showType(loFrame_) << ")" << LogIO::POST;
254 os << "- Reference time: " << me.get(Unit(String("d"))).getValue()
255 << " [day]" << LogIO::POST;
256 os << "- Reference direction: [" << dirvec[0] << ", " << dirvec[1]
257 << "] (" << md.getRefString() << ") " << LogIO::POST;
258
259 //Print summary (signal)
260 os << "[Signal side band]" << LogIO::POST;
261 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqid << ")"
262 << LogIO::POST;
263 os << "- Reference value: " << refval << " [Hz] ("
264 << MFrequency::showType(tabframe) << ") = "
265 << sigrefval << " [Hz] (" << MFrequency::showType(loFrame_)
266 << ")" << LogIO::POST;
267 os << "- Reference pixel: " << refpix << LogIO::POST;
268 os << "- Increment: " << increment << " [Hz]" << LogIO::POST;
269
270 // Calculate image band incr and refval in loFrame_
271 Double imgincr = -increment;
272 Double imgrefval = 2 * lo1Freq_ - sigrefval;
273 Double imgrefval_tab = imgrefval;
274 // Convert imgrefval back to table base frame
275 if (tabframe != loFrame_)
276 imgrefval = tobframe(Quantum<Double>(imgrefval, "Hz")).get("Hz").getValue();
277 // Set new frequencies table
278 uInt fIDnew = freqTab_.addEntry(refpix, imgrefval, imgincr);
279 // Update FREQ_ID in table.
280 freqIdVec = fIDnew;
281
282 // Print summary (Image side band)
283 os << "[Image side band]" << LogIO::POST;
284 os << "- IFNO: " << imgTab_p->getIF(0) << " (FREQ_ID = " << freqIdVec(0)
285 << ")" << LogIO::POST;
286 os << "- Reference value: " << imgrefval << " [Hz] ("
287 << MFrequency::showType(tabframe) << ") = "
288 << imgrefval_tab << " [Hz] (" << MFrequency::showType(loFrame_)
289 << ")" << LogIO::POST;
290 os << "- Reference pixel: " << refpix << LogIO::POST;
291 os << "- Increment: " << imgincr << " [Hz]" << LogIO::POST;
292};
293
294Bool STSideBandSep::checkFile(const string name, string type)
295{
296 File file(name);
297 if (!file.exists()){
298 return false;
299 } else if (type.empty()) {
300 return true;
301 } else {
302 // Check for file type
303 switch (tolower(type[0])) {
304 case 'f':
305 return file.isRegular(True);
306 case 'd':
307 return file.isDirectory(True);
308 case 's':
309 return file.isSymLink();
310 default:
311 throw AipsError("Invalid file type. Available types are 'file', 'directory', and 'symlink'.");
312 }
313 }
314};
315
316bool STSideBandSep::getLo1FromAsdm(const string asdmname,
317 const double refval,
318 const double refpix,
319 const double increment,
320 const int nChan)
321{
322 // Check for relevant tables.
323 string spwname = asdmname + "/SpectralWindow.xml";
324 string recname = asdmname + "/Receiver.xml";
325 if (!checkFile(spwname) || !checkFile(recname)) {
326 throw(AipsError("Could not find subtables in ASDM"));
327 }
328
329 return false;
330
331};
332
333
334bool STSideBandSep::getLo1FromScanTab(CountedPtr< Scantable > &scantab,
335 const double refval,
336 const double refpix,
337 const double increment,
338 const int nChan)
339{
340 LogIO os(LogOrigin("STSideBandSep","getLo1FromScanTab()", WHERE));
341 // Check for relevant tables.
342 const TableRecord &rec = scantab->table().keywordSet() ;
343 String spwname, recname;
344 if (rec.isDefined("ASDM_SPECTRALWINDOW") && rec.isDefined("ASDM_RECEIVER")){
345 spwname = rec.asString("ASDM_SPECTRALWINDOW");
346 recname = rec.asString("ASDM_RECEIVER");
347 }
348 else {
349 // keywords are not there
350 os << LogIO::WARN
351 << "Could not find necessary table names in scantable header."
352 << LogIO::POST;
353 return false;
354 }
355 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
356 throw(AipsError("Could not find relevant subtables in MS"));
357 }
358 // Get root MS name
359 string msname;
360 const String recsuff = "/ASDM_RECEIVER";
361 String::size_type pos;
362 pos = recname.size()-recsuff.size();
363 if (recname.substr(pos) == recsuff)
364 msname = recname.substr(0, pos);
365 else
366 throw(AipsError("Internal error in parsing table name from a scantable keyword."));
367
368 if (!checkFile(msname))
369 throw(AipsError("Internal error in parsing MS name from a scantable keyword."));
370
371 return getLo1FromAsisTab(msname, refval, refpix, increment, nChan);
372
373};
374
375bool STSideBandSep::getLo1FromAsisTab(const string msname,
376 const double refval,
377 const double refpix,
378 const double increment,
379 const int nChan)
380{
381 LogIO os(LogOrigin("STSideBandSep","getLo1FromAsisTab()", WHERE));
382 os << "Searching an LO1 frequency in '" << msname << "'" << LogIO::POST;
383 // Check for relevant tables.
384 const string spwname = msname + "/ASDM_SPECTRALWINDOW";
385 const string recname = msname + "/ASDM_RECEIVER";
386 if (!checkFile(spwname,"directory") || !checkFile(recname,"directory")) {
387 throw(AipsError("Could not find relevant tables in MS"));
388 }
389
390 Table spwtab_ = Table(spwname);
391 String asdmSpw;
392 ROTableRow spwrow(spwtab_);
393 const Double rtol = 0.01;
394 for (uInt idx = 0; idx < spwtab_.nrow(); idx++){
395 const TableRecord& rec = spwrow.get(idx);
396 // Compare nchan
397 if (rec.asInt("numChan") != (Int) nChan)
398 continue;
399 // Compare increment
400 Double asdminc;
401 Array<Double> incarr = rec.asArrayDouble("chanWidthArray");
402 if (incarr.size() > 0)
403 asdminc = incarr(IPosition(1, (uInt) refpix));
404 else
405 asdminc = rec.asDouble("chanWidth");
406 if (abs(asdminc - abs(increment)) > rtol * abs(increment))
407 continue;
408 // Compare refval
409 Double asdmrefv;
410 Array<Double> refvarr = rec.asArrayDouble("chanFreqArray");
411 if (refvarr.size() > 0){
412 const uInt iref = (uInt) refpix;
413 const Double ratio = refpix - (Double) iref;
414 asdmrefv = refvarr(IPosition(1, iref))*(1.-ratio)
415 + refvarr(IPosition(1,iref+1))*ratio;
416 }
417 else {
418 const Double ch0 = rec.asDouble("chanFreqStart");
419 const Double chstep = rec.asDouble("chanFreqStep");
420 asdmrefv = ch0 + chstep * refpix;
421 }
422 if (abs(asdmrefv - refval) < 0.5*abs(asdminc)){
423 asdmSpw = rec.asString("spectralWindowId");
424 break;
425 }
426 }
427
428 if (asdmSpw.empty()){
429 os << LogIO::WARN << "Could not find relevant SPW ID in " << spwname << LogIO::POST;
430 return false;
431 }
432 else {
433 os << asdmSpw << " in " << spwname
434 << " matches the freqeuncies of signal side band." << LogIO::POST;
435 }
436
437 Table rectab_ = Table(recname);
438 ROTableRow recrow(rectab_);
439 for (uInt idx = 0; idx < rectab_.nrow(); idx++){
440 const TableRecord& rec = recrow.get(idx);
441 if (rec.asString("spectralWindowId") == asdmSpw){
442 const Array<Double> loarr = rec.asArrayDouble("freqLO");
443 lo1Freq_ = loarr(IPosition(1,0));
444 os << "Found LO1 Frequency in " << recname << ": "
445 << lo1Freq_ << " [Hz]" << LogIO::POST;
446 return true;
447 }
448 }
449 os << LogIO::WARN << "Could not find " << asdmSpw << " in " << recname
450 << LogIO::POST;
451 return false;
452};
453
454// String STSideBandSep::()
455// {
456
457// };
458
459} //namespace asap
Note: See TracBrowser for help on using the repository browser.