Top 10 Rules of Thumb for Galaxy Fitting

(Actually, there are 11....)

1. To save from making subtle yet serious mistakes, it is usually better work with data that have pixel values in [counts] rather than count rate [counts/sec] or other units, e.g. μJy, ergs/sec, etc..

• Image units that are close to the basic counting units (e.g. counts or electrons, as opposed to count rate) help one to know when the sky value is correct or incorrect. Getting the sky right can be a tricky business for a number of reasons: flatfielding errors, neighboring contamination, and a mismatch between the model profile and the data, etc.. A small error in measuring the sky by even a few electrons may affect the fit results significantly. When the pixel values of an image are divided by the exposure time, say hundreds or thousands of seconds, and you make an error by a few counts, the apparent degree of the error is reduced by that arbitrary factor and may go unnoticed. For example, 1.203×10-2 versus 1.20×10-2 [counts/sec] may represent a significantly large difference, but you may not know it either because the software you're using may not show you enough decimal places or, more likely, you may think the number is too small to matter without thinking about it in a more proper context. With units in [counts] it is much easier to detect errors visually by moving the mouse cursor around the image. It is also easier to quickly compare the sky level with uncertainties in the mean of the sky, i.e. RMS/√Npix, which is also a very small number when the image is divided by the exposure time. Moreover, in counts or electrons, it is fairly easy to compute basic statistics in one's head to immediately know what's reasonable or not. As an explicit example, consider whether the difference between 1.203e-2 versus 1.20-2 electrons/second given 12500 second exposure time, and over 1e4 pixels is statistically significant by doing the calculation in your head. Then compare that with when the units are in counts: 150.375 versus 150.00 electrons over 1e4 pixels. In both examples the significance of the difference is about 3σ, but it is much easier to know that in the latter situation. When the image data units are in [counts/sec] it is easy to make mistakes that are otherwise trivial to avoid when the units are in [counts].

• For GALFIT to create a sigma (σ) image internally, the image ADU has to be in units of [counts] rather than [counts/second], such that:

GAIN [e-/count] × ADU [count] × NCOMBINE = [electrons],

where GAIN is the detector gain parameter from the image header, for a single image. NCOMBINE is the number of images combined. Feeding GALFIT an image in [counts/sec] and asking it to create σ-image internally, bad and strange results may happen, even if the residuals may look fine, because the relative scaling of noise between the background and the data could be wrong.

If one wants to work with images in [counts/sec] or other units, it is not advisable to let GALFIT create the sigma image internally. See here and here for more details.

• Sometimes images may have units of [counts/sec] but the standard header keywords BUNIT doesn't reflect that, EXPTIME is not one second, and GAIN still has units of [counts/ADU]. This is confusing because it is common to expect GAIN × ADU to have units of [electrons]. For statistics and to detect errors, it may be necessary to know how to properly convert from image ADUs back to electrons.
2. Make sure the background sky value is correct.

• As discussed in this link, the sky is an important parameter that the user should make sure to get right. A wrong sky value is the most common cause for systematic errors. GALFIT can fit the sky as a free parameter, but it is necessary for users to confirm that the sky is reasonable, not significantly affected by a mismatch in the profiles between the data and the model, and not have flatfielding problems.

• SExtractor sky is not the best value to use for galaxy fitting. While the sky determined by SExtractor is fine for many purposes, that value often over-predicts the background level for galaxy fitting because there is no way for SExtractor to know when a galaxy profile ``ends.'' So it determines the sky by deciding when the gradient flattens out enough compared to the background noise. Overestimating the intrinsic sky level suppresses the Sérsic index of a model, causing the measured size to appear smaller than actual and the luminosity to appear fainter. This would affect high Sérsic index objects more than low Sérsic index objects, because high index models have more extended wings. The proper way to determine the sky is to use a technique which can extrapolate the flux out to infinity if there are no other objects nearby. Model fitting is the only way to properly estimate the sky level for highly extended sources.

• It is almost always better to work with images where the pixel units are in [counts] instead of count rate [counts per second]. As mentioned above in Item 1 but worth repeating: for catching fitting errors, it is hard to overemphasize how much more intuitive it is to perform analysis on images where the pixels are in counts rather than counts/sec or other units. If you have a data image in counts/sec and also have a sigma image (also in units of counts/sec) to go along with it, you can simply multiply both images by the exposure time.

• If one must work with images in [counts/sec] or other units for some reason, make sure to compare the sky GALFIT estimates with not only the background mean value determined some other way, but also with the uncertainty in the mean, which is not just the sky RMS but the sky RMS/√Npix, where Npix is the number of pixels used to estimate the mean sky value.

• One way to make sure the sky is fitted well is to try fitting a galaxy and its neighbors with multiple components, to remove as much of the residuals as possible. If all the components behave (i.e. fairly low Sérsic indices below 4, size not crazy, axis ratio reasonable, the components make sense) and the sky value estimates don't change when you introduce more and more complexity into the model, then it is likely that the sky estimate is pretty robust. If the sky value does change significantly from one fit to the next, then the value from the more complex analysis is generally more reliable, if all the components make sense. For this to work well, the fitting region has to also include a lot of sky area.

3. Make sure the sigma image is correct.

• The reason for using a sigma image and what it is are discussed here. A correct sigma image should look like the data, qualitatively, except the pixel data values will be scaled like the square root of the pixel value. Please make sure to inspect the sigma image visually if you are providing one. If you are not providing a sigma image, please follow the instructions here for how to let GALFIT create a sigma image.

• If you want GALFIT to calculate the sigma image, make sure that the data have units of [counts] in ADUs, not [counts/second] or other units, such that:

GAIN × ADU × NCOMBINE = [electrons].

Otherwise, GALFIT will not know how to calculate a proper sigma image. Here, again, GAIN is that for a single image, and NCOMBINE is the number of images combined to create the image.

• The units for the sigma image and the data image should agree, e.g. one shouldn't be in units of [counts] and the other in [MJy/sr], or [counts/second], etc..

• The exposure time map that many data reduction software provide (e.g. MULTIDRIZZLE), sometimes with the file extension ".wht," is not a sigma image.

4. Use a good PSF for convolution.

• The use of the convolution PSF, and how to make one, are explained here.

• Make sure the PSF has a very high signal-to-noise compared to the data, is sky subtracted, is better than Nyquist sampled (FWHM > 2 pixels), and is large enough in radius (diameter > 30 times FWHM).

• If you have doubts about whether the PSF image size is too small or large, try testing the fit by using different size PSFs. If the answer doesn't change between using a smaller PSF image and a large one, then that is good: either one will do, but the smaller one will save time if you have to fit the data several times. If the answer does change, the result obtained using a larger PSF image is always more accurate unless the PSF has problems at large radius, e.g. neighboring contamination.

5. Use object masking and simultaneous fitting to remove neighboring contamination.

• Every object in the image should either be fitted simultaneously or masked out. The way to mask out objects in an image is either to use SExtractor segmentation image or to do so manually, as discussed here.

• GALFIT can take SExtractor segmentation image directly as a bad pixel mask: all non-zero pixels in the segmentation image are ignored by GALFIT in the fit. However, be sure to "zero out" regions in the mask that you want GALFIT to fit or else GALFIT won't see it.

6. Use the convolution box correctly/wisely.

• The reason why the convolution box is needed is discussed here.

• Objects with high concentrations often require larger convolution boxes than objects with flatter centers, because more of the flux is distributed from the center into the wings. Convolving an image of the sky does nothing.

• When in doubt, fitting with larger convolution box is always more accurate than with a smaller box, if there is a difference between the two. If there is no difference, choosing the smaller convolution box will save you time if you have to fit the data several times. Figure out what is large enough by trial and error. Or, start with a small convolution box until one is happy with a result, then run a last pass by using a large convolution box.

7. Start out simple, build up complexity.

• GALFIT has no problem fitting multiple, non-overlapping, objects simultaneously from the start, even in an automated way. However, if the goal is to do detailed decomposition on an object using multiple components, it is generally a good idea to start from the simplest and most basic model, then build up complexity as the situation requires. For instance, when fitting a quasar host galaxy (i.e. a galaxy with a strong central point source), it is usually better to fit the quasar component first, then take the result and add a galaxy component beneath. There are several reasons for doing so. One is that it saves time and is a smarter way to go than to throw in two components blindly from the start. Otherwise, if multiple components fitting the same region of space have bad initial guesses, there is a chance that one or more of the components may get moved out of the way by the dominant one in the first few iterations.

8. Be on the look-out for small object sizes/axis ratio, and any parameters enclosed in *...*.

• When the size (effective radius, scale length, etc.) is too small (i.e. r < 0.5 pixel), or when the axis ratio is too small ( q < 0.1), the user should beware. GALFIT will start to converge inefficiently in such situations or not at all, because the gradient of the models starts to become unresolved. When this happens one should check to see what in the image is causing this behavior (e.g. neighboring contamination, a nucleated core). If this behavior persists, take the best fit result and hold the parameters fixed to the minimum radius/axis ratio and restart the fit.

• Usually, when the size of a component starts to be < 0.5 pixels, it is generally advisable to replace that component with a PSF function instead of trying to fit it using other functions.

• As of Version 3.0.1 (July 15, 2010), GALFIT will warn you if a parameter gets into problematic regimes by enclosing the value in *...*, e.g. *0.09*. GALFIT will also place an "ERROR" header keyword into the model image slice [2]. Any parameter in *...* ought to be reset to some other values, held fixed, or replaced with other components, and then one should rerun GALFIT. If there is any parameter enclosed in *...*, one cannot assume that other parameters without *...* have correctly converged; chances are good they have not.

9. Look at the model images.

• This may sound obvious, but it often goes neglected. When the signal-to-noise of an object is low, the residual image may look perfectly fine, but the model may be affected by a neighboring object. Or, the component one thinks is being fitted might be fitting something else instead. Be sure to inspect the model every time.

• Note that running GALFIT by doing ``galfit -o3 filename'' will create an image block, ``subcomps.fits,'' where each model component is produced separately in its own image slice. When fitting multiple components, especially when using high order modes, it is a good idea to not only check the overall model, but also the individual model components.

10. Do not abuse high order parameters (Fourier modes, C0, bending modes, spiral structures, truncations).

• Just like Rule #6, high order parameters should be used incrementally rather than from the first attempt at fitting a galaxy. High-order parameters (diskiness/boxiness, Fourier modes, spiral rotation functions, truncations) are modifications to the best fitting ellipse models. As such, they should be used only once you have obtained the best fitting, classical, ellipsoid models. If one throws in these parameters before obtaining a good classical fit, the large initial mismatch between the data and the model guess will look like large amplitude perturbations, causing high-order parameters to go hay-wire and possibly not converge at all. The relevant issue here is not one about degeneracy. Even though there is degeneracy, the main issue is about the concept of using high order parameters as perturbations to something. If the goal is to quantify deviations from a best fitting ellipse, rather than from any ellipse, the best fitting ellipse has to be obtained first.

11. Avoid Using Constraint Files.

• You may use a constraint file when you don't want GALFIT to freely search the parameter space. However, this feature is used much more frequently than recommended. Numerically, it is not easy to find a global minimum when regions of parameter space are fenced off, which may introduce new local minima which the code used to be able to ``go around." A lot of research has been dedicated to addressing this issue in the computation world. The Levenberg-Marquardt algorithm GALFIT uses is one of the simplest implementation and may not suffice to deal with parameter constraints optimally. However, more to the point...,

• If one has to use parameter constraints, often times there is something abnormal about the analysis situation. For instance, sometimes the reason the Sérsic index grows unrealistically large is due to the presence of a central point source, due to there being a neighboring galaxy, or due to the presence of another galaxy component not being accounted. While it is possible to put a constraint on fitting parameters, it might only lead to completely meaningless results rather than one which is closest to underlying truth. It is therefore better to figure out what is causing the problem than to let parameter constraints dictate the solution.

• There is one exception to the rule about using constraint files: constraining parameters of different components using hard constraints (see EXAMPLE.CONSTRAINTS in GALFIT package) is a perfectly fine thing to do because the process is mathematically well defined.

• It is also perfectly ok to hold components fixed to values using the toggle flag in the main GALFIT menu file.