import sys import os import casac import string import inspect import pdb from parameter_check import * casapath=os.environ['AIPSPATH'] import asap as sd os.environ['AIPSPATH']=casapath def sdfit(sdfile=None,telescope=None,fluxunit=None,specunit=None,frame=None,doppler=None,scanlist=None,field=None,iflist=None,fitmode=None,maskline=None,invertmask=None,nfit=None,fitfile=None,plotlevel=None,thresh=None,min_nchan=None,avg_limit=None,box_size=None,edge=None,xstat=None): """ASAP SD line-fitting task (prototype): Task Version 2007-03-08 STM Keyword arguments: sdfile -- name of input SD dataset default: none - must input file name example: 'mysd.asap' See sdcal for allowed formats. telescope -- the telescope name or characteristics options: (str) name or (list) list of gain info default: '' (none set) example: telescope='GBT' for the GBT telescope='AT' for one of the default scopes known to ASAP. telescope=[104.9,0.43] diameter(m), ap.eff. telescope=[0.743] gain in Jy/K telescope='FIX' to change default fluxunit see description below fluxunit -- units for line flux options: (str) 'K','Jy','' default: '' (keep current fluxunit) WARNING: For GBT data, see description below. specunit -- units for spectral axis options: (str) 'channel','km/s','GHz','MHz','kHz','Hz','' default: '' (keep current specunit) 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] fitmode -- mode for fitting options: (str) 'list','auto' 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 autopars (see below) 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 lines to fit in in maskline region 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. fitfile -- name of output file for fit results default: no output fit file example: 'mysd.fit' plotlevel -- control for plotting of results options: (int) 0=none, 1=some, 2=more default: 0 (no plotting) example: plotlevel=1 plots fit and residual no hardcopy available for fitter WARNING: be careful plotting OTF data with lots of fields AUTOPARS: the following parameters are used for fitmode='auto' ONLY ------------------------------------------------------------------- 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 consequtive 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. ------------------------------------------------------------------- xstat -- RETURN ONLY: a Python dictionary of line statistics keys: 'peak','cent','fwhm' 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. Furthermore, it assumes that any selection of scans, IFs, polarizations, and time and channel averaging/smoothing has also already been done (in sdcal) as there are no controls for these. Note that you can run sdcal with calmode='none' and do selection, writing out a new scantable. Note that multiple scans and IFs can in principle be handled, but we recommend that you use scanlist, field, and iflist to give a single selection for each fit. 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. Does not handle multiple polarizations. 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. Therefore, if you wish to change the fluxunit (see below), then you need to tell it what to do. If you set telescope='AT' it will use its internal defaults. If you set telescope 'GBT', it will use an approximate aperture efficiency conversion. See sdcal description for information on fluxunit conversion. """ a=inspect.stack() stacklevel=0 for k in range(len(a)): if (string.find(a[k][1], 'ipython console') > 0): stacklevel=k break myf=sys._getframe(stacklevel).f_globals myf['__last_taskname']='sdfit' tb=myf['tb'] im = myf['im'] qa = myf['qa'] me = myf['me'] sm = myf['sm'] ia = myf['ia'] cl = myf['cl'] pl = myf['pl'] rg = myf['rg'] default=myf['default'] plotxy=myf['plotxy'] ### ### #Handle globals or user over-ride of arguments # function_signature_defaults=dict(zip(sdfit.func_code.co_varnames,sdfit.func_defaults)) for item in function_signature_defaults.iteritems(): key,val = item keyVal = eval(key) if (keyVal == None): #user hasn't set it - use global/default pass else: #user has set it - use over-ride # sys._getframe(1).f_globals[key]=keyVal myf[key]=keyVal ### # #Add type/menu/range error checking here # if type(mask)==str: mask=[mask] # # LIST PARAMETER NAMES HERE arg_names=['sdfile','telescope','fluxunit','specunit','frame','doppler',\ 'scanlist','field','iflist','fitmode','maskline','invertmask',\ 'nfit','fitfile','plotlevel',\ 'thresh','min_nchan','avg_limit','box_size','edge','xstat'] # # Assign parameter values arg_values = [] for ptst in arg_names: stst = myf[ptst] if ( type(stst) == str ): # Strip off leading and trailing whitespace stst = stst.strip() myf[ptst] = stst sex = ptst+"='"+stst+"'" else: sex = ptst+"="+str(stst) exec(sex) arg_values = arg_values + [stst] # # LIST ALLOWED PARAMETER TYPES HERE arg_types=[str,(str,list),str,str,str,str,(list,int),str,(int,list),str,\ list,bool,(int,list),str,int,(float,int),(int,float),\ (int,float),(float,int),(list,int),dict] #parameter_printvalues(arg_names,arg_values,arg_types) # # Check allowed types and values try: parameter_checktype(arg_names,arg_values,arg_types) parameter_checkmenu('fitmode',fitmode,['list','auto']) parameter_checkmenu('fluxunit',fluxunit,['K','k','Jy','JY','jy','']) parameter_checkmenu('specunit',specunit,['channel','km/s','GHz','MHz','kHz','Hz','']) parameter_checkmenu('frame',frame,['REST','TOPO','LSRD','LSRK','BARY','GEO','GALACTO','LGROUP','CMB','']) parameter_checkmenu('doppler',doppler,['RADIO','OPTICAL','Z','BETA','GAMMA','radio','optical','z','']) except TypeError, e: print "sdfit -- TypeError: ", e return except ValueError, e: print "sdfit -- OptionError: ", e return ### ### ### Now the actual task code ### try: # pdb.set_trace() #load the data without averaging s=sd.scantable(sdfile,False) # determine current fluxunit fluxunit_now = s.get_fluxunit() print "Current fluxunit = "+fluxunit_now # set flux unit string (be more permissive than ASAP) if ( fluxunit == 'k' ): fluxunit = 'K' elif ( fluxunit == 'JY' or fluxunit == 'jy' ): fluxunit = 'Jy' # fix the fluxunit if necessary # No way to query scantable to find what telescope is # so rely on user input if ( telescope == 'FIX' or telescope == 'fix' ): if ( fluxunit != '' ): if ( fluxunit == fluxunit_now ): print "No need to change default fluxunits" else: s.set_fluxunit(fluxunit) print "Reset default fluxunit to "+fluxunit fluxunit_now = s.get_fluxunit() else: print "Warning - no fluxunit for set_fluxunit" elif ( fluxunit=='' or fluxunit==fluxunit_now ): print "No need to convert fluxunits" elif ( telescope == 'GBT' or telescope == 'gbt' ): print "Convert fluxunit to "+fluxunit # THIS IS THE CHEESY PART # Calculate ap.eff eta at rest freq # Use Ruze law # eta=eta_0*exp(-(4pi*eps/lambda)**2) # with print "Using GBT parameters" eps = 0.390 # mm eta_0 = 0.71 # at infinite wavelength # Ideally would use a freq in center of # band, but rest freq is what I have rf = s.get_restfreqs()[0]*1.0e-9 # GHz eta = eta_0*pl.exp(-0.001757*(eps*rf)**2) print "Calculated ap.eff. eta = %5.3f " % (eta) print "At rest frequency %5.3f GHz" % (rf) D = 104.9 # 100m x 110m print "Assume phys.diam D = %5.1f m" % (D) s.convert_flux(eta=eta,d=D) print "Successfully converted fluxunit to "+fluxunit elif ( telescope == 'AT' or telescope == 'at' ): s.convert_flux() elif ( type(telescope) == list ): # User input telescope params if ( len(telescope) > 1 ): D = telescope[0] eta = telescope[1] print "Use phys.diam D = %5.1f m" % (D) print "Use ap.eff. eta = %5.3f " % (eta) s.convert_flux(eta=eta,d=D) elif ( len(telescope) > 0 ): jypk = telescope[0] print "Use gain = %6.4f Jy/K " % (jypk) s.convert_flux(jyperk=jypk) else: print "Empty telescope list" else: # Unknown telescope type print "Unknown telescope - cannot convert" # set spectral axis unit if ( specunit != '' ): s.set_unit(specunit) specunit_now = specunit else: specunit_now = s.get_unit() # reset frame and doppler if needed if ( frame != '' ): s.set_freqframe(frame) else: print 'Using current frequency frame' if ( doppler != '' ): if ( doppler == 'radio' ): ddoppler = 'RADIO' elif ( doppler == 'optical' ): ddoppler = 'OPTICAL' elif ( doppler == 'z' ): ddoppler = 'Z' else: ddoppler = doppler s.set_doppler(ddoppler) else: print 'Using current doppler convention' # Prepare a selection sel=sd.selector() # Scan selection if ( type(scanlist) == list ): # is a list scans = scanlist else: # is a single int, make into list scans = [ scanlist ] if ( len(scans) > 0 ): sel.set_scans(scans) # Select source names if ( field != '' ): sel.set_name(field) # NOTE: currently can only select one # set of names this way, will probably # need to do a set_query eventually # IF selection if ( type(iflist) == list ): # is a list ifs = iflist else: # is a single int, make into list ifs = [ iflist ] if ( len(ifs) > 0 ): # Do any IF selection sel.set_ifs(ifs) # Apply the selection (if any) s.set_selection(sel) del sel # Make line region masks and list of line regions if ( fitmode == 'list' ): # Assume the user has given a list of lines # e.g. maskline=[[3900,4300]] for a single line if ( len(maskline) > 0 ): # There is a user-supplied channel mask for lines if ( not invertmask ): linemask=s.create_mask(maskline) domask = True doguess = True # Make sure this is a list-of-lists # (e.g. [[1,10],[20,30]]) if ( type(maskline[0]) == list ): ll = len(maskline) nlines = len(maskline) linelist = maskline else: # Not a list, turn into one nlines = len(maskline)/2 linelist = [] for i in range(nlines): lo = maskline[2*i] up = maskline[2*i+1] linelist = linelist + [[lo,up]] else: # invert regions linemask=s.create_mask(maskline,invert=True) domask = True nlines = 1 linelist=[] doguess = False else: # Use whole region nlines = 1 linelist=[] doguess = True domask = False print "Identified ",nlines," regions for fitting" if ( doguess ): "Will use these as starting guesses" else: "No starting guesses available" else: # Fit mode AUTO print "Trying AUTO mode" if ( len(maskline) > 0 ): # There is a user-supplied channel mask for lines linemask=s.create_mask(maskline,invert=invertmask) domask = True else: domask = False # A single region nlines = 1 doguess = False # If we have line regions, get starting guesses linemax=[] lineeqw=[] linecen=[] if ( doguess ): # For each line get guess of max, cen and estimated # equivalent width (sum/max) if( len(linelist) > 0): for x in linelist: msk = s.create_mask(x) maxl = s.stats('max',msk) suml = s.stats('sum',msk) if ( maxl[0] != 0.0 ): eqw = suml[0]/maxl[0] else: eqw = 0.0 linemax = linemax + [maxl[0]] lineeqw = lineeqw + [eqw] cen = 0.5*(x[0] + x[1]) linecen = linecen + [cen] del msk else: # For what its worth, do stats on unmasked region if ( domask ): maxl = s.stats('max',linemask) suml = s.stats('sum',linemask) else: maxl = s.stats('max') suml = s.stats('sum') if ( maxl[0] != 0.0 ): eqw = suml[0]/maxl[0] else: eqw = 0.0 linemax = linemax + [maxl[0]] lineeqw = lineeqw + [eqw] # I dont know a better way to specify the center of the spectrum cen = s.nchan/2 linecen = linecen + [cen] # Now the line fitting if (nlines > 0): if ( fitmode == 'list' ): # Get number of things to fit from nfit list if ( type(nfit) == list ): numfit = len(nfit) nlist = nfit else: numfit = 1 nlist = [nfit] # Drop extra over nlines numfit = min(numfit,nlines) ncomps = 0 comps = [] for i in range(numfit): ncomps = ncomps + nlist[i] comps += [nlist[i]] else: # Auto mode - one comp per line region # Overwriting user-supplied nfit numfit = nlines ncomps = nlines comps = [] for i in range(nlines): comps += [1] print "Will fit ",ncomps," components in ",numfit," regions" if (numfit > 0): # Fit the line using numfit gaussians # Assume the nfit list matches maskline f=sd.fitter() f.set_function(gauss=ncomps) if ( domask ): f.set_scan(s,linemask) else: f.set_scan(s) if ( doguess ): # Guesses using max, cen, and eqw # NOTE: should there be user options here? n = 0 for i in range(numfit): if ( comps[i] == 1 ): # use guess maxl = linemax[i] fwhm = 0.7*abs( lineeqw[i] ) cenl = linecen[i] f.set_gauss_parameters(maxl, cenl, fwhm,\ component=n) n = n + 1 else: # cannot guess for multiple comps yet n = n + comps[i] else: # No guesses print "Fitting lines without starting guess" # Now fit f.fit() fstat = f.get_parameters() # Check for convergence goodfit = ( len(fstat['errors']) > 0 ) if ( goodfit ): # Plot residuals if ( abs(plotlevel) > 0): # each IF is separate panel, pols stacked f.plot(residual=True) # Retrieve fit parameters xstat['nfit'] = ncomps xstat['peak'] = [] xstat['cent'] = [] xstat['fwhm'] = [] for i in range(ncomps): fstat = f.get_parameters(i) xstat['peak'] = xstat['peak'] \ + [[fstat['params'][0],\ fstat['errors'][0]]] xstat['cent'] = xstat['cent'] \ + [[fstat['params'][1],\ fstat['errors'][1]]] xstat['fwhm'] = xstat['fwhm'] \ + [[fstat['params'][2],\ fstat['errors'][2]]] # Store fit if ( fitfile != '' ): f.store_fit(fitfile) else: # Did not converge print "ERROR: Fit failed to converge" xstat['nfit'] = -ncomps # Clean up del f # Assign variables for return myf['xstat']=xstat # Final clean up del s if (domask): del linemask # DONE except TypeError, e: print "sdfit -- TypeError: ", e return except ValueError, e: print "sdfit -- OptionError: ", e return except Exception, instance: print '***Error***',instance return saveinputs=myf['saveinputs'] saveinputs('sdfit','sdfit.last') ###End of sdfit def sdfit_defaults(): a=inspect.stack() stacklevel=0 for k in range(len(a)): if (string.find(a[k][1], 'ipython console') > 0): stacklevel=k break myf=sys._getframe(stacklevel).f_globals myf['sdfile']='' myf['telescope']='' myf['fluxunit']='' myf['specunit']='' myf['frame']='' myf['doppler']='' myf['scanlist']=[] myf['field']='' myf['iflist']=[] myf['fitmode']='auto' myf['maskline']=[] myf['invertmask']=False myf['nfit']=0 myf['fitfile']='' myf['plotlevel']=0 myf['thresh']=5 myf['min_nchan']=3 myf['avg_limit']=4 myf['box_size']=0.2 myf['edge']=[0] myf['xstat']={} def sdfit_description(key='sdfit'): desc={'sdfit': 'ASAP SD line fitting task\n', 'sdfile': 'name of input SD dataset', 'telescope': 'name or param of telescope for flux conversion', 'fluxunit': 'units for line flux (K,Jy,def=current)', 'specunit': 'units for spectral axis (def=current)', 'frame': 'frequency reference frame, e.g. LSRK (def=current)', 'doppler': 'doppler convention, e.g. RADIO (def=current)', 'scanlist': 'list of scans to use (e.g. [1,2,3,4])', 'field': 'string for selection by source name (e.g. FLS3a*)', 'iflist': 'list of IF ids to select (e.g. [0,1])', 'fitmode': 'SD line fitting mode (auto,list)', 'maskline': 'list of regions for line fitting (e.g. [[2500,3000],[3500,4500]])', 'invertmask': 'invert maskline regions (True,False)', 'nfit': 'list of number of Gaussians to fit in each line region (e.g. [2,1])', 'fitfile': 'output file name for fit results', 'plotlevel': 'control for plotting of results (0=none,1+=some)', 'thresh': '(Auto) S/N ratio limit for linefinder', 'min_nchan': '(Auto) min number of consecutive chans above threshold', 'avg_limit': '(Auto) max number of chans to average for broad lines', 'box_size': '(Auto) running mean box size as fraction of spectrum', 'edge': '(Auto) channels to drop from beginning and end', 'xstat': '(RETURN ONLY) Python dictionary with line fitting output', } return desc[key] ###End of task