source: trunk/src/PKSFiller.cpp@ 2580

Last change on this file since 2580 was 2268, checked in by Malte Marquarding, 13 years ago

Ignore flagrow for rpfits data for now

File size: 10.4 KB
RevLine 
[1795]1//
2// C++ Implementation: PKSFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2010
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>
[1916]22#include <casa/Utilities/DataType.h>
[1816]23#include <casa/Logging/LogIO.h>
[1795]24
[1904]25#include <casa/Containers/Record.h>
[1795]26#include <measures/Measures/MDirection.h>
27#include <measures/Measures/MeasConvert.h>
28
29#include <atnf/PKSIO/PKSrecord.h>
30#include <atnf/PKSIO/PKSreader.h>
31
32#ifdef HAS_ALMA
33 #include <casa/System/ProgressMeter.h>
34#endif
35#include <casa/Logging/LogIO.h>
36
37#include <time.h>
[1916]38#include <sstream>
[1795]39
40#include "STDefs.h"
41#include "STAttr.h"
42
43#include "PKSFiller.h"
44#include "STHeader.h"
45
46using namespace casa;
47
48namespace asap {
49
50
51PKSFiller::PKSFiller( CountedPtr< Scantable > stbl ) :
52 FillerBase(stbl),
[1815]53 reader_(0)
[1795]54{
55 setReferenceRegex(".*(e|w|_R)$");
56}
57
58PKSFiller::~PKSFiller()
59{
60 close();
61}
62
[1904]63bool PKSFiller::open( const std::string& filename, const Record& rec)
[1795]64{
65 Bool haveBase, haveSpectra;
66
67 String inName(filename);
68 Path path(inName);
69 inName = path.expandedName();
70
71 File file(inName);
72 filename_ = inName;
73
74 // Create reader and fill in values for arguments
75 String format;
76 Vector<Bool> beams, ifs;
77 Vector<uInt> nchans,npols;
78
[1880]79 String antenna("");
[1916]80 Bool getPt = False;
[1815]81
[1916]82 // parsing MS options
83 if ( rec.isDefined( "ms" ) ) {
84 Record msrec = rec.asRecord( "ms" ) ;
85 //msrec.print( cout ) ;
86 if ( msrec.isDefined( "getpt" ) ) {
87 getPt = msrec.asBool( "getpt" ) ;
88 }
89 if ( msrec.isDefined( "antenna" ) ) {
90 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
91 Int antInt = msrec.asInt( "antenna" ) ;
92 ostringstream oss ;
93 oss << antInt ;
94 antenna = String( oss ) ;
95 }
96 else {
97 antenna = msrec.asString( "antenna" ) ;
98 }
99 }
100 }
101
[1815]102 reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
103 nchans, npols, haveXPol_, haveBase, haveSpectra);
104 if (reader_.null() == True) {
[1795]105 return False;
106 }
107
108 if (!haveSpectra) {
109 reader_ = 0;
110 throw(AipsError("No spectral data in file."));
111 }
[1816]112
113 LogIO os( LogOrigin( "PKSFiller", "open()", WHERE ) ) ;
[1795]114 nBeam_ = beams.nelements();
115 nIF_ = ifs.nelements();
116 // Get basic parameters.
117 if ( anyEQ(haveXPol_, True) ) {
118 os << "Cross polarization present" << LogIO::POST;
119 for (uInt i=0; i< npols.nelements();++i) {
120 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
121 }
122 }
[1815]123 STHeader header;
124 header.nchan = max(nchans);
125 header.npol = max(npols);
126 header.nbeam = nBeam_;
[1795]127
[1815]128 Int status = reader_->getHeader(header.observer, header.project,
129 header.antennaname, header.antennaposition,
130 header.obstype,
131 header.fluxunit,
132 header.equinox,
133 header.freqref,
134 header.utc, header.reffreq,
135 header.bandwidth);
[1795]136
137 if (status) {
138 reader_ = 0;
139 throw(AipsError("Failed to get header."));
140 }
[1816]141 os << "Found " << header.antennaname << " data." << LogIO::POST;
[1815]142 if ((header.obstype).matches("*SW*")) {
[1795]143 // need robust way here - probably read ahead of next timestamp
144 os << "Header indicates frequency switched observation.\n"
145 << "setting # of IFs = 1 " << LogIO::POST;
146
147 nIF_ = 1;
[1815]148 header.obstype = String("fswitch");
[1795]149 }
150 // Determine Telescope and set brightness unit
151 Bool throwIt = False;
[1815]152 Instrument inst = STAttr::convertInstrument(header.antennaname, throwIt);
[1795]153
154 if (inst==ATMOPRA || inst==TIDBINBILLA) {
[1815]155 header.fluxunit = "K";
[1795]156 } else {
157 // downcase for use with Quanta
[1815]158 if (header.fluxunit == "JY") {
159 header.fluxunit = "Jy";
[1795]160 }
161 }
162 STAttr stattr;
[1815]163 header.poltype = stattr.feedPolType(inst);
164 header.nif = nIF_;
165 header.epoch = "UTC";
166 // *** header.frequnit = "Hz"
[1795]167 // Apply selection criteria.
168 Vector<Int> ref;
169
170 Vector<Int> start(nIF_, 1);
171 Vector<Int> end(nIF_, 0);
172 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
[1815]173 setHeader(header);
[1795]174 //For MS, add the location of POINTING of the input MS so one get
175 //pointing data from there, if necessary.
176 //Also find nrow in MS
177 nInDataRow = 0;
178 if (format == "MS2") {
179 Path datapath(inName);
180 String ptTabPath = datapath.absoluteName();
181 Table inMS(ptTabPath);
182 nInDataRow = inMS.nrow();
183 ptTabPath.append("/POINTING");
184 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
[1815]185 if (header.antennaname.matches("GBT")) {
[1795]186 String GOTabPath = datapath.absoluteName();
187 GOTabPath.append("/GBT_GO");
188 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
189 }
190 }
[1815]191 String freqFrame = header.freqref;
[1795]192 //translate frequency reference frame back to
193 //MS style (as PKSMS2reader converts the original frame
194 //in FITS standard style)
195 if (freqFrame == "TOPOCENT") {
196 freqFrame = "TOPO";
197 } else if (freqFrame == "GEOCENER") {
198 freqFrame = "GEO";
199 } else if (freqFrame == "BARYCENT") {
200 freqFrame = "BARY";
201 } else if (freqFrame == "GALACTOC") {
202 freqFrame = "GALACTO";
203 } else if (freqFrame == "LOCALGRP") {
204 freqFrame = "LGROUP";
205 } else if (freqFrame == "CMBDIPOL") {
206 freqFrame = "CMB";
207 } else if (freqFrame == "SOURCE") {
208 freqFrame = "REST";
209 }
210 // set both "FRAME" and "BASEFRAME"
211 table_->frequencies().setFrame(freqFrame, false);
[1815]212 table_->frequencies().setFrame(freqFrame, true);
[1795]213
214 return true;
215}
216
217void PKSFiller::close( )
218{
[1815]219 if (reader_.null() != False) {
220 reader_ = 0;
[1795]221 }
222 table_ = 0;
223}
224
225void PKSFiller::fill( )
226{
227 int status = 0;
[1816]228 LogIO os( LogOrigin( "PKSFiller", "fill()", WHERE ) ) ;
[1795]229
230/**
231 Int beamNo, IFno, refBeam, scanNo, cycleNo;
232 Float azimuth, elevation, focusAxi, focusRot, focusTan,
233 humidity, parAngle, pressure, temperature, windAz, windSpeed;
234 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
235 String fieldName, srcName, tcalTime, obsType;
236 Vector<Float> calFctr, sigma, tcal, tsys;
237 Matrix<Float> baseLin, baseSub;
238 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
239 Matrix<Float> spectra;
240 Matrix<uChar> flagtra;
241 Complex xCalFctr;
242 Vector<Complex> xPol;
243**/
244
245 Double min = 0.0;
246 Double max = nInDataRow;
247#ifdef HAS_ALMA
248 ProgressMeter fillpm(min, max, "Data importing progress");
249#endif
250 PKSrecord pksrec;
251 pksrec.srcType=-1;
252 int n = 0;
253 bool isGBTFITS = false ;
[1815]254 if ((table_->getAntennaName().find( "GBT" ) != String::npos)
[1795]255 && File(filename_).isRegular()) {
256 FILE *fp = fopen( filename_.c_str(), "r" ) ;
257 fseek( fp, 640, SEEK_SET ) ;
258 char buf[81] ;
[2163]259 size_t tmp = fread( buf, 80, 1, fp ) ;
260 (void) tmp;
[1795]261 buf[80] = '\0' ;
262 if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
263 isGBTFITS = true ;
264 }
265 fclose( fp ) ;
266 }
267
268 while ( status == 0 ) {
269
270 status = reader_->read(pksrec);
271 if ( status != 0 ) break;
272 n += 1;
273 Regex filterrx(".*[SL|PA]$");
274 Regex obsrx("^AT.+");
[1815]275 if ( table_->getAntennaName().matches(obsrx) &&
[1795]276 pksrec.obsType.matches(filterrx)) {
277 //cerr << "ignoring paddle scan" << endl;
278 continue;
279 }
280
281 Vector<Double> sratedbl(pksrec.scanRate.nelements());
282 convertArray(sratedbl, pksrec.scanRate);
[1816]283 setScanRate(sratedbl);
[1795]284
285 Regex rx(getReferenceRegex());
286 Regex rx2("_S$");
287 Int match = pksrec.srcName.matches(rx);
288 std::string srcname;
[1872]289 Int srctype = Int(SrcType::NOTYPE);
[1795]290 if (match) {
291 srcname = pksrec.srcName;
[1872]292 srctype = Int(SrcType::PSOFF);
[1795]293 } else {
294 srcname = pksrec.srcName.before(rx2);
[1872]295 srctype = Int(SrcType::PSON);
[1795]296 }
[1872]297 if ( pksrec.srcType != Int(SrcType::NOTYPE)) {
[1795]298 srctype = pksrec.srcType ;
299 }
300 setTime(pksrec.mjd, pksrec.interval);
301 setSource(srcname, srctype, pksrec.fieldName,
302 pksrec.srcDir, pksrec.srcPM, pksrec.srcVel);
303
304 if (nBeam_ > 1 ) {
305 setReferenceBeam(pksrec.refBeam-1);
306 }
[1805]307
[1795]308 setMolecule(pksrec.restFreq);
309 setTcal(pksrec.tcalTime, pksrec.tcal);
310 setWeather(pksrec.temperature, pksrec.pressure,
311 pksrec.humidity, pksrec.windSpeed,
312 pksrec.windAz);
313 setFocus(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
314 pksrec.focusRot);
315 setDirection(pksrec.direction, pksrec.azimuth, pksrec.elevation);
316
[1805]317 setFrequency(Double(pksrec.spectra.nrow()/2),
318 pksrec.refFreq, pksrec.freqInc);
319
[1816]320 // Note: this is only one value for all polarisations!
[2268]321
322 //setFlagrow(pksrec.flagrow);
[1795]323 // Turn the (nchan,npol) matrix and possible complex xPol vector
324 // into 2-4 rows in the scantable
325 Vector<Float> tsysvec(1);
326 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
327 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
328 uInt polno =0;
329 for ( uInt i=0; i< npol; ++i ) {
330 tsysvec = pksrec.tsys(i);
331 if (isGBTFITS) {
332 polno = pksrec.polNo ;
333 } else {
334 polno = i;
335 }
[1828]336 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
337 pksrec.beamNo-1);
[1795]338 setSpectrum(pksrec.spectra.column(i), pksrec.flagged.column(i), tsysvec);
339 commitRow();
340 }
341 if ( haveXPol_[0] ) {
342 // no tsys given for xpol, so emulate it
343 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
344 // add real part of cross pol
345 polno = 2;
346 Vector<Float> r(real(pksrec.xPol));
347 // make up flags from linears
348 /// @fixme this has to be a bitwise or of both pols
349 /// pksrec.flagged.column(0) | pksrec.flagged.column(1);
350
[1831]351 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
352 pksrec.beamNo-1);
[1795]353 setSpectrum(r, pksrec.flagged.column(0), tsysvec);
354 commitRow();
355
356 // ad imaginary part of cross pol
357 polno = 3;
358 Vector<Float> im(imag(pksrec.xPol));
[1830]359 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
360 pksrec.beamNo-1);
[1795]361 setSpectrum(im, pksrec.flagged.column(0), tsysvec);
362 commitRow();
363 }
364#ifdef HAS_ALMA
365 fillpm._update(n);
366#endif
367 }
368 if (status > 0) {
369 close();
370 throw(AipsError("Reading error occured, data possibly corrupted."));
371 }
372#ifdef HAS_ALMA
373 fillpm.done();
374#endif
375}
376
377}//namespace asap
Note: See TracBrowser for help on using the repository browser.