source: trunk/src/STFiller.cpp@ 1082

Last change on this file since 1082 was 1081, checked in by mar637, 18 years ago

using obsType to skip paddle scans

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