Changeset 1391 for trunk


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
Files:
23 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/asapfitter.py

    r1306 r1391  
    2525        self._p = None
    2626        self._selection = None
     27        self.uselinear = False
    2728
    2829    def set_data(self, xdat, ydat, mask=None):
     
    7475        Set the function to be fit.
    7576        Parameters:
    76             poly:    use a polynomial of the order given
     77            poly:    use a polynomial of the order given with nonlinear least squares fit
     78            lpoly:   use polynomial of the order given with linear least squares fit
    7779            gauss:   fit the number of gaussian specified
    7880        Example:
    7981            fitter.set_function(gauss=2) # will fit two gaussians
    80             fitter.set_function(poly=3)  # will fit a 3rd order polynomial
     82            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
     83            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
    8184        """
    8285        #default poly order 0
     
    8689            n = kwargs.get('poly')
    8790            self.components = [n]
     91            self.uselinear = False
     92        elif kwargs.has_key('lpoly'):
     93            self.fitfunc = 'poly'
     94            n = kwargs.get('lpoly')
     95            self.components = [n]
     96            self.uselinear = True
    8897        elif kwargs.has_key('gauss'):
    8998            n = kwargs.get('gauss')
     
    91100            self.fitfuncs = [ 'gauss' for i in range(n) ]
    92101            self.components = [ 3 for i in range(n) ]
     102            self.uselinear = False
    93103        else:
    94104            msg = "Invalid function type."
     
    146156            if len(fxdpar) and fxdpar.count(0) == 0:
    147157                 raise RuntimeError,"No point fitting, if all parameters are fixed."
    148             converged = self.fitter.fit()
     158            if self.uselinear:
     159                converged = self.fitter.lfit()
     160            else:
     161                converged = self.fitter.fit()
    149162            if not converged:
    150163                raise RuntimeError,"Fit didn't converge."
     
    501514                             array(self.data._getmask(self._fittedrow),
    502515                                   copy=False))
    503            
     516                             
    504517            ylab = self.data._get_ordinate_label()
    505518
  • trunk/python/asapmath.py

    r1362 r1391  
    22from asap import rcParams
    33from asap import print_log
     4from asap import selector
    45
    56def average_time(*args, **kwargs):
     
    113114    return s
    114115
     116def dototalpower(calon, caloff, tcalval=0.0):
     117    """
     118    Do calibration for CAL on,off signals.
     119    Adopted from GBTIDL dototalpower
     120    Parameters:
     121        calon:         the 'cal on' subintegration
     122        caloff:        the 'cal off' subintegration
     123        tcalval:       user supplied Tcal value
     124    """
     125    varlist = vars()
     126    from asap._asap import stmath
     127    stm = stmath()
     128    stm._setinsitu(False)
     129    s = scantable(stm._dototalpower(calon, caloff, tcalval))
     130    s._add_history("dototalpower",varlist)
     131    print_log()
     132    return s
     133
     134def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
     135    """
     136    Calculate a quotient (sig-ref/ref * Tsys)
     137    Adopted from GBTIDL dosigref
     138    Parameters:
     139        sig:         on source data
     140        ref:         reference data
     141        smooth:      width of box car smoothing for reference
     142        tsysval:     user specified Tsys (scalar only)
     143        tauval:      user specified Tau (required if tsysval is set)
     144    """
     145    varlist = vars()
     146    from asap._asap import stmath
     147    stm = stmath()
     148    stm._setinsitu(False)
     149    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
     150    s._add_history("dosigref",varlist)
     151    print_log()
     152    return s
     153
     154def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     155    """
     156    Calibrate GBT position switched data
     157    Adopted from GBTIDL getps
     158    Currently calps identify the scans as position switched data if they
     159    contain '_ps' in the source name. The data must contains 'CAL' signal
     160    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
     161    need to be present in the source name field.
     162    (GBT MS data reading process to scantable automatically append these
     163    id names to the source names)
     164
     165    Parameters:
     166        scantab:       scantable
     167        scannos:       list of scan numbers
     168        smooth:        optional box smoothing order for the reference
     169                       (default is 1 = no smoothing)
     170        tsysval:       optional user specified Tsys (default is 0.0,
     171                       use Tsys in the data)
     172        tauval:        optional user specified Tau
     173        tcalval:       optional user specified Tcal (default is 0.0,
     174                       use Tcal value in the data)
     175    """
     176    varlist = vars()
     177    # check for the appropriate data
     178    s = scantab.get_scan('*_ps*')
     179    if s is None:
     180        msg = "The input data appear to contain no position-switch mode data."
     181        if rcParams['verbose']:
     182            print msg
     183            return
     184        else:
     185            raise TypeError(msg)
     186    ssub = s.get_scan(scannos)
     187    if ssub is None:
     188        msg = "No data was found with given scan numbers!"
     189        if rcParams['verbose']:
     190            print msg
     191            return
     192        else:
     193            raise TypeError(msg)
     194    ssubon = ssub.get_scan('*calon')
     195    ssuboff = ssub.get_scan('*[^calon]')
     196    if ssubon.nrow() != ssuboff.nrow():
     197        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
     198        if rcParams['verbose']:
     199            print msg
     200            return
     201        else:
     202            raise TypeError(msg)
     203    cals = dototalpower(ssubon, ssuboff, tcalval)
     204    sig = cals.get_scan('*ps')
     205    ref = cals.get_scan('*psr')
     206    if sig.nscan() != ref.nscan():
     207        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
     208        if rcParams['verbose']:
     209            print msg
     210            return
     211        else:
     212            raise TypeError(msg)
     213
     214    #for user supplied Tsys
     215    if tsysval>0.0:
     216        if tauval<=0.0:
     217            msg = "Need to supply a valid tau to use the supplied Tsys"
     218            if rcParams['verbose']:
     219                print msg
     220                return
     221            else:
     222                raise TypeError(msg)
     223        else:
     224            sig.recalc_azel()
     225            ref.recalc_azel()
     226            #msg = "Use of user specified Tsys is not fully implemented yet."
     227            #if rcParams['verbose']:
     228            #    print msg
     229            #    return
     230            #else:
     231            #    raise TypeError(msg)
     232            # use get_elevation to get elevation and
     233            # calculate a scaling factor using the formula
     234            # -> tsys use to dosigref
     235
     236    #ress = dosigref(sig, ref, smooth, tsysval)
     237    ress = dosigref(sig, ref, smooth, tsysval, tauval)
     238    ress._add_history("calps", varlist)
     239    print_log()
     240    return ress
     241
     242def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     243    """
     244    Do full (but a pair of scans at time) processing of GBT Nod data
     245    calibration.
     246    Adopted from  GBTIDL's getnod
     247    Parameters:
     248        scantab:     scantable
     249        scannos:     a pair of scan numbers, or the first scan number of the pair
     250        smooth:      box car smoothing order
     251        tsysval:     optional user specified Tsys value
     252        tauval:      optional user specified tau value (not implemented yet)
     253        tcalval:     optional user specified Tcal value
     254    """
     255    varlist = vars()
     256    from asap._asap import stmath
     257    stm = stmath()
     258    stm._setinsitu(False)
     259
     260    # check for the appropriate data
     261    s = scantab.get_scan('*_nod*')
     262    if s is None:
     263        msg = "The input data appear to contain no Nod observing mode data."
     264        if rcParams['verbose']:
     265            print msg
     266            return
     267        else:
     268            raise TypeError(msg)
     269
     270    # need check correspondance of each beam with sig-ref ...
     271    # check for timestamps, scan numbers, subscan id (not available in
     272    # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
     273    # and beam 1 - ref...)
     274    # First scan number of paired scans or list of pairs of
     275    # scan numbers (has to have even number of pairs.)
     276
     277    #data splitting
     278    scan1no = scan2no = 0
     279
     280    if len(scannos)==1:
     281        scan1no = scannos[0]
     282        scan2no = scannos[0]+1
     283        pairScans = [scan1no, scan2no]
     284    else:
     285        #if len(scannos)>2:
     286        #    msg = "calnod can only process a pair of nod scans at time."
     287        #    if rcParams['verbose']:
     288        #        print msg
     289        #        return
     290        #    else:
     291        #        raise TypeError(msg)
     292        #
     293        #if len(scannos)==2:
     294        #    scan1no = scannos[0]
     295        #    scan2no = scannos[1]
     296        pairScans = list(scannos)
     297
     298    if tsysval>0.0:
     299        if tauval<=0.0:
     300            msg = "Need to supply a valid tau to use the supplied Tsys"
     301            if rcParams['verbose']:
     302                print msg
     303                return
     304            else:
     305                raise TypeError(msg)
     306        else:
     307            scantab.recalc_azel()
     308    resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
     309    resspec._add_history("calnod",varlist)
     310    print_log()
     311    return resspec
     312
     313def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     314    """
     315    Calibrate GBT frequency switched data.
     316    Adopted from GBTIDL getfs.
     317    Currently calfs identify the scans as frequency switched data if they
     318    contain '_fs' in the source name. The data must contains 'CAL' signal
     319    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
     320    need to be present in the source name field.
     321    (GBT MS data reading via scantable automatically append these
     322    id names to the source names)
     323
     324    Parameters:
     325        scantab:       scantable
     326        scannos:       list of scan numbers
     327        smooth:        optional box smoothing order for the reference
     328                       (default is 1 = no smoothing)
     329        tsysval:       optional user specified Tsys (default is 0.0,
     330                       use Tsys in the data)
     331        tauval:        optional user specified Tau
     332    """
     333    varlist = vars()
     334    from asap._asap import stmath
     335    stm = stmath()
     336    stm._setinsitu(False)
     337
     338#    check = scantab.get_scan('*_fs*')
     339#    if check is None:
     340#        msg = "The input data appear to contain no Nod observing mode data."
     341#        if rcParams['verbose']:
     342#            print msg
     343#            return
     344#        else:
     345#            raise TypeError(msg)
     346    s = scantab.get_scan(scannos)
     347    del scantab
     348
     349    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
     350    resspec._add_history("calfs",varlist)
     351    print_log()
     352    return resspec
     353
    115354def simple_math(left, right, op='add', tsys=True):
    116355    """
  • trunk/python/asapplotter.py

    r1358 r1391  
    741741        return userlabel or d[mode]
    742742
     743    def plotazel(self, scan=None):
     744        """
     745        plot azimuth and elevation  versus time of a scantable
     746        """
     747        import pylab as PL
     748        from matplotlib.dates import DateFormatter, timezone, HourLocator, MinuteLocator, DayLocator
     749        from matplotlib.ticker import MultipleLocator
     750        from matplotlib.numerix import array, pi
     751        self._data = scan
     752        dates = self._data.get_time()
     753        t = PL.date2num(dates)
     754        tz = timezone('UTC')
     755        PL.cla()
     756        PL.ioff()
     757        PL.clf()
     758        tdel = max(t) - min(t)
     759        ax = PL.subplot(2,1,1)
     760        el = array(self._data.get_elevation())*180./pi
     761        PL.ylabel('El [deg.]')
     762        dstr = dates[0].strftime('%Y/%m/%d')
     763        if tdel > 1.0:
     764            dstr2 = dates[len(dates)-1].strftime('%Y/%m/%d')
     765            dstr = dstr + " - " + dstr2
     766            majloc = DayLocator()
     767            minloc = HourLocator(range(0,23,12))
     768            timefmt = DateFormatter("%b%d")
     769        else:
     770            timefmt = DateFormatter('%H')
     771            majloc = HourLocator()
     772            minloc = MinuteLocator(20)
     773        PL.title(dstr)
     774        PL.plot_date(t,el,'b,', tz=tz)
     775        #ax.grid(True)
     776        ax.yaxis.grid(True)
     777        yloc = MultipleLocator(30)
     778        ax.set_ylim(0,90)
     779        ax.xaxis.set_major_formatter(timefmt)
     780        ax.xaxis.set_major_locator(majloc)
     781        ax.xaxis.set_minor_locator(minloc)
     782        ax.yaxis.set_major_locator(yloc)
     783        if tdel > 1.0:
     784            labels = ax.get_xticklabels()
     785        #    PL.setp(labels, fontsize=10, rotation=45)
     786            PL.setp(labels, fontsize=10)
     787        # Az plot
     788        az = array(self._data.get_azimuth())*180./pi
     789        if min(az) < 0:
     790            for irow in range(len(az)):
     791                if az[irow] < 0: az[irow] += 360.0
     792
     793        ax = PL.subplot(2,1,2)
     794        PL.xlabel('Time (UT)')
     795        PL.ylabel('Az [deg.]')
     796        PL.plot_date(t,az,'b,', tz=tz)
     797        ax.set_ylim(0,360)
     798        #ax.grid(True)
     799        ax.yaxis.grid(True)
     800        #hfmt = DateFormatter('%H')
     801        #hloc = HourLocator()
     802        yloc = MultipleLocator(60)
     803        ax.xaxis.set_major_formatter(timefmt)
     804        ax.xaxis.set_major_locator(majloc)
     805        ax.xaxis.set_minor_locator(minloc)
     806        ax.yaxis.set_major_locator(yloc)
     807        if tdel > 1.0:
     808            labels = ax.get_xticklabels()
     809            PL.setp(labels, fontsize=10)
     810        PL.ion()
     811        PL.draw()
     812
     813    def plotpointing(self, scan=None):
     814        """
     815        plot telescope pointings
     816        """
     817        import pylab as PL
     818        from matplotlib.dates import DateFormatter, timezone
     819        from matplotlib.ticker import MultipleLocator
     820        from matplotlib.numerix import array, pi, zeros
     821        self._data = scan
     822        dir = array(self._data.get_directionval()).transpose()
     823        ra = dir[0]*180./pi
     824        dec = dir[1]*180./pi
     825        PL.cla()
     826        PL.ioff()
     827        PL.clf()
     828        ax = PL.axes([0.1,0.1,0.8,0.8])
     829        ax = PL.axes([0.1,0.1,0.8,0.8])
     830        ax.set_aspect('equal')
     831        PL.plot(ra,dec, 'b,')
     832        PL.xlabel('RA [deg.]')
     833        PL.ylabel('Declination [deg.]')
     834        PL.title('Telescope pointings')
     835        [xmin,xmax,ymin,ymax] = PL.axis()
     836        PL.axis([xmax,xmin,ymin,ymax])
     837        PL.ion()
     838        PL.draw()
     839
  • trunk/python/scantable.py

    r1373 r1391  
    516516        """
    517517        return self._get_column(self._getdirection, row)
     518
     519    def get_directionval(self, row=-1):
     520        """
     521        Get a list of Positions on the sky (direction) for the observations.
     522        Return a float for each integration in the scantable.
     523        Parameters:
     524            row:    row no of integration. Default -1 return all rows
     525        Example:
     526            none
     527        """
     528        return self._get_column(self._getdirectionvec, row)
     529
    518530
    519531    def set_unit(self, unit='channel'):
     
    12151227
    12161228
    1217     def poly_baseline(self, mask=None, order=0, plot=False, insitu=None):
     1229    def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
    12181230        """
    12191231        Return a scan which has been baselined (all rows) by a polynomial.
     
    12241236                        indivual fit has to be approved, by typing 'y'
    12251237                        or 'n'
     1238            uselin:     use linear polynomial fit
    12261239            insitu:     if False a new scantable is returned.
    12271240                        Otherwise, the scaling is done in-situ
     
    12401253            f = fitter()
    12411254            f.set_scan(self, mask)
    1242             f.set_function(poly=order)
     1255            #f.set_function(poly=order)
     1256            if uselin:
     1257                f.set_function(lpoly=order)
     1258            else:
     1259                f.set_function(poly=order)
    12431260            s = f.auto_fit(insitu, plot=plot)
    12441261            s._add_history("poly_baseline", varlist)
  • 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.