source: trunk/src/STFiller.cpp@ 1605

Last change on this file since 1605 was 1586, checked in by Malte Marquarding, 15 years ago

Ticket #165: have removed the hard-coding of parallactifying the data. NOTE THIS breaks the Table structure as I have moved the PARANGLE column from the main table into the FOCUS table. We need to have a new release. Also one needs to explicitly tell the scantable via rc or member function to enable parallactifying

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1//
2// C++ Implementation: STFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006-2007
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iostream.h>
13#include <casa/iomanip.h>
14
15#include <casa/Exceptions.h>
16#include <casa/OS/Path.h>
17#include <casa/OS/File.h>
18#include <casa/Quanta/Unit.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/ArrayLogical.h>
21#include <casa/Utilities/Regex.h>
22
23#include <casa/Containers/RecordField.h>
24
25#include <tables/Tables/TableRow.h>
26
27#include <atnf/PKSIO/PKSrecord.h>
28#include <atnf/PKSIO/PKSreader.h>
29#ifdef HAS_ALMA
30 #include <casa/System/ProgressMeter.h>
31#endif
32
33#include "STDefs.h"
34#include "STAttr.h"
35
36#include "STFiller.h"
37#include "STHeader.h"
38
39using namespace casa;
40
41namespace asap {
42
43STFiller::STFiller() :
44 reader_(0),
45 header_(0),
46 table_(0),
47 refRx_(".*(e|w|_R)$")
48{
49}
50
51STFiller::STFiller( CountedPtr< Scantable > stbl ) :
52 reader_(0),
53 header_(0),
54 table_(stbl),
55 refRx_(".*(e|w|_R)$")
56{
57}
58
59STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
60 reader_(0),
61 header_(0),
62 table_(0),
63 refRx_(".*(e|w|_R)$")
64{
65 open(filename, whichIF, whichBeam);
66}
67
68STFiller::~STFiller()
69{
70 close();
71}
72
73void STFiller::open( const std::string& filename, int whichIF, int whichBeam )
74{
75 if (table_.null()) {
76 table_ = new Scantable();
77 }
78 if (reader_) { delete reader_; reader_ = 0; }
79 Bool haveBase, haveSpectra;
80
81 String inName(filename);
82 Path path(inName);
83 inName = path.expandedName();
84
85 File file(inName);
86 if ( !file.exists() ) {
87 throw(AipsError("File does not exist"));
88 }
89 filename_ = inName;
90
91 // Create reader and fill in values for arguments
92 String format;
93 Vector<Bool> beams, ifs;
94 Vector<uInt> nchans,npols;
95 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
96 nchans, npols, haveXPol_,haveBase, haveSpectra
97 )) == 0 ) {
98 throw(AipsError("Creation of PKSreader failed"));
99 }
100 if (!haveSpectra) {
101 delete reader_;
102 reader_ = 0;
103 throw(AipsError("No spectral data in file."));
104 return;
105 }
106 nBeam_ = beams.nelements();
107 nIF_ = ifs.nelements();
108 // Get basic parameters.
109 if ( anyEQ(haveXPol_, True) ) {
110 pushLog("Cross polarization present");
111 for (uInt i=0; i< npols.nelements();++i) {
112 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
113 }
114 }
115 if (header_) delete header_;
116 header_ = new STHeader();
117 header_->nchan = max(nchans);
118 header_->npol = max(npols);
119 header_->nbeam = nBeam_;
120
121 Int status = reader_->getHeader(header_->observer, header_->project,
122 header_->antennaname, header_->antennaposition,
123 header_->obstype,
124 header_->fluxunit,
125 header_->equinox,
126 header_->freqref,
127 header_->utc, header_->reffreq,
128 header_->bandwidth);
129
130 if (status) {
131 delete reader_;
132 reader_ = 0;
133 delete header_;
134 header_ = 0;
135 throw(AipsError("Failed to get header."));
136 }
137 if ((header_->obstype).matches("*SW*")) {
138 // need robust way here - probably read ahead of next timestamp
139 pushLog("Header indicates frequency switched observation.\n"
140 "setting # of IFs = 1 ");
141 nIF_ = 1;
142 header_->obstype = String("fswitch");
143 }
144 // Determine Telescope and set brightness unit
145
146
147 Bool throwIt = False;
148 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
149
150 if (inst==ATMOPRA || inst==TIDBINBILLA) {
151 header_->fluxunit = "K";
152 } else {
153 // downcase for use with Quanta
154 if (header_->fluxunit == "JY") {
155 header_->fluxunit = "Jy";
156 }
157 }
158 STAttr stattr;
159 header_->poltype = stattr.feedPolType(inst);
160 header_->nif = nIF_;
161 header_->epoch = "UTC";
162 // *** header_->frequnit = "Hz"
163 // Apply selection criteria.
164 Vector<Int> ref;
165 ifOffset_ = 0;
166 if (whichIF>=0) {
167 if (whichIF>=0 && whichIF<nIF_) {
168 ifs = False;
169 ifs(whichIF) = True;
170 header_->nif = 1;
171 nIF_ = 1;
172 ifOffset_ = whichIF;
173 } else {
174 delete reader_;
175 reader_ = 0;
176 delete header_;
177 header_ = 0;
178 throw(AipsError("Illegal IF selection"));
179 }
180 }
181 beamOffset_ = 0;
182 if (whichBeam>=0) {
183 if (whichBeam>=0 && whichBeam<nBeam_) {
184 beams = False;
185 beams(whichBeam) = True;
186 header_->nbeam = 1;
187 nBeam_ = 1;
188 beamOffset_ = whichBeam;
189 } else {
190 delete reader_;
191 reader_ = 0;
192 delete header_;
193 header_ = 0;
194 throw(AipsError("Illegal Beam selection"));
195 }
196 }
197 Vector<Int> start(nIF_, 1);
198 Vector<Int> end(nIF_, 0);
199 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0]);
200 table_->setHeader(*header_);
201 //For MS, add the location of POINTING of the input MS so one get
202 //pointing data from there, if necessary.
203 //Also find nrow in MS
204 nInDataRow = 0;
205 if (format == "MS2") {
206 Path datapath(inName);
207 String ptTabPath = datapath.absoluteName();
208 Table inMS(ptTabPath);
209 nInDataRow = inMS.nrow();
210 ptTabPath.append("/POINTING");
211 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
212 if ((header_->antennaname).matches("GBT")) {
213 String GOTabPath = datapath.absoluteName();
214 GOTabPath.append("/GBT_GO");
215 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
216 }
217 }
218 //table_->focus().setParallactify(true);
219}
220
221void STFiller::close( )
222{
223 delete reader_;reader_=0;
224 delete header_;header_=0;
225 table_ = 0;
226}
227
228int asap::STFiller::read( )
229{
230 int status = 0;
231
232 Double min = 0.0;
233 Double max = nInDataRow;
234#ifdef HAS_ALMA
235 ProgressMeter fillpm(min, max, "Data importing progress");
236#endif
237 PKSrecord pksrec;
238 int n = 0;
239 while ( status == 0 ) {
240 status = reader_->read(pksrec);
241 if ( status != 0 ) break;
242 n += 1;
243
244 Regex filterrx(".*[SL|PA]$");
245 Regex obsrx("^AT.+");
246 if ( header_->antennaname.matches(obsrx) &&
247 pksrec.obsType.matches(filterrx)) {
248 //cerr << "ignoring paddle scan" << endl;
249 continue;
250 }
251 TableRow row(table_->table());
252 TableRecord& rec = row.record();
253 // fields that don't get used and are just passed through asap
254 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
255 // MRC changed type from double to float
256 Vector<Double> sratedbl(pksrec.scanRate.nelements());
257 convertArray(sratedbl, pksrec.scanRate);
258 *srateCol = sratedbl;
259 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
260 *spmCol = pksrec.srcPM;
261 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
262 *sdirCol = pksrec.srcDir;
263 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
264 *svelCol = pksrec.srcVel;
265 // the real stuff
266 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
267 *fitCol = -1;
268 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
269 *scanoCol = pksrec.scanNo-1;
270 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
271 *cyclenoCol = pksrec.cycleNo-1;
272 RecordFieldPtr<Double> mjdCol(rec, "TIME");
273 *mjdCol = pksrec.mjd;
274 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
275 *intCol = pksrec.interval;
276 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
277 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
278 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
279 *fieldnCol = pksrec.fieldName;
280 // try to auto-identify if it is on or off.
281 Regex rx(refRx_);
282 Regex rx2("_S$");
283 Int match = pksrec.srcName.matches(rx);
284 if (match) {
285 *srcnCol = pksrec.srcName;
286 } else {
287 *srcnCol = pksrec.srcName.before(rx2);
288 }
289 //*srcnCol = pksrec.srcName;//.before(rx2);
290 *srctCol = match;
291 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
292 *beamCol = pksrec.beamNo-beamOffset_-1;
293 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
294 Int rb = -1;
295 if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
296 *rbCol = rb;
297 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
298 *ifCol = pksrec.IFno-ifOffset_- 1;
299 uInt id;
300 /// @todo this has to change when nchan isn't global anymore
301 id = table_->frequencies().addEntry(Double(header_->nchan/2),
302 pksrec.refFreq, pksrec.freqInc);
303 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
304 *mfreqidCol = id;
305
306 id = table_->molecules().addEntry(pksrec.restFreq);
307 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
308 *molidCol = id;
309
310 id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
311 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
312 *mcalidCol = id;
313 id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
314 pksrec.humidity, pksrec.windSpeed,
315 pksrec.windAz);
316 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
317 *mweatheridCol = id;
318
319 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
320 id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
321 pksrec.focusTan, pksrec.focusRot);
322 *mfocusidCol = id;
323 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
324 *dirCol = pksrec.direction;
325 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
326 *azCol = pksrec.azimuth;
327 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
328 *elCol = pksrec.elevation;
329
330 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
331 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
332 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
333
334 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
335 // Turn the (nchan,npol) matrix and possible complex xPol vector
336 // into 2-4 rows in the scantable
337 Vector<Float> tsysvec(1);
338 // Why is pksrec.spectra.ncolumn() == 3 for haveXPol_ == True
339 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
340 for ( uInt i=0; i< npol; ++i ) {
341 tsysvec = pksrec.tsys(i);
342 *tsysCol = tsysvec;
343 *polnoCol = i;
344
345 *specCol = pksrec.spectra.column(i);
346 *flagCol = pksrec.flagged.column(i);
347 table_->table().addRow();
348 row.put(table_->table().nrow()-1, rec);
349 }
350 if ( haveXPol_[0] ) {
351 // no tsys given for xpol, so emulate it
352 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
353 *tsysCol = tsysvec;
354 // add real part of cross pol
355 *polnoCol = 2;
356 Vector<Float> r(real(pksrec.xPol));
357 *specCol = r;
358 // make up flags from linears
359 /// @fixme this has to be a bitwise or of both pols
360 *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
361 table_->table().addRow();
362 row.put(table_->table().nrow()-1, rec);
363 // ad imaginary part of cross pol
364 *polnoCol = 3;
365 Vector<Float> im(imag(pksrec.xPol));
366 *specCol = im;
367 table_->table().addRow();
368 row.put(table_->table().nrow()-1, rec);
369 }
370#ifdef HAS_ALMA
371 fillpm._update(n);
372#endif
373 }
374 if (status > 0) {
375 close();
376 throw(AipsError("Reading error occured, data possibly corrupted."));
377 }
378#ifdef HAS_ALMA
379 fillpm.done();
380#endif
381 return status;
382}
383
384}//namespace asap
Note: See TracBrowser for help on using the repository browser.