Program: XGAUPROF
#begin section gettingstarted
GETTING STARTED
===============
XGAUPROF fits parameters of a number of Gauss related
functions to data from a GIPSY set or to data from file.
The functions include also fitting of skewness and
kurtosis parameters and Lorentzian line widths.
Enter the name of a GIPSY set in the 'Setname' input field.
Do not specify any axes or ranges yet.
If the set exists, you get a 'Profile axis' menu to the
right of the 'Setname' field. Select the axis that represents
the X-axis of your profile(s). A default grid range will
appear. Enter a box. If you press in an empty box
field, a default will appear. The first profile in the box
will be plotted and if you see a green curve (initial guess
curve) you can push the button to fit the profile.
#end section gettingstarted
#begin section buttons
DESCRIPTION OF THE BUTTONS AND INPUT FIELDS
===========================================
MAIN GUI
Usually you will start the program without specifying any
keywords. To get results you will have to enter the name
of a GIPSY image first (if you want to import data from other
sources, see description at ).
Then the 'Profile axis' menu will become active. This menu
contains the names of all the axes in your set. After this,
the default limits in grid coordinates of your profiles is
displayed in the input field 'Grid range'. This value can
be changed at any time to exclude unsuitable data. The input
field 'Phys. range' corresponds to the grid range. Here you
can change the limits in physical coordinates.
Pressing in the 'Box' input field will generate a default
box. Note that the default generates the maximum possible box,
so if you start with a big cube, then it can take a while
before all data is read. Finally you get a default in the
'Position' input field. Usually this is the first position
in the box and you immediately will see the corresponding
profile in the plot at the left of the gui.
You can go trough all the selected profiles with the
'<' and '>' buttons. If you have initial guesses (usually
plotted in green) then you can try to make a fit by pressing
the button. Results will appear in the panels just
below the input fields.
The output of results will be discussed later.
It is possible to exclude data by clicking with the left
hand mouse button in the neighbourhood of the point that you
want to exclude. In the window there is also
an option to exclude points by giving the range of points
that you want to include!
You can set a clip or range to the Y data in the Clip
input field in the same window. Values
below the clip (1 value entered) are marked in the plot
and get weight zero in the fit. If you entered two values,
then all values outside this range are marked and get
weight zero in the fit.
The button is designed to allow
the storage of the fit results on disk. It is activated only
after the selection of a proper file name in the menu.
The button redraws data and initial estimates.
This button starts a menu with the options:
1) 'save plot'
This option results in sending your current plot
to a printer or a PostScript file. In a sub menu you
will be asked to give a file name and to select a 'Plotter'.
If the file already exists, you will be warned. Writing
is started after pressing the button.
The file name that you entered, is displayed in the
gui and the GIPSY log.
2) 'save profile data'
Save current profile data to a file on disk. This can be
handy if you want to import your GIPSY data in a package
that reads ASCII files.
Be sure that you made a fit before saving the data, otherwise
the fitted amplitudes and the residuals are dummies.
Here is part of an example output file:
! X GRID X PHYSICAL AMPLITUDE FITTED AMPL RESIDUAL
-49 816.783 -0.229097 2.63399e-47 -0.229097
-48 820.927 -0.245038 7.42847e-46 -0.245038
.......
.......
3) 'Fit output file (short)'
Open a file to which the fit results can be sent. If
you entered a valid name, the button in the main gui
is enabled. Pressing this button results in writing the
current fit data to the ASCII file until the program is
finished. The layout of the file
is such that all output of one fit is written onto one
line. A header will inform you about the meaning of the
numbers. If a fit could not be made, a message is given.
4) 'Fit output file (long)'
Same as above, but the results are written in a so called
log format.
5) 'Events input file'
If you want to process a couple of profiles, it is sometimes
easier to store events in a file. There are two methods.
The common method is using the Ggi event player which is started
by clicking the right mouse button in the background of the
main gui.
However, this application has also it's own event player and
this one works as follows: Create an ASCII file on
disk and write the following lines in that file:
setname=real
profile=freq
box=-10 -10 10 10
GUI_FILE=2
PROMPTER_OK=y
PROMPTER_OK=y
fit=y
fstore=y
next=y
fit=y
fstore=y
Give it the name events.txt and enter this name in the
input field under , 'Events input file'.
The text will be read and events are generated immediately.
The events are just keywords with a value. The problem
is to find which keyword corresponds to a button.
If you are running XGAUPROF, type on the Hermes command
line: XGAUPROF MONITOR=Y
Press a button and you will see a log of the keyword(s)
that is/are involved. However you have to distinguish
yourself which events are generated by you and which are
generated by the system as a reaction on your action.
Note that this mechanism is different from the so called
GIPSY defaults file.
Pressing will pop up a window for importing
data from OTHER sources than a GIPSY image. Both reading
from text file and reading from a GIPSY table are supported.
The syntax for the rows is displayed in the help text
of the 'Rows' input field.
Example: Read all x and y values from file 'perfil.dat':
Filename: perfil.dat
Column X: 1
Y: 2
Errors :
Rows : :
(Pressing in the rows field has the same effect
as entering the colon.)
Note that in order to get initial guesses for the least
squares fit, your Y values must not be too small
(> 1e-8).
For the least squares fit routine it doesn't matter how
the X-values in your data are distributed, however for
the initial estimates routine we need equally spaced
X-values. Therefore XGAUPROF will eliminate duplicates
in X and uses a spline interpolation to create a dummy
profile with equidistant X-values and without blanks.
If a spline interpolation could not be made (e.g. all
Y values blank except one) then no initial estimates
could be found.
Note that the fit routine works with the raw data.
The use of these inputs is straightforward.
If you fix the amplitude range, then each profile in the
box will get the same range in Y-values (=amplitudes).
The button clears the plot and redraws data and
estimates.
'amp. %'
Set the percentage of the maximum at which the user wants
the width. This option works only for standard gaussians
because only for this function there is a simple relation
between fitted dispersion and the width at arbitrary height.
Turn valleys into peaks with this button. The original data
is multiplied by -1. In order to restore the original data you
have to push the button again.
'tolerance'
This window contains two parts. The first part contains
parameters which control the least squares fit routine.
The tolerance is a stop criterion. its default is zero,
which implies that the program determines the minimum value.
If the fit routine cannot improve the reduced chi-squared
by 'tolerance', it stops and returns the current fit
parameters.
'lambda'
Lambda is the mixing parameter. it determines the initial
weight of steepest descent method relative to the Taylor
method. Lambda should be a small value (i.e. 0.01).
It can only be zero when the partial derivatives are
independent of the parameters which is not the case in
this program. In practice, this values will almost never
be changed.
'max. Its'
The least squares routine fails to produce results if it
needs more than this number of iterations to achieve
convergence.
'Include range'
Enter the range in physical coordinates for which you want
to include the data. Other points are masked and get zero
weight in the least squares fit. If you enter an empty input
then the entire data range is taken into account and the
masked points are reset.
'Clip'
You can set a clip or range to the Y data. Values
below the clip (1 value entered) are marked in the plot
and get weight zero in the fit. If you entered two values,
then all values outside this range are marked and get
weight zero in the fit. The clip is applied to all the
profiles not only the current.
'Min. ampl.'
Initial estimates with an amplitude smaller than this value
are discarded. See also 'Q'.
'Min. disp.'
Initial estimates with a dispersion smaller than this value
are discarded. The input has the same units as the X-values
in the plot. Usually these are physical units and not
grids. See also 'Q'.
'Q'
This is called the 'smoothing factor'. 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 stimated. The maximum amplitude
is derived from the observed profile.
The default is a series of Q's (1 2 3 4 & 5). If you enter
more than one Q then the Q with the best (i.e. lowest chi-
squared) initial estimates is used for the least squares fit.
The estimate filters 'Min. ampl.' and Min. disp.' apply
only to this 'best' Q.
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 any dispersion results in the same reduced chi-squared
and therefore there is no convergence.
'Median filter'
If a profile contains spurious elements, you can use a
three pixel median filter to smooth the data a bit.
In many cases this is less work than excluding these points
by hand.
There are three options in the weights menu. First, the
default, is uniform weight, i.e. each non blank point in
the profile gets weight 1. The reduced chi-square that is
calculated is corrected for the current rms of that profile
afterwards.
Second option is an automatic weighting. The program calculates
a second profile that is a three pixel median filter of the
original profile. The distance in Y between these two profiles
is a measure for a weight of each point. With this method,
spurious elements get less weight in the fit, which improves
the fit results for some profiles.
The third option allows you to use errors from an
ASCII file on disk as specified in the menu.
If the uncertainty of a data point is given by E then the
corresponding weight is equal to 1/(E*E) normalized to the
average of all weighting factors (Bevington, Data Reduction and
Error Analysis, Ch 10).
This window contains input parmaters and fit results.
Gaussian parameters:
The param window allows you to control the initial estimates
as proposed by the program. If the program cannot find any
estimates, perhaps you can and this is the place to do that.
This is also the window where you tell the fit routine
to set gaussian parameters to fixed or free. The more you fix,
the better the fit.
Background:
The background is taken into account in the fit. It has
three parameters. A constant, a linear and a quadratic term.
In the default situation only the constant term is a
free parameter.
Slider:
Each parameter has a button. If you push this button,
you activate the slider for this parameter. While moving the
slider, you can watch the parameter value changing in the
value field and in the plot.
rms:
Initial estimates can only be determined if a reasonable
value for the rms of the noise is given. If you have a value
in mind, push the button to the right of the rms
input field and enter your value. Otherwise, a default is
calculated. The profile data is sorted and the distance
between the quartiles is taken a a measure of rms of the
noise. Usually this value is higher than expected, but
it works well to find initial estimates.
Below right in this window there are two buttons and
that have the same function as the buttons with the same
labels in the main gui.
A menu is presented with four options:
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.
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.
4) Voigt. Fit Lorentzian and Doppler widths together
with center and amplitude.
The mathematics behind all this is worked out
in a separate document.
Select the number of gaussians that you want to fit.
Show this document. Push the button again to get rid of
the help window.
OUTPUT
------
The different functions generate different output, but there
are some quantities that are the same for all. First there
is output in the status area; the number of iteration to
reach convergence, the total flux, reduced chi-square, the
standard deviation of the residuals and
the current value of the fit tolerance are displayed.
In the so called output panel, the fit parameters are listed.
These are the parameters that characterize a distribution
like the line strength (or area) of the profile, the
mean X value (which is only the center if the distribution
is symmetric like the standard Gauss and the Voigt profile),
Further the maximum amplitude is given and also the
values for the background are listed. For the gauss-Hermite
series, the skewness and kurtosis are listed. The real function
parameters of the Gauss-Hermite functions are not listed in
the gui log, but can be found in the GIPSY log. These 'real'
function parameters are transformed into properties of a
distributions. An extra parameter gives the position
of the maximum which is, as mentioned before, for the
GH-series not the same as the value of the mean X.
A description of the functions and their parameters can be
found in a separate document. Here is a summary for the
GH-series.
The GH-series with h3 and h4:
f = a.Exp[-1/2.g^2] *
{ 1 + h3[c1.g+c3.g^3] + h4[c5+c2.g^2+c4.g^4] }
and:
g = (x-b)/c
c1..c5 have fixed values.
The output in short format to file disk has a table header
and columns with results. The first line can look like:
! GAUSS fit parameter table created: Tue Jul 6 16:25:04 1999
!RCHI2 |SIG.RES |AREA | err | MEAN | err | DISP. | err | ......
8.97364 1.04838 1030.28 141.687 2237.89 ......
8.97364 1.04838 1030.28 141.687 2237.89 .....
.....
etc. etc.
The meaning of the titles:
RCHI2: reduced chi-squared. This value is scaled by the
rms of the noise that you entered. If you did not
enter a rms, the value 1.0 is substituted.
SIG.RES: Standard deviation of the residuals.
AREA: Total flux under fitted function.
err: Error on this parameter.
MEAN: Mean of gaussian of first component.
DISP: Dispersion
FW at 50.0 pct: Full width at certain percentage of maximum.
AMPL. Amplitude
Z0/1/2 Background terms constant, linear and quadratic.
RA PARAM Position of profile in grids. These are the axis
names of the box.
INITIAL ESTIMATES
-----------------
In a separate document we describe the method and give
a reference to the article of U.J. Schwarz on which the
estimate method is based.
LEAST SQUARES FIT
-----------------
The least squares fit minimizes the reduced chi-square
which relates a model profile to the real data. If the
decrement in the value of the reduced chi-squared is less
than a given tolerance, then the iteration stops and
the fitted parameters are returned.
Maximum accuracy can be obtained by setting the tolerance to
0 which is also the default value.
BLANKS
------
Blanks in the profile get weight zero in the fit. For
the estimate routine, blanks are replaced by the current
zero level. In the plot window the blanks are indicated by
character 'b'.
#end section buttons
#begin section keywords
KEYWORDS
========
CHANGE THE SIZE OF YOUR GUI AT STARTUP
--------------------------------------
keyword: GUISIZE=
arguments: size in X in pixels
size in Y in pixels
height of residual window in pixels
For building events files you need to know the keywords used by
the application. Here is a (incomplete) list with possible keywords
with some default values. You can use the Ggi monitor function to
find the keywords corresponding to certain events. The monitor
tool is started when you click the right mouse button in the
background of the main gui.
AMPPERC= 50
CHARSIZE= 1.10
COLCON= 8
COLEST= 3
COLFIT= 2
COLFRAME= 1
COLMED= 4
CONNECT= NO
CRITAMP= 0
CRITDISP= 0
DRAWZERO= YES
ESTFILTER= NO
FITFILE= profilefit.dat
FUNCTION= 0
GGIOPT= FIXCOL
GUI_FILE=
HCFILENAME= profilefit.ps
HELP=
LAB= 0.01
LWIDTH= 2
LWIDTHCUR= 2
MAXITS= 200
NCOMP= 1
PARAM=
PLOTDENS= 200
PLOTEST= YES
PLOTMED= NO
PLOTRMS= NO
PLOTTER= 0
PROFFILE= profiledata.dat
Q= 1 2 3 4 5
RMS= 1.0
SUBDIV= 1 1
TOLERANCE= 0
VIEWPORT= 0.15 0.1 0.95 0.9
WEIGHTS= 0
Z0_EST= 0
Z1_EST= 0
Z1_FIX= YES
Z2_EST= 0
Z2_FIX= YES
#end section keywords
#begin section about
Purpose: Fit parameters of a standard Gauss function,
a Gauss-Hermite series or a Voigt lineshape to data from
GIPSY sets or ASCII files
Category: ANALYSIS, PROFILES, PLOTTING
File: xgauprof.src
Author: M.G.R. Vogelaar (GUI: J.P. Terlouw)
Keywords: See section keywords
Version: 2.1
Example:
Notes:
Related document: xgauprof.ps
Updates: Feb 22, 1999: VOG, Document created.
Jun 15, 2001: VOG, Help updated, connect lines added.
Apr 15, 2009: VOG, Removed unused macro NINT and inserted
string.h in import.c and param.c
Jun 28, 2010: JPT, Changed log font.
#end section about