Reducing Palomar 200-inch LFC Data:
Part2 - Photometry

Roy Gal, 30 June 2004

For some reason the EQUINOX header keyword is not set properly in the image headers of mosaicked LFC images. This is needed for SExtractor to output correct coordinates. Just do:
cl> hedit sc1604_i.fits EQUINOX 2000.0 add+ ver-
and similarly for any images you will be running through SExtractor.

The first thing we need to do is obtain the photometric zero points, color terms, and airmass terms to calibrate our data. The equation we will use to calibrate our objects is:
mcalib = minst + ZP + B*airmass + C*color
Here, ZP is the zero point, B is the coefficient for the airmass term (where everything is corrected to airmass=1), and C is the coefficient for the color term. The color terms we use are (r'-i') for the r' and i' filters, and (i'-z') for the z' filter.

I perform my calibrations and photometric measurements using data in counts/sec!! This means that all standards AND all target data will need to be normalized to 1 second exposures !!

Photometering the Standards

We already normalized the standards, so we are ready to run SExtractor on them. I won't go into the details of SExtractor here, but you SHOULD NOT USE IT BLINDLY. Read the official documentation as well as SExtractor For Dummies.

To run SExtractor, we need a configuration file and an output parameter file. Examples of each are provided at the links below:
Configuration file named sex.config.stds
Parameter file named sex.param.stds

If you download these into the directory where you will be running SExtractor, that's all you need. We create a script to run SExtractor on all the standard star images:

% awk '{print "sex -c sex.config.stds -CATALOG_NAME mos"$1".cat -WEIGHT_IMAGE weight_"$1".fits -CHECKIMAGE_NAME mos_"$1"_ap.fits mos_"$1"_norm.fits"}' stdnums.list > gosex.stds
% chmod +x gosext.stds
% ./gosex.stds
I am going to use the IRAF task fitparams to fit the photometric calibration. For this, we will need a few input files:
  1. The observation files, listing each standard, the filter, airmass, instrumental magnitude, and magnitude error. We need one file for each filter pair we have a color term for. For instance, if we have r',i',z' observations, we need TWO observation files - one with r',i', and one with i',z'. The format of this file is:
    sa104_428       r  1.571 15.860 0.02
    *               i  1.523 15.778 0.02
    sa104_428       r  1.458 15.785 0.02
    * i  1.506 15.784 0.02
    sa104_428       r  1.236 15.778 0.02
    * i  1.214 15.719 0.02
    sa104_428       r  1.241 15.760 0.02
    * i  1.217 15.728 0.02
    pg1047+003a     r  1.202 16.768 0.02
    * i  1.197 16.762 0.02
    pg1047+003a     r  1.200 16.722 0.02
    * i  1.198 16.758 0.02
    pg1047+003a     r  1.198 16.734 0.02
    * i  1.197 16.762 0.02
    There is one entry for each pair of observations. The first column is the object name, which must EXACTLY match the entries in the catalog file, described below. The second column is the filter name, followed by the airmass, the instrumental magnitude, and the magnitude error. The airmass can be found in each image's header using
    cl> imhead mos_703.fits lo+ | grep AIRMASS
    AIRMASS =                1.571  / Airmass at start
    The magnitudes are MAG_BEST and MAGERR_BEST from the SExtractor output file. To find your object in this file, look in your image, find the approximate x,y of your standard(s), and find the correspond entry. Note that the entries are sortedroughly in ascending y coordinate. The observation file should be kept in your local working directory.
  2. A catalog file, with the known magnitudes of the objects. These also must contain filter pairs. You can dowmload my catalog files for ri (sdss_ri.stds) and iz (sdss_iz.stds). Note that these contain just a few stars. You can edit these files to include more standards from the SDSS Standards paper (Smith et al. 2002).

    I usually place this file in a pers subdirectory in my iraf directory (in my case /Users/gal/Work/iraf/pers/); I add the line

    set     pers            = "home$pers/"
    to my (or as well.
  3. A configuration file, which tells fitparams what all the columns in the other files mean, and what equation to fit. Once again, you need one config file for each filter pair, and place them in your pers directory. Example files ri (ristd.config) and iz (izstd.config) can be downloaded.

Once you have the observations, catalog, and configuration files in place, you can run (in IRAF):

cl> noao
cl> digiphot
cl> photcal
cl> fitparams obs=ri.stdobs catalog=pers$sdss_ri.stds con=pers$ristd.config par=ri.coefs
cl> fitparams obs=iz.stdobs catalog=pers$sdss_iz.stds con=pers$izstd.config par=iz.coefs
You will get a plot like this showing the standard magnitude on the x-axis and the difference between this and the fit on the y-axis. The top of the screen will also show you the number of points, number of rejected points, and the RMS. Hitting "k" in this window will switch to airmass vs. residual; this is a very useful plot to delete outliers. If it looks like one observation is way out of line, put your cursor over it and hit "d". Then hit "f" to refit! If you are happy with the fit, hit "q" to quit. You are then asked if you want to save the results, say "y", and then "next". You now fit the second magnitude of the pair. When you are happy with this fit, hit "q", "y", "n", "y". The results are in the file ri.coefs (or iz.coefs). Click here for an example results file. Note that there are separate sections for the two different filters you are fitting! We use the R and I sections from the ri.coefs file, and the Z section from the iz.coefs file. The lines of interest are the three numbers below the line
        values  3
and will have numbers like
The three numbers are the zero point ZP (relative to 30.0), the airmass coefficient B, and the color coefficient C. Write these down somewhere!

Photometering Target Observations

Now that we know are calibration coefficients, we are ready to set up SExtractor for our deep combined target observations. We will be using SExtractor in dual-image mode, where the multi-filter ultradeep observation is used for object detection, and the individual filter images are used for measurement.

This scheme proves to be a little tricky. The problem is that the image sizes of the ultradeep image and the individual filter images are NOT the same, and there are integer pixel offsets between them. So, what we need to do is create versions of the individual filter images that are the same size as the ultradeep image, with objects at the same coordinates. This requires care, so do this carefully and check your work!

The first thing to do is get the image size of the ultradeep image. This will be the largest of all the images since the largest number of dithers were used to make it:

cl> imhead sc1604_1_deep.fits
sc1604_1_deep.fits[9009,8787][real]: sc1604_1
So this image is 9009x8787 pixels. Next, we make totally blank images for the r,i,z images that are this same size. While we are doing this, we might as well make "blanks" for the weight images as well:
cl> imarith sc1604_1_deep * 0. sc1604_1_r_blank
cl> imarith sc1604_1_deep * 0. sc1604_1_i_blank
cl> imarith sc1604_1_deep * 0. sc1604_1_z_blank
cl> imarith sc1604_1_deep * 0. sc1604_1_r_weight_blank
cl> imarith sc1604_1_deep * 0. sc1604_1_i_weight_blank
cl> imarith sc1604_1_deep * 0. sc1604_1_z_weight_blank
Now, we need to find the offsets between the ultradeep and individual filter images. I've done this in a very hokey way, but it seems to work. You could use tools such as geotran, or write a script to do it. Some display tools like Starlink's Gaia can do this interactively as well.

First, display the ultradeep image. We want to make sure it is not displayed with any binning, so we need to display less than 8192 pixels in each dimension:

ms> set stdimage=imt8192
ms> mscdispl sc1604_1_deep[1000:8000,1000:8000] 1
ms> mscexam log=deep.log keeplog+
Zoom in near the center of the frame, and using the ds9 zoom window in the top right, find the center of a nice, round star, not too bright. Hit "x" to print out the coordinate. Do this for 5 or so stars in the center of the image. I like to print a hardcopy of the image area I am using and mark/number the stars I am using. When you are done with mscexam, hit "q" to quit.

Now, display a single filter image and identify the same objects:

ms> mscdispl sc1604_r[1000:8000,1000:8000] 1
ms> mscexam log=r.log keeplog+
We can paste the two log files together and use awk to find out the offsets:
% paste deep.log r.log > foo
% awk 'NR >1 {print $1-$4,$2-$5}' foo
The x,y offsets should be the same for each object; some might be different by one or two pixels since this is centroiding by eye. But the correct offset should be obvious. Write down these offsets for the r,i, and z images.

Let's say you got an offset of 195,22 (x,y) for the r image from your awk command above. Now, we need to copy the r image into the "blank" we made, at the right location! For this, we need to know the size of the original image.

cl> imhead sc1604_1_r
sc1604_1_r[8804,8765][real]: sc1604_1
cl> imcopy sc1604_1_r.fits[1:8804,1:8765] sc1604_1_r_blank.fits[195:8998,22:8786]
cl> mscdispl sc1604_1_deep[1000:8000,1000:8000] 1
cl> mscdispl sc1604_1_r_blank[1000:8000,1000:8000] 1

Blink the two frames to make sure you got the offsets right. Sometimes I find I am off by one pixel in x or y or both. If so, delete the output frame, redo imarith to make a new blank frame, and redo the imcopy with adjusted x,y offsets. Repeat this procedure, checking the offsets, for the i and z frames. They will all have different offsets!

We will also need to make weight maps for each frame that are the same size. We've already created the blanks. We are going to use the exposure map outputs from mscstack to make the weight maps. These maps will tell us the noise properties of each pixel based on the number of exposures that went into each one. We imcopy the exposure masks into their corresponding blanks just as for each image:

cl> imcopy sc1604_1_r_exp.fits[1:8804,1:8765] sc1604_1_r_weight_blank.fits[195:8998,22:8786]
By running mscstat on the exposure image, we can find out the maximum total exposure time in each pixel. We will want to normalize the weight maps by this exposure time so that the maximum weight is 1:
cl> mscstat sc1604_1_r_weight_blank.fits
#               IMAGE      NPIX      MEAN    STDDEV       MIN       MAX
 sc1604_1_r_weight_blank.fits  77167060     2276.     1508.        0.     3600.
cl> mscarith sc1604_1_r_weight_blank / 3600. sc1604_1_r_weight
cl> imdel sc1604_1_r_weight_blank
Finally, we need to normalize the deep exposures, dividing by the exposure time of a SINGLE exposure (usually 450s for our data). This is because we have averaged in mscstack, not summed. NOTE: If the exposures that went into mscstack are NOT all the same length, you need to divide by the average exposure time.
cl> mscarith sc1604_1_r_blank / 450. sc1604_1_r_shift_norm
Remember, you need to normalize the image and weight map for the ultradeep image as well as the separate filter images!

At this point, we have the necessary images to run SExtractor.

IMPORTANT NOTE: If your data are not all photometric, we will also need to run SExtractor on some photometric data as well! In that case, you will need to go through the shifting procedure discussed above for the individual photometric frames as well as the deep frames!!!

So, first we run SExtractor on the ultradeep (r+i+z) image to detect objects. Example config and param files are available:

To run SExtractor:

% sex -c sc1_deep_sex.config ultradeepimage.fits

Next, we will run SExtractor on the remaining images (combined images in individual filters and single photometric exposures). We will use SExtractor in two-image mode, where detection and setting of apertures is done on the ultradeep image, but all measurements come from the individual filter images. To do this, we need a slightly different config file, wherein we specify the weight maps to use for each image.

Example config and param files for running SExtractor on r-band combined image, using ultradeep image for detection:
target.param (same as for deep image)

NOTE: When using this file for different images, you need to change the following:
GAIN (this should be the gain times the total integration time in seconds if your image is in counts/sec. The gain is ~2 for the LFC CCDs)
and, if necessary the detection thresholds.

You can use the same default.param file as for your ultradeep image.

To run SExtractor in two image mode:

% sex -c sc1_r_sex.config ultradeepimage.fits,singlefilterimage.fits