source: trunk/src/STFiller.cpp @ 1439

Last change on this file since 1439 was 1438, checked in by Malte Marquarding, 16 years ago

allow ALMA specific build

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