source: trunk/src/STFiller.cpp@ 1446

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

allow ALMA specific build

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