Gauss-Hermite Pixel Fitting Software
***************************
PART I: GENERAL INFORMATION
***************************
Author
------
* Roeland P. van der Marel,
1992-1994 : algorithm and program development
1994-1997 : minor modifications and installation on Sun Ultra 170
address : Space Telescope Science Institute
Research Programs Office (RPO)
3700 San Martin Drive
Baltimore, MD 21218
Tel : (+1) 410 338 4931
Fax : (+1) 410 338 4596
e-mail : marel@stsci.edu
homepage : http://sol.stsci.edu/~marel/
Summary of the method
---------------------
The line-of-sight velocity profile of the stars in a galaxy is
inferred from an observed galaxy spectrum, by fitting the spectrum
with the convolution of a template spectrum and a velocity profile
that is parameterized as a Gauss-Hermite series. The best-fitting
parameters are determined by chi-squared minimization in Pixel space.
References
----------
The method and an application are discussed in the following paper:
- van der Marel 1994, MNRAS, 270, 271
which can be retrieved form my homepage.
Acknowledgments
---------------
If you have found this software useful for your research, I would
appreciate an acknowledgment to `use of the Gauss-Hermite Pixel
Fitting Software developed by R.P. van der Marel'.
Bug Reports
-----------
Please send bug reports and important comments or questions to
marel@stsci.edu
Caveats
-------
This software comes without guarantee and on an `as is' basis. It is
not being actively maintained or upgraded.
All original testing of the code was done on a UNIX Sparc 2
workstation with the SunOS operating system in the period
1992-1994. In 1997 I ensured that the code would run on a Sun Ultra
170 workstation with the solaris operating system. However, little
detailed testing of the code was performed on this operating system.
It is not guaranteed that the software will work without problems on
other platforms, with different operating systems and different
compilers. However, the code is close to standard fortran, so if any,
only minor revisions will be necessary.
The programs read and write files in IRAF (.imh) format. Use of the
IRAF format requires linking of the code to standard IRAF libraries.
I will consider all requests for use of this software. However, I do
not allow people who have received the software from me to distribute
it to third parties. Anyone who uses this software is encouraged to
send me an email with their address, so that I can send the reports of
bugs and revisions.
The determination of stellar kinematics, and in particular
line-of-sight velocity profiles, from galaxy spectra can be very
tricky. Proper care should be given to details such as the continuum
subtraction, low-frequency filtering and template matching. Using this
software as a black-box in a quick-and-dirty manner will almost
certainly lead to erroneous results.
Directory structure
-------------------
The directory structure of the software is as follows:
prog : source of the main programs
bin : executables of the main programs
lib : object libraries
sub : subroutines used by the programs
iraf : subroutines for reading IRAF files
examp : examples of the use of the programs
Installation
------------
To install the software edit the file setup.com. All site dependent
environment variables should be made to point at the correct directories
and libraries. Then issue the UNIX command
prompt> source install.com
This will construct all executables and put them in the directory bin.
Once the programs have been installed and the executables created, the
environment variables in setup.com are not needed anymore. However, they
are required every time any of the programs is recompiled. Users who
will be modifying the programs regularly may wish to add the command
source ....path..../setup.com
to their .cshrc file (where ....path.... should be the path to the directory
that contains the file setup.com).
Main Programs
-------------
There are 5 main programs:
preptemp.f : Prepares the template spectrum. A 1D template spectrum is
read. It is convolved with Gaussian VPs
of various dispersions and the results are written to
a new 2D output frame
prepgal.f : Prepares the galaxy spectrum. A mask file is written for
use in the VP fitting. The radius corresponding to each row
is calculated and written to file. The galaxy spectrum is
rebinned spatially (upon request), and written to a new
2D frame. Another 2D frame is written which contains the
formal noise error in each pixel.
pixfitgau.f: The best-fitting Gaussian VPs are determined to the
rebinned galaxy spectrum, using the 2D template frame
produced by preptemp.f. The parameters of the Gaussians
are written to file, as well as the parameters of the
best-fitting continuum in the fit. Optionally, 2D frames
are written containing the fit to each row, the residuals
of the fit, and the weighted residuals (i.e., the residual
divided by the formal noise).
pixfitGH.f : The best-fitting Gauss-Hermite moments are determined to
the rebinned galaxy spectrum, while keeping the
Gaussian fit parameters fixed to the values determined
by pixfitgau.f. The results are written to file.
Optionally, 2D frames are written containing the fit to
each row, the residuals of the fit, and the weighted
residuals.
pixfitemis.f : A residuals frame created with either pixfitgau.f or
pixfitGH.f is read. This program then fits single or double
Gaussians to emission lines that show up in
the residuals. The results are written to file.
Optionally, 2D frames are written containing the fit to
each row, the residuals of the fit, and the weighted
residuals. Also optionally, Ascii files with the data and
the fits are written (so as to allow easy plotting with
e.g. mongo).
The programs should be run in the order in which they are listed. Each
program reads output of one or more of the previous programs. Obviously,
pixfitGH.f and pixfitemis.f need not be run if one is not interested in
velocity profile shapes or emission lines.
All programs read and write all their in- and output data in the IRAF .imh
format. FITS files can be transformed to this format using the task
dataio.rfits in IRAF.
Auxiliary program
-----------------
There is one small auxiliary program:
rmscalc.f : This program reads the (.gau1) output file of pixfitgau.f,
and calculates the RMS velocity SQRT[V^2+sig^2]
and its error. Results are written to a new (.gau3)
output file. The systemic velocity that is to be assumed
can either be given by hand, or be calculated by the
program as the straight mean of all the velocities in
the .gau1 file.
Examples
--------
The directory examp contains example galaxy and star spectra, as well as
examples of the use of the programs. To run the examples issue the UNIX command
prompt> source examples.com
To understand the examples read the comments in the file examples.com.
Known Bugs:
-----------
Currently none.
*********************************
PART II: Discussion of the method
*********************************
Introduction
------------
This manual describes the use of the Gauss-Hermite Pixel Fitting
software which extracts line strengths, radial velocities, velocity
dispersions, and the Gauss-Hermite velocity profile coefficients from
galaxy spectra.
1. Method
---------
The basic assumption is that the galaxy spectrum is the convolution of
a suitable template spectrum and a certain line-of-sight velocity
profile (VP).
The template spectrum should be of similar type as the averaged light
of the galaxy population. In practice K0 III templates work reasonably
well for elliptical galaxies. To obtain accurate results for the VP
deviations from a Gaussian it is often better to adopt a more general
approach. For example, one might adopt a weighted mix of stellar
spectra of different spectral types that optimizes the fit to the
galaxy spectrum (see e.g.: Rix & White 1992; van der Marel 1994).
The velocity profile is expanded in a so-called Gauss-Hermite series
(van der Marel & Franx 1993; see also Gerhard 1993 for a related but
different approach). The lowest order term of the series is a
Gaussian. The higher order terms are orthogonal to this Gaussian. The
lowest order anti-symmetric deviations from a Gaussian are quantified
by the Gauss-Hermite moments h3 and h5, the lowest order symmetric
deviations by h4 and h6. The Gauss-Hermite expansion is useful for
various reasons: (i) it is a straightforward extension of the
conventional assumption of Gaussian VPs; (ii) it yields a set of well
defined uncorrelated parameters; (iii) these parameters can be
measured with accuracy from real data with reliable error bars. Other
than that, there is nothing magical about the expansion. The
program was not designed to determine VPs in an unparameterized manner.
The basic problem (where o stands for convolution)
Galaxy = Template o VP
can be solved in many different ways to obtain the best fitting VP
(e.g., Franx & Illingworth 1988; Bender 1990; Rix & White 1992;
van der Marel & Franx 1993; Kuijken & Merrifield 1993; Saha & Williams 1994;
van der Marel 1994). In the Pixel Fitting method we minimize the
chi_squared defined by:
chi_squared = SUM_pixels [ Gal - (Temp o VP) ]^2
This is is done in a brute force manner. The convolutions (Temp o VP) are
actually evaluated, and simply compared to the data in a chi^2 sense.
Results obtained with the Gauss-Hermite Pixel Fitting Software have
been compared to the results of other methods, in particular to that
of Rix & White (1992) [which also minimizes chi_squared in pixel
space] and of van der Marel & Franx (1993) [which minimizes
chi_squared in Fourier space]. Agreement was found to be excellent.
This appears to indicate that: (i) the software contains no major
bugs; and (ii) that the Gauss-Hermite moments can be determined
accurately from data of sufficient quality.
The different methods do appear to return somewhat different formal
error bars on the fitted parameters in cases where the chi^2 of the
best fit is significantly larger than the number of degrees of freedom
(the number of pixels that is fitted). This occurs often in practice
as a result of template mismatching. Error bars returned by the Pixel
Fitting Software can be a factor 2 to 3 smaller than those returned by
the Fourier Fitting Software of van der Marel & Franx (1993), and a
factor 1 to 2 smaller than those returned by the method of Rix & White
(1992). In cases where there are no systematic errors (only
Poisson Noise), the formal errors of the different methods all agree.
Template spectra cannot generally be expected to fit the galaxy
continuum. Fitting the continuum requires a full stellar population
synthesis. To date, no method has been developed to simultaneously do
a stellar population synthesis and a stellar kinematical
analysis. Hence, all existing methods for VP analysis essentially
solve
Galaxy = [Template o VP] + arbitrary continuum terms
as a fudge to correct for this.
In our Gauss-Hermite Fourier Fitting Method, a sum of continuum terms
is fitted to the data simultaneously with the broadened template. The
continuum terms are Legendre polynomials up to a given order, to be
specified by the user.
References:
-----------
Bender R., 1990, A&A, 229, 441
Franx M., Illingworth G.D., 1988, ApJ, 327, L55
Gerhard O.E., 1993, MNRAS, 265, 213
Kuijken K., Merrifield M.R., 1993, MNRAS, 264, 712
Rix H.W., White S.D.M., 1992, MNRAS, 254, 389
Saha P., Williams T.B., 1994, AJ, 107, 1295
van der Marel R.P., Franx M., 1993, ApJ, 407, 525
van der Marel R.P., 1994, MNRAS, 270, 271
***********************************************
PART III: Description of program in- and output
***********************************************
0. Input formats
----------------
To successfully run the programs one needs a 1D-template spectrum and a
2D galaxy spectrum, both in IRAF .imh format. The galaxy frame must have
the spectra along its rows. It is assumed that the spectra are completely
reduced, sky-subtracted, and most importantly: rebinned
logarithmically in wavelength. The rebinning and number of pixels etc.
should be the same for the galaxy and template spectra.
!!! If the input spectra do not conform to these conventions,
!!! the software will produce erroneous or no output.
The spectra should be present in the directory from which the
executables are run. No additional files are required. The programs ask
for input on the command line.
1. preptemp.f
-------------
The program pretemp.f reads an input template spectrum. It convolves
this spectrum with Gaussian VPs, and writes the results to a 2D
output frame. The program asks the user to specify the following
parameters:
velscale : the number of km/s per pixel in the input spectrum.
velshift : the mean velocity of the Gaussian VPs. This should be
an estimate of the velocity difference between the
galaxy spectrum that will be analyzed and the input
template spectrum.
sigmin, sigstep, Nsig: these parameters specify the array of velocity
dispersions that will be used. The output frame will
have Nsig rows. Each row i will contain the convolution of the
template spectrum with a VP of dispersion sig(i), where:
sig(i) = sigmin + [sigstep*(i-1)]
The range sigmin to {sigmin+[sigstep*(Nsig-1)]} should
bracket the region within which the velocity dispersions of
the galaxy spectra fall.
iextend : When doing the convolution it is necessary to know the
behavior of the input spectrum outside the range where
it is known. The parameter iextend specifies how the
extrapolation is done.
1 : set the spectrum to zero
2 : set the spectrum to the boundary value
3 : extend the spectrum using the gradient between the
first and last column.
ires : When doing the convolution the program specifies the VP
on a finer grid in velocity than the input spectrum. This
is necessary to obtain proper results. The pixel size
of the VP array is set equal to
velscale / ires
where ires should be an integer power of two. One should
choose ires such that (velscale/ires) << sigmin.
name : The name of the 2D output frame
An intrinsic parameter of the program (which need not be specified by the
user) is sigmany. VPs are set to zero beyond sigmany times the Gaussian
velocity dispersion. The default value is sigmany = 6.0.
On output the program also writes a file name.log, containing a
listing of the parameter settings that were used. This file is read by
subsequent programs.
2. prepgal.f
------------
This program reads a 2D input frame with galaxy spectra. It then prepares
for the VP fitting to be done by subsequent programs. It first writes
a mask file which contains a listing of the pixels that should be
masked in the fitting. To do so, the user has to specify
a number of parameters:
sky ? : The program asks whether the last row of the frame contains the
sky spectrum that was subtracted from each of the other rows.
If this is not the case, the program assumes the sky value
to be zero. If a sky spectrum is present, the program
asks for two more parameters:
threshold : all pixels in which sky > threshold will be masked. This
is useful in cases where strong sky lines did not properly
subtract.
neighmask : the first neighmask neighbors of pixels with sky > threshold
are also masked.
The user is also asked to specify any additional regions that need to be
masked. This is useful, e.g., in cases where the galaxy spectra contain
emission lines.
Then the user is asked for:
ncolfirst, ncollast : all subsequent VP fitting will only use the
columns ncolfirst to ncollast in the galaxy spectrum.
One generally doesn't want these to be the true first and
last column of the spectrum, because the value of the
broadened and shifted template cannot be properly calculated
close to the edges of the spectrum.
The program proceeds by calculating the formal noise error in each pixel.
For this it needs:
readnoise : the read-out noise per pixel in the input spectrum (in e/pix).
gain : the gain of the input spectrum (in e/ADU).
The noise of the sky spectrum is also taken into account in the calculation.
Then the program assigns a radius to each row of the input frame
according to: radius = r0 + (scale*(irow-nrowgcen))
The parameters in this expression must be set by the user.
scale : the number of arcsec/pixel in the input spectrum
nrowgcen : the row in the input spectrum that is closest to the
galaxy center
r0 : the radius corresponding to row nrowgcen
If the user does not want to assign a radius to each row he can specify
this. In this case the program will automatically set:
scale = 1.0 , nrowgcen = 0 , r0 = 0.0
The user is then asked for:
nrowfirst,nrowlast : only rows between these extremes will be used in the
VP fitting.
A kinematical analysis of a galaxy spectrum requires sufficient S/N.
The program therefore offers an option to rebin all rows between
nrowfirst,nrowlast in the spatial direction so as to achieve some S/N.
It therefore asks for:
stonmin : the minimum average S/N per pixel in each row of the
rebinned frame.
The program first rebins the rows closest to nrowgcen until it has
sufficient S/N. It then proceeds by rebinning the row numbers larger
than nrowgcen, and subsequently those with row numbers smaller than
nrowgcen.
The program asks for the name of the 2D output frame. It writes two 2D
output frames:
name.imh : frame with rebinned galaxy spectra
nameE.imh : frame with the formal noise error in each pixel
The program also writes:
name.mask : the information on the mask to apply in the fitting
name.rad : information on the rebinning and radial scale
name.log : listing of parameter settings that were used
These files are read by the subsequent VP fitting programs.
3. pixfitgau.f
--------------
This program asks for the names of the 2D output frames generated by
preptemp.f and prepgal.f. It then determines the best-fitting Gaussian
VP to each row of the rebinned galaxy spectrum. It minimizes chi^2 in
pixel space, leaving out all pixels that were masked by prepgal.f.
The program asks for:
npolmax : the maximum polynomial order used in the continuum fitting.
The results are written to files name.gau*, where name is the name of the
frame with the rebinned galaxy spectra.
name.gau1 : lists in table format:
i : input row number
rad : the corresponding radius
chisq : the chi^2 of the best fit
gam, dgam : the best-fitting line-strength and its error
v, sigvel : the best-fitting mean velocity and its error
sig, dsig : the best-fitting dispersion and its error
name.gau2 : lists in table format the parameters of the continuum fit.
These results are read by the program pixfitGH.f, but
will not generally be of interest to the user.
One more file is written:
name.gaulog : listing of parameter settings that were used.
These output files are read by pixfitGH.f.
Upon request, the program can write 2D frames with the best fit to
each row, the residuals, and the weighted residuals (i.e., the
residual divided by the formal noise). If constructed, these frames
are called: nameGF.imh, nameGR.imh and nameGW.imh. They only contain
the columns ncolfirst to ncollast as defined in the mask file (because
outside this region the fit is meaningless), and hence generally have
a somewhat smaller format than the input file name.imh. These frames
are not used by subsequent programs. They are written solely for the
purpose of visualization.
Users who wish to experiment with the software might find it useful to
write a little program that reads on input the mask file and the file
with weighted residuals, and on output writes a new mask file in which
all pixels with very large residuals are masked. By running the VP
fitting programs a second time with the new mask, somewhat more
accurate results may be expected.
4. pixfitGH.f
-------------
This program asks for the name of the 2D frame analyzed by
pixfitgau.f. It then reads a variety of output files of previous
programs, and determines the best-fitting Gauss-Hermite moments, while
keeping the Gaussian fit parameters fixed to those obtained by
pixfitgau.f. It minimizes chi^2 in pixel space, leaving out all pixels
that were masked by prepgal.f.
The program asks for:
maxGH : the maximum Gauss-Hermite order to be determined
i12zero : switch that determines whether h1 and h2 are to be kept fixed
to zero during the fit.
The results are written to name.GH1, where name is the name of the
frame with the rebinned galaxy spectra. It lists in tabular format
i : input row number
rad : the corresponding radius
h0, dh0 : the zeroth order Gauss-Hermite moment and
its formal error.
etc.
One more file is written:
name.GHlog : listing of parameter settings that were used.
Upon request, the program can write 2D frames with the best fit to
each row, the residuals, and the weighted residuals (i.e., the
residual divided by the formal noise). If constructed, these frames
are called: nameHF.imh, nameHR.imh and nameHW.imh. They only contain
the columns ncolfirst to ncollast as defined in the mask file (because
outside this region the fit is meaningless), and hence generally have
a somewhat smaller format than the input file name.imh. These frames
are written solely for the purpose of visualization.
If one runs the program with maxGH=2 and i12zero=0, one generally
finds values of h1 and h2 very close to zero, indicating that
pixfitgau.f has properly found the best-fitting Gaussian. However, if
one runs the program with maxGH > 2 and i12zero=0 one can find values
of h1 and h2 of order a few percent. This is because the
Gauss-Hermite moments are uncorrelated only in an idealized
situation. The orthogonality of the Gauss-Hermite functions is broken
when one includes continuum fitting, masking, etc. On the other hand,
the Gauss-Hermite moments are always close to being uncorrelated, and
the results for h3, h4, etc. seldom depend strongly on the value of
i12zero.
5. pixfitemis.f
---------------
This program asks for the name of a 2D frame that has been analyzed by
pixfitgau.f or pixfitGH.f. It then reads the residuals file created by
either of these programs (which one of the two must be specified by
the user). The program then fits either single or double Gaussians to
emission lines (or more generally: peaks) that show up in the
residuals. It does so by minimizing chi^2 in pixel space. The errors
in the individual pixel values are taken into account.
The program asks for a variety of input information:
pcen : the pixel number (need not be an integer) in the input
RESIDUALS file that corresponds to emission at the
systemic velocity.
irow1, irow2 : the first and last row in the input residuals frame
to which Gaussians are to be fitted.
vmax : the fits are done over the region -vmax to vmax (with
respect to the systemic velocity given by the user).
Pixels outside this region are ignored.
iverb : determines whether verbose output should (1) or should
not (0) be printed on the screen while doing the fits.
ngauss : the number of Gaussians that is to be fitted to each emission
line. Note that it is often useful to fit two Gaussians,
even if there is only one line. This allows, e.g.,
emission line asymmetries to be determined.
The program requires initial guesses for its parameters. The parameters are:
a1, (a4) : Gaussian amplitude(s)
a2, (a5) : Gaussian mean velocity(ies) [in km/s]
a3, (a6) : Gaussian dispersion(s) x SQRT(2) [in km/s]
There are two possibilities. Either an input file must be provided which lists
for each row between irow1 and irow2 the following information:
irow, a1, a2, a3, a4, a5, a6
Or alternatively the required (3 or 6) initial guesses may be given on
the command line. In this case the user can specify whether these
initial guesses should be used for each row, or just for the first
row. In the latter case, the result for each row will be used as
initial guess for the next row.
The results are written files name.GH1 and name.GH2, where name is
the name of the input residuals frame. All output velocities are in km/s.
name.GH1 lists in tabular format:
i : input row number
rad : the corresponding radius
counts, dcounts : the total number of counts in the best fit, and
corresponding formal error
v, dv : the mean velocity of the best fit, and its formal
RMS, dRMS : the RMS velocity of the best fit, and its formal error
sig, disg : the dispersion of the best fit, and its formal error
skew, dskew : the skewness of the best fit, and its formal error
chisq : the chi_squared of the best fit. To interpret this number it
should be compared to the number of pixels that was actually
fitted to, which is printed on the screen.
name.GH2 lists in tabular format:
i : input row number
rad : the corresponding radius
a1, da1
....
a6, da6 : the parameters of the best fitting double or single Gaussian
function, with its formal errors.
Upon request, the program can write 2D frames with the best fit to
each row, the residuals, and the weighted residuals (i.e., the
residual divided by the formal noise). These frames contain all
columns of the input residuals frame, not just the ones that were
fitted to. If constructed, these frames are called: nameEF.imh,
nameER.imh and nameEW.imh.
Also upon request, the program can write Ascii files with the data for
each row as a function of velocity (name.em3) and the fits for each
row as a function of velocity (name.em4). These files allow easy
plotting with, e.g., mongo.
Finally, one more file is written:
name.emlog : lists the parameter settings that were used.
6. User Beware: calculation of redshift
---------------------------------------
To calculate the mean velocity of the best fit (v), the program
determines the shift in logarithmically rebinned spectral space. The
result is multiplied by the speed of light (c). This yields the
correct answer for nearby objects, for which this program was
originally developed. However, for distant objects one is generally
interested in the redshift z. Given the definition of v, the redshift
follows from:
z = EXP(v/c) - 1
Note that this is not the same as the relativistic Doppler equation.
7. Error messages -- Array sizes
--------------------------------
The parameters listed in paramval.h in the directory /prog determine
the sizes of the arrays that are reserved by the programs. When the
input data is too large for these arrays an error message is
generated, and the program aborts. In this case one should increase
the sizes of the parameters in paramval.h, and reinstall the programs.
Alternatively, your machine might have too little memory available to
reserve all the required space. In this case one might wish to decrease
the sizes of the parameters in paramval.h, and reinstall the programs.
Roeland van der Marel.
Return to my home page.
Last modified September 19, 2004.
Roeland van der Marel,
marel@stsci.edu.
Copyright © 1998 The Association of Universities for
Research in Astronomy, Inc. All Rights Reserved.