source: branches/newfiller/src/PKSFiller.cpp@ 1998

Last change on this file since 1998 was 1816, checked in by Malte Marquarding, 14 years ago

added setScanrate. some more logging. added setFlagrow

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