The QUEST Data Processing Software Pipeline

Prior to starting the PHOTOMETRY program, each night's drift-scan data are separated into smaller files on disk. Each column of data, typically spanning ∼120° in R.A., is saved on disk as separate images (frames) of 640 × 2400 pixels each. This is the pixel size of the CCD arrays (plus 40 overscan pixels) and corresponds to 0.15° × 0.58° on the sky. A typical night's data consist of 28 columns, with four passbands per column and 200 drift-scan frames per passband. Each column is processed independently. Within each column, the frames of each passband are also processed independently in the initial stages (preprocessing and object detection). To get CCD coordinate transformations and astrometric solutions, scan-wide solutions are obtained for each passband. For seeing and photometry, the frames are once again processed independently, but information from the CCD coordinate and astrometric transformations is used to tie the photometry across passbands.

The first stage of the PHOTOMETRY program is to bias-subtract, dark-subtract, and flat-field the data on a frame-by-frame basis for each of the 112 CCDs in the camera. For each row of each frame, the first correction is to subtract the overall bias level of the CCD readout amplifier. This level is determined by median-averaging the overscan pixels in each row of the frame. Doing this on a row-by-row basis corrects for occasional low-level (a few counts) intermittent changes in the bias level that sometimes occur during image readout. The second correction is to subtract the fixed, column-dependent deviations in the dark level that are caused by nonnegligible dark currents and the presence of bad pixels in some of the CCDs (bad pixels appear as bad columns in drift scans). This column-dependent correction (the dark field) is determined from a dark drift scan. Finally, variations in the response of the CCD are corrected by dividing by a flat field that is derived from a twilight drift scan. The dark and twilight flat fields required for these corrections are taken nightly, but we create only a single set of correction fields per lunation. For this, we use the darks and flats taken on the first dark night (Johnson filters) and the first gray night (SDSS filters).

To make the dark correction field, a single raw dark frame is first corrected by subtracting the overscan bias level. The dark field is then divided into 240 subframes of 10 rows each. These subframes are median-averaged together, pixel by pixel, to yield a single average subframe. This subframe is then reduced to a single-row image by computing the average signal in each column. This final result is the dark correction field that is subtracted from each row of the drift scan to be corrected.

To make the flat correction field, there is an initial selection of the appropriate twilight frames to use for the correction. Ideally the sky level in counts should be above 5,000 so that the signal-to-noise ratio is high, yet lower than 10,000 to avoid levels that are close to the full well of some of the CCDs. During twilight, we take a long drift scan (∼20 frames), starting when the sky is bright and ending when it is dark to ensure that at least three consecutive frames are taken in each filter that meet this criterion. These chosen frames are then processed, like the darks, to remove the overscan bias. Then the dark correction, previously determined from the dark field, is subtracted from each row of each frame to remove column-dependent biases and dark currents. Next a one-dimensional, row-dependent fit to the sky level is determined for each frame. Dividing each frame by this fit removes the sky background that varies monotonically with time during the twilight scan. The three frames are then median-averaged together. This average frame is further divided into 240 subframes, 10 rows each, which are median-averaged together to yield a single subframe. A final averaging of the subframe rows together and normalizing of the average level to unity yields the flat-field correction.

To check that no mistakes have been made in the creation of the dark and flat fields (owing, for example, to processing errors, lights left on in the dome, scattered moon light, or electronics problems) we compare newly created darks and flats to those from the previous lunation. Fractional differences smaller than a few percent are acceptable. If large deviations are found, however, a new set of darks and flats are created using dark and twilight scans from a subsequent night. Occasionally, we have used night-sky drift scans for flat fields when we have not been able to obtain twilight flats (owing to weather or equipment problems). The twilight and night-sky scans yield equally valid flat fields if the night-sky background is high enough (this depends on the passband). For drift scans, fringing is not an issue because fringe effects are averaged out in the scanning processes.

Object detection is carried out on one frame at a time, separately for each of the four filters. It is a two-pass procedure. The first pass is a simple, threshold-based routine. In the second pass, these initial detections are refined by PSF fitting. This fitting procedure more carefully separates closely spaced stars (separation ) which would otherwise be counted as single objects. Since the sky level across a whole frame can vary (owing to twilight, moon rise, extinction changes, etc.), each frame is subdivided into thousands of small regions and the local-sky level for each region is calculated as a two-stage clipped average of the pixel values.

In the first pass, the detection algorithm steps through the pixels in an image in a raster fashion, assembling objects with pixels where the flux exceeds 2.5 standard deviations above the local sky background. An object is completed when it is surrounded on all sides by below-threshold pixels. Objects are kept only if they consist of at least four adjacent pixels above threshold. This criterion voids the detection of single-pixel noise artifacts such as cosmic ray hits. Because the CCDs are thinned to ∼15 μm thickness and because exposure times for drift scans are only a few minutes, cosmic ray hits exceeding a few pixels in length rarely appear. At this stage of the pipeline, detections in one filter are not tied to those in another (this occurs later; see § 3.7).

During the raster scan of each row of a given CCD, a map of bad columns for that CCD is consulted as each pixel is evaluated. When an object is found that overlaps or abuts a bad column, the boundary for that object is terminated at the bad column. No further pixels from the row are associated with the object. There is no attempt to smooth over the bad column. At this state in the pipeline, there is no special treatment for objects abutting bad columns or any edges of the frame. There is also no effort made to merge the objects at the start or end row of a frame with the previous or following frame in the scan sequence. Later, when fluxes are measured, flags are set to indicate which objects are at frame boundaries or connect to bad columns, and there are additional checks to remove bad detections.

In the second pass of the detection routine, the first step is to obtain an accurate model of the PSF for each frame (this accounts for changes in the PSF over time owing to variable seeing conditions). The PSF is viewed as a two-dimensional distribution of light impacting an imaging device, which the atmosphere and the instrument produce from an idealized point source such as a distant star. We do not attempt to describe this as an analytic function but instead use the actual measurements of stars on the image that have typical shapes and good statistics. The two-dimensional flux distribution varies slightly depending on the location within a pixel where the object is centered; i.e., the PSF is slightly different for a star centered at the center of a pixel than for a star centered near the corner of a pixel. To take this into account, the central pixel is divided into 4 × 4 segments. Selected good-quality stars are sorted into 16 classes, each class consisting of images centered on one of the 4 × 4 segments of the central pixel. For each class the pixel fluxes of many images are summed and then normalized to a total unit flux. This produces 16 different PSF tables. In the PSF fitting of a particular object, the segment that holds the central pixel of the object is first evaluated. This determines the appropriate PSF table to use.

The PSF fitting procedure for each object is an iterative one. The fluxes in the individual pixels of the object (determined in the first pass) are examined. If there is a simple maximum the object is fit to a single PSF. If there is more than one local maximum the object is fit to a number of separate PSFs corresponding to the number of local maxima. Both the fluxes and positions of the PSFs are simultaneously fit to each object. Then the fluxes expected from the fitted PSFs are subtracted from the object image. If there is a significant positive residual, the residual flux is considered as a separate object and a new PSF is added centered on the residual flux. This allows the detection algorithm to detect multiple images or faint stars previously obscured by a close bright star. The procedure is then repeated for a maximum number of iterations specified by the user, typically three fitting iterations in total. In each iteration, the best fit is determined by a search for the χ² minimum. The search is implemented using Powell's multidimensional methods with Brent's inverse parabolic algorithm for each successive line minimization (Press et al. 1986)

To obtain the transformation from CCD pixel coordinates to astrometric positions, objects detected in each filter are matched to stars in the USNO A-2.0 catalog (Monet 1988). All the objects detected with a given CCD during an entire drift scan are pooled together into a single catalog so that a single average transformation, linear in pixel coordinates, can be found for the entire scan. Within a given catalog, the pixel coordinate corresponding to an object's image row number is indexed cumulatively from the starting row of the drift-scan image. A linear solution is possible in the decl. direction because each CCD covers only a small angle in declination. The solution is necessarily linear in the R.A. direction because any nonlinear distortions are averaged in the scanning direction. Once the average solution is found, residual deviations from the average solution (owing to small clocking errors, any small telescope motions during the scan, and the effect of line drops discussed above) are tracked as a function of right ascension and used as second-order corrections to the average solution.

To derive an astrometric transformation, the first step is to make a crude estimate of the pixel positions of the USNO catalog stars, given the known telescope pointing, the known pixel scale of the drift-scan images, and the catalog positions. A triangle-matching algorithm is then used with selected bright objects to associate the catalog stars to the detected stars, and then the linear transformation is computed. We use the matching algorithm distributed by NOAO with their IRAF image processing software (see FOCAS routine at http://iraf.noao.edu). This is described in detail by Valdes et al. (1995). In our pipeline, the matching is done one frame at a time, each time using the solution from the previous frame to predict the pixel positions of the catalog stars in the next frame. Here we choose a frame size of 2,400 rows because it is convenient and usually yields a sufficient number of stars to yield a good solution, but the choice is not critical. Once all frames have been solved, the average solution is determined for the whole scan. If there are any frames that do not yield good solutions, they are skipped and the average solution is assumed to hold for these frames.

It should be noted that the USNO catalog is in J2000 coordinates, which is not the natural coordinate system that is observed on any given night. The USNO coordinates are precessed to the epoch of the observations before the matching is attempted. Once the final transformations are obtained, the astrometric coordinates of the detected objects are precessed back to J2000 coordinates. The SLALIB software library (Wallace 1994) is used to make these transformations.

To judge the precision of the astrometry of this procedure we compared the right ascension (R.A.) and declination (decl.) from the QUEST survey with the astrometric standards of Stone et al. (1999). We found that the mean values of USNO and the Stone standards disagree by about 0.2''. To bring our data into agreement with the Stone standards we applied a final calibration correction of typically -0.1^'' in R.A. and -0.2^'' in decl. After this final correction (which was applied in the ANALYSIS program described below) the comparison with the astrometric standards is shown in Figure 3. We see from this comparison that the astrometric precision of the drift-scan data is 0.1''. However, over the entire area of our survey, the systematic error is limited by that of the USNO catalog that we use for calibration (∼0.2^''; see Assafin et al. 2001).

Zoom In Zoom Out Reset image size

Bright objects should be present in all four passbands in each of the 28 columns. Fainter objects or those with extremely blue or red color may be bright enough for detection in some of the filters but below the threshold for detection in others. So that we can measure the flux or an upper limit to the flux in the passband where a detected object is below the detection threshold, and so that detections of a given object in multiple passbands can be identified as the same object, we calculate linear transformations between CCDs that relate their pixel coordinates. As in the astrometric solutions (see § 3.4), all the detections by a given CCD in a given drift scan are pooled into a single catalog, ignoring frame boundaries, and pixel rows are indexed from the starting row of the drift scan.

The process starts with the selection of a lead CCD. Although the choice is not critical, we pick the CCD for which the corresponding passband is most likely to yield detections of any given object (i.e., the passband with the faintest magnitude limit). For the Johnson filters it is the R filter, and for the SDSS filters, it is the r^' filter. Next, given the astrometric solutions previously determined for each passband (see § 3.4), the objects detected in each of the nonlead CCDs are matched to the objects from the lead CCD. Once matched, a Powell minimization algorithm (Press et al. 1986) is used to determine the linear equations transforming the nonlead pixel coordinates to the lead-CCD pixel coordinates. As in the case of the astrometric solution described in § 3.4, a linear transformation is expected, but we track the residuals to the solution as a function of R.A. and use these as second-order corrections to the transformation. The accuracy of the final CCD-to-CCD transformations has been measured to be 0.08 pixels in both the x and y directions on the CCDs. Note that this is more accurate than the astrometric solutions discussed above (see § 3.4) because no transformation to a catalog is involved, other than to match the stars on two different CCDs. The final transformations are used to match the objects detected in the four different passbands to each other and to register the apertures and PSF fits used for photometry (see §§ 3.7 and 3.8).

This routine takes the list of detected bright objects and calculates the average seeing for each frame in each filter by performing a clipped mean of the measured FWHM values (determined from the PSF fitting). These frame-wise computed values of FWHM, as well as the number of objects per frame, are stored for later use in the aperture photometry routine.

Before making any photometric measurements, the first step is to merge the detection list of the different passbands. Typically the CCDs with different filters have widely varying fluxes for a given object. For example, in the Johnson set, R and I usually have the highest flux, B less, and U the least. Often objects are detected in the reddest filters and not detected in the U filter. To obtain a flux measurement in a passband where the object is too faint for detection but its position is known from its detection in another passband, the list of objects detected in all four passbands in a scan are merged. Using the CCD-to-CCD coordinate transformations (see § 3.5), the pixel coordinates of these merged objects are transformed into the coordinate system of each separate CCD. For each object, a flux is measured using an aperture centered at the transformed location. Even if there is no detectable signal for a given passband, a measurement is made. In cases where an object is detected in multiple passbands including the lead CCD, the position is taken to be the one determined from the lead CCD, as this is likely to give the best estimate of the centroid.

Aperture photometry is performed in a standard way with four different apertures: one with a one-FWHM-radius circle aperture, one with a two-FWHM radius, and two with fixed apertures with 3- and 10-pixel radii. The one-FWHM aperture is optimal for photometry, as the contribution of the sky noise is minimal. Fluxes are measured in the other apertures to allow for later discrimination of extended or blended objects from point sources (see § 5). The choices for these larger apertures were not tied to any specific survey goal, but are in the size range useful for characterizing the shapes of extended sources. For each aperture the flux is taken to be the total flux in the aperture minus the sky flux (the number of pixels times the average sky background per pixel). The average sky background is obtained as the clipped mean sky level in an annulus (typically of inside and outside radii of 20 and 25 pixels respectively) around the object.

The error on the flux, σ_F, is calculated by adding in quadrature the errors in the object flux within the aperture, the error in the sky flux within the aperture, and the error in the estimate of the sky background:

where F is the object flux in counts (ADC units), G is the gain in counts per electron, N_obj is the number of pixels inside the aperture, N_sky is the number of pixels in the sky annulus, and σ_sky is the sky noise. We do not separately account for read noise because the measured sky noise already includes the contribution from read noise. Gains are checked every few months or any time adjustments have been made to the electronics. The flux and its error are then converted to the instrumental magnitude, M_I, and its error, σ_M, using the following formulas:

For the one-FWHM aperture, the ellipticity, orientation, and size of the minor and major axes of the image, the FWHM, the sky level, the maximum pixel value and the number of pixels above threshold are recorded. Error flags are set to indicate problems such as saturated pixels, badly measured sky background, and negative flux. There is a nonlinear response of the CCDs, which only becomes significant for U- and B-band observations on dark nights when the sky background is very low, but it is corrected for in the calibration of the instrumental magnitudes (see § 6.2).

The QUEST pipeline measures PSF as well as aperture photometry for all objects. As in the aperture photometry, the list of object positions is provided by the detection routine and is the merged list from all passbands. Much of the PSF photometry code is very similar to that used in the PSF detection routines (second pass of object detection). The main difference is that object positions are not allowed to vary in the PSF fitting. The positions determined at detection are frozen. Only fluxes are allowed to vary. As in the detection stage, a set of empirical PSFs is constructed and fitted to each object. Objects within 13 pixels of each other are gathered into a group, and their fluxes are simultaneously fit using the appropriate PSFs. A 13-pixel separation is chosen because it is much larger than the FWHM, ensuring that any one group will not contain flux from a neighboring group. Simultaneous fitting of all the group members ensures proper handling of objects whose PSFs overlap. The flux errors are obtained by evaluating the flux range, for a given object, for which the resulting reduced χ² varies by unity. This PSF routine is quite robust. It fails to give a meaningful flux <1% of the time.

The errors obtained from the χ² fits are essentially statistical errors. To get a measure of systematic errors in this procedure, we compared the results of repeated measurements of the same object in the same filter carried out on different nights. It was found that a systematic error of 0.015 mag had to be added in quadrature to the statistical errors to bring the errors into agreement with the rms of the repeated measurements. No specific reason for this systematic error could be found (or else we would have corrected the problem). Likely causes are variable extinctions from night to night and frame to frame that are not completely taken out by the extinction corrections, and variations in the PSF with position in each field.

If all nights were perfectly photometric, and all observations were made at the same air mass, we would not need this correction. On nights with a significant cloud cover or other adverse weather conditions, the dome does not open and no observations are made. However, even on nights when data are recorded, there is often a significant amount of extinction. In this survey we also use the nonphotometric nights with extinctions <0.5 mag. We thus have to make a correction for the extinction. The only large area survey that covers all of the area of the Palomar QUEST survey is the USNO survey (Monet 1988). We started out using the USNO catalog for the extinction correction but found that the photometry of the USNO catalog was too variable to be useful for this purpose (exceeding 10%). We therefore adopted the following procedure using only our own survey data, making use of the fact that we covered each area of sky in the survey multiple times over different nights.

The procedure starts by taking each drift-scan strip (recall that these are 4.6° wide in declination) and dividing it into 1° segments in R.A. (about 4 minutes of scan time). In each passband all the observations in a given R.A. segment are collected and the brightness of a sample of stars is compared. We then assume that the observation with the highest flux for the selected stars is a good approximation to photometric and assign an extinction correction to all of the other observations to bring them up to the level of the brightest night. No other requirements are made for the brightest night.

The distribution of the extinction corrections (the number of 1° segments vs. their extinction, as defined above) is shown in Figure 4. We find that ∼80% of the data are extincted by <0.2 mag. In other words, 80% of the data yield photometry consistent to ∼20% under repeated observation. If the extinction correction is larger than 0.5 mag, the data are flagged as bad by the ANALYSIS program and not used for science data analysis. About 10% of the data are lost due to this cut. These extinction corrections are fed into the ANALYSIS program described above and are used to correct the instrumental magnitudes. The precision of these corrections is discussed in § 6.3 below.

Zoom In Zoom Out Reset image size

The Sloan Digital Sky Survey (SDSS) Data Release 4 (Adelman-McCarthy et al. 2006) is used to do the photometric calibration of the QUEST survey, using areas of sky where the two surveys overlap. Several clear nights of QUEST data, both with Johnson and with SDSS filters, are selected (these are nights with the least extinction based on comparison with at least four other nights covering the same area). Objects in QUEST are matched to those in SDSS (using objects classified as stars by SDSS) and the QUEST instrumental magnitudes (after being corrected for extinction as described above) are compared to the SDSS-quoted magnitudes. The calibration consists of calculating the zero-point offset, i.e., the magnitude that has to be added to the QUEST instrumental magnitudes to bring them into agreement with SDSS. This is straightforward for the SDSS filters because those are the filters used by SDSS. To calibrate the data with Johnson filters, an algorithm (Fukugita et al. 1996) was used to convert the magnitudes of the SDSS catalog to Johnson magnitudes with sufficient precision to compare to QUEST data with Johnson filters.

This calibration procedure was carried out for all 112 CCDs, separately with the Johnson and SDSS filter sets, so there were a total of 224 calibrations (for each CCD, there are only two filter possibilities). Furthermore, the calibration of a single CCD was allowed to vary as a function of the column across the CCD, the magnitude of the object, and the sky brightness. The full calibration of the camera thus consisted of 112 CCDs  × 2 filters  × 10 regions across a CCD  × 6 magni-tude bins × 6 sky brightness bins = 80,640 calibration parameters. As discussed in Baltay et al. (2007), several of the CCDs in the full array are not functional, and so the number of calibration parameters is slightly smaller. Thanks to the substantial overlap area with the SDSS survey we could use ∼10⁶ stars to evaluate all of these parameters. These parameters were fed into the ANALYSIS Program described above and were used to transform the instrumental magnitudes (already corrected for extinction) to obtain the final calibrated magnitudes.

Extensive checks were carried out to verify the quality of the photometric calibration. In one of the most stringent tests, we looked at the absolute photometric magnitude calibration by comparing with the SDSS data using a 1500 deg² region where the two surveys overlap. We compared each of the eight filters separately. The comparison for the SDSS r^' filter is shown in Figure 5. We use in this comparison the entire QUEST data set from September 2003 to September 2006, with the standard quality cuts applied in the ANALYSIS program (see § 7), including both photometric and non-photometric nights. The comparison, Figure 5a, has a standard deviation σ = 0.11 mag. This includes both the SDSS and QUEST statistical errors as well as the calibration errors.

Zoom In Zoom Out Reset image size

To separate the statistical and calibration errors, we plot the distribution of the magnitude difference between SDSS and QUEST divided by the total error, where the total error is the SDSS and QUEST errors on each object added in quadrature to the estimated calibration error. We vary the calibration error estimate and choose the one that gives σ = 1 for this distribution. For the r^' filter, the best estimate of the overall systematic calibration error was 0.08 mag, as shown in Figure 5b. The results for the estimated calibration error for all of the other passbands including both photometric and non-photometric nights vary between 0.06 and 0.12 mag, except for the U and B filters where the errors are larger. There is also typically 15% of the data in a non-Gaussian tail.

This rough overall photometric calibration error is for the entire survey, including non-photometric nights and with no quality cuts. When we restrict the comparison to photometric nights only, the calibration errors are considerably smaller. Furthermore this procedure uses a single set of magnitude zero points for all nights, with the relative extinction corrections as described in § 6.1. No attempt was made at this point to determine color terms in the corrections, though the ANALYSIS program has procedures to do this. In any particular science analysis where better photometric precision is required, calibration procedures appropriate for that particular analysis will have to be used. For the general goals outlined in § 1, the precision indicated by Figure 5 is adequate.

For example, to check the systematic errors of the relative magnitude calibration we took frames with repeated measurements of the same area of the sky spread out over the three years of the survey. We normalized the different epochs to each other using a sample of 10 or more field stars, and plotted the distribution of the deviation of the magnitudes of the individual measurements from the average magnitude. An example plot for the Johnson R filter, using data for bright stars with six or more measurements at different epochs, is shown in Figure 6a. The standard deviation of this distribution is σ = 0.03 mag. This error includes both the statistical and the systematic components of the error. To separate the two we used a procedure similar to that described above. We divided each deviation from the average by the error, where the error was taken as both the statistical and estimated systematic errors added in quadrature, and varied the estimated error. We needed a systematic error of 0.02 mag to produce a distribution with σ = 1 as shown in Figure 6b. Thus, our best estimate for the systematic error on relative photometry is 2%. It is also gratifying that a single Gaussian with σ = 1 is a good fit to the data with very little non-Gaussian tail. By careful selection and further calibration of the data, with cuts to eliminate the lower-quality observations, we obtain a systematic error of 0.7% (Bauer 2008).

Zoom In Zoom Out Reset image size