.help apvariance Aug90 noao.twodspec.apextract

.ce
Variance Weighted and Cleaned Extractions


There are two types of aperture extraction (estimating the background
subtracted flux across a fixed width aperture at each image line or
column) in the APEXTRACT package.  One is a simple sum of pixel values
across an aperture.  It is selected by specifying "none" for the
\fIweights\fR parameter.  The second type weights each pixel in the sum
by it's estimated variance based on a spectrum model and detector noise
parameters.  This type of extraction is selected by specifying
"variance" for the weighting parameter.  These two extractions are
defined by the following equations.

.nf
	none:	S = sum { I - B }
    variance:	S = sum { (P**2 / V) (I - B) / P } / sum { P**2 / V }
.fi

S is the one dimensional spectrum flux at a particular wavelength (line
or column along the dispersion axis).  The sum is over all pixels at
that wavelength within the aperture limits.  If the aperture endpoints
occupy only a fraction of a pixel then the pixel value above the
background is multiplied by the fraction.  I is the pixel value and B
is the estimated background at that pixel (see \fBapbackground\fR), P
is estimated normalized profile value for that pixel (see
\fBapprofile\fR), and V is the estimated variance of the pixel based on
the noise model described below.  Note that the quantity (I-B)/P is an
independent estimate of the total flux from one pixel since the
integral of P is one and it is these estimates that are variance
weighted.

Variance weighting is often called "optimal" extraction since it
produces the best unbiased signal-to-noise estimate of the flux in the
two dimensional profile.  The theory and application of this type of
weighting has been described in several papers.  The ones which were
closely examined and used as a model for the algorithms in this
software are "An Optimal Extraction Algorithm for CCD Spectroscopy",
PASP 98, 609, 1986, by Keith Horne and "The Extraction of Highly
Distorted Spectra", PASP 100, 1032, 1989, by Tom Marsh.

The noise model for the image data used in the variance weighting,
cleaning, and profile fitting consists of a constant gaussian noise and
a photon count dependent poisson noise.  The signal is related to the
number of photons detected in a pixel by a \fRgain\fR parameter given
as the number of photons per data number.  The gaussian noise is given
by a \fIreadnoise\fR parameter which is a defined as a sigma in
photons.  The poisson noise is approximated as gaussian with sigma
given by the number of photons.

Some additional effects which should be considered in principle, and
which are possibly important in practice, are that the variance
estimate should be based on the actual number of photons detected before
correction for pixel sensitivity; i.e. before flat field correction.
Furthermore the uncertainty in the flat field should also be included
in the weighting.  However, the profile must be determined free of
sensitivity effects including rapid larger scale variations such as
fringing.  Thus, ideally one should input the unflat-fielded
observation and the flat field data and carry out the extractions with
the above points in mind.  However, due to the complexity often
involved in basic CCD reductions and special steps required for
producing spectroscopic flat fields this level of sophistication is not
provided by the current package.  The package does provide, however,
for propagation of an approximate uncertainty in the background
estimate when using background subtraction.

The noise model is described by the following equations.

.nf
    (1) V = max (VMIN, (R**2 + I + VB) / G**2)
	    max (VMIN, (R**2 + S * P + B + VB) / G**2)

    (2) VB = 0.                 if (B = 0)
	   = B / (N - 1)        if (B > 0)

    (3) VMIN = 1 / G**2         if (R = 0)
	       R**2 / G**2      if (R > 0)
.fi

V is the desired variance of a pixel to use for variance weighting.  R
is the photon read out noise specified by the parameter \fIreadnoise\fR
and G is the photon per data value gain specified by the parameter
\fIgain\fR.  There are two forms to (1).  The first is used in the
initial pass of estimating the spectrum flux S and the actual pixel
value I (which includes any background) is used for the poisson term.
The other form is used in a second pass (and further passes if
cleaning) using the estimated data value based on the normalized
profile P scaled to the estimated total flux plus the estimated
background B; i.e. I estimated = S * P + B.

The background variance VB is computed using the poisson noise model
based on the estimated background counts.  If no background subtraction
is done then both B and VB are set to zero.  If a background is
determined the background is either an average or function fit to
pixels in defined background regions.  If a fit is used B need not be a
constant.   Because the background estimate is based on a finite number of
pixels, the poisson variance estimate is divided by the number N (minus
one) of pixels used in determining the background.  The number of
pixels used includes any box car smoothing.  Thus, the larger the
number of background pixels the smaller the background noise
contribution to the variance weighting.  This method is only
approximate since no correction is made for the number of degrees of
freedom and correlations when using the fitting method of background
estimation.

VMIN is a minimum variance need to avoid generating zero or negative
variances from the data.  The definition of VMIN is such that if a zero
read out noise is specified (which is certainly possible such as with
photon counting detectors) then a minimum of 1 photon is imposed.
Otherwise the minimum is set by the read out noise even if the poisson
count part is (unphysically) negative.

One deviation from the linear photon response mode which is considered
is saturation.   A data level specified by the parameter
\fIsaturation\fR is used to exclude data from the profile fitting.
During extraction the saturated pixels are not treated any differently
than unsaturated pixels except that dispersion points with saturated
pixels are flagged by reversing the sign of the final estimated sigma;
the sigma output is enabled with the \fIextras\fR parameter.  Exclusion
of saturated pixels from the extraction, as is done with deviant
pixels, was tried but this resulted in higher noise in the spectrum.

If removal of cosmic rays and other deviant pixels is desired, called
cleaning and selected with a \fIclean\fR parameter, they are
iteratively rejected based on the estimated variance and excluded from
the weighted sum.  Note that a cleaned extraction is always variance
weighted regardless of the value of the \fIweights\fR parameter.  This
makes sense since the detector noise parameters must be specified and
the spectrum profile computed, so all of the computational effort must
be done anyway, and the variance weighting is as good or superior to a
simple unweighted extraction.

The detection and removal of deviant pixels is straightforward.  Based
on the noise model described earlier, pixels deviating by more than a
specified number of sigma (square root of the variance) above or below
the model are removed from the weighted sum.  A new spectrum estimate
is made and the rejection is repeated.  The rejections are made one at
a time starting with the most deviant and up to half the pixels in the
aperture may be rejected.  The total number of rejected pixels in the
spectrum is recorded in the logfile and a profile plot of data and
model profile is recorded in the plotfile.

As a final step when computing a weighted/cleaned spectrum the total
fluxes from the weighted spectrum and the simple unweighted spectrum
(excluding any deviant and saturated pixels) are computed and a
"bias" factor of the ratio of the two fluxes is multiplied into
the weighted spectrum and the sigma estimate.  This makes the total
fluxes the same.  The bias factor is recorded in the logfile
if one is kept.  Also a check is made for unusual bias factors.
If the two fluxes disagree by more than a factor of two a warning
is given on the standard output and the logfile with the individual
total fluxes as well as the bias factor.  If the bias factor is
negative a warning is also given and no bias factor is applied.
.ih
SEE ALSO
apbackground approfiles apall apsum
.endhelp
