sdfit.py

Go to the documentation of this file.

#
# This file was generated using xslt from its XML file
#
# Copyright 2009, Associated Universities Inc., Washington DC
#
import sys
import os
from  casac import *
import string
from taskinit import casalog
#from taskmanager import tm
import task_sdfit
def sdfit(infile='', antenna=0, fluxunit='', telescopeparm='', specunit='', restfreq='', frame='', doppler='', scanlist=[], field='', iflist=[], pollist=[], fitfunc='gauss', fitmode='auto', maskline=[], invertmask=False, nfit=[], thresh=5.0, min_nchan=3, avg_limit=4, box_size=0.2, edge=[0], outfile='', overwrite=False, plotlevel=0):

        """ASAP SD task: fit a spectral line
        Keyword arguments:
        infile -- name of input SD dataset
                default: none - must input file name
                example: 'mysd.asap'
                         See sdcal for allowed formats.
        antenna -- antenna name or id (only effective for MS input). 
        fluxunit -- units for line flux
                options: (str) 'K','Jy',''
                default: '' (keep current fluxunit)
                WARNING: For GBT data, see description below.
            >>> fluxunit expandable parameter
                 telescopeparm -- the telescope characteristics
                        options: (str) name or (list) list of gain info
                        default: '' (none set)
                        example: if telescopeparm='', it tries to get the telescope
                                 name from the data.
                                 Full antenna parameters (diameter,ap.eff.) known
                                 to ASAP are
                                 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA', 'DSS-43',
                                 'CEDUNA','HOBART'. For GBT, it fixes default fluxunit
                                 to 'K' first then convert to a new fluxunit.
                                 telescopeparm=[104.9,0.43] diameter(m), ap.eff.
                                 telescopeparm=[0.743] gain in Jy/K
                                 telescopeparm='FIX' to change default fluxunit
                                 see description below

        specunit -- units for spectral axis
                options: (str) 'channel','km/s','GHz','MHz','kHz','Hz',''
                default: '' (=current)
                example: this will be the units for maskline
            >>> specunit expandable parameters
                 restfreq -- rest frequency
                         available type includes float, int, string, list of float, 
                         list of int, list of string, and list of dictionary. the 
                         default unit of restfreq in case of float, int, or string 
                         without unit is Hz. string input can be a value only 
                         (treated as Hz) or a value followed by unit for which 'GHz',
                         'MHz','kHz',and 'Hz' are available. 
                         a list can be used to set different rest frequencies for 
                         each IF. the length of list input must be nIF. dictionary 
                         input should be a pair of molecule name and frequency with 
                         keys of 'name' and 'value', respectively. values in the 
                         dictionary input follows the same manner as for single 
                         float or string input. 
                         example: 345.796
                                  '1420MHz'
                                  [345.8, 347.0, 356.7]
                                  ['345.8MHz', '347.0MHz', '356.7MHz']
                                  [{'name':'CO','value':345}]
        frame -- frequency frame for spectral axis
                options: (str) 'LSRK','REST','TOPO','LSRD','BARY',
                         'GEO','GALACTO','LGROUP','CMB'
                default: currently set frame in scantable
                WARNING: frame='REST' not yet implemented
        doppler -- doppler mode
                options: (str) 'RADIO','OPTICAL','Z','BETA','GAMMA'
                default: currently set doppler in scantable
        scanlist -- list of scan numbers to process
                default: [] (use all scans)
                example: [21,22,23,24]
        field -- selection string for selecting scans by name
                default: '' (no name selection)
                example: 'FLS3a*'
                this selection is in addition to scanlist and iflist
        iflist -- list of IF id numbers to select
                default: [] (use all IFs)
                example: [15]
        pollist -- list of polarization id numbers to select
                default: [] (use all polarizations)
                example: [1]
        fitfunc -- function for fitting
                options: (str) 'gauss','lorentz'
                default: 'gauss'
        fitmode -- mode for fitting
                options: (str) 'list','auto','interact'
                default: 'auto'
                example: 'list' will use maskline to define regions to
                                fit for lines with nfit in each
                         'auto' will use the linefinder to fit for lines
                                using the following parameters
                         'interact' allows adding and deleting mask 
                                regions by drawing rectangles on the plot 
                                with mouse. Draw a rectangle with LEFT-mouse 
                                to ADD the region to the mask and with RIGHT-mouse 
                                to DELETE the region. 

            >>> fitmode expandable parameters             
                 thresh -- S/N threshold for linefinder
                         default: 5
                         example: a single channel S/N ratio above which the channel is
                                  considered to be a detection
                 min_nchan -- minimum number of consecutive channels for linefinder
                         default: 3
                         example: minimum number of consecutive channels required to pass threshold
                 avg_limit -- channel averaging for broad lines
                         default: 4
                         example: a number of consecutive channels not greater than
                                  this parameter can be averaged to search for broad lines
                 box_size -- running mean box size
                         default: 0.2
                         example: a running mean box size specified as a fraction
                                  of the total spectrum length
                 edge -- channels to drop at beginning and end of spectrum
                         default: 0
                         example: [1000] drops 1000 channels at beginning AND end
                                  [1000,500] drops 1000 from beginning and 500 from end

                 Note: For bad baselines threshold should be increased,
                 and avg_limit decreased (or even switched off completely by
                 setting this parameter to 1) to avoid detecting baseline
                 undulations instead of real lines.

        maskline -- list of mask regions to INCLUDE in LINE fitting
                default: all
                example: maskline=[[3900,4300]] for a single region, or
                         maskline=[[3900,4300],[5000,5400]] for two, etc.
        invertmask -- invert mask (EXCLUDE masklist instead)
                options: (bool) True, False
                default: False
                example: invertmask=True, then will make one region that is
                         the exclusion of the maskline regions
        nfit -- list of number of gaussian/lorentzian lines to fit in in maskline region (ignored when fitmode='auto')
                default: 0 (no fitting)
                example: nfit=[1] for single line in single region,
                         nfit=[2] for two lines in single region,
                         nfit=[1,1] for single lines in each of two regions, etc.
        outfile -- name of output file for fit results
                default: no output fit file
                example: 'mysd.fit'
        overwrite -- overwrite the outfile if already exists
                options: (bool) True, False
                default: False
        plotlevel -- control for plotting of results
                options: (int) 0=none, 1=some, 2=more
                default: 0 (no plotting)
                example: plotlevel=1 plots fit
                         plotlevel=2 plots fit and residual 
                         no hardcopy available for fitter
                WARNING: be careful plotting OTF data with lots of fields

        -------------------------------------------------------------------
        Returns a Python dictionary of line statistics
                keys:    'peak','cent','fwhm','nfit'
                example: each value is a list of lists with one list of
                         2 entries [fitvalue,error] per component.
                         e.g. xstat['peak']=[[234.9, 4.8],[234.2, 5.3]]
                         for 2 components.


        DESCRIPTION:

        Task sdfit is a basic line-fitter for single-dish spectra.
        It assumes that the spectra have been calibrated in sdcal
        or sdreduce.

        Furthermore, it assumes that any selection of scans, IFs,
        polarizations, and time and channel averaging/smoothing has
        also already been done (in other sd tasks) as there are no controls
        for these.  Note that you can use sdsave to do selection, writing
        out a new scantable.

        Note that multiple scans, IFs, and polarizations can in principle 
        be handled, but we recommend that you use scanlist, field, iflist, 
        and pollist to give a single selection for each fit.

        Currently, you can choose Gaussian or Lorentzian profile as a 
        fitting model. 

        For complicated spectra, sdfit does not do a good job of
        "auto-guessing" the starting model for the fit.  We recommend
        you use sd.fitter in the toolkit which has more options, such
        as fixing components in the fit and supplying starting guesses
        by hand.

        WARNING: sdfit will currently return the fit for the first
        row in the scantable.

        ASAP recognizes the data of the "AT" telescopes, but currently
        does not know about the GBT or any other telescope. This task
        does know about GBT. Telescope name is obtained from the data.
        If you wish to change the fluxunit (see below), and telescopeparm='',
        for the AT telescopes it will use internal telescope parameters for
        flux conversion. For GBT, it will use an approximate aperture
        efficiency conversion.  If you give telescopeparm a list, then 
        if the list has a single float it is assumed to be the gain in Jy/K, 
        if two or more elements they are assumed to be telescope diameter (m) 
        and aperture efficiency respectively.

        WARNING for the GBT raw SDFITS format data as input:
        SDtasks are able to handle GBT raw SDFITS format data since the 
        data filler is available. However, the functionality is not well 
        tested yet, so that there may be unknown bugs.  

  
        """
        if type(scanlist)==int: scanlist=[scanlist]
        if type(iflist)==int: iflist=[iflist]
        if type(pollist)==int: pollist=[pollist]
        if type(nfit)==int: nfit=[nfit]
        if type(edge)==int: edge=[edge]

#
#    The following is work around to avoid a bug with current python translation
#
        mytmp = {}

        mytmp['infile'] = infile
        mytmp['antenna'] = antenna
        mytmp['fluxunit'] = fluxunit
        mytmp['telescopeparm'] = telescopeparm
        mytmp['specunit'] = specunit
        mytmp['restfreq'] = restfreq
        mytmp['frame'] = frame
        mytmp['doppler'] = doppler
        mytmp['scanlist'] = scanlist
        mytmp['field'] = field
        mytmp['iflist'] = iflist
        mytmp['pollist'] = pollist
        mytmp['fitfunc'] = fitfunc
        mytmp['fitmode'] = fitmode
        mytmp['maskline'] = maskline
        mytmp['invertmask'] = invertmask
        mytmp['nfit'] = nfit
        mytmp['thresh'] = thresh
        mytmp['min_nchan'] = min_nchan
        mytmp['avg_limit'] = avg_limit
        mytmp['box_size'] = box_size
        mytmp['edge'] = edge
        mytmp['outfile'] = outfile
        mytmp['overwrite'] = overwrite
        mytmp['plotlevel'] = plotlevel
        pathname='file:///'+os.environ.get('CASAPATH').split()[0]+'/share/xml/'
        trec = casac.utils().torecord(pathname+'sdfit.xml')

        casalog.origin('sdfit')
        if trec.has_key('sdfit') and casac.utils().verify(mytmp, trec['sdfit']) :
            result = task_sdfit.sdfit(infile, antenna, fluxunit, telescopeparm, specunit, restfreq, frame, doppler, scanlist, field, iflist, pollist, fitfunc, fitmode, maskline, invertmask, nfit, thresh, min_nchan, avg_limit, box_size, edge, outfile, overwrite, plotlevel)

        else :
          result = False
        return result