COPYRIGHT (c) 1999 Kapteyn Astronomical Institute University of Groningen, The Netherlands All Rights Reserved. Program: XGAUFIT Purpose: XGAUFIT estimates parameters for multi component one dimensional gauss functions using data in selected profiles. It can use these estimates to fit the parameters to the data in a least squares fit. Output is send to subsets in a set. Category: ANALYSIS, PROFILES, FITTING File: xgaufit.src Author: M.G.R. Vogelaar (GUI: J.P. Terlouw) Keywords: This is a list of keywords from an example run of the application. It should be obvious to which gui elements these keywords correspond. AMPPERC= 50 AXUNITS= KM/S BOXAXES= RA DEC CRITAMP= 0 CRITDISP= 0 DATAUNITS= W.U. FUNCTION= 0 GAUSET=fittedprofiles GRIDRANGE= -30 30 LAB= 0.01 MAXITS= 200 NCOMP= 1 OUTSET=fittedparams PARAM= PHYSRANGE= 552.668 1548.98 PROFILE= freq PROFINDX= Q= 1 2 3 4 5 RMS= 1 SETNAME= real TOLERANCE= 0.001 VSYS= 1050 Z0_EST= 0 Z0_FIX= YES Z1_EST= 0 Z1_FIX= YES Z2_EST= 0 Z2_FIX= YES Buttons etc.: Notes about the Graphical User Interface (GUI): User input in this application is coupled to the Hermes keyword mechanism. Keywords can still be specified on the Hermes command line and also the syntax for input of numbers and text is the same as for keywords. Also, Emacs-like text editing in input fields is possible and we added mouse functionality. E.g. the right mouse button will erase text in an input field. With the middle button you can paste text stored in your mouse buffer. , and ctrl-U Wipe current text from input field ctrl-A Go to beginning of field ctrl-E Go to end of field ctrl-K Wipe from cursor position to end of line (ctrl-Y does NOT yank back) Go to next input field. Control keys in the Help window: ctrl-S Search forward ctrl-R Search backwards ctrl-V Scroll text one page down (same as Page-Up) ctrl-P Scroll text one line up ctrl-Z Scroll text one line down Notes: 1) Do not close pop-up windows using the cross symbol on the upper windows bar. This will crash your application. Instead use the close button or repeat pressing the button in the main gui with which you opened your pop-up window. 2) On some systems, a menu item is un-selectable if the Caps-Lock is on! MAIN WINDOW: Setname: input field. Enter name of set which contains your profile data. Profile axis: Menu List of possible profile directions. Grid range: Input field. Default is the entire profile length in grids. If you change one of the values, the corresponding physical range is altered too. Phys range: Input field. Default is the entire profile length in physical coordinates. If you change one of the values, the corresponding grid range is altered too. Box: Input field. The limits in your set between which profiles are examined. A box consist of two vectors. The first gives the low limits and the second gives the high limits. The dimension of a vector is the set dimension - 1. If you press you will get a default box which includes all profiles in the set. Vsys: Input field Value used in sorting when you want to sort in distance to a position along the profile axis. Usually this will be the systemic velocity of a galaxy. Parameter output: Input field Name of the output set on disk which holds all the fitted parameters in subsets. Size of the output set is determined by the program. Fitted data output: Input field Name of the output set on disk which holds the fitted data. The output set is copied from the input set. Initially all data is set to blank. You can get a residual cube by subtracting this set from the input set e.g. with task SUB. GO/STOP: Button Start calculating initial estimates and fitting the parameters to the data in the profiles. If you want to abort this process, press this button again and wait a while for the program to respond. BUTTONS in button bar: FILE ==== Menu options: 1) Get estimates from these subsets instead of calculating them. You must give the name of the set and a subset specification. E.g. myestset param 1 2 3. You need as much subsets as components x 3. 2) exit program FILTERS ======= Amp. range: Input field No values means no filtering. One value filters out all amplitudes below this value and two values filter out all amplitudes below the first and all amplitudes above the second value. Vel. range: Input field Disp.range: Input field See above. These last input fields need values in physical coordinates. FIT OPTIONS =========== Tolerance: Input field Stop criterion for least squares fit (see description below). Lambda: Input field Mixing parameter (see description below) Max. its: Input field If the stop criterion is not reached after this number of least squares iterations, then abort the fitting of the current profile. Write blanks in the output set in the corresponding position. RMS: Input field Set the value with which the value of chi-squared is scaled to get a value for the reduced chi-squared that is related to the noise in your profile. The default is 1. Min. ampl: Input field Discard initial estimates with an amplitude lower than this value. Min. disp: Input field Reject initial estimates with a dispersion lower than this value. Q: Input field Smoothing factors (see description below) PARAM ===== Set fit parameters to fixed or free in the fit. If you want to set a parameter to fixed, you have to enter a value for this parameter. For the background terms, an initial estimate must be given in the input fields. The defaults for all the background terms are zero. FUNCTION ======== Menu with the following functions (see description below) Gauss Gauss Hermite (h3) Gauss Hermite (h3,h4) Voigt NCOMP ===== Menu with possible number of components per function. SORT ==== Menu with sort options: None Peak, Disp, Vsys Select sort option if you have selected more than one component to fit and want to sort these components. See also see description below. HELP ==== Pop up this dc1 document Description: BASIC INPUT The decomposition of profiles into Gaussian components can be a powerful method in the study of the kinematics of neutral hydrogen. The name of the set for which you want to study your profiles is entered in the main gui of the program. The input field is called 'Setname'. If a set is known, then the profile direction is selected from a menu labeled 'Profile axis'. Usually this will be the FREQ or VEL axis of your data set. Note that you can select a profile direction as long as your input set has more than one axis so that input is not restricted to a three dimensional data cube. The length of your profile axis must be greater than 1 and less than 4096. The length of a profile is set in the input fields 'Grid range' and 'Phys. range'. A default is displayed after the name of a set and the profile axis are known. You can change these values either in grids or physical coordinates. If you change one field you see that the values in the other field will change also. HINT: Inspect with XGAUPROF whether you can exclude unwanted data by decreasing the grid range. However if your grid range is too small you can run into trouble while fitting base lines. The 'Box' field is blank. Pressing results in a default box which is the entire data range in grids. Note that the box can be given in any coordinates as you were used to in the non-gui GIPSY tasks. In general, the program will work on profiles in any direction and instead of velocity you could think of any kind of physical coordinate. The program automatically inserts the correct units in the header of the output set. PREPARING THE FIT Estimates calculated by the program ----------------------------------- You can fit the parameters of four different functions to the profile data. The results are written to a GIPSY output set on disk, each parameter in a different subset per profile position. But the least squares fit needs a good starting point for these parameters and an important part of this program is dedicated to finding good estimates by fitting second degree polynomials to the profiles (See: Schwarz, 1968, Bull.Astr.Inst. Netherlands, Volume 19, 405-413). An important role here play the so called smoothing parameter 'Q'. The initial estimate routine uses 2Q+1 points around a maximum to fit a second order polynomial. The coefficient of the second-order term of the polynomial is an approximation of the second derivative of the observed profile. With a moments method the center and dispersion are estimated. The maximum amplitude is derived from the observed profile. In practice there is not one smoothing factor that gives the best (in terms of chi-squared) estimates for all profiles. Therefore you can enter in the 'FIT OPTIONS' menu a series of smoothing factors ('Q'), which are all used in all profiles to get the best estimates. You can exclude estimates by entering values for a minimum amplitude ('Min.ampl') and the minimum dispersion ('Min.disp'), but usually you leave these fields blank and set a filter on the least squares fit results. Note that there is one special case build in the program. If a profile contains blanks or zero-level values except for one, then the estimate routine fails to return an estimate. One can treat this peak as a peak with known amplitude and center, and a full width at half maximum of one pixel (which is automatically converted to a dispersion). However, the least squares fit routine will not fit such profiles because it converges to zero dispersion. Estimates from an input set --------------------------- Alternatively it is also possible to read estimates from a set. This set could be an output set of a previous run of XGAUFIT. If not, then you must be sure that the box of your estimates set matches the box of your input set and that the data has the same physical coordinates as the set with the profiles. Not only the set name must be entered, but also the subsets, e.g.: myestimates param 1 2 3 With this option it is possible to separate the data with which you create estimates, and create fits. For instance you can create estimates with clipped or conditional transferred data and use these estimates to fit parameters using the original profile data in a second step. Such methods avoid the noisy parts of your data and result in 'cleaner' output. THE FIT PROCEDURE ----------------- The actual fitting is a least-squares fit of a function to a set of data points. The method used is described in: Marquardt, J.Soc.Ind.Appl.Math. 11, 431 (1963). It is a mixture of the steepest descent method and the Taylor method. Weights ------- The fit routine uses equal weights for all data that is not blank. Blank values get weight zero. Tolerance --------- The fit is controlled by three parameters which are set in the 'FIT OPTIONS' window. Most important parameter is the stopping criterion 'Tolerance'. Fitting is stopped when successive iterations fail to produce a decrement in reduced chi-squared less than TOLERANCE= The value cannot be less than a certain minimum as set by the system. This means that maximum accuracy can be obtained with TOLERANCE=0.0 Usually only signal is present in a small region in the profile while the chi square is calculated over the entire baseline. In this case the Chi square primarily represents the fitting of the baseline instead of the gauss and then a very tiny decrease of the Chi square might very well be significant for the fitting of the gauss to the signal. So one should be very careful not setting TOLERANCE= too high. If TOLERANCE=0 then the minimization is stopped at machine precision. Increasing the tolerance results in a shorter process time. If the tolerance cannot be reached within the maximum number of iterations, the fit failed and blanks are written in the output set. HINT: In most cases the default tolerance 0.001 gives fast performance and good fit results. Mixing factor ------------ Lambda in the 'FIT OPTIONS' window is the mixing parameter in the least squares fit routine. Lambda determines the initial weight of steepest descent method relative to the Taylor method. It should be a small value (i.e. 0.001). Maximum number of iterations ---------------------------- Fitting is also controlled by the maximum number of iterations 'Max. Its.' in the 'FIT OPTIONS' window. If the maximum number of iterations is exceeded, then the fit failed and blanks are written in the output set. Fixed or free parameters ------------------------ In the 'PARAM' window you can alter the status of the fit parameters. Default, all function related parameters are 'free' parameters (i.e. to be fitted) and the background parameters (a constant, a linear and a quadratic term) are fixed at a value 0.0. Error on fitted parameters -------------------------- The errors on the fitted parameters are calculated as real errors by means of the covariance matrix in the fit and the average deviation between data and fit over the entire profile. However, if the region with signal is short compared to the full length of the profile, the average deviation between data and fit is primarily determined along the part of the baseline where no signal is present and then it becomes almost equal to the noise in the profile and the real errors determined in this way virtually become formal errors. Note that if Hanning smoothing has been applied, the subsets are two by two correlated. The errors on the parameters are then a factor sqrt(2) higher. MULTI COMPONENTS AND SORTING ---------------------------- If you fitted parameters of multicomponent functions, then it could be useful to sort the fitted values. Usually you will sort components in INCREASING order of distance to a user given value on the profile axis (often the systemic velocity) as given in Vsys in the main window. The sort menu is found in the upper right area of the main gui. Other sorting options are: sort components in DECREASING order of peak value and sort components in DECREASING order of dispersion FILTERS ------- The program has three filters for the output. You can enter ranges for amplitude, velocity (or whatever is on the profile axis) and dispersion. Fitted parameters outside these ranges are filtered out, i.e. they are set to blank in the output. If no values are entered then nothing will be filtered. If only one value is entered for one of the filters, then everything BELOW this value is filtered out. In the gui log you can monitor how many profile fits are rejected by these filters. Warning: The velocity filter is straightforward for standard Gauss functions. If one fits parameters of the Gauss-Hermite series, then the filter applies to parameter b which is only the mean (velocity) of the profile distribution if the parameters h3 and h4 are zero. FUNCTIONS --------- The function menu provides the choice of four functions of which you can fit the parameters to the data in your profiles. 1) A standard Gaussian. Parameters are amplitude, center, and dispersion, 2) Gauss-Hermite polynomial (h3). Parameters are a,b,c, (which are NOT the same as amplitude, center and dispersion) and h3. The parameters are translated to the more familiar amplitude, center, dispersion and skewness. All parameters mentioned above are saved in the output set. 3) Gauss-Hermite polynomial (h3, h4). Same as above, but an extra parameter h4 is included. The value for this h4 parameter is translated to kurtosis. Note that if one of h3 and h4 is not zero, the mean of the distribution is not the position of the maximum. The program calculates this parameter separately. (Reference; Marel, P. van der, Franx, M., A new method for the identification of non-gaussian line profiles in elliptical galaxies. A.J., 407 525-539, 1993 April 20). 4) Voigt. Fit Lorentzian and Doppler widths together with center and amplitude. The mathematics behind all this is worked out in a separate document: (http://www.astro.rug.nl/~gipsy/xgauprof/index.html) OUTPUT ------ A name is sufficient to specify the output set. If the set exists, you are asked to confirm to write it over. Before writing any data, an existing set will always be deleted first to avoid complaints about incompatible sizes if you change boxes or grid ranges and do not alter the name of the output set. The program creates the new set on disk and copies relevant header items. It replaces the profile axis by a new axis called PARAM-XGAUFIT. Its length depends on the maximum number of components to be fitted and the selected function. It is also possible to save the estimates only (Est.only button in 'FIT OPTIONS' window). The Hermes log displays the structure of your output. This information, together with relevant units, is also written to the comment fields header of the output set so that you can always retrieve subset information e.g.: HISTORY ITEM=C MODE=R INSET=myout or: FIXHED INSET=myout ITEM= which shows header items at set level or FIXHED INSET=myout p ITEM= to show units , minimum, maximum and number of blanks at subset level. Not only the fitted parameters are saved, but also the errors on these parameters, the background parameters, total flux, number of iterations, reduced Chi-squared of the fit, the standard deviation of the residuals and the smoothing factor 'Q' used by the estimate routine. After finishing the fitting process, the program writes a description of the subsets to the comment fields of the header of the output set. The units, minimum and maximum values and the number of blanks are written on subset level in the FITS keywords BUNIT, DATAMAX, DATAMIN and NBLANK. Output of the fitted profiles ----------------------------- There is an option to reconstruct your data using the fitted parameters. These 'fitted profiles' are written to a set which is a copy of the input set, i.e. the header is a copy and the data is set to 'blank'. The size of the output is limited by the values for the box and the profile range (Grid range or Phys. range). The header keywords for the minimum and maximum data value and the number of blanks are NOT updated. This update can be done with task MNMX. Each fitted profile is written to this output set. If the fit was not successful a blank appears in all the profile subsets of the output set. A message is appended to the HISTORY of the output set. The message contains the task name and time/date. Updates: Jul 20, 1999: VOG, Document created. Jul 27, 2000: VOG, g->Yfit was not initialized to NULL. This seemed only to cause problems on PC's(?) Nov 6, 2001: JPT, Wait for processing finished before quitting. Apr 15, 2009: VOG, Removed macro NINT Jun 28, 2010: JPT, Changed log font.