Changeset 1391 for trunk/src


Ignore:
Timestamp:
07/30/07 11:59:36 (17 years ago)
Author:
Malte Marquarding
Message:

merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.

Location:
trunk/src
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STAttr.cpp

    r1345 r1391  
    7373   Float D = 1.0;
    7474   switch (inst) {
     75      case ALMA:
     76        {
     77           D = 12.0;
     78        }
     79        break;
     80
    7581      case ATMOPRA:
    7682        {
    … …  
    9298        {
    9399           D = 30.0;
     100        }
     101        break;
     102      case GBT:
     103        {
     104           D = 104.9;
    94105        }
    95106        break;
    … …  
    338349  if (t==String("DSS-43")) {
    339350    inst = TIDBINBILLA;
     351  } else if (t==String("ALMA")) {
     352    inst = ALMA;
    340353  } else if (t==String("ATPKSMB")) {
    341354    inst = ATPKSMB;
    … …  
    346359  } else if (t==String("CEDUNA")) {
    347360    inst = CEDUNA;
     361  } else if (t==String("GBT")) {
     362    inst = GBT;
    348363  } else if (t==String("HOBART")) {
    349364    inst = HOBART;

  • trunk/src/STDefs.h

    r1103 r1391  
    4141
    4242  enum Instrument {UNKNOWNINST=0,
     43                   ALMA,
    4344                   ATPKSMB,
    4445                   ATPKSHOH,
    … …  
    4647                   TIDBINBILLA,
    4748                   CEDUNA,
     49                   GBT,
    4850                   HOBART,
    4951                   N_INSTRUMENTS};

  • trunk/src/STFiller.cpp

    r1188 r1391  
    2626
    2727#include <atnf/PKSIO/PKSreader.h>
     28//#include <casa/System/ProgressMeter.h>
     29
    2830
    2931#include "STDefs.h"
    … …  
    134136  }
    135137  // Determine Telescope and set brightness unit
     138
    136139
    137140  Bool throwIt = False;
    … …  
    184187  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
    185188  table_->setHeader(*header_);
     189  //For MS, add the location of POINTING of the input MS so one get
     190  //pointing data from there, if necessary.
     191  //Also find nrow in MS
     192  nInDataRow = 0;
     193  if (format == "MS2") {
     194    Path datapath(inName);
     195    String ptTabPath = datapath.absoluteName();
     196    Table inMS(ptTabPath);
     197    nInDataRow = inMS.nrow();
     198    ptTabPath.append("/POINTING");
     199    table_->table().rwKeywordSet().define("POINTING", ptTabPath);
     200    if ((header_->antennaname).matches("GBT")) {
     201      String GOTabPath = datapath.absoluteName();
     202      GOTabPath.append("/GBT_GO");
     203      table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
     204    }
     205  }
     206     
    186207}
    187208
    … …  
    209230  Complex         xCalFctr;
    210231  Vector<Complex> xPol;
     232  Double min = 0.0;
     233  Double max = nInDataRow;
     234  //ProgressMeter fillpm(min, max, "Data importing progress");
     235  int n = 0;
    211236  while ( status == 0 ) {
    212237    status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
    … …  
    220245                          spectra, flagtra, xCalFctr, xPol);
    221246    if ( status != 0 ) break;
     247    n += 1;
     248   
    222249    Regex filterrx(".*[SL|PA]$");
    223250    Regex obsrx("^AT.+");
    … …  
    250277    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
    251278    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
     279    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
     280    *fieldnCol = fieldName;
    252281    // try to auto-identify if it is on or off.
    253282    Regex rx(".*[e|w|_R]$");
    … …  
    340369      row.put(table_->table().nrow()-1, rec);
    341370    }
     371    //fillpm._update(n);
    342372  }
    343373  if (status > 0) {
    … …  
    345375    throw(AipsError("Reading error occured, data possibly corrupted."));
    346376  }
     377  //fillpm.done();
    347378  return status;
    348379}

  • trunk/src/STFiller.h

    r1353 r1391  
    9999  casa::String filename_;
    100100  casa::CountedPtr< Scantable > table_;
    101   casa::Int nIF_, nBeam_, nPol_, nChan_;
     101  casa::Int nIF_, nBeam_, nPol_, nChan_, nInDataRow;
    102102  casa::uInt ifOffset_, beamOffset_;
    103103  casa::Vector<casa::Bool> haveXPol_;

  • trunk/src/STFitter.cpp

    r1232 r1391  
    329329}
    330330
     331bool Fitter::lfit() {
     332  LinearFit<Float> fitter;
     333  CompoundFunction<Float> func;
     334
     335  uInt n = funcs_.nelements();
     336  for (uInt i=0; i<n; ++i) {
     337    func.addFunction(*funcs_[i]);
     338  }
     339
     340  fitter.setFunction(func);
     341  //fitter.setMaxIter(50+n*10);
     342  // Convergence criterium
     343  //fitter.setCriteria(0.001);
     344
     345  // Fit
     346  Vector<Float> sigma(x_.nelements());
     347  sigma = 1.0;
     348
     349  parameters_.resize();
     350  parameters_ = fitter.fit(x_, y_, sigma, &m_);
     351  std::vector<float> ps;
     352  parameters_.tovector(ps);
     353  setParameters(ps);
     354
     355  error_.resize();
     356  error_ = fitter.errors();
     357
     358  chisquared_ = fitter.getChi2();
     359
     360  residual_.resize();
     361  residual_ =  y_;
     362  fitter.residual(residual_,x_);
     363  // use fitter.residual(model=True) to get the model
     364  thefit_.resize(x_.nelements());
     365  fitter.residual(thefit_,x_,True);
     366  return true;
     367}
    331368
    332369std::vector<float> Fitter::evaluate(int whichComp) const

  • trunk/src/STFitter.h

    r1028 r1391  
    2727//#                        AUSTRALIA
    2828//#
    29 //# $Id:$
     29//# $Id$
    3030//#---------------------------------------------------------------------------
    3131#ifndef STFITTER_H
    … …  
    6464  void reset();
    6565  bool fit();
     66  // Fit via linear method
     67  bool lfit();
    6668  bool computeEstimate();
    6769

  • trunk/src/STMath.cpp

    r1384 r1391  
    4848#include "STAttr.h"
    4949#include "STMath.h"
     50#include "STSelector.h"
    5051
    5152using namespace casa;
    … …  
    536537}
    537538
     539// dototalpower (migration of GBTIDL procedure dototalpower.pro)
     540// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
     541// do it for each cycles in a specific scan.
     542CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
     543                                              const CountedPtr< Scantable >& caloff, Float tcal )
     544{
     545if ( ! calon->conformant(*caloff) ) {
     546    throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
     547  }
     548  setInsitu(false);
     549  CountedPtr< Scantable > out = getScantable(caloff, false);
     550  Table& tout = out->table();
     551  const Table& tcon = calon->table();
     552  Vector<Float> tcalout;
     553  Vector<Float> tcalout2;  //debug
     554
     555  if ( tout.nrow() != tcon.nrow() ) {
     556    throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
     557  }
     558  // iteration by scanno or cycle no.
     559  TableIterator sit(tout, "SCANNO");
     560  TableIterator s2it(tcon, "SCANNO");
     561  while ( !sit.pastEnd() ) {
     562    Table toff = sit.table();
     563    TableRow row(toff);
     564    Table t = s2it.table();
     565    ScalarColumn<Double> outintCol(toff, "INTERVAL");
     566    ArrayColumn<Float> outspecCol(toff, "SPECTRA");
     567    ArrayColumn<Float> outtsysCol(toff, "TSYS");
     568    ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
     569    ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
     570    ROScalarColumn<uInt> outpolCol(toff, "POLNO");
     571    ROScalarColumn<Double> onintCol(t, "INTERVAL");
     572    ROArrayColumn<Float> onspecCol(t, "SPECTRA");
     573    ROArrayColumn<Float> ontsysCol(t, "TSYS");
     574    ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
     575    //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
     576
     577    for (uInt i=0; i < toff.nrow(); ++i) {
     578      //skip these checks -> assumes the data order are the same between the cal on off pairs
     579      //
     580      Vector<Float> specCalon, specCaloff;
     581      // to store scalar (mean) tsys
     582      Vector<Float> tsysout(1);
     583      uInt tcalId, polno;
     584      Double offint, onint;
     585      outpolCol.get(i, polno);
     586      outspecCol.get(i, specCaloff);
     587      onspecCol.get(i, specCalon);
     588      Vector<uChar> flagCaloff, flagCalon;
     589      outflagCol.get(i, flagCaloff);
     590      onflagCol.get(i, flagCalon);
     591      outtcalIdCol.get(i, tcalId);
     592      outintCol.get(i, offint);
     593      onintCol.get(i, onint);
     594      // caluculate mean Tsys
     595      uInt nchan = specCaloff.nelements();
     596      // percentage of edge cut off
     597      uInt pc = 10;
     598      uInt bchan = nchan/pc;
     599      uInt echan = nchan-bchan;
     600
     601      Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
     602      Vector<Float> testsubsp = specCaloff(chansl);
     603      MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
     604      MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
     605      MaskedArray<Float> spdiff = spon-spoff;
     606      uInt noff = spoff.nelementsValid();
     607      //uInt non = spon.nelementsValid();
     608      uInt ndiff = spdiff.nelementsValid();
     609      Float meantsys;
     610
     611/**
     612      Double subspec, subdiff;
     613      uInt usednchan;
     614      subspec = 0;
     615      subdiff = 0;
     616      usednchan = 0;
     617      for(uInt k=(bchan-1); k<echan; k++) {
     618        subspec += specCaloff[k];
     619        subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
     620        ++usednchan;
     621      }
     622**/
     623      // get tcal if input tcal <= 0
     624      String tcalt;
     625      Float tcalUsed;
     626      tcalUsed = tcal;
     627      if ( tcal <= 0.0 ) {
     628        caloff->tcal().getEntry(tcalt, tcalout, tcalId);
     629        if (polno<=3) {
     630          tcalUsed = tcalout[polno];
     631        }
     632        else {
     633          tcalUsed = tcalout[0];
     634        }
     635      }
     636
     637      Float meanoff;
     638      Float meandiff;
     639      if (noff && ndiff) {
     640         //Debug
     641         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
     642         meanoff = sum(spoff)/noff;
     643         meandiff = sum(spdiff)/ndiff;
     644         meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
     645      }
     646      else {
     647         meantsys=1;
     648      }
     649
     650      tsysout[0] = Float(meantsys);
     651      MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
     652      MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
     653      MaskedArray<Float> sig =   Float(0.5) * (mcaloff + mcalon);
     654      //uInt ncaloff = mcaloff.nelementsValid();
     655      //uInt ncalon = mcalon.nelementsValid();
     656
     657      outintCol.put(i, offint+onint);
     658      outspecCol.put(i, sig.getArray());
     659      outflagCol.put(i, flagsFromMA(sig));
     660      outtsysCol.put(i, tsysout);
     661    }
     662    ++sit;
     663    ++s2it;
     664  }
     665  return out;
     666}
     667
     668//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
     669// observatiions.
     670// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
     671//        no smoothing).
     672// output: resultant scantable [= (sig-ref/ref)*tsys]
     673CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
     674                                          const CountedPtr < Scantable >& ref,
     675                                          int smoothref,
     676                                          casa::Float tsysv,
     677                                          casa::Float tau )
     678{
     679if ( ! ref->conformant(*sig) ) {
     680    throw(AipsError("'sig' and 'ref' scantables are not conformant."));
     681  }
     682  setInsitu(false);
     683  CountedPtr< Scantable > out = getScantable(sig, false);
     684  CountedPtr< Scantable > smref;
     685  if ( smoothref > 1 ) {
     686    float fsmoothref = static_cast<float>(smoothref);
     687    std::string inkernel = "boxcar";
     688    smref = smooth(ref, inkernel, fsmoothref );
     689    ostringstream oss;
     690    oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
     691    pushLog(String(oss));
     692  }
     693  else {
     694    smref = ref;
     695  }
     696  Table& tout = out->table();
     697  const Table& tref = smref->table();
     698  if ( tout.nrow() != tref.nrow() ) {
     699    throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
     700  }
     701  // iteration by scanno? or cycle no.
     702  TableIterator sit(tout, "SCANNO");
     703  TableIterator s2it(tref, "SCANNO");
     704  while ( !sit.pastEnd() ) {
     705    Table ton = sit.table();
     706    Table t = s2it.table();
     707    ScalarColumn<Double> outintCol(ton, "INTERVAL");
     708    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
     709    ArrayColumn<Float> outtsysCol(ton, "TSYS");
     710    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
     711    ArrayColumn<Float> refspecCol(t, "SPECTRA");
     712    ROScalarColumn<Double> refintCol(t, "INTERVAL");
     713    ROArrayColumn<Float> reftsysCol(t, "TSYS");
     714    ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
     715    ROScalarColumn<Float> refelevCol(t, "ELEVATION");
     716    for (uInt i=0; i < ton.nrow(); ++i) {
     717
     718      Double onint, refint;
     719      Vector<Float> specon, specref;
     720      // to store scalar (mean) tsys
     721      Vector<Float> tsysref;
     722      outintCol.get(i, onint);
     723      refintCol.get(i, refint);
     724      outspecCol.get(i, specon);
     725      refspecCol.get(i, specref);
     726      Vector<uChar> flagref, flagon;
     727      outflagCol.get(i, flagon);
     728      refflagCol.get(i, flagref);
     729      reftsysCol.get(i, tsysref);
     730
     731      Float tsysrefscalar;
     732      if ( tsysv > 0.0 ) {
     733        ostringstream oss;
     734        Float elev;
     735        refelevCol.get(i, elev);
     736        oss << "user specified Tsys = " << tsysv;
     737        // do recalc elevation if EL = 0
     738        if ( elev == 0 ) {
     739          throw(AipsError("EL=0, elevation data is missing."));
     740        } else {
     741          if ( tau <= 0.0 ) {
     742            throw(AipsError("Valid tau is not supplied."));
     743          } else {
     744            tsysrefscalar = tsysv * exp(tau/elev);
     745          }
     746        }
     747        oss << ", corrected (for El) tsys= "<<tsysrefscalar;
     748        pushLog(String(oss));
     749      }
     750      else {
     751        tsysrefscalar = tsysref[0];
     752      }
     753      //get quotient spectrum
     754      MaskedArray<Float> mref = maskedArray(specref, flagref);
     755      MaskedArray<Float> mon = maskedArray(specon, flagon);
     756      MaskedArray<Float> specres =   tsysrefscalar*((mon - mref)/mref);
     757      Double resint = onint*refint*smoothref/(onint+refint*smoothref);
     758
     759      //Debug
     760      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
     761      // fill the result, replay signal tsys by reference tsys
     762      outintCol.put(i, resint);
     763      outspecCol.put(i, specres.getArray());
     764      outflagCol.put(i, flagsFromMA(specres));
     765      outtsysCol.put(i, tsysref);
     766    }
     767    ++sit;
     768    ++s2it;
     769  }
     770  return out;
     771}
     772
     773CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
     774                                     const std::vector<int>& scans,
     775                                     int smoothref,
     776                                     casa::Float tsysv,
     777                                     casa::Float tau,
     778                                     casa::Float tcal )
     779
     780{
     781  setInsitu(false);
     782  STSelector sel;
     783  std::vector<int> scan1, scan2, beams;
     784  std::vector< vector<int> > scanpair;
     785  std::vector<string> calstate;
     786  String msg;
     787
     788  CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
     789  CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
     790
     791  std::vector< CountedPtr< Scantable > > sctables;
     792  sctables.push_back(s1b1on);
     793  sctables.push_back(s1b1off);
     794  sctables.push_back(s1b2on);
     795  sctables.push_back(s1b2off);
     796  sctables.push_back(s2b1on);
     797  sctables.push_back(s2b1off);
     798  sctables.push_back(s2b2on);
     799  sctables.push_back(s2b2off);
     800
     801  //check scanlist
     802  int n=s->checkScanInfo(scans);
     803  if (n==1) {
     804     throw(AipsError("Incorrect scan pairs. "));
     805  }
     806
     807  // Assume scans contain only a pair of consecutive scan numbers.
     808  // It is assumed that first beam, b1,  is on target.
     809  // There is no check if the first beam is on or not.
     810  if ( scans.size()==1 ) {
     811    scan1.push_back(scans[0]);
     812    scan2.push_back(scans[0]+1);
     813  } else if ( scans.size()==2 ) {
     814   scan1.push_back(scans[0]);
     815   scan2.push_back(scans[1]);
     816  } else {
     817    if ( scans.size()%2 == 0 ) {
     818      for (uInt i=0; i<scans.size(); i++) {
     819        if (i%2 == 0) {
     820          scan1.push_back(scans[i]);
     821        }
     822        else {
     823          scan2.push_back(scans[i]);
     824        }
     825      }
     826    } else {
     827      throw(AipsError("Odd numbers of scans, cannot form pairs."));
     828    }
     829  }
     830  scanpair.push_back(scan1);
     831  scanpair.push_back(scan2);
     832  calstate.push_back("*calon");
     833  calstate.push_back("*[^calon]");
     834  CountedPtr< Scantable > ws = getScantable(s, false);
     835  uInt l=0;
     836  while ( l < sctables.size() ) {
     837    for (uInt i=0; i < 2; i++) {
     838      for (uInt j=0; j < 2; j++) {
     839        for (uInt k=0; k < 2; k++) {
     840          sel.reset();
     841          sel.setScans(scanpair[i]);
     842          sel.setName(calstate[k]);
     843          beams.clear();
     844          beams.push_back(j);
     845          sel.setBeams(beams);
     846          ws->setSelection(sel);
     847          sctables[l]= getScantable(ws, false);
     848          l++;
     849        }
     850      }
     851    }
     852  }
     853
     854  // replace here by splitData or getData functionality
     855  CountedPtr< Scantable > sig1;
     856  CountedPtr< Scantable > ref1;
     857  CountedPtr< Scantable > sig2;
     858  CountedPtr< Scantable > ref2;
     859  CountedPtr< Scantable > calb1;
     860  CountedPtr< Scantable > calb2;
     861
     862  msg=String("Processing dototalpower for subset of the data");
     863  ostringstream oss1;
     864  oss1 << msg  << endl;
     865  pushLog(String(oss1));
     866  // Debug for IRC CS data
     867  //float tcal1=7.0;
     868  //float tcal2=4.0;
     869  sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
     870  ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
     871  ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
     872  sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
     873
     874  // correction of user-specified tsys for elevation here
     875
     876  // dosigref calibration
     877  msg=String("Processing dosigref for subset of the data");
     878  ostringstream oss2;
     879  oss2 << msg  << endl;
     880  pushLog(String(oss2));
     881  calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
     882  calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
     883
     884  // iteration by scanno or cycle no.
     885  Table& tcalb1 = calb1->table();
     886  Table& tcalb2 = calb2->table();
     887  TableIterator sit(tcalb1, "SCANNO");
     888  TableIterator s2it(tcalb2, "SCANNO");
     889  while ( !sit.pastEnd() ) {
     890    Table t1 = sit.table();
     891    Table t2= s2it.table();
     892    ArrayColumn<Float> outspecCol(t1, "SPECTRA");
     893    ArrayColumn<Float> outtsysCol(t1, "TSYS");
     894    ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
     895    ScalarColumn<Double> outintCol(t1, "INTERVAL");
     896    ArrayColumn<Float> t2specCol(t2, "SPECTRA");
     897    ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
     898    ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
     899    ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
     900    for (uInt i=0; i < t1.nrow(); ++i) {
     901      Vector<Float> spec1, spec2;
     902      // to store scalar (mean) tsys
     903      Vector<Float> tsys1, tsys2;
     904      Vector<uChar> flag1, flag2;
     905      Double tint1, tint2;
     906      outspecCol.get(i, spec1);
     907      t2specCol.get(i, spec2);
     908      outflagCol.get(i, flag1);
     909      t2flagCol.get(i, flag2);
     910      outtsysCol.get(i, tsys1);
     911      t2tsysCol.get(i, tsys2);
     912      outintCol.get(i, tint1);
     913      t2intCol.get(i, tint2);
     914      // average
     915      // assume scalar tsys for weights
     916      Float wt1, wt2, tsyssq1, tsyssq2;
     917      tsyssq1 = tsys1[0]*tsys1[0];
     918      tsyssq2 = tsys2[0]*tsys2[0];
     919      wt1 = Float(tint1)/tsyssq1;
     920      wt2 = Float(tint2)/tsyssq2;
     921      Float invsumwt=1/(wt1+wt2);
     922      MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
     923      MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
     924      MaskedArray<Float> avspec =  invsumwt * (wt1*mspec1 + wt2*mspec2);
     925      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
     926      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
     927      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
     928      Array<Float> avtsys =  tsys1;
     929
     930      outspecCol.put(i, avspec.getArray());
     931      outflagCol.put(i, flagsFromMA(avspec));
     932      outtsysCol.put(i, avtsys);
     933    }
     934    ++sit;
     935    ++s2it;
     936  }
     937  return calb1;
     938}
     939
     940//GBTIDL version of frequency switched data calibration
     941CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
     942                                      const std::vector<int>& scans,
     943                                      int smoothref,
     944                                      casa::Float tsysv,
     945                                      casa::Float tau,
     946                                      casa::Float tcal )
     947{
     948
     949
     950  STSelector sel;
     951  CountedPtr< Scantable > ws = getScantable(s, false);
     952  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
     953  CountedPtr< Scantable > calsig, calref, out;
     954
     955  //split the data
     956  sel.setName("*_fs");
     957  ws->setSelection(sel);
     958  sig = getScantable(ws,false);
     959  sel.reset();
     960  sel.setName("*_fs_calon");
     961  ws->setSelection(sel);
     962  sigwcal = getScantable(ws,false);
     963  sel.reset();
     964  sel.setName("*_fsr");
     965  ws->setSelection(sel);
     966  ref = getScantable(ws,false);
     967  sel.reset();
     968  sel.setName("*_fsr_calon");
     969  ws->setSelection(sel);
     970  refwcal = getScantable(ws,false);
     971
     972  calsig = dototalpower(sigwcal, sig, tcal=tcal);
     973  calref = dototalpower(refwcal, ref, tcal=tcal);
     974
     975  out=dosigref(calsig,calref,smoothref,tsysv,tau);
     976
     977  return out;
     978}
     979
    538980
    539981CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
    … …  
    12501692  vec = 0;
    12511693  pols->table_.rwKeywordSet().define("nPol", Int(1));
    1252   pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
     1694  //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
     1695  pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
    12531696  std::vector<CountedPtr<Scantable> > vpols;
    12541697  vpols.push_back(pols);

  • trunk/src/STMath.h

    r1374 r1391  
    132132                                        bool preserve = true );
    133133
     134  /**
     135    * Calibrate total power scans (translated from GBTIDL)
     136    * @param calon uncalibrated Scantable with CAL noise signal
     137    * @param caloff uncalibrated Scantable with no CAL signal
     138    * @param tcal optional scalar Tcal, CAL temperature (K)
     139    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
     140    * (spectrum - average of the two CAL on and off spectra;
     141    * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
     142    */
     143  casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
     144                                            const casa::CountedPtr<Scantable>& caloff,
     145                                            casa::Float tcal=1.0 );
     146
     147  /**
     148    * Combine signal and reference scans (translated from GBTIDL)
     149    * @param sig Scantable which contains signal scans
     150    * @param ref Scantable which contains reference scans
     151    * @param smoothref optional Boxcar smooth width of the reference scans
     152    * default: no smoothing (=1)
     153    * @param tsysv optional scalar Tsys value at the zenith, required to
     154    * set tau, as well
     155    * @param tau optional scalar Tau value
     156    * @return casa::CountedPtr<Scantable> which holds combined scans
     157    * (spectrum = (sig-ref)/ref * Tsys )
     158    */
     159  casa::CountedPtr<Scantable> dosigref( const casa::CountedPtr<Scantable>& sig,
     160                                        const casa::CountedPtr<Scantable>& ref,
     161                                        int smoothref=1,
     162                                        casa::Float tsysv=0.0,
     163                                        casa::Float tau=0.0 );
     164
     165 /**
     166    * Calibrate GBT Nod scan pairs (translated from GBTIDL)
     167    * @param s Scantable which contains Nod scans
     168    * @param scans Vector of scan numbers
     169    * @param smoothref optional Boxcar smooth width of the reference scans
     170    * @param tsysv optional scalar Tsys value at the zenith, required to
     171    * set tau, as well
     172    * @param tau optional scalar Tau value
     173    * @param tcal optional scalar Tcal, CAL temperature (K)
     174    * @return casa::CountedPtr<Scantable> which holds calibrated scans
     175    */
     176  casa::CountedPtr<Scantable> donod( const casa::CountedPtr<Scantable>& s,
     177                                     const std::vector<int>& scans,
     178                                     int smoothref=1,
     179                                     casa::Float tsysv=0.0,
     180                                     casa::Float tau=0.0,
     181                                     casa::Float tcal=0.0 );
     182
     183  /**
     184    * Calibrate frequency switched scans (translated from GBTIDL)
     185    * @param s Scantable which contains frequency switched  scans
     186    * @param scans Vector of scan numbers
     187    * @param smoothref optional Boxcar smooth width of the reference scans
     188    * @param tsysv optional scalar Tsys value at the zenith, required to
     189    * set tau, as well
     190    * @param tau optional scalar Tau value
     191    * @param tcal optional scalar Tcal, CAL temperature (K)
     192    * @return casa::CountedPtr<Scantable> which holds calibrated scans
     193    */
     194  casa::CountedPtr<Scantable> dofs( const casa::CountedPtr<Scantable>& s,
     195                                    const std::vector<int>& scans,
     196                                    int smoothref=1,
     197                                    casa::Float tsysv=0.0,
     198                                    casa::Float tau=0.0,
     199                                    casa::Float tcal=0.0 );
     200
     201
    134202  casa::CountedPtr<Scantable>
    135203    freqSwitch( const casa::CountedPtr<Scantable>& in );

  • trunk/src/STMathWrapper.h

    r1353 r1391  
    9191                                               preserve ) ); }
    9292
     93  ScantableWrapper dototalpower( const ScantableWrapper& calon,
     94                             const ScantableWrapper& caloff, casa::Float tcal= 0 )
     95  { return ScantableWrapper( STMath::dototalpower( calon.getCP(), caloff.getCP(), tcal ) ); }
     96
     97  ScantableWrapper dosigref( const ScantableWrapper& sig,
     98                             const ScantableWrapper& ref,
     99                             int smoothref = 0, casa::Float tsysv=0.0, casa::Float tau=0.0)
     100  { return ScantableWrapper( STMath::dosigref( sig.getCP(), ref.getCP(), smoothref, tsysv, tau ) ); }
     101
     102  ScantableWrapper donod( const ScantableWrapper& s,
     103                          const std::vector<int>& scans,
     104                          int smoothref = 0,
     105                          casa::Float tsysv=0.0, casa::Float tau=0.0, casa::Float tcal=0.0 )
     106  { return ScantableWrapper( STMath::donod( s.getCP(), scans, smoothref, tsysv, tau, tcal ) ); }
     107
     108  ScantableWrapper dofs( const ScantableWrapper& s,
     109                         const std::vector<int>& scans,
     110                         int smoothref = 0,
     111                         casa::Float tsysv=0.0, casa::Float tau=0.0, casa::Float tcal=0.0 )
     112  { return ScantableWrapper( STMath::dofs( s.getCP(), scans, smoothref, tsysv, tau, tcal ) ); }
     113
    93114  ScantableWrapper
    94115    freqSwitch( const ScantableWrapper& in )

  • trunk/src/STSelector.cpp

    r1004 r1391  
    8383void asap::STSelector::setName( const std::string& sname )
    8484{
    85   std::string sql = "SELECT from $1 WHERE SRCNAME == pattern('"+sname+"')";
     85  std::string sql = "SELECT FROM $1 WHERE SRCNAME == pattern('"+sname+"')";
    8686  setTaQL(sql);
    8787}

  • trunk/src/STTcal.cpp

    r856 r1391  
    6868}
    6969
     70/*** rewrite this for handling of GBT data
    7071uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
    7172{
    … …  
    9091  return resultid;
    9192}
     93***/
     94
     95uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
     96{
     97  // test if this already exists
     98  // TT - different Tcal values for each polarization, feed, and
     99  // data description. So there may be multiple entries for the same
     100  // time stamp.
     101  uInt resultid;
     102  uInt rno = table_.nrow();
     103  //table_.addRow();
     104  // get last assigned tcal_id and increment
     105  if ( rno == 0 ) {
     106    resultid = 0;
     107  }
     108  else {
     109    idCol_.get(rno-1, resultid);
     110    resultid++;
     111  }
     112  table_.addRow();
     113  tcalCol_.put(rno, cal);
     114  timeCol_.put(rno, time);
     115  idCol_.put(rno, resultid);
     116  return resultid;
     117}
     118
    92119
    93120void STTcal::getEntry( String& time, Vector<Float>& tcal, uInt id )

  • trunk/src/STWriter.cpp

    r1390 r1391  
    116116
    117117  // Extract the header from the table.
     118  // this is a little different from what I have done
     119  // before. Need to check with the Offline User Test data
    118120  STHeader hdr = in->getHeader();
    119121  //const Int nPol  = hdr.npol;
    … …  
    123125  Vector<uInt> nPol(nIF),nChan(nIF);
    124126  Vector<Bool> havexpol(nIF);
     127  String fluxUnit = hdr.fluxunit;
     128
    125129  nPol = 0;nChan = 0; havexpol = False;
    126130  for (uint i=0;i<ifs.size();++i) {
    … …  
    264268  pushLog(String(oss));
    265269  writer_->close();
    266 
     270  //if MS2 delete POINTING table exists and copy the one in the keyword
     271  if ( format_ == "MS2" ) {
     272    replacePtTab(table, filename);
     273  }
    267274  return 0;
    268275}
    … …  
    314321}
    315322
    316 
    317 }
     323// For writing MS data, if there is the reference to
     324// original pointing table it replace it by it.
     325void STWriter::replacePtTab (const Table& tab, const std::string& fname)
     326{
     327  String oldPtTabName = fname;
     328  oldPtTabName.append("/POINTING");
     329  if ( tab.keywordSet().isDefined("POINTING") ) {
     330    String PointingTab = tab.keywordSet().asString("POINTING");
     331    if ( Table::isReadable(PointingTab) ) {
     332      Table newPtTab(PointingTab, Table::Old);
     333      newPtTab.copy(oldPtTabName, Table::New);
     334      ostringstream oss;
     335      oss << "STWriter: copied  " <<PointingTab  << " to " << fname;
     336      pushLog(String(oss));
     337    }
     338  }
     339}
     340
     341}

  • trunk/src/STWriter.h

    r1353 r1391  
    8181                      casa::Vector<casa::Complex>& xpol,
    8282                      const casa::Table& tab);
     83
     84  void replacePtTab(const casa::Table& tab, const std::string& fname);
     85
    8386  std::string     format_;
    8487  PKSwriter* writer_;

  • trunk/src/Scantable.cpp

    r1375 r1391  
    10281028}
    10291029
     1030std::string asap::Scantable::getAntennaName() const
     1031{
     1032  String out;
     1033  table_.keywordSet().get("AntennaName", out);
     1034  return out;
     1035}
     1036
     1037int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
     1038{
     1039  String tbpath;
     1040  int ret = 0;
     1041  if ( table_.keywordSet().isDefined("GBT_GO") ) {
     1042    table_.keywordSet().get("GBT_GO", tbpath);
     1043    Table t(tbpath,Table::Old);
     1044    // check each scan if other scan of the pair exist
     1045    int nscan = scanlist.size();
     1046    for (int i = 0; i < nscan; i++) {
     1047      Table subt = t( t.col("SCAN") == scanlist[i]+1 );
     1048      if (subt.nrow()==0) {
     1049        cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
     1050        ret = 1;
     1051        break;
     1052      }
     1053      ROTableRow row(subt);
     1054      const TableRecord& rec = row.get(0);
     1055      int scan1seqn = rec.asuInt("PROCSEQN");
     1056      int laston1 = rec.asuInt("LASTON");
     1057      if ( rec.asuInt("PROCSIZE")==2 ) {
     1058        if ( i < nscan-1 ) {
     1059          Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
     1060          if ( subt2.nrow() == 0) {
     1061            cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
     1062            ret = 1;
     1063            break;
     1064          }
     1065          ROTableRow row2(subt2);
     1066          const TableRecord& rec2 = row2.get(0);
     1067          int scan2seqn = rec2.asuInt("PROCSEQN");
     1068          int laston2 = rec2.asuInt("LASTON");
     1069          if (scan1seqn == 1 && scan2seqn == 2) {
     1070            if (laston1 == laston2) {
     1071              cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1072              i +=1;
     1073            }
     1074            else {
     1075              cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1076            }
     1077          }
     1078          else if (scan1seqn==2 && scan2seqn == 1) {
     1079            if (laston1 == laston2) {
     1080              cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
     1081              ret = 1;
     1082              break;
     1083            }
     1084          }
     1085          else {
     1086            cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
     1087            ret = 1;
     1088            break;
     1089          }
     1090        }
     1091      }
     1092      else {
     1093        cerr<<"The scan does not appear to be standard obsevation."<<endl;
     1094      }
     1095    //if ( i >= nscan ) break;
     1096    }
     1097  }
     1098  else {
     1099    cerr<<"No reference to GBT_GO table."<<endl;
     1100    ret = 1;
     1101  }
     1102  return ret;
     1103}
     1104
     1105std::vector<double>  asap::Scantable::getDirectionVector(int whichrow) const
     1106{
     1107  Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
     1108  std::vector<double> dir;
     1109  Dir.tovector(dir);
     1110  return dir;
     1111}
     1112
    10301113}
    10311114 //namespace asap

  • trunk/src/Scantable.h

    r1385 r1391  
    381381    { STFitEntry fe; fitTable_.getEntry(fe, mfitidCol_(row)); return fe; }
    382382
     383  //Added by TT
     384  /**
     385   * Get the antenna name
     386   * @return antenna name string
     387   */
     388  std::string getAntennaName() const;
     389
     390  /**
     391   * For GBT MS data only. check a scan list
     392   * against the information found in GBT_GO table for
     393   * scan number orders to get correct pairs.
     394   * @param[in] scan list
     395   * @return status
     396   */
     397  int checkScanInfo(const std::vector<int>& scanlist) const;
     398
     399  /**
     400   * Get the direction as a vector, for a specific row
     401   * @param[in] whichrow the row numbyyer
     402   * @return the direction in a vector
     403   */
     404  std::vector<double> getDirectionVector(int whichrow) const;
     405
    383406private:
    384407

  • trunk/src/ScantableWrapper.h

    r1385 r1391  
    199199    { return table_->columnNames(); }
    200200
     201  std::string getAntennaName() const
     202    { return table_->getAntennaName(); }
     203
     204  int checkScanInfo(const vector<int>& scanlist) const
     205    { return table_->checkScanInfo(scanlist); }
     206
     207  std::vector<double>  getDirectionVector(int whichrow) const
     208    { return table_->getDirectionVector(whichrow); }
     209
    201210private:
    202211  casa::CountedPtr<Scantable> table_;

  • trunk/src/python_Fitter.cpp

    r894 r1391  
    5555        .def("reset", &Fitter::reset)
    5656        .def("fit", &Fitter::fit)
     57        .def("lfit", &Fitter::lfit)
    5758        .def("evaluate", &Fitter::evaluate)
    5859      ;

  • trunk/src/python_STMath.cpp

    r1308 r1391  
    5252        .def("_auto_quotient", &STMathWrapper::autoQuotient)
    5353        .def("_quotient", &STMathWrapper::quotient)
     54        .def("_dototalpower", &STMathWrapper::dototalpower)
     55        .def("_dosigref", &STMathWrapper::dosigref)
     56        .def("_donod", &STMathWrapper::donod)
     57        .def("_dofs", &STMathWrapper::dofs)
    5458        .def("_stats", &STMathWrapper::statistic)
    5559        .def("_freqswitch", &STMathWrapper::freqSwitch)

  • trunk/src/python_Scantable.cpp

    r1360 r1391  
    102102    .def("_getdirection", &ScantableWrapper::getDirectionString,
    103103         (boost::python::arg("whichrow")=0) )
     104    .def("get_antennaname", &ScantableWrapper::getAntennaName)
    104105    .def("_flag", &ScantableWrapper::flag)
    105106    .def("_save",  &ScantableWrapper::makePersistent)
    … …  
    121122    .def("_recalcazel", &ScantableWrapper::calculateAZEL)
    122123    .def("_setsourcetype", &ScantableWrapper::setSourceType)
     124    .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
    123125  ;
    124126};
Note: See TracChangeset for help on using the changeset viewer.