MGFIT N1 N2 N3 Infile Outfile THRLO THRHI NSPLIT FIT RESET SET
Fits multiple gaussian components to spectra
allowing many different constraints on parameters.
*** NB. MGFIT has a reject cycle loop. If this is set, then on exit the
spectra will have masked points. This allows one to examine
which parts have been masked, but you may want to keep an
unaltered copy of the spectra as well.
mgfit uses Marquardt's method of minimisation. This involves a mixture of
steepest descent and using a second-derivative matrix to get to the minimum.
In difficult cases (many variables) you will see it coming up with the same
value of chi**2 several times in a row as it moves more towards steepest
descent having found that the quadratic method is not doing well. This is
normal behaviour. You should however try to give it as reasonable a start
as possible in such cases because it is not guaranteed to work. Identical
chi**2 are a sign of degeneracy or near degeneracy in the model variables
being fit (see below) so you may want to check whether you can re-specify
your model.
This method is pretty poor at finding the global minimum, and you may have
to hold its hand a lot on its way to it. It has the advantage that once
there, one can get uncertainties from it.
Parameters:
N1, N2 -- Range of slots to fit.
N3 -- First slot to store fits. Useful for comparing model
with data, 0 to ignore.
Infile -- Input file defining the model to be used; see below.
Outfile -- Output file with fitted parameters.
THRLO, THRHI -- Lower and upper sigma reject thresholds. Pixels will
be rejected and the fit repeated until no more are rejected.
Thresholds measured in units of RMS of fit. See note above.
NSPLIT -- Number to split each exposure in as a means of accounting
for binary motion during an exposure. In this case the
model will consist of a trapezoidal integration over the
exposure. Set = 1 for no account for smearing. This slows
done the fits since every gaussian and derivatives have to
computed NSPLIT times over.
FIT -- There are two fit methods:
'I' -- The variables are fitted separately for
each spectrum.
'A' -- The variable parameters are fitted to all
the spectra at once.
RESET -- Yes to reset to the starting values provided in the input
file, No to use the previous spectrum's final fitted
values.
SET -- yes to set bad pixel mask (if you know that some
parts will not be fitted, it is better to reject
them a priori than let the routine try to do it for
you.)
Input files:
Because of the complexity of the possible models, file input is used.
This section now describes the format required. It does so with a
series of examples of increasing complexity. If you are new to this you
may want to adopt a similar approach to the fitting.
Fitting a single gaussian plus a polynomial.
-------------------------------
poly: 6500. $const $grad
gaussian: 6562.76 $off $fwhm 2.
$const = 1.
$grad = 0.
$off = 0.
$fwhm = 25.
------------------------------
WARNING: make sure you include the '.'s in the numbers. I don't know
whether this is a feature or a bug, but I can't be bothered to fix it.
The polynomial is indicated by starting "poly:". This must be exactly
this with no leading blanks. It is then defined by a pivot wavelength
and coefficients. In this case 6500 is a fixed value whereas $const
$grad are variables (the constant and gradient terms of the poly) or
at least potentially variable -- see below for how they can be fixed.
The gaussian is similarly specified starting with "gaussian:". In this
case the parameters are (1) the central wavelength (angstroms),
(2) a velocity offset (km/s), (3) a full width half max (Angstroms) and
(4) the peak height in mJy. In the example given, only the offset and
FWHM are variable.
The variable names are case sensitive and must always begin with a $.
Each variable must be given an initial value as in "$fwhm = 25.".
The initial value lines can have no leading blanks. The names can be
8 characters at most, excluding the $.
FIXING "VARIABLES"
One often wants to switch which parameters vary or not. This can be done
by adding an "F" for fixed at the end of its initialisation line, e.g.:
$fwhm = 25. F
In fact in complicated cases it is likely that you would want most of the
variables to be held fixed in this manner while you optimise subsets of them.
This should be emphasized: the routine cannot magically find a solution
if you give it too poor a start, and in complex cases you will need to "hold
its hand".
"VARIABLES" VERSUS "PARAMETER DEFINITIONS"
In the example above each parameter defining the model is associated
with just one variable. However, it is useful to distinguish between
"variables" and "parameter definitions". For example consider the
following:
-------------------------------
poly: 6500. $const $grad
gaussian: 6562.76 $off $fwhm $area/$fwhm
$const = 1.
$grad = 0.
$off = 0.
$fwhm = 25.
$area = 20.
------------------------------
The last item of the gaussian represents its height and is here
"$area/$fwhm$". $area is the variable that will be fitted here
and the height of a gaussian is proportional to its area divided
by its width as specified. "$area/$fwhm" is a "parameter definition"
while $area and $fwhm are the "variables". Parameter definitions
tell mgfit how to compute the value of a parameter from the available
variables and can consist of series of terms joined by the *, /, + or
-. Brackets are not supported at the moment and there must be no
blanks. Thus "2.5*$fwhm+10.-$height*$area" is legal, while
"2.*($fwhm+10.)" is not.
The differentiation between the variables and the parameters required
to specify the gaussian permits great flexibility in fitting. Here is
another example:
----------------------------------------
gaussian: 6562 $off 6562*$vfwhm/$c $haght
gaussian: 4861 $off 4861*$vfwhm/$c $hbhgt
$off = 0.
$vfwhm = 1000.
$hahgt = 2.
$hbhgt = 1.
$c = 2.9979e5 F
----------------------------------------
This fits one gaussian to Halpha and another to Hbeta. They are
forced to have the same width in velocity (units of km/s, $vfwhm)
by the constructs 6562*$vfwhm/$c etc. They are allowed different
heights but have the same velocity offset from their central
wavelengths ($off). The use of one variable in many different
components is a good way of reducing the chances of degeneracy
in the fits.
DEGENERACY
With flexibility comes the room for error and it is very easy to
set up a file with degeneracy in the variables which may cause
mgfit to bomb out, although not crash molly (I hope). e.g.
---------------------------------
gaussian: 6562 $off1+$off2 10. 1.
$off1 = 0.
$off2 = 10.
---------------------------------
is clearly a degenerate case. In general it is not easy to check
for these beforehand and so I have not tried to. It is up to you
to be sensible when creating the fit file. What happens during
minimisation is that a matrix with near zero determinant is found
and you will get a message about a routine called GAUSSJ failing.
This will cause mgfit to either move on to the next spectrum or
stop altogether. If this happens, check your model. It may also
be a result of starting too far from the solution; you can use the
no variable option (see next) to check for this.
NO VARIABLES
A file like this:
---------------------------------
gaussian: 6562 0. 10. $height
$height = 1. F
---------------------------------
has no variables and no fit will be performed. However, if you have
specified an output slot, it will compute a model and store it there.
This is a useful way to develop initial models in hard cases where if
your initial guess is too far off, the routine fails to find a minimum.
It also allows you to compute your model for arbitrary spectra.
SPECTRUM-TO-SPECTRUM VARIATIONS: ACCESSING HEADER PARAMETERS
You may want to fit a model to all spectra at once but account for
some sort of variation from spectrum-to-spectrum. For instance suppose
you were fitting a line gaussian in shape but with a height that
you suspected varied sinusoidally with time. Then a file like so could
be what you want:
---------------------------------
gaussian: 6562 0. 10. $const+$hnorm*soff
$const = 1.
$hnorm = 1.
---------------------------------
The new element here is "soff" (NB it has no $). This will be interpreted
as a header parameter that must be of type double. It is up to you
to have set this for every spectrum. That is, every spectrum will need
to have a double precision header item called 'soff' with a value set.
Note that case DOES matter here and you cannot have blanks in its name.
In this case soff could have been set to the sine of the HJDs for example.
RESOLUTION PARAMETER
As of June 1999, molly now recognises a parameter representing the
FWHM resolution of a spectrum as in:
---------------------------------
gaussian: 6562 0. $fwhm $height
resolution: FWHM
$fwhm = 2.
$height = 1.
---------------------------------
The meaning of this is as follows: suppose you have several spectra of
the same object, taken with different resolutions. Then by specifying
a resolution, in this case as a header parameter FWHM, the program
tries to find an "underlying model" representing the lines. i.e. it
will fit a single $fwhm and $height, blurr them by the appropriate
resolution and then fit to the data. The resolution is added in
in quadrature so that the fitted FWHM = SQRT(FWHM**2+$fwhm**2) and the
fitted height = $height*FWHM/SQRT(FWHM**2+$fwhm**2).
Another use for the resolution parameter is to provide a lower limit to
the fitted widths. This can be useful in preventing narrow spikes being
fitted to cosmic rays for instance.
The resolution parameter only has any meaning if you have some gaussian
components.
ORBITAL FITS
mgfit allows you to specify orbital parameters for each gaussian.
Thus:
---------------------------------
gaussian: 6562 $gamma 10. 1. $kx $ky 2450000 $period
$gamma = 0.
$kx = 100.
$ky = 150.
$period = 0.1
---------------------------------
will fit a fixed gaussian moving on an orbit of the form:
V = gamma - Kx*cos(2*PI*phase) + Ky*sin(2*PI*phase)
with the phase computed from the ephemeris which can be fitted
just like the other parameters. (Note here that the zero point
of the ephemeris is fixed here to prevent degeneracy. See above.)
---------------------------------
gaussian: 6562 $gamma 10. 1. $kx $ky $kx2 $ky2 2450000 $period
$gamma = 0.
$kx = 100.
$ky = 150.
$kx2 = 10.
$ky2 = -5.
$period = 0.1
---------------------------------
will also fit a second harmonic:
V = gamma - Kx*cos(2*PI*phase) + Ky*sin(2*PI*phase)
- Kx2*cos(4*PI*phase) + Ky2*sin(4*PI*phase)
MISCELLANY
The central wavelength of the gaussians is also a possible variable,
thus:
gaussian: $wave 0. 10. 1.
is perfectly legal. The only thing that cannot be fitted is the
polynomial pivot wavelength; I just felt it was unlikely ever to be
much use to fit this one.
Here is a way of specifying a multigaussian fit for a fixed ephemeris:
gaussian: 6562 $off1 $fwhm1 $height1 $kx1 $ky1 ZeroO PeriodO
gaussian: 6562 $off2 $fwhm2 $height2 $kx2 $ky2 ZeroO PeriodO
etc, where ZeroO and PeriodO is the orbital ephemeris set using
"Phase"
Any line that cannot be translated into any of the above is
regarded as a comment.
PROBLEMS
Mind you have a carriage return at the end of the last line of
the input file otherwise it will go unrecognised.
This command belongs to the classes: fitting , measurement