CFHT queue update, Part 4: so…close

So. Close.

The most recent data has come back from CFHT, and we ALMOST got all our pointings. We are at 95.5% completed on our observing run. MegaCam will be back on the ‘scope fo Dec 13-20, so we may get the last pointing then, but we certainly have lots of data! Enough to keep us busy for a bit.

CFHT queue update, Part 3: almost there

On the eve of leaving for San Diego, I thought I’d drop in an update on the CFHT program (search ‘rogerson’) we’re currently running. This is the 3rd data release (see DR1 and DR2):

we’re now very close to finishing the entire Stripe 82. Note this is only for quasars between 1.73 < z < 2.51, hence the patchy work on the CFHT coverage.

As mentioned above, I’m heading down to San Diego for a run on the Hale 5m telescope at the Palomar Observatory. Stay tuned to this blog for daily updates on that.

Scripting the Prescription

I’ve found it very difficult to keep the photometry steps in my head, and having very limited knowledge of IRAF and its scripting abilities, I’ve generated a Python script. This script (below) does two things: 1. spits out the needed IRAF commands allowing for a quick copy/paste into the IRAF terminal, 2. easy documentation of photometry files for use in time series analysis.

 

——Python to IRAF script—–

author: Jesse Rogerson

date: 9 November 2012

 

import sys
import numpy as np
import pyfits
import math

print ‘This Script is to prompt the user on the steps to calculate’
print ‘the g,r, and i magnitudes from the CFHT data.’
print ‘To use properly, simply follow the instructions and prompts’
print ‘throughout.’
obj=raw_input(‘Name of object:’)
print ”
print ”
print ‘STEP1’
print ‘Unpack and display the 3 fits files to be used (g,r,i)’
print ‘STEP2’
print ‘Use ds9 to determine the coordinates of the quasar to be measured’
print ‘and 5 comparison stars for aperture correction’
print ‘place the comparison coords in ‘+obj+’.coo.1′
print ‘place quasar coords in ‘+obj+’.coo.2′
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP3’
print ‘INPUT FILES:’
gfile=raw_input(‘Name of g filter fits file:’)
rfile=raw_input(‘Name of r filter fits file:’)
ifile=raw_input(‘Name of i filter fits file:’)
ccdNum=raw_input(‘CCDnumber:’)
compCoo=obj+’.coo.1′
objCoo=obj+’.coo.2′
print ‘Name of comparison star coordinatefile:’+compCoo
print ‘Name of quasar coordinate file:’+objCoo
print ”
print ”
print ‘STEP4’
print ‘In IRAF, load noao.digiphot.apphot, set pars files’
print ‘datapars,photpars,fitskypars’
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP5’
#photometer the fits files given the coordinate files
#output files will be *.fits.mag.#
print ‘Photometer commands’
print ‘Run the following commands in IRAF:’
print ‘phot’+’ ‘+gfile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+gfile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
print ‘phot’+’ ‘+rfile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+rfile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
print ‘phot’+’ ‘+ifile+'[ccd’+ccdNum+’] coo=’+compCoo+’ interac-‘
print ‘phot’+’ ‘+ifile+'[ccd’+ccdNum+’] coo=’+objCoo+’ interac-‘
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP6’
#create aperture correction files *.apcor.1
print ‘load noao.digiphot.photcal’
print ‘Run the following commands in IRAF’
print ‘watch for convergence….’
print ”
print ”
print ‘mkapf ‘+gfile+’.mag.1 9 ‘+gfile+’.apcor.1 small=1 large=9′
print ‘mkapf ‘+rfile+’.mag.1 9 ‘+rfile+’.apcor.1 small=1 large=9′
print ‘mkapf ‘+ifile+’.mag.1 9 ‘+ifile+’.apcor.1 small=1 large=9′
var=raw_input(‘continue? (y/n)’)
print ”
print ”
print ‘STEP7’
print ‘Extract quasar magnitude and apply aperture correction’
print ‘Run the following commands in IRAF’
print ‘load noao.digiphot.ptools’
#g:
print ‘txdump ‘+gfile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+gfile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+gfile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+gfile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+gfile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
#r:
print ‘txdump ‘+rfile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+rfile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+rfile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+rfile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+rfile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
#i:
print ‘txdump ‘+ifile+’.mag.2 mag yes | scan (x)’
print ‘txdump ‘+ifile+’.mag.2 merr yes | scan (y)’
print ‘fields ‘+ifile+’.apcor.1 3 | sort num+ rev+ | scan (z)’
print ‘y=sqrt(y**2.+z**2.)’
print ‘fields ‘+ifile+’.apcor.1 2 | sort num+ rev+ | scan (z)’
print ‘x=x+z’
print ‘printf(“%7.3f %6.3f %8.3f”,x,y,z, >> “‘+ifile+’.mea”)’
var=raw_input(‘continue? (y/n)’)
print ”
print ”

print ‘Extracting Header Info for calculations’

files=(gfile,rfile,ifile)
#files=(‘1572573p.fits’,’1572590p.fits’,’1572602p.fits’)
outfile=open(path+obj+’.tbl’,’w’)
#outfile=open(path+’002845.77+010648.3.tbl’,’w’)

outfile.write(‘#inst_mag merr AIRMASS PHOTO_C’
+’PHOTO_CS PHOTO_X PHOT_K MJD \n’)

for f in files:
hdulist=pyfits.open(f)
airmass=hdulist[0].header[‘AIRMASS’]
pc=hdulist[0].header[‘PHOT_C’]
pcs=hdulist[0].header[‘PHOT_CS’]
px=hdulist[0].header[‘PHOT_X’]
pk=hdulist[0].header[‘PHOT_K’]
mjd=hdulist[0].header[‘MJDATE’]
hdulist.close()
mag=np.genfromtxt(path+f+’.mea’,dtype=str)
outfile.write(mag[0]+’ ‘+mag[1]+’ ‘+str(airmass)+’ ‘+str(pc)+’ ‘+str(pcs)
+’ ‘+str(px)+’ ‘+str(pk)+’ ‘+str(mjd)+’\n’)

outfile.close()

dT=np.genfromtxt(path+obj+’.tbl’,skip_header=1,dtype=float)
g=dT[0,:]
r=dT[1,:]
i=dT[2,:]

#R = r + r_c + r_K*(X_r-1)
#Rerr= sqrt(r_err**2.+r_cs**2.)
R=r[0]+r[3]+(r[6]*(r[2]-1.))
Rerr=math.sqrt(r[1]**2.+r[4]**2.)

#G = [ g + g_c + g_K*(X_g-1) – g_X*R ] / [1-g_X]
#Gerr = sqrt(g_err**2.+g_cs**2.)
G=(g[0]+g[3]+(g[6]*(g[2]-1.)) – (g[5]*R))/(1.-g[5])
Gerr=math.sqrt(g[1]**2.+g[4]**2.)

#I = [ i + i_c + i_K*(X_i-1) + i_X*R ] / [1+i_X)
#Ierr = sqrt(i_err**2.+i_cs**2.)
I=(i[0]+i[3]+(i[6]*(i[2]-1.)) + (i[5]*R))/(1.+i[5])
Ierr=math.sqrt(i[1]**2.+i[4]**2.)

print G, Gerr
print R, Rerr
print I, Ierr

outfile=open(path+’BPV/CFHTprop/PHOT/’+obj+’.CFHTphot.dat’,’w’)
outfile.write(g[7]+’ ‘+str(G)+’ ‘+str(Gerr)+’ ‘+str(R)
+’ ‘+str(Rerr)+’ ‘+str(I)+’ ‘+str(Ierr))
outfile.close()
#G-R; err = sqrt(Gerr**2.+Rerr**2.) #G-R & err
#R-I; err = sqrt(Rerr**2.+Ierr**2.) #R-I & err

Photometry, first steps

As the CFHT pointings (see my queue update) are not designed to hit one, but many targets, I was forced to use the pointing information to determine which *.fits file a given Stripe82 object was located. Further to the problem, the MegaPrime camera as 36 CCDs in its almost 1 square degree field of view. In order to determine on which of these 36 CCD a given Stripe 82 quasar would land, I used the following formulae:

Parameters:

CRVAL1 – RA (in deg) of pointing centre
CRVAL2 – Dec (in deg) of pointing centre

RA/Dec – quasar RA/Dec

CRPIX1, CRPIX2, CD1_1, CD2_2 – constants given by the CCD layout.

Calculation:
x=CRPIX1 + (RA – CRVAL1)/CD1_1
y=CRPIX2 + (Dec- CRVAL2)/CD2_2

If I calculate the (x,y) (in pix) using the constants from all 36 CCDs, only one CCD will satisfy: 0 < x < 2048,0 < y < 4612. That is the CCD on which the quasar is located. Labeled *.fits[ccd##] x y

With chip number and (x,y) coordinates in hand, I can begin photometry. Using Phot, in $/noao/digiphot/apphot, requires the use of 4 parameter files: datapars, centerpars, photpars, fitskypars. Notable parameters: datapars.ccdread, gain, airmass, filter, exptime, photpars.apertur, zmag, fitskypars.annulus,dannulus

photpars.aperture allows you to choose a number of apertures from which to sample the object, showing the flattening of photons gained with ever increasing apertures.

1572573p.fits[ccd08]
1690.088 3079.221 (estimated x=1606, 2963)

aper     sum(counts)     area(pixels)     flux(counts)     mag     merr
2.00     7505.262         12.87262         4447.409         20.171 0.017
3.00     14969.63         28.55213         8187.152         19.508 0.013
4.00     23573.49         50.56833         11561.13         19.133 0.012
5.00     33010.21         78.8496           14279.73         18.904 0.011
6.00     43157.03         113.218           16262.43         18.763 0.011
7.00     54404.07         154.346           17739.63         18.669 0.012
9.00     79897.76         254.7832         19374.78         18.573 0.014
12.00  128163.3          452.8017         20601.61         18.506 0.017
15.00   188872.4         706.9811         20931.22         18.489 0.021

It’s clear the magnitude measured here is levelling off, the best mag_err is at an aperture of 6 pixels. The question now is: how do I convert this magnitude to an an actual magnitude. Looking at the CFHT image header I’m given:

PHOT_C = 26.4850 / Elixir zero point – measured for camera run
PHOT_CS = 0.0098 / Elixir zero point – scatter
PHOT_NS = 56 / Elixir zero point – N stars
PHOT_NM = 18 / Elixir zero point – N images
PHOT_C0 = 26.4600 / Elixir zero point – nominal
PHOT_X = 0.1480 / Elixir zero point – color term
PHOT_K = -0.1500 / Elixir zero point – airmass term
PHOT_C1 = ‘g_SDSS ‘ / Elixir zero point – color 1
PHOT_C2 = ‘r_SDSS ‘ / Elixir zero point – color 2
COMMENT Photometric Analysis is incomplete for this image.
COMMENT MAG_SAT and MAG_LIM cannot be determined.
COMMENT Formula for Photometry, based on keywords given in this header:
COMMENT m = -2.5*log(DN) + 2.5*log(EXPTIME)
COMMENT M = m + PHOT_C + PHOT_K*(AIRMASS – 1) + PHOT_X*(PHOT_C1 – PHOT_C2)

It seems te formulae at the bottom are the key.

 

 

Data Overview: CFHT MegaCam

contacts
qsoteam@cfht.hawaii.edu – program questions
dads@cfht.hawaii.edu – data questions

Data downloaded from CFHT is NOT raw. It is detrended by Elixir to correct for dark, bias, and flat-fielding. All images that have been run through this reduction is suffixed with a ‘p’ indicating they have been processed. Each observing group is given a grade for quality, 1 through 5, where 1 is the best.
RAW MegaPrime data have WCS (World Coordinate System) information loaded, but after processing, there is a precision of 0.5 to 1 arcsec across the CCD. This indicates, i can search for my objects by RA/DEC, to identify them, and measure them. Precision is loaded in the FITS header as ‘CERROR.’

quote:
‘The goal of the Elixir calibration is to provide the photometric equation parameter allowing the measurement of magnitudes in the SDSS system.’

Note:
Dec value of each pointing will either be: 0.0 deg, +0.77 deg, or -0.77 deg.
RA value of each pointing will be any step between 0 and 60, 310 and 360

each pointing RA and Dec is the CENTRE of the image, therefore the bounding of the images are:
top right = Ra+0.5 deg,Dec+0.5 deg
top left = Ra-0.5 deg, Dec+0.5 deg
bottom right = Ra+0.5 deg, Dec-0.5 deg
bottom left = Ra-0.5 deg, Dec-0.5 deg

Note that CFHT provides data in *.fz compression. Can use fpack, and funpack to work with these.

The relevant URL is:
http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fpack/
The Megacam page about image compression is:
http://www3.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/cfht/cfitsio_comp.html