source: trunk/src/STFiller.cpp@ 1114

Last change on this file since 1114 was 1110, checked in by mar637, 18 years ago

added detection of telesope for paddle scan rejection as this was rejecting all GBT scans

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