source: trunk/src/STFiller.cpp@ 1500

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