Changeset 1387 for branches

Timestamp:

07/27/07 02:00:22 (17 years ago)

Author:

TakTsutsumi

Message:

merged from NRAO version of ASAP 2.1 with ALMA specific modifications

Location:

branches/alma/src

Files:

• STAttr.cpp (modified) (4 diffs)
• STFiller.cpp (modified) (9 diffs)
• STFitter.cpp (modified) (1 diff)
• STMath.cpp (modified) (3 diffs)
• STSelector.cpp (modified) (1 diff)
• STTcal.cpp (modified) (2 diffs)
• STWriter.cpp (modified) (5 diffs)
• Scantable.cpp (modified) (1 diff)
• Templates.cpp (added)
• python_Fitter.cpp (modified) (1 diff)
• python_STMath.cpp (modified) (1 diff)
• python_Scantable.cpp (modified) (2 diffs)

branches/alma/src/STAttr.cpp

@@ -73 +73 @@
    Float D = 1.0;
    switch (inst) {
+      case ALMA:
+        {
+           D = 12.0;
+        }
+        break;
+
       case ATMOPRA:
         {
@@ -92 +98 @@
         {
            D = 30.0;
+        }
+        break;
+      case GBT:
+        {
+           D = 104.9;
         }
         break;
@@ -338 +349 @@
   if (t==String("DSS-43")) {
     inst = TIDBINBILLA;
+  } else if (t==String("ALMA")) {
+    inst = ALMA;
   } else if (t==String("ATPKSMB")) {
     inst = ATPKSMB;
@@ -346 +359 @@
   } else if (t==String("CEDUNA")) {
     inst = CEDUNA;
+  } else if (t==String("GBT")) {
+    inst = GBT;
   } else if (t==String("HOBART")) {
     inst = HOBART;

branches/alma/src/STFiller.cpp

@@ -26 +26 @@
 
 #include <atnf/PKSIO/PKSreader.h>
+#include <casa/System/ProgressMeter.h>
+
 
 #include "STDefs.h"
@@ -117 +119 @@
                                   header_->freqref,
                                   header_->utc, header_->reffreq,
-                                  header_->bandwidth);
+                                  header_->bandwidth,
+                                  header_->fluxunit);
 
   if (status) {
@@ -135 +138 @@
   // Determine Telescope and set brightness unit
 
+
   Bool throwIt = False;
   Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
-  header_->fluxunit = "Jy";
+  //header_->fluxunit = "Jy";
   if (inst==ATMOPRA || inst==TIDBINBILLA) {
      header_->fluxunit = "K";
@@ -184 +188 @@
   reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
   table_->setHeader(*header_);
+  //For MS, add the location of POINTING of the input MS so one get
+  //pointing data from there, if necessary.
+  //Also find nrow in MS 
+  nInDataRow = 0;
+  if (format == "MS2") {
+    Path datapath(inName); 
+    String ptTabPath = datapath.absoluteName();
+    Table inMS(ptTabPath);
+    nInDataRow = inMS.nrow();
+    ptTabPath.append("/POINTING");
+    table_->table().rwKeywordSet().define("POINTING", ptTabPath);
+    if ((header_->antennaname).matches("GBT")) {
+      String GOTabPath = datapath.absoluteName();
+      GOTabPath.append("/GBT_GO");
+      table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
+    }
+  }
+     
 }
 
@@ -209 +231 @@
   Complex         xCalFctr;
   Vector<Complex> xPol;
+  Double min = 0.0;
+  Double max = nInDataRow;
+  ProgressMeter fillpm(min, max, "Data importing progress");
+  int n = 0;
   while ( status == 0 ) {
     status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
@@ -220 +246 @@
                           spectra, flagtra, xCalFctr, xPol);
     if ( status != 0 ) break;
+    n += 1;
+    
     Regex filterrx(".*[SL|PA]$");
     Regex obsrx("^AT.+");
@@ -250 +278 @@
     RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
     RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
+    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
+    *fieldnCol = fieldName;
     // try to auto-identify if it is on or off.
     Regex rx(".*[e|w|_R]$");
@@ -340 +370 @@
       row.put(table_->table().nrow()-1, rec);
     }
+    fillpm._update(n);
   }
   if (status > 0) {
@@ -345 +376 @@
     throw(AipsError("Reading error occured, data possibly corrupted."));
   }
+  fillpm.done();
   return status;
 }

branches/alma/src/STFitter.cpp

@@ -329 +329 @@
 }
 
+bool Fitter::lfit() {
+  LinearFit<Float> fitter;
+  CompoundFunction<Float> func;
+
+  uInt n = funcs_.nelements();
+  for (uInt i=0; i<n; ++i) {
+    func.addFunction(*funcs_[i]);
+  }
+
+  fitter.setFunction(func);
+  //fitter.setMaxIter(50+n*10);
+  // Convergence criterium
+  //fitter.setCriteria(0.001);
+
+  // Fit
+  Vector<Float> sigma(x_.nelements());
+  sigma = 1.0;
+
+  parameters_.resize();
+  parameters_ = fitter.fit(x_, y_, sigma, &m_);
+  std::vector<float> ps;
+  parameters_.tovector(ps);
+  setParameters(ps);
+
+  error_.resize();
+  error_ = fitter.errors();
+
+  chisquared_ = fitter.getChi2();
+
+  residual_.resize();
+  residual_ =  y_;
+  fitter.residual(residual_,x_);
+  // use fitter.residual(model=True) to get the model
+  thefit_.resize(x_.nelements());
+  fitter.residual(thefit_,x_,True);
+  return true;
+}
 
 std::vector<float> Fitter::evaluate(int whichComp) const

branches/alma/src/STMath.cpp

@@ -48 +48 @@
 #include "STAttr.h"
 #include "STMath.h"
+#include "STSelector.h"
 
 using namespace casa;
@@ -536 +537 @@
 }
 
+// dototalpower (migration of GBTIDL procedure dototalpower.pro)
+// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
+// do it for each cycles in a specific scan.
+CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
+                                              const CountedPtr< Scantable >& caloff, Float tcal )
+{
+if ( ! calon->conformant(*caloff) ) {
+    throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
+  }
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(caloff, false);
+  Table& tout = out->table();
+  const Table& tcon = calon->table();
+  Vector<Float> tcalout;
+  Vector<Float> tcalout2;  //debug
+
+  if ( tout.nrow() != tcon.nrow() ) {
+    throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
+  }
+  // iteration by scanno or cycle no.
+  TableIterator sit(tout, "SCANNO");
+  TableIterator s2it(tcon, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table toff = sit.table();
+    TableRow row(toff);
+    Table t = s2it.table();
+    ScalarColumn<Double> outintCol(toff, "INTERVAL");
+    ArrayColumn<Float> outspecCol(toff, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(toff, "TSYS");
+    ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
+    ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
+    ROScalarColumn<uInt> outpolCol(toff, "POLNO");
+    ROScalarColumn<Double> onintCol(t, "INTERVAL");
+    ROArrayColumn<Float> onspecCol(t, "SPECTRA");
+    ROArrayColumn<Float> ontsysCol(t, "TSYS");
+    ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
+    //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
+
+    for (uInt i=0; i < toff.nrow(); ++i) {
+      //skip these checks -> assumes the data order are the same between the cal on off pairs
+      //
+      Vector<Float> specCalon, specCaloff;
+      // to store scalar (mean) tsys
+      Vector<Float> tsysout(1);
+      uInt tcalId, polno;
+      Double offint, onint;
+      outpolCol.get(i, polno);
+      outspecCol.get(i, specCaloff);
+      onspecCol.get(i, specCalon);
+      Vector<uChar> flagCaloff, flagCalon;
+      outflagCol.get(i, flagCaloff);
+      onflagCol.get(i, flagCalon);
+      outtcalIdCol.get(i, tcalId);
+      outintCol.get(i, offint);
+      onintCol.get(i, onint);
+      // caluculate mean Tsys
+      uInt nchan = specCaloff.nelements();
+      // percentage of edge cut off
+      uInt pc = 10;
+      uInt bchan = nchan/pc;
+      uInt echan = nchan-bchan;
+
+      Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
+      Vector<Float> testsubsp = specCaloff(chansl);
+      MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
+      MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
+      MaskedArray<Float> spdiff = spon-spoff;
+      uInt noff = spoff.nelementsValid();
+      //uInt non = spon.nelementsValid();
+      uInt ndiff = spdiff.nelementsValid();
+      Float meantsys;
+
+/**
+      Double subspec, subdiff;
+      uInt usednchan;
+      subspec = 0;
+      subdiff = 0;
+      usednchan = 0;
+      for(uInt k=(bchan-1); k<echan; k++) {
+        subspec += specCaloff[k];
+        subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
+        ++usednchan;
+      }
+**/
+      // get tcal if input tcal <= 0
+      String tcalt;
+      Float tcalUsed;
+      tcalUsed = tcal;
+      if ( tcal <= 0.0 ) {
+        caloff->tcal().getEntry(tcalt, tcalout, tcalId);
+        if (polno<=3) {
+          tcalUsed = tcalout[polno];
+        }
+        else {
+          tcalUsed = tcalout[0];
+        }
+      }
+
+      Float meanoff;
+      Float meandiff;
+      if (noff && ndiff) {
+         //Debug
+         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
+         meanoff = sum(spoff)/noff;
+         meandiff = sum(spdiff)/ndiff;
+         meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
+      }
+      else {
+         meantsys=1;
+      }
+
+      tsysout[0] = Float(meantsys);
+      MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
+      MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
+      MaskedArray<Float> sig =   Float(0.5) * (mcaloff + mcalon);
+      //uInt ncaloff = mcaloff.nelementsValid();
+      //uInt ncalon = mcalon.nelementsValid();
+
+      outintCol.put(i, offint+onint);
+      outspecCol.put(i, sig.getArray());
+      outflagCol.put(i, flagsFromMA(sig));
+      outtsysCol.put(i, tsysout);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return out;
+}
+
+//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
+// observatiions.
+// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
+//        no smoothing).
+// output: resultant scantable [= (sig-ref/ref)*tsys]
+CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
+                                          const CountedPtr < Scantable >& ref,
+                                          int smoothref,
+                                          casa::Float tsysv,
+                                          casa::Float tau )
+{
+if ( ! ref->conformant(*sig) ) {
+    throw(AipsError("'sig' and 'ref' scantables are not conformant."));
+  }
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(sig, false);
+  CountedPtr< Scantable > smref;
+  if ( smoothref > 1 ) {
+    float fsmoothref = static_cast<float>(smoothref);
+    std::string inkernel = "boxcar";
+    smref = smooth(ref, inkernel, fsmoothref );
+    ostringstream oss;
+    oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
+    pushLog(String(oss));
+  }
+  else {
+    smref = ref;
+  }
+  Table& tout = out->table();
+  const Table& tref = smref->table();
+  if ( tout.nrow() != tref.nrow() ) {
+    throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
+  }
+  // iteration by scanno? or cycle no.
+  TableIterator sit(tout, "SCANNO");
+  TableIterator s2it(tref, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table ton = sit.table();
+    Table t = s2it.table();
+    ScalarColumn<Double> outintCol(ton, "INTERVAL");
+    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(ton, "TSYS");
+    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
+    ArrayColumn<Float> refspecCol(t, "SPECTRA");
+    ROScalarColumn<Double> refintCol(t, "INTERVAL");
+    ROArrayColumn<Float> reftsysCol(t, "TSYS");
+    ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
+    ROScalarColumn<Float> refelevCol(t, "ELEVATION");
+    for (uInt i=0; i < ton.nrow(); ++i) {
+
+      Double onint, refint;
+      Vector<Float> specon, specref;
+      // to store scalar (mean) tsys
+      Vector<Float> tsysref;
+      outintCol.get(i, onint);
+      refintCol.get(i, refint);
+      outspecCol.get(i, specon);
+      refspecCol.get(i, specref);
+      Vector<uChar> flagref, flagon;
+      outflagCol.get(i, flagon);
+      refflagCol.get(i, flagref);
+      reftsysCol.get(i, tsysref);
+
+      Float tsysrefscalar;
+      if ( tsysv > 0.0 ) {
+        ostringstream oss;
+        Float elev;
+        refelevCol.get(i, elev);
+        oss << "user specified Tsys = " << tsysv;
+        // do recalc elevation if EL = 0
+        if ( elev == 0 ) {
+          throw(AipsError("EL=0, elevation data is missing."));
+        } else {
+          if ( tau <= 0.0 ) {
+            throw(AipsError("Valid tau is not supplied."));
+          } else {
+            tsysrefscalar = tsysv * exp(tau/elev);
+          }
+        }
+        oss << ", corrected (for El) tsys= "<<tsysrefscalar;
+        pushLog(String(oss));
+      }
+      else {
+        tsysrefscalar = tsysref[0];
+      }
+      //get quotient spectrum
+      MaskedArray<Float> mref = maskedArray(specref, flagref);
+      MaskedArray<Float> mon = maskedArray(specon, flagon);
+      MaskedArray<Float> specres =   tsysrefscalar*((mon - mref)/mref);
+      Double resint = onint*refint*smoothref/(onint+refint*smoothref);
+
+      //Debug
+      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
+      // fill the result, replay signal tsys by reference tsys
+      outintCol.put(i, resint);
+      outspecCol.put(i, specres.getArray());
+      outflagCol.put(i, flagsFromMA(specres));
+      outtsysCol.put(i, tsysref);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
+                                     const std::vector<int>& scans,
+                                     int smoothref,
+                                     casa::Float tsysv,
+                                     casa::Float tau,
+                                     casa::Float tcal )
+
+{
+  setInsitu(false);
+  STSelector sel;
+  std::vector<int> scan1, scan2, beams;
+  std::vector< vector<int> > scanpair;
+  std::vector<string> calstate;
+  String msg;
+
+  CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
+  CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
+
+  std::vector< CountedPtr< Scantable > > sctables;
+  sctables.push_back(s1b1on);
+  sctables.push_back(s1b1off);
+  sctables.push_back(s1b2on);
+  sctables.push_back(s1b2off);
+  sctables.push_back(s2b1on);
+  sctables.push_back(s2b1off);
+  sctables.push_back(s2b2on);
+  sctables.push_back(s2b2off);
+
+  //check scanlist
+  int n=s->checkScanInfo(scans);
+  if (n==1) {
+     throw(AipsError("Incorrect scan pairs. "));
+  }
+
+  // Assume scans contain only a pair of consecutive scan numbers.
+  // It is assumed that first beam, b1,  is on target.
+  // There is no check if the first beam is on or not.
+  if ( scans.size()==1 ) {
+    scan1.push_back(scans[0]);
+    scan2.push_back(scans[0]+1);
+  } else if ( scans.size()==2 ) {
+   scan1.push_back(scans[0]);
+   scan2.push_back(scans[1]);
+  } else {
+    if ( scans.size()%2 == 0 ) {
+      for (uInt i=0; i<scans.size(); i++) {
+        if (i%2 == 0) {
+          scan1.push_back(scans[i]);
+        }
+        else {
+          scan2.push_back(scans[i]);
+        }
+      }
+    } else {
+      throw(AipsError("Odd numbers of scans, cannot form pairs."));
+    }
+  }
+  scanpair.push_back(scan1);
+  scanpair.push_back(scan2);
+  calstate.push_back("*calon");
+  calstate.push_back("*[^calon]");
+  CountedPtr< Scantable > ws = getScantable(s, false);
+  uInt l=0;
+  while ( l < sctables.size() ) {
+    for (uInt i=0; i < 2; i++) {
+      for (uInt j=0; j < 2; j++) {
+        for (uInt k=0; k < 2; k++) {
+          sel.reset();
+          sel.setScans(scanpair[i]);
+          sel.setName(calstate[k]);
+          beams.clear();
+          beams.push_back(j);
+          sel.setBeams(beams);
+          ws->setSelection(sel);
+          sctables[l]= getScantable(ws, false);
+          l++;
+        }
+      }
+    }
+  }
+
+  // replace here by splitData or getData functionality
+  CountedPtr< Scantable > sig1;
+  CountedPtr< Scantable > ref1;
+  CountedPtr< Scantable > sig2;
+  CountedPtr< Scantable > ref2;
+  CountedPtr< Scantable > calb1;
+  CountedPtr< Scantable > calb2;
+
+  msg=String("Processing dototalpower for subset of the data");
+  ostringstream oss1;
+  oss1 << msg  << endl;
+  pushLog(String(oss1));
+  // Debug for IRC CS data
+  //float tcal1=7.0;
+  //float tcal2=4.0;
+  sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
+  ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
+  ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
+  sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
+
+  // correction of user-specified tsys for elevation here
+
+  // dosigref calibration
+  msg=String("Processing dosigref for subset of the data");
+  ostringstream oss2;
+  oss2 << msg  << endl;
+  pushLog(String(oss2));
+  calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
+  calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
+
+  // iteration by scanno or cycle no.
+  Table& tcalb1 = calb1->table();
+  Table& tcalb2 = calb2->table();
+  TableIterator sit(tcalb1, "SCANNO");
+  TableIterator s2it(tcalb2, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table t1 = sit.table();
+    Table t2= s2it.table();
+    ArrayColumn<Float> outspecCol(t1, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(t1, "TSYS");
+    ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
+    ScalarColumn<Double> outintCol(t1, "INTERVAL");
+    ArrayColumn<Float> t2specCol(t2, "SPECTRA");
+    ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
+    ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
+    ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
+    for (uInt i=0; i < t1.nrow(); ++i) {
+      Vector<Float> spec1, spec2;
+      // to store scalar (mean) tsys
+      Vector<Float> tsys1, tsys2;
+      Vector<uChar> flag1, flag2;
+      Double tint1, tint2;
+      outspecCol.get(i, spec1);
+      t2specCol.get(i, spec2);
+      outflagCol.get(i, flag1);
+      t2flagCol.get(i, flag2);
+      outtsysCol.get(i, tsys1);
+      t2tsysCol.get(i, tsys2);
+      outintCol.get(i, tint1);
+      t2intCol.get(i, tint2);
+      // average
+      // assume scalar tsys for weights
+      Float wt1, wt2, tsyssq1, tsyssq2;
+      tsyssq1 = tsys1[0]*tsys1[0];
+      tsyssq2 = tsys2[0]*tsys2[0];
+      wt1 = Float(tint1)/tsyssq1;
+      wt2 = Float(tint2)/tsyssq2;
+      Float invsumwt=1/(wt1+wt2);
+      MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
+      MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
+      MaskedArray<Float> avspec =  invsumwt * (wt1*mspec1 + wt2*mspec2);
+      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
+      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
+      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
+      Array<Float> avtsys =  tsys1;
+
+      outspecCol.put(i, avspec.getArray());
+      outflagCol.put(i, flagsFromMA(avspec));
+      outtsysCol.put(i, avtsys);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return calb1;
+}
+
+//GBTIDL version of frequency switched data calibration
+CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
+                                      const std::vector<int>& scans,
+                                      int smoothref,
+                                      casa::Float tsysv,
+                                      casa::Float tau,
+                                      casa::Float tcal )
+{
+
+
+  STSelector sel;
+  CountedPtr< Scantable > ws = getScantable(s, false);
+  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
+  CountedPtr< Scantable > calsig, calref, out;
+
+  //split the data
+  sel.setName("*_fs");
+  ws->setSelection(sel);
+  sig = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fs_calon");
+  ws->setSelection(sel);
+  sigwcal = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fsr");
+  ws->setSelection(sel);
+  ref = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fsr_calon");
+  ws->setSelection(sel);
+  refwcal = getScantable(ws,false);
+
+  calsig = dototalpower(sigwcal, sig, tcal=tcal);
+  calref = dototalpower(refwcal, ref, tcal=tcal);
+
+  out=dosigref(calsig,calref,smoothref,tsysv,tau);
+
+  return out;
+}
+
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
@@ -1250 +1692 @@
   vec = 0;
   pols->table_.rwKeywordSet().define("nPol", Int(1));
-  pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
+  //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
+  pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
   std::vector<CountedPtr<Scantable> > vpols;
   vpols.push_back(pols);

branches/alma/src/STSelector.cpp

@@ -83 +83 @@
 void asap::STSelector::setName( const std::string& sname )
 {
-  std::string sql = "SELECT from $1 WHERE SRCNAME == pattern('"+sname+"')";
+  std::string sql = "SELECT FROM $1 WHERE SRCNAME == pattern('"+sname+"')";
   setTaQL(sql);
 }

branches/alma/src/STTcal.cpp

@@ -68 +68 @@
 }
 
+/*** rewrite this for handling of GBT data 
 uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
 {
@@ -90 +91 @@
   return resultid;
 }
+***/
+
+uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
+{
+  // test if this already exists
+  // TT - different Tcal values for each polarization, feed, and
+  // data description. So there may be multiple entries for the same
+  // time stamp.
+  uInt resultid;
+  uInt rno = table_.nrow();
+  //table_.addRow();
+  // get last assigned tcal_id and increment
+  if ( rno == 0 ) {
+    resultid = 0;
+  }
+  else {
+    idCol_.get(rno-1, resultid);
+    resultid++;
+  }
+  table_.addRow();
+  tcalCol_.put(rno, cal);
+  timeCol_.put(rno, time);
+  idCol_.put(rno, resultid);
+  return resultid;
+}
+
 
 void STTcal::getEntry( String& time, Vector<Float>& tcal, uInt id )

branches/alma/src/STWriter.cpp

@@ -116 +116 @@
 
   // Extract the header from the table.
+  // this is a little different from what I have done
+  // before. Need to check with the Offline User Test data
   STHeader hdr = in->getHeader();
   //const Int nPol  = hdr.npol;
@@ -123 +125 @@
   Vector<uInt> nPol(nIF),nChan(nIF);
   Vector<Bool> havexpol(nIF);
+  String fluxUnit = hdr.fluxunit;
+
   nPol = 0;nChan = 0; havexpol = False;
   for (uint i=0;i<ifs.size();++i) {
@@ -137 +141 @@
                                hdr.antennaname, hdr.antennaposition,
                                hdr.obstype, hdr.equinox, hdr.freqref,
-                               nChan, nPol, havexpol, False);
+                               nChan, nPol, havexpol, False, fluxUnit);
   if ( status ) {
     throw(AipsError("Failed to create output file"));
@@ -265 +269 @@
   pushLog(String(oss));
   writer_->close();
-
+  //if MS2 delete POINTING table exists and copy the one in the keyword
+  if ( format_ == "MS2" ) {
+    replacePtTab(table, filename);
+  }
   return 0;
 }
@@ -315 +322 @@
 }
 
-
-}
+// For writing MS data, if there is the reference to
+// original pointing table it replace it by it.
+void STWriter::replacePtTab (const Table& tab, const std::string& fname)
+{
+  String oldPtTabName = fname;
+  oldPtTabName.append("/POINTING");
+  if ( tab.keywordSet().isDefined("POINTING") ) {
+    String PointingTab = tab.keywordSet().asString("POINTING");
+    if ( Table::isReadable(PointingTab) ) {
+      Table newPtTab(PointingTab, Table::Old);
+      newPtTab.copy(oldPtTabName, Table::New);
+      ostringstream oss;
+      oss << "STWriter: copied  " <<PointingTab  << " to " << fname;
+      pushLog(String(oss));
+    }
+  }
+}
+
+}

branches/alma/src/Scantable.cpp

@@ -1028 +1028 @@
 }
 
+std::string asap::Scantable::getAntennaName() const
+{
+  String out;
+  table_.keywordSet().get("AntennaName", out);
+  return out;
+}
+
+int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
+{
+  String tbpath;
+  int ret = 0;
+  if ( table_.keywordSet().isDefined("GBT_GO") ) {
+    table_.keywordSet().get("GBT_GO", tbpath);
+    Table t(tbpath,Table::Old);
+    // check each scan if other scan of the pair exist
+    int nscan = scanlist.size();
+    for (int i = 0; i < nscan; i++) {
+      Table subt = t( t.col("SCAN") == scanlist[i]+1 );
+      if (subt.nrow()==0) {
+        cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
+        ret = 1;
+        break;
+      }
+      ROTableRow row(subt);
+      const TableRecord& rec = row.get(0);
+      int scan1seqn = rec.asuInt("PROCSEQN");
+      int laston1 = rec.asuInt("LASTON");
+      if ( rec.asuInt("PROCSIZE")==2 ) {
+        if ( i < nscan-1 ) {
+          Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
+          if ( subt2.nrow() == 0) {
+            cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
+            ret = 1;
+            break;
+          }
+          ROTableRow row2(subt2);
+          const TableRecord& rec2 = row2.get(0);
+          int scan2seqn = rec2.asuInt("PROCSEQN");
+          int laston2 = rec2.asuInt("LASTON");
+          if (scan1seqn == 1 && scan2seqn == 2) {
+            if (laston1 == laston2) {
+              cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              i +=1;
+            }
+            else {
+              cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+            }
+          }
+          else if (scan1seqn==2 && scan2seqn == 1) {
+            if (laston1 == laston2) {
+              cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
+              ret = 1;
+              break;
+            }
+          }
+          else {
+            cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
+            ret = 1;
+            break;
+          }
+        }
+      }
+      else {
+        cerr<<"The scan does not appear to be standard obsevation."<<endl;
+      }
+    //if ( i >= nscan ) break;
+    }
+  }
+  else {
+    cerr<<"No reference to GBT_GO table."<<endl;
+    ret = 1;
+  }
+  return ret;
+}
+
+std::vector<double>  asap::Scantable::getDirectionVector(int whichrow) const
+{
+  Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
+  std::vector<double> dir;
+  Dir.tovector(dir);
+  return dir;
+}
+
 }
  //namespace asap

branches/alma/src/python_Fitter.cpp

@@ -55 +55 @@
         .def("reset", &Fitter::reset)
         .def("fit", &Fitter::fit)
+        .def("lfit", &Fitter::lfit)
         .def("evaluate", &Fitter::evaluate)
       ;

branches/alma/src/python_STMath.cpp

@@ -52 +52 @@
         .def("_auto_quotient", &STMathWrapper::autoQuotient)
         .def("_quotient", &STMathWrapper::quotient)
+        .def("_dototalpower", &STMathWrapper::dototalpower)
+        .def("_dosigref", &STMathWrapper::dosigref)
+        .def("_donod", &STMathWrapper::donod)
+        .def("_dofs", &STMathWrapper::dofs)
         .def("_stats", &STMathWrapper::statistic)
         .def("_freqswitch", &STMathWrapper::freqSwitch)

branches/alma/src/python_Scantable.cpp

@@ -102 +102 @@
     .def("_getdirection", &ScantableWrapper::getDirectionString,
          (boost::python::arg("whichrow")=0) )
+    .def("get_antennaname", &ScantableWrapper::getAntennaName)
     .def("_flag", &ScantableWrapper::flag)
     .def("_save",  &ScantableWrapper::makePersistent)
@@ -121 +122 @@
     .def("_recalcazel", &ScantableWrapper::calculateAZEL)
     .def("_setsourcetype", &ScantableWrapper::setSourceType)
+    .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
   ;
 };