source: trunk/src/STFiller.cpp@ 1490

Last change on this file since 1490 was 1467, checked in by Malte Marquarding, 16 years ago

fix for cast from double to float in scanrate field

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