#!/usr/bin/env python import os from pyraf import iraf import pyfits from numpy import * from pylab import * from time import time iraf.noao(_doprint=0) iraf.digiphot(_doprint=0) iraf.daophot(_doprint=0) iraf.daophot.cache = 'yes' iraf.delete('*.?') iraf.delete('*psf.fits,*sys.fits,*.sys,*geodb,*.?.fits,*.??.fits') # CROWED-FIELD PHOTOMETRY PIPELINE (CPP) FOR q957g FRAMES imlist = 'imlist' # list of frames (input) wfile = 'q957g_BDR' # basic data release (output) fin = open(imlist,'r') time1 = time() N = 200 # number of frames (total number of available frames = 200) Nstar = 7 # number of stars (total number of available stars = 19): X, G, F, H, E, D, R ... ghr = 0.028 # galaxy-to-H ratio in the g band size = 32 # half-size of subframes (windows) for IMFITFITS scale = 3.59234 # CCD properties: scale (pixel/arcsecond), gain (e-/count) and readnoise (e-/pixel) gain = 2.796 readnoise = 7. t = zeros(N,float) qA = zeros(N,float) qB = zeros(N,float) chi2 = zeros(N,float) gal = zeros(N,float) stars = zeros((Nstar,N),float) #starserr = zeros((Nstar,N),float) ellipticity = zeros(N,float) fwhm = zeros(N,float) PA = zeros(N,float) for i in range(0,N): filename = fin.readline() filename = filename[:-1] hdr = pyfits.getheader(filename+'.fits') seeing = hdr['seeing'] t[i] = hdr['mjd'] - 50000 # Parameters for the initial estimation of stellar background and instrumental flux (and the improvement of stellar position) iraf.datapars.noise = 'poisson' # (a) "electronic" noise model and gain factor iraf.datapars.readnoise = readnoise iraf.datapars.epadu = gain iraf.centerpars.calgorithm = 'centroid' # (b) centre iraf.centerpars.cbox = 2*seeing iraf.fitskypars.salgorithm = 'median' # (c) background iraf.fitskypars.annulus = 4*seeing+1 iraf.fitskypars.dannulus = 2*seeing iraf.photpars.weighting = 'constant' # (d) flux iraf.photpars.apertur = 2*seeing # IRAF initial coordinates of the stars in *.coo.1 (improved coordinates and IRAF initial fluxes in *.mag.1) # Procedure: (a) geomap: *.xy.1 (T) -> *.geodb (TL), (b) geoxytran: SDSS19.dat (input) + *.geodb (TL) + *.xy.1 (T) -> *.coo.1 (output) iraf.immatch.geomap('xy1/'+filename+'.xy.1',filename+'.geodb','1','1024','1','1024',fitgeo='rotate',inter='no',verbose='no') iraf.immatch.geoxytran('SDSS19.dat',filename+'.coo.1',filename+'.geodb','xy1/'+filename+'.xy.1',direct='forward') iraf.phot(filename,filename+'.coo.1',filename+'.mag.1','',verify='no',verbose='no') # IRAF initial coordinates of the quasar images in *.coo.2 (improved coordinates and IRAF initial fluxes in *.mag.2) iraf.centerpars.cbox = 11 # photometry in a crowed field (3" apertures and external background) iraf.fitskypars.annulus = 45 iraf.fitskypars.dannulus = 11 iraf.photpars.apertur = 11 iraf.immatch.geoxytran('reference2',filename+'.coo.2',filename+'.geodb','xy1/'+filename+'.xy.1',direct='forward') iraf.phot(filename,filename+'.coo.2',filename+'.mag.2','',verify='no',verbose='no') s = iraf.pdump(filename+'.mag.2','xc','id=1',Stdout=1) # reading the improved quasar coordinates xqA = float(s[0]) s = iraf.pdump(filename+'.mag.2','yc','id=1',Stdout=1) yqA = float(s[0]) s = iraf.pdump(filename+'.mag.2','xc','id=2',Stdout=1) xqB = float(s[0]) s = iraf.pdump(filename+'.mag.2','yc','id=2',Stdout=1) yqB = float(s[0]) # Coordinates of the centre of the lens system (xc,yc) and subframe with quasar images (*.sys.fits) xc = (xqA + xqB)/2 yc = (yqA + yqB)/2 section = '['+str(int(xc)-size+1)+':'+str(int(xc)+size)+','+str(int(yc)-size+1)+':'+str(int(yc)+size)+']' iraf.imcopy(filename+section,filename+'.sys',verbose='no') fluxqA = iraf.pdump(filename+'.mag.2','flux','id=1',Stdout=1) # reading the quasar fluxes fluxqB = iraf.pdump(filename+'.mag.2','flux','id=2',Stdout=1) # IRAF background around the quasar images (around the central coordinates in *.coo.3) fcoo3 = open(filename+'.coo.3','w') fcoo3.write(str(xc)+' '+str(yc)+'\n') fcoo3.close() iraf.phot(filename,filename+'.coo.3',filename+'.mag.3','',verify='no',verbose='no') msky = iraf.pdump(filename+'.mag.3','msky','id=1',Stdout=1) # PSF (H) star: reading its improved coordinates (xPSF,yPSF), associated background (skyPSF) and flux (fluxPSF), cutting to the PSF # subframe (*.psf.fits) and removing the local background s = iraf.pdump(filename+'.mag.1','xc','id=4',Stdout=1) # id = 4 => star H (PSF) xPSF = float(s[0]) s = iraf.pdump(filename+'.mag.1','yc','id=4',Stdout=1) yPSF = float(s[0]) skyPSF = iraf.pdump(filename+'.mag.1','msky','id=4',Stdout=1) fluxPSF = iraf.pdump(filename+'.mag.1','flux','id=4',Stdout=1) section = '['+str(int(xPSF)-size+1)+':'+str(int(xPSF)+size)+','+str(int(yPSF)-size+1)+':'+str(int(yPSF)+size)+']' iraf.imcopy(filename+section,filename+'.psf',verbose='no') iraf.imarith(filename+'.psf','-',float(skyPSF[0]),filename+'.psf') # Estimating the initial value of the rotation: rotation(TL) + rotation(reference frame) fgeodb = open(filename+'.geodb') for j in range(0,20): s = fgeodb.readline() if 'xrotation' in s: words = s.split() angle = float(words[1]) angle = angle + 3.1564 if angle > 180.: angle = angle - 360 if angle < -180.: angle = angle + 360 # IMFITFITS stars photometry for ind in range(Nstar): s = iraf.pdump(filename+'.mag.1','xc','id='+str(ind+1),Stdout=1) # IRAF improved coordinates x = float(s[0]) s = iraf.pdump(filename+'.mag.1','yc','id='+str(ind+1),Stdout=1) y = float(s[0]) if (int(x)-size)>0 and(int(y)-size)>0 and(int(x)+size)<=1024 and(int(y)+size)<=1024: # rejecting stars close to CCD edges sky = iraf.pdump(filename+'.mag.1','msky','id='+str(ind+1),Stdout=1) # IRAF backgrounds flux = iraf.pdump(filename+'.mag.1','flux','id='+str(ind+1),Stdout=1) # IRAF instrumental fluxes section = '['+str(int(x)-size+1)+':'+str(int(x)+size)+','+str(int(y)-size+1) +':'+str(int(y)+size)+']' iraf.imcopy(filename+section,filename+'.'+str(ind),verbose='no') # stellar subframes (*.ind.fits) fini = open(filename+'.ini.'+str(ind),'w') # ini files (*.ini.ind) fini.write('background '+sky[0]+'\n') fini.write('psfblur 0\n') fini.write('ngalaxy 0\n') fini.write('nptsrc 1\n') fini.write('ndevauc 0\n') fini.write('ptsrc[0].intens '+flux[0]+'\n') fini.write('ptsrc[0].x 0'+'\n') fini.write('ptsrc[0].y 0'+'\n') fini.write('rotation '+str(angle)+'\n') fini.write('scale '+str(scale)+'\n') fini.write('sy_over_sx 1.'+'\n') x0 = size-1+x-int(x) # positions in subframes y0 = size-1+y-int(y) fini.write('x0 '+str(x0)+'\n') fini.write('y0 '+str(y0)+'\n') fini.close() ffree = open(filename+'.free.'+str(ind),'w') # free files (*.free.ind): 4 free pars sky0 = str(float(sky[0])*0.1) ffree.write('background '+sky0+'\n') flux0 = str(float(flux[0])*0.1) ffree.write('ptsrc[0].intens '+flux0+'\n') ffree.write('x0 '+str(1.)+'\n') ffree.write('y0 '+str(1.)+'\n') ffree.close() s0 = 'imfitfits '+filename+'.free.'+str(ind)+' '+filename+'.'+str(ind)+' '+filename+'.psf '+filename+'.mod.'+str(ind)+ ' < '+filename+'.ini.'+str(ind)+' | adjcoords > '+filename+'.out.'+str(ind) # running IMFITFITS -> *.out.ind & *.mod.ind.fits os.system(s0) f1 = open(filename+'.out.'+str(ind),'r') for j in range(40): s = f1.readline() if 'normalized' in s: words = s.split() norm = float(words[4]) # reading the PSF normalization factor if 'ptsrc[0].intens' in s: break f1.close() words = s.split() stars[ind,i] = float(words[1]) # reading the IMFITFITS final instrumental fluxes # IMFITFITS quasar images photometry # Constraints on ptsrc[1].x, ptsrc[1].y, devauc[0].x, devauc[0].y, devauc[0].reff, devauc[0].bovera and devauc[0].angle. The constraints # were obtained by # (a) Bernstein et al. 1997, ApJ 483, L79-L82: optical astrometry and surface photometry of the lens galaxy (HST frames) # (b) Keeton, Kochanek & Falco 1998, ApJ 509, 561-578: optical structure of the lens galaxy (Reff, e = 1 - b/a and PA from a de Vaucouleurs # profile fit to HST frames) # (c) CASTLES Web site: astrometry in agreement with Ref. (a) fini = open(filename+'.ini.sys','w') # ini file (*.ini.sys) fini.write('background '+msky[0]+'\n') fini.write('psfblur 0\n') fini.write('ngalaxy 0\n') fini.write('nptsrc 2\n') fini.write('ndevauc 1\n') fini.write('ptsrc[0].intens '+fluxqA[0]+'\n') fini.write('ptsrc[0].x 0\n') fini.write('ptsrc[0].y 0\n') fini.write('ptsrc[1].intens '+fluxqB[0]+'\n') fini.write('ptsrc[1].x '+str(-1.229)+'\n') # constraint 1 (C1) fini.write('ptsrc[1].y '+str(-6.048)+'\n') # C2 fluxgal = stars[3,i]*ghr fini.write('devauc[0].intens '+str(fluxgal)+'\n') # additional constraint on the lens galaxy flux fini.write('devauc[0].x -1.406\n') # C3 fini.write('devauc[0].y -5.027\n') # C4 fini.write('devauc[0].reff 4.63\n') # C5 fini.write('devauc[0].bovera 0.79\n') # C6 fini.write('devauc[0].angle 49\n') # C7 fini.write('rotation '+str(angle)+'\n') fini.write('scale '+str(scale)+'\n') fini.write('sy_over_sx 1.'+'\n') x0 = size-1+xc-int(xc) +(xqA-xqB)/2 # position of A in the system subframe y0 = size-1+yc-int(yc) +(yqA-yqB)/2 fini.write('x0 '+str(x0)+'\n') fini.write('y0 '+str(y0)+'\n') fini.close() ffree = open(filename+'.free.sys','w') # free file (*.free.sys): 6 free pars msky0 = str(float(msky[0])*0.1) ffree.write('background '+msky0+'\n') fluxqA0 = str(float(fluxqA[0])*0.1) fluxqB0 = str(float(fluxqB[0])*0.1) # fluxgal0 = str(fluxgal*0.1) ffree.write('ptsrc[0].intens '+fluxqA0+'\n') ffree.write('ptsrc[1].intens '+fluxqB0+'\n') # ffree.write('devauc[0].intens '+fluxgal0+'\n') ffree.write('rotation '+str(1.)+'\n') ffree.write('x0 '+str(1.)+'\n') ffree.write('y0 '+str(1.)+'\n') ffree.close() s0 = 'imfitfits '+filename+'.free.sys '+filename+'.sys '+filename+'.psf '+filename+'.mod.sys'+ ' < '+filename+'.ini.sys | adjcoords > '+filename+'.out.sys' # running IMFITFITS -> *.out.sys & *.mod.sys.fits os.system(s0) iraf.imarith(filename+'.sys','-',filename+'.mod.sys',filename+'.res.sys') # obtaining the residual frame (*.res.sys.fits) f1 = open(filename+'.out.sys','r') for j in range(40): s = f1.readline() if 'background' in s: words = s.split() background = float(words[1]) # reading the IMFITFITS final background (QSOs) if 'rotation' in s: words = s.split() rotation = float(words[1]) # reading the IMFITFITS final rotation if 'x0' in s: words = s.split() x0 = float(words[1]) # reading the IMFITFITS final position of A in the system subfranme if 'y0' in s: words = s.split() y0 = float(words[1]) if 'normalized' in s: words = s.split() norm = float(words[4]) # reading (again) the PSF normalization factor if 'ptsrc[0].intens' in s: words = s.split() qA[i] = float(words[1]) # reading the IMFITFITS final instrumental flux of A if 'ptsrc[1].intens' in s: words = s.split() qB[i] = float(words[1]) # reading the IMFITFITS final instrumental flux of B if 'devauc[0].intens' in s: words = s.split() fgal = float(words[1]) gal[i] = fgal/stars[3,i] # reading the galaxy-to-H ratio in the g band break f1.close() # Quality of the fit to the subframe of interest (lens system): chi2 iraf.imarith(filename+'.res.sys','*',filename+'.res.sys',filename+'.res.1.sys') # *.res.1.sys.fits contains Ni = (Fi - F*i)^2 iraf.imarith(filename+'.sys','/',gain,filename+'.1.sys') # *.1.sys.fits contains Fi/gain iraf.imarith(filename+'.1.sys','+',(readnoise/gain)**2,filename+'.2.sys') # *.2.sys.fits contains Di = Fi/gain + (rn/gain)^2 iraf.imdiv(filename+'.res.1.sys',filename+'.2.sys',filename+'.res.2.sys',rescale='norescale') # *.res.2.sys.fits contains Ni/Di im = pyfits.getdata(filename+'.res.2.sys.fits') chi2[i] = im.sum()/4090. # Quality of the frame stddev = hdr['stddev'] # (a) shape of images (H star) iraf.starfind(filename+'.3',filename+'.obj.1',seeing/2,5*stddev,roundhi=1.0,verbose='no') fobj1 = open(filename+'.obj.1','r') for j in range(17): s = fobj1.readline() words = s.split() fwhm[i] = 2*float(words[4])/scale # FWHM (arcsec) ellipticity[i] = float(words[5]) # ellipticity = 1 - b/a PA[i] = float(words[6]) # PA (deg) fin.close() fin = open(imlist,'r') fw = open(wfile,'w') for i in range(0,N): filename = fin.readline() filename = filename[:-1] hdr = pyfits.getheader(filename+'.fits') seeing = hdr['seeing'] fwhmhead = seeing/scale fwhmdev = fwhmhead - fwhm[i] # deviation between FWHM in header and measured FWHM (arcsec) image = iraf.pdump(filename+'.mag.1','image','id=1',Stdout=1) fluxA = iraf.pdump(filename+'.mag.2','flux','id=1',Stdout=1) area = iraf.pdump(filename+'.mag.2','area','id=1',Stdout=1) msky = iraf.pdump(filename+'.mag.2','msky','id=1',Stdout=1) fluxA = float(fluxA[0]) area = float(area[0]) msky = float(msky[0]) SNR = fluxA*gain/sqrt(fluxA*gain+area*(msky*gain+readnoise*readnoise)) # (b) signal-to-noise ratio (A quasar image) # Relative instrumental magnitudes A - E and B - E (checking the pipeline) relA = -2.5*log10(qA[i]/stars[4,i]) relB = -2.5*log10(qB[i]/stars[4,i]) # Making the q957g_BDR file: (1) frame number, (2) frame label, (3) MJD - 50000 (days), (4) FWHM ("), (5) ellipticity, (6) SNR, (7) flux of A (counts), (8) flux of B (counts), (9) chi2, (10-16) fluxes of X, G, F, H, E, D, R (counts), and (17-18) A - E and B - E (mag) print i+1, image[0][0:5], str(t[i])[0:9], str(fwhm[i])[0:4], str(ellipticity[i])[0:5], str(int(SNR)), qA[i], qB[i], str(chi2[i])[0:4], for j in range(Nstar-1): print stars[j,i], print stars[Nstar-1,i], print str(relA)[0:5], str(relB)[0:5] print >> fw, i+1, image[0][0:5], str(t[i])[0:9], str(fwhm[i])[0:4], str(ellipticity[i])[0:5], str(int(SNR)), qA[i], qB[i], str(chi2[i])[0:4], for j in range(Nstar-1): print >> fw, stars[j,i], print >> fw, stars[Nstar-1,i], print >> fw, str(relA)[0:5], str(relB)[0:5] fin.close() fw.close() # Runtime time2 = time() dtime = time2 - time1 hours = int(dtime/3600) minutes = int((dtime-hours*3600)/60) secs = int(dtime-hours*3600-minutes*60) print hours, 'h', minutes, 'min', secs, 'sec' # Figure of instrumental magnitudes A and B (calibrated from the SDSS magnitude of the H star) norm = stars[3,:] qA = -2.5*log10(qA/norm) + 15.1164 qB = -2.5*log10(qB/norm) + 15.1164 subplot(211, xlim=(3350,4250), ylim=(17.45,17.05), autoscale_on=False) plot(t, qA, 'ro') title('A image') subplot(212, xlim=(3350,4250), ylim=(17.45,17.05), autoscale_on=False) plot(t, qB, 'go') title('B image') xlabel('MJD-50000') show()