Constructing the PSF with IRAF/DAOPHOT

Concise instructions

Kari Nilsson, 13.10.2020

In order to construct the PSF, three tasks need to be set up correctly and run. In addition, five more tasks are "embedded" into these three, so the total number of tasks to set up is 8.

Some settings depend on the nature of the data (FWHM, crowdiness of the field, etc.), so do not just copy the values. Think what they mean.

The tasks are:

phot ( datapars, centerpars, fitskypars and photpars embedded
pstsel (datapars, daopars embedded)
psf ( datapars, daopars embedded ).
seepsf (converts the PSF image into a "normal" image)

To make the PSF photometry and subtract the PSF from the image you need :

group (a task relevant to crowded field photometry only, but needs to ve run anyway)
nstar (performs PSF fitting photmetry, also a preparatory step for substar)
substar (subtracts the PSF from stars).

But before going into the details, first a few words about choosing the PSF model

1. Choosing the model and PSF stars

There are two types of PSF models one can use :

analytical PSF
empirical PSF

Analytical PSF means that the PSF is represented by a mathematical function, like the Gaussian or the Moffat finction. The advantage is that the model is noise-free, which means it can be scaled freely without introducing additional noise. The drawback is at true PSF shapes are very complex and not very well represented by a an axialsymmetrical function (and all IRAF analytical PSFs fall into this category).

An empirical PSF is constructed from a star/stars in the image. In an ideal case, there is an unsaturated star close to the target and this star is at least a magnitude brighter that the target. The PSF star needs to be close, because the shape of the PSF changes over the field of view (FOV) due to the design of the system and misalignment of the optical elements.

An empirical PSF needs to be brighter than the target to avoid adding too much noise into the fitting processs. The empirical PSF will contain random noise, since it is derived directly from the image. When such a PSF is subtracted from the target, it is first scaled to the brightness of the target, and the noise is scaled with the same factor. If the PSF star is brighter than the target, this scaling factor is < 1, meaning that the noise is also reduced.

The noise can also be reduced by combining several stars into the PSF model. This is easy to in IRAF/DAOPHOT. However, choosing stars over the whole FOV means that the PSF may not represent the PSF at the target location very well due to aberrations.

There is really no single guideline how to select PSF stars, it all depends on the situatuion. Testing various options and choosing the best one is the best way to proceed. Here are some guidelines, though:

In the best case there is an unsaturated and isolated star close to the target, which is at least one magnitude brighter than the target. Derive an empirical PSF from that star.
Most likely the situation is not as ideal as above. Consider next using several stars, choosing them as close to the target as possible. The number of stars should be high enough to reduce the noise in the PSF image below the noise in the target image, preferably at least by a factor of ~ 2.
If there is no way to make an empirical PSF with reasonable S/N, you need to use an analytcal PSF. Use function moffat25 rather than a Gaussian, since the former usually represents the PSF better, especially further away from the center.
Using an analytical PSF will produce nasty residuals, since the real PSFs are not axially symmetric. Whether this is a problem depends on the circumstances. For instance, consider the task of fitting a host galaxy + AGN nucleus model to imaging data. If the host galaxy is bright and well resolved, the PSF quality has less impact compared to a situation where the host galaxy is barely resolved. In the latter case even tiny errors in the PSF will lead to a large error in the derived host galaxy parameters.
The author of these instructions has software to fit more complex, non-axiallysymmetric PSF models to stellar images. Basically the wavefront is modeled with 15 first Zernike polynomials and the focal plane image is derived. The polynomials are adjusted to find the best-fit model. This approach can produce very complex, non-axialsymmetric, PSF models. But mind you, it is not very easy to use, and convergence is sloooow.

It also happens often that the final emipirical PSF is not completely "clean", i.e. there are faint stars overlapping the PSF image. These must be edited away by e.g. replacing the affected areas with corresponfing areas on the opposite side of the PSF. Ask for details if you need this.

2. Startup

Start IRAF & ds9 in the usual manner. Remember to do the startup xterm or xgterm in order to enable the plots. When you have the IRAF cl prompt, type

  cl> noao
  cl> digiphot
  cl> daophot

3. phot

3.1 Setting up phot

  cl> epar phot

Set up the phot as follows. Note that * = parameter changed from default, otherwise keep the default. After editing, crtl-D exits.

*       image = "0905sum"       Input image(s)
*      coords = ""              Input coordinate list(s) (default: image.coo.?)
       output = "default"       Output photometry file(s) (default: image.mag.?
      skyfile = ""              Input sky value file(s)
    (plotfile = "")             Output plot metacode file
    (datapars = "")             Data dependent parameters
  (centerpars = "")             Centering parameters
  (fitskypars = "")             Sky fitting parameters
    (photpars = "")             Photometry parameters
*(interactive = yes)            Interactive mode?
    (radplots = no)             Plot the radial profiles?
   (icommands = "")             Image cursor: [x y wcs] key [cmd]
   (gcommands = "")             Graphics cursor: [x y wcs] key [cmd]
       (wcsin = )_.wcsin)       The input coordinate system (logical,tv,physica
      (wcsout = )_.wcsout)      The output coordinate system (logical,tv,physic
       (cache = )_.cache)       Cache the input image pixels in memory?
      (verify = )_.verify)      Verify critical phot parameters?
      (update = )_.update)      Update critical phot parameters?
     (verbose = )_.verbose)     Print phot messages?
    (graphics = )_.graphics)    Graphics device
     (display = )_.display)     Display device
        (mode = "ql")

It also

3.2 datapars


  cl> datapars
  
       (scale = 1.)             Image scale in units per pixel
*    (fwhmpsf = 4.)             FWHM of the PSF in scale units
    (emission = yes)            Features are positive?
       (sigma = 0.)             Standard deviation of background in counts
     (datamin = INDEF)          Minimum good data value
     (datamax = INDEF)          Maximum good data value
       (noise = "poisson")      Noise model
     (ccdread = "")             CCD readout noise image header keyword
        (gain = "")             CCD gain image header keyword
*  (readnoise = 3.5)            CCD readout noise in electrons
*      (epadu = 0.15)           Gain in electrons per count
    (exposure = "")             Exposure time image header keyword
     (airmass = "")             Airmass image header keyword
      (filter = "")             Filter image header keyword
     (obstime = "")             Time of observation image header keyword
       (itime = 1.)             Exposure time
    (xairmass = INDEF)          Airmass
     (ifilter = "INDEF")        Filter
       (otime = "INDEF")        Time of observation
        (mode = "ql")

Note that gain and readout noise are their effective values. They are used to predict the noise in a pixel, so if there are wrong, noise estimates and sometimes the PSF will be wrong.

3.3 centerpars

  
  cl> centerpars

* (calgorithm = "centroid")     Centering algorithm
*       (cbox = 7.)             Centering box width in scale units
  (cthreshold = 0.)             Centering threshold in sigma above background
  (minsnratio = 1.)             Minimum signal-to-noise ratio for centering alg
    (cmaxiter = 10)             Maximum iterations for centering algorithm
*   (maxshift = 5.)             Maximum center shift in scale units
       (clean = no)             Symmetry clean before centering
      (rclean = 1.)             Cleaning radius in scale units
       (rclip = 2.)             Clipping radius in scale units
      (kclean = 3.)             K-sigma rejection criterion in skysigma
    (mkcenter = no)             Mark the computed center
        (mode = "ql")

I usually set cbox = 2 times FWHM.

3.4 fitskypars

  cl> fitskypars

* (salgorithm = "mode")         Sky fitting algorithm
*    (annulus = 40.)            Inner radius of sky annulus in scale units
*   (dannulus = 10.)            Width of sky annulus in scale units
    (skyvalue = 0.)             User sky value
    (smaxiter = 10)             Maximum number of sky fitting iterations
     (sloclip = 0.)             Lower clipping factor in percent
     (shiclip = 0.)             Upper clipping factor in percent
    (snreject = 50)             Maximum number of sky fitting rejection iterati
   (sloreject = 3.)             Lower K-sigma rejection limit in sky sigma
   (shireject = 3.)             Upper K-sigma rejection limit in sky sigma
       (khist = 3.)             Half width of histogram in sky sigma
*    (binsize = 0.001)          Binsize of histogram in sky sigma
      (smooth = no)             Boxcar smooth the histogram
       (rgrow = 0.)             Region growing radius in scale units
       (mksky = no)             Mark sky annuli on the display
        (mode = "ql")

The background is measured from a ring centered on the target and with radius annulus and width dannulus. Obviously, then annulus has to be large enough so that you are measuring empty sky, not the wings of the target. It doesn't matter if there are a few faint objects the sky ring (the "mode" algorithm of IRAF takes care of the contamination to a large extent), but generally you want to measure the sky from a region devoid of any targets.

3.5 photpars

  cl> photpars
  
   (weighting = "constant")     Photometric weighting scheme
*  (apertures = "30")           List of aperture radii in scale units
        (zmag = 25.)            Zero point of magnitude scale
     (mkapert = no)             Draw apertures on the display
        (mode = "ql")

3.6 Running phot

  cl> phot

Iraf will ask you to confirm a bunch of parameters, just hit enter until the questions stop. You may get an error

    Warning: Graphics overlay not available for display device.

but you can ignore this.

Place the cursor over the star you want to measure (= a star you want to use for the PSF) an press space bar to measure it. Repeat for all stars. When yo are done, press q in DS9 and again q in IRAF window.

The result is a file called *.mag.?, where * is the image name and ? a running number increased automatically each time you run phot (i.e. phot doesn't overwrite files).

4. pstsel

4.1 daopars


  cl> daopars
  
*   (function = "moffat25")     Form of analytic component of psf model
    (varorder = 0)              Order of empirical component of psf model
      (nclean = 0)              Number of cleaning iterations for computing psf
   (saturated = no)             Use wings of saturated stars in psf model compu
    (matchrad = 3.)             Object matching radius in scale units
*     (psfrad = 40.)            Radius of psf model in scale units
*     (fitrad = 6.)             Fitting radius in scale units
    (recenter = yes)            Recenter stars during fit?
      (fitsky = no)             Recompute group sky value during fit?
    (groupsky = yes)            Use group rather than individual sky values?
    (sannulus = 0.)             Inner radius of sky fitting annulus in scale un
   (wsannulus = 11.)            Width of sky fitting annulus in scale units
     (flaterr = 0.75)           Flat field error in percent
     (proferr = 5.)             Profile error in percent
     (maxiter = 50)             Maximum number of fitting iterations
     (clipexp = 6)              Bad data clipping exponent
   (cliprange = 2.5)            Bad data clipping range in sigma
    (mergerad = INDEF)          Critical object merging radius in scale units
 (critsnratio = 1.)             Critical S/N ratio for group membership
    (maxnstar = 10000)          Maximum number of stars to fit
    (maxgroup = 60)             Maximum number of stars to fit per group
        (mode = "ql")

psfrad determines the size of the PSF. fitrad should be set to approx 1.5 times FWHM.

4.2 pstsel

  cl> epar pstsel

*       image = "0905sum"       Image for which to build psf star list
*    photfile = "default"       Photometry file (default: image.mag.?)
      pstfile = "default"       Output psf star list file (default: image.pst.?
      maxnpsf = 25              Maximum number of psf stars
     (mkstars = no)             Mark deleted and accepted psf stars
    (plotfile = "")             Output plot metacode file
    (datapars = "")             Data dependent parameters
     (daopars = "")             Psf fitting parameters
 (interactive = no)             Select psf stars interactively?
    (plottype = "mesh")         Default plot type (mesh|contour|radial)
   (icommands = "")             Image cursor: [x y wcs] key [cmd]
   (gcommands = "")             Graphics cursor: [x y wcs] key [cmd]
       (wcsin = )_.wcsin)       The input coordinate system (logical,tv,physica
      (wcsout = )_.wcsout)      The output coordinate system (logical,tv,physic
       (cache = )_.cache)       Cache the input image pixels in memory?
      (verify = )_.verify)      Verify critical pstselect parameters?
      (update = )_.update)      Update critical pstselect parameters?
     (verbose = )_.verbose)     Print pstselect messages?
    (graphics = )_.graphics)    Graphics device
     (display = )_.display)     Image display device
        (mode = "ql")

4.3 Running pstsel

  cl> pstsel

Again, IRAF asks to confirm a number of things, just hit enter. When the task has run, you have a file called *.pst.*

5. psf

5.1 setting up psf

  cl> epar psf  

*       image = "0905sum"       Input image(s) for which to build PSF
     photfile = "default"       Input photometry file(s) (default: image.mag.?)
*     pstfile = "default"       Input psf star list(s) (default: image.pst.?)
     psfimage = "default"       Output PSF image(s) (default: image.psf.?)
     opstfile = "default"       Output PSF star list(s) (default: image.pst.?)
    groupfile = "default"       Output PSF star group file(s) (default: image.p
    (plotfile = "")             Output plot metacode file
    (datapars = "")             Data dependent parameters
*    (daopars = "")             Psf fitting parameters
   (matchbyid = yes)            Match psf star list to photometry file(s) by id
*(interactive = no)             Compute the psf interactively?
     (mkstars = no)             Mark deleted and accepted psf stars?
   (showplots = no)             Show plots of PSF stars?
    (plottype = "mesh")         Default plot type (mesh|contour|radial)
   (icommands = "")             Image cursor: [x y wcs] key [cmd]
   (gcommands = "")             Graphics cursor: [x y wcs] key [cmd]
       (wcsin = )_.wcsin)       The input coordinate system (logical,tv,physica
      (wcsout = )_.wcsout)      The output coordinate system (logical,tv,physic
       (cache = )_.cache)       Cache the input image pixels in memory?
      (verify = )_.verify)      Verify critical psf parameters?
      (update = )_.update)      Update critical psf parameters?
     (verbose = )_.verbose)     Print psf messages?
    (graphics = )_.graphics)    Graphics device
     (display = )_.display)     Display device
        (mode = "ql")

5.2 Running psf

  cl> psf

Confirm again multiple times. When finished, you will have the PSF image *.psf.?.fits.

Displaying this image will show something which doesn't quite look like a star. This is due to the way IRAF stores the PSF. The analytical part of the PSF is stored in the PSF image header. The residuals wrt the analytical part are stored as pixel data and shown when you display the image.

Note also that the PSF image is supersampled with 2 times higher resolution than the original image. IRAF can use this format, but any external program cannot necessarily do it. You can use seepsf to convert the "funny" IRAF format PSF into a normal stellar image.