source: trunk/src/STFiller.cpp @ 1140

Last change on this file since 1140 was 1139, checked in by mar637, 18 years ago

temporary hack for import of methanol multibeam data

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.7 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    }
[1139]168  } else {
169    // hack Multibeam Methanol until pksreader is patched
170     if ( (header_->obstype) == "MX" && header_->nbeam  == 7 ) {
171        pushLog("Guessing this is Methanol Multibeam Data .\n"
172                "Only importing first IF...");
173        ifs = False;
174        ifs[0] = True;
175        header_->nif = 1;
176        nIF_ = 1;
177        ifOffset_ = 0;
178     }
[332]179  }
[405]180
[332]181  beamOffset_ = 0;
182  if (whichBeam>=0) {
[805]183    if (whichBeam>=0 && whichBeam<nBeam_) {
[1060]184      beams = False;
185      beams(whichBeam) = True;
[805]186      header_->nbeam = 1;
187      nBeam_ = 1;
188      beamOffset_ = whichBeam;
189    } else {
190      delete reader_;
191      reader_ = 0;
192      delete header_;
193      header_ = 0;
194      throw(AipsError("Illegal Beam selection"));
195    }
[332]196  }
197  Vector<Int> start(nIF_, 1);
198  Vector<Int> end(nIF_, 0);
[1060]199  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
[846]200  table_->setHeader(*header_);
[805]201}
[405]202
[805]203void STFiller::close( )
204{
205  delete reader_;reader_=0;
206  delete header_;header_=0;
207  table_ = 0;
[2]208}
209
[805]210int asap::STFiller::read( )
211{
[87]212  int status = 0;
213
[17]214  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
[87]215  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
[2]216    humidity, parAngle, pressure, temperature, windAz, windSpeed;
[73]217  Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
[1072]218  String          fieldName, srcName, tcalTime, obsType;
[2]219  Vector<Float>   calFctr, sigma, tcal, tsys;
220  Matrix<Float>   baseLin, baseSub;
221  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2);
222  Matrix<Float>   spectra;
223  Matrix<uChar>   flagtra;
224  Complex         xCalFctr;
225  Vector<Complex> xPol;
[805]226  while ( status == 0 ) {
227    status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
[1072]228                          srcName, srcDir, srcPM, srcVel, obsType, IFno,
229                          refFreq, bandwidth, freqInc, restFreq, tcal, tcalTime,
[805]230                          azimuth, elevation, parAngle, focusAxi,
231                          focusTan, focusRot, temperature, pressure,
232                          humidity, windSpeed, windAz, refBeam,
233                          beamNo, direction, scanRate,
234                          tsys, sigma, calFctr, baseLin, baseSub,
235                          spectra, flagtra, xCalFctr, xPol);
236    if ( status != 0 ) break;
[1081]237    Regex filterrx(".*[SL|PA]$");
[1110]238    Regex obsrx("^AT.+");
239    if ( header_->antennaname.matches(obsrx) && obsType.matches(filterrx)) {
240        //cerr << "ignoring paddle scan" << endl;
[1081]241        continue;
242    }
[805]243    TableRow row(table_->table());
244    TableRecord& rec = row.record();
[999]245    // fields that don't get used and are just passed through asap
246    RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
247    *srateCol = scanRate;
248    RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
249    *spmCol = srcPM;
250    RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
251    *sdirCol = srcDir;
252    RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
253    *svelCol = srcVel;
254    // the real stuff
[972]255    RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
256    *fitCol = -1;
[805]257    RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
258    *scanoCol = scanNo-1;
259    RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
260    *cyclenoCol = cycleNo-1;
261    RecordFieldPtr<Double> mjdCol(rec, "TIME");
262    *mjdCol = mjd;
263    RecordFieldPtr<Double> intCol(rec, "INTERVAL");
264    *intCol = interval;
265    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
266    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
267    // try to auto-identify if it is on or off.
[922]268    Regex rx(".*[e|w|_R]$");
[942]269    Regex rx2("_S$");
[805]270    Int match = srcName.matches(rx);
271    if (match) {
272      *srcnCol = srcName;
273    } else {
274      *srcnCol = srcName.before(rx2);
275    }
276    //*srcnCol = srcName;//.before(rx2);
277    *srctCol = match;
278    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
279    *beamCol = beamNo-beamOffset_-1;
280    RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
281    Int rb = -1;
282    if (nBeam_ > 1 ) rb = refBeam-1;
283    *rbCol = rb;
284    RecordFieldPtr<uInt> ifCol(rec, "IFNO");
285    *ifCol = IFno-ifOffset_- 1;
286    uInt id;
287    /// @todo this has to change when nchan isn't global anymore
288    id = table_->frequencies().addEntry(Double(header_->nchan/2),
289                                            refFreq, freqInc);
290    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
291    *mfreqidCol = id;
[25]292
[805]293    id = table_->molecules().addEntry(restFreq);
294    RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
295    *molidCol = id;
[414]296
[805]297    id = table_->tcal().addEntry(tcalTime, tcal);
298    RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
299    *mcalidCol = id;
300    id = table_->weather().addEntry(temperature, pressure, humidity,
301                                    windSpeed, windAz);
302    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
303    *mweatheridCol = id;
304    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
305    id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
306    *mfocusidCol = id;
307    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
308    *dirCol = direction;
[922]309    RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
[805]310    *azCol = azimuth;
[922]311    RecordFieldPtr<Float> elCol(rec, "ELEVATION");
[805]312    *elCol = elevation;
[405]313
[805]314    RecordFieldPtr<Float> parCol(rec, "PARANGLE");
315    *parCol = parAngle;
[332]316
[805]317    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
318    RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
319    RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
[405]320
[805]321    RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
322    // Turn the (nchan,npol) matrix and possible complex xPol vector
323    // into 2-4 rows in the scantable
324    Vector<Float> tsysvec(1);
325    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
326    uInt npol = (spectra.ncolumn()==1 ? 1: 2);
327    for ( uInt i=0; i< npol; ++i ) {
328      tsysvec = tsys(i);
329      *tsysCol = tsysvec;
330      *polnoCol = i;
[405]331
[805]332      *specCol = spectra.column(i);
333      *flagCol = flagtra.column(i);
334      table_->table().addRow();
335      row.put(table_->table().nrow()-1, rec);
[2]336    }
[1060]337    if ( haveXPol_[0] ) {
[805]338      // no tsys given for xpol, so emulate it
339      tsysvec = sqrt(tsys[0]*tsys[1]);
340      *tsysCol = tsysvec;
341      // add real part of cross pol
342      *polnoCol = 2;
343      Vector<Float> r(real(xPol));
344      *specCol = r;
345      // make up flags from linears
346      /// @fixme this has to be a bitwise or of both pols
347      *flagCol = flagtra.column(0);// | flagtra.column(1);
348      table_->table().addRow();
349      row.put(table_->table().nrow()-1, rec);
350      // ad imaginary part of cross pol
351      *polnoCol = 3;
352      Vector<Float> im(imag(xPol));
353      *specCol = im;
354      table_->table().addRow();
355      row.put(table_->table().nrow()-1, rec);
[2]356    }
[805]357  }
358  if (status > 0) {
359    close();
360    throw(AipsError("Reading error occured, data possibly corrupted."));
361  }
[2]362  return status;
363}
[805]364
365}//namespace asap
Note: See TracBrowser for help on using the repository browser.