source: branches/alma/src/STFiller.cpp @ 1451

Last change on this file since 1451 was 1446, checked in by TakTsutsumi, 16 years ago

Merged recent updates (since 2007) from nrao-asap

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