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 !!
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: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 ....etc.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 startThe 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.
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 login.cl (or loginuser.cl) as well.
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.coefsYou 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
-3.16986
-0.2333422
0.1187139
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!
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_1So 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_blankNow, 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
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_blankFinally, 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_normRemember, 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:
sc1_deep_sex.config
target.param
To run SExtractor:
% sex -c sc1_deep_sex.config ultradeepimage.fits
Example config and param files for running SExtractor on r-band combined image, using ultradeep image for detection:
sc1_r_sex.config
target.param (same as for deep image)
NOTE: When using this file for different images, you need to change the following:
CATALOG_NAME
MAG_ZEROPOINT
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)
WEIGHT_IMAGE
CHECKIMAGE_NAME
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