#! /usr/bin/env python # 2016-Jul-06 shaw@noao.edu import copy import sys import numpy as np from pyraf import iraf from pyraf.iraf import gemini, gemtools, gmos import fileSelect as fs def gmos_lsns_proc(): ''' GMOS Data Reduction Cookbook companion script to the chapter: "Reduction of Longslit Nod & Shuffle Spectra with IRAF" PyRAF script to: Process GMOS spectra for XXX, in program YY. The names for the relevant header keywords and their expected values are described in the DRC chapter entitled "Supplementary Material" Perform the following starting in the parent work directory: cd /path/to/work_directory Place the fileSelect.py module in your work directory. Now execute this script from the unix prompt: python gmos_lsns_proc.py ''' print ("### Begin Processing GMOS/LS_NS Spectra ###") print ("###") print ("=== Creating MasterCals ===") # This whole example depends upon first having built an sqlite3 database of metadata: # cd ./raw # python obslog.py obsLog.sqlite3 dbFile='raw/obsLog.sqlite3' # From the work_directory: # Create a query dictionary of essential parameter=value pairs. # Select bias exposures within a few days of the Oct epoch observations: qd = {'Oct': {'use_me':1, 'Instrument':'GMOS-N','CcdBin':'2 1','RoI':'Full', 'Disperser':'R400+_%','CentWave':740.0,'AperMask':'NS0.75arcsec', 'Object':'%.w%', 'DateObs':'2003-10-20:2003-11-02'} } # Make additional copies for other epochs of observations; # this will yield ~50 bias exposures per epoch. qd['Nov'] = copy.deepcopy(qd['Oct']) qd['Nov']['DateObs'] = '2003-11-22:2003-11-25' qd['Dec'] = copy.deepcopy(qd['Oct']) qd['Dec']['DateObs'] = '2003-11-17:2003-12-26' qd['NovDec'] = copy.deepcopy(qd['Oct']) qd['NovDec']['DateObs'] = '2003-11-22:2003-12-26' qd['Jan'] = copy.deepcopy(qd['Oct']) qd['Jan']['DateObs'] = '2004-01-20:2004-02-16' qd['All'] = copy.deepcopy(qd['Oct']) qd['All']['DateObs'] = '2003-11-22:2004-02-16' print (" --Creating Bias MasterCals--") # Use primarily the default task parameters. gemtools.gemextn.unlearn() # Disarm a bug in gbias gmos.gbias.unlearn() biasPars = { 'logfile':'biasLog.txt','rawpath':'./raw/','fl_vardq':'yes','verbose':'no' } epochs = ['Oct','NovDec','Jan'] for e in epochs: # The following SQL generates the list of bias files to process. SQL = fs.createQuery('bias', qd[e]) biasFiles = fs.fileListQuery(dbFile, SQL, qd[e]) # The str.join() funciton is needed to transform a python list into a # comma-separated string of file names that IRAF can understand. if len(biasFiles) > 1: gmos.gbias(','.join(str(x) for x in biasFiles), 'MCbias'+e, **biasPars) # Clean up iraf.imdel('gN200*.fits,gsN200*.fits') print (" --Creating Dark MasterCal--") # Use primarily the default task parameters. #gemtools.gemextn.unlearn() # Disarm a bug in gbias gmos.gnsdark.unlearn() darkPars = { 'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','bias':'MCbiasNovDec', 'rawpath':'./raw/','fl_vardq':'yes','fl_inter':'no', 'logfile':'gnsdarkLog.txt','verbose':'no' } darkFiles = fs.fileListQuery(dbFile, fs.createQuery('dark', qd['All']), qd['All']) gmos.gnsdark(','.join(str(x) for x in darkFiles), 'MCdark', **darkPars) # Fudge the number of N&S cycles to match most science files. iraf.hedit('MCdark'+'.fits[0]','NODCOUNT',10,add='no', update='yes',verify='no') print (" -- Creating Wavelength Calibration MasterCals --") # Reduce the Arcs to perform bias correction, and mosaic the extensions. arcPars = { 'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_gscrrej':'no', 'fl_dark':'no','fl_flat':'no','fl_gmosaic':'yes','fl_fixpix':'no', 'fl_gsappwave':'yes','fl_oversize':'no','fl_vardq':'no','fl_fulldq':'no', 'rawpath':'./raw','fl_inter':'no','logfile':'gsreduceLog.txt','verbose':'no' } # Wavelength solution might be improved if done interactively. gmos.gswavelength.unlearn() wavePars = { 'coordlist':'gmos$data/CuAr_GMOS.dat','fwidth':6,'nsum':50,'function':'chebyshev', 'order':4,'fl_inter':'no','logfile':'gswaveLog.txt','verbose':'no' } for e in epochs: arcFiles = fs.fileListQuery(dbFile, fs.createQuery('arc', qd[e]), qd[e]) if len(arcFiles) > 0: gmos.gsreduce(','.join(str(x) for x in arcFiles), bias='MCbias'+e, **arcPars) gmos.gswavelength(','.join('gs'+str(x) for x in arcFiles), **wavePars) iraf.imdel('gN200*.fits') print (" -- Creating GCAL Spectral Flat-Fields for Science Targets --") # Set the task parameters. gmos.gireduce.unlearn() gmos.gsflat.unlearn() sciFlatFlags = { 'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_dark':'no', 'fl_fixpix':'no','fl_oversize':'no','fl_vardq':'yes','fl_fulldq':'no', 'rawpath':'./raw','fl_inter':'no','fl_detec':'yes', 'fl_double':'yes','nshuffle':1536, 'function':'spline3','order':'5,5,17', 'logfile':'gsflatLog.txt','verbose':'no' } # Response curve fitting should be done interactively. # At the prompt, select image line=2300 for the fit; hit the second time. for e in epochs: flatFiles = fs.fileListQuery(dbFile, fs.createQuery('gcalFlat', qd[e]), qd[e]) if len(flatFiles) > 0: gmos.gsflat(','.join(str(x) for x in flatFiles), 'MCflat'+e, bias='MCbias'+e, **sciFlatFlags) # Mask emission lines. file = 'MCflat' + e + '.fits[sci,2]' if e == 'Oct': iraf.imreplace(file+'[830:852,*]', 1.) else: iraf.imreplace(file+'[840:865,*]', 1.) print ("=== Processing Standard Stars ===") # Associate target names (in the headers) with the magic "iraf name", # the Arc exposure, and matching flat-fields. stdPars = { 'BD+28d4211':{'arc':'gsN20031028S0152','irafName':'bd284211', 'flat':['N20031028S0153','N20031028S0164']}, 'Hiltner600':{'arc':'gsN20040130S0126','irafName':'hilt600', 'flat':['N20040128S0048','N20040129S0115']} } print (" -- Creating GCAL Spectral Flat-Fields for Standard Stars --") stdFlatFlags = { 'fl_over':'no','fl_trim':'no','fl_bias':'no','fl_dark':'no', 'fl_fixpix':'no','fl_oversize':'no','fl_vardq':'yes','fl_fulldq':'no', 'rawpath':'./','fl_inter':'yes','fl_detec':'yes', 'function':'spline3','order':'5,5,17', 'logfile':'gsflatLog.txt','verbose':'no' } for targ,p in stdPars.iteritems(): flatFiles = p['flat'] outFile = targ + '_flat' gmos.gsflat(','.join('gs'+str(x) for x in flatFiles), outFile, **stdFlatFlags) # Mask emission lines. ext = outFile + '.fits[sci,2]' if e == 'Oct': iraf.imreplace(ext+'[830:852,*]', 1.) else: iraf.imreplace(ext+'[840:865,*]', 1.) iraf.imdel('gN200*.fits,gsN200*.fits') print (" -- Basic Processing --") qd['Std'] = {'use_me':1, 'Instrument':'GMOS-N','CcdBin':'2 1','RoI':'Full', 'Disperser':'R400+_%','CentWave':740.0,'AperMask':'NS0.75arcsec', 'Object':'BD+28d4211', 'DateObs':'*', } # Use primarily the default task parameters. gmos.gsreduce.unlearn() stdFlags = { 'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_gscrrej':'no', 'fl_dark':'no','fl_flat':'yes','fl_gmosaic':'yes','fl_fixpix':'yes', 'fl_gsappwave':'yes','fl_oversize':'no', 'fl_vardq':'no','fl_fulldq':'no','rawpath':'./raw', 'fl_inter':'no','logfile':'gsreduceLog.txt','verbose':'no' } for e in epochs: bias = 'MCbias'+e for targ in stdPars.keys(): qd[e]['Object'] = targ flat = targ + '_flat' stdFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qd[e]), qd[e]) if len(stdFiles) > 0: gmos.gsreduce(','.join(str(x) for x in stdFiles), bias=bias, flat=flat, **stdFlags) # Clean up iraf.imdel('gN200*.fits') print (" -- Wavelength calibrate, extract, and combine --") # The paired Std exposures are offset 10 arcsec from each other: # e.g., XOFFSET = -1.25 (both); 'YOFFSET= 0' for 0143, and '=10.' for 0144. # Set task parameters. gmos.gstransform.unlearn() transPars = { 'fl_vardq':'no','interptype':'linear','fl_flux':'yes', 'logfile':'gstransLog.txt' } gmos.gsextract.unlearn() extrPars = { 'apwidth':3.,'fl_inter':'no','find':'yes', 'trace':'yes','tfunction':'spline3','torder':'5','tnsum':20, 'background':'fit','bfunction':'chebyshev','border':2, 'fl_vardq':'no','logfile':'gsextrLog.txt' } gemtools.gemcombine.unlearn() combPars = { 'combine':'average','reject':'none','logfile':'combLog.txt','verbose':'no' } for targ,p in stdPars.iteritems(): qs = qd['Std'] qs['Object'] = targ stdFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qs), qs) gmos.gstransform (','.join('gs'+str(x) for x in stdFiles), wavtraname=p['arc'], **transPars) gmos.gsextract(','.join('tgs'+str(x) for x in stdFiles), **extrPars) gemtools.gemcombine(','.join('etgs'+str(x) for x in stdFiles), targ, **combPars) iraf.imdel('gN200*.fits,gsN200*.fits') print (" -- Interactive sensitivity calibration --") stars = str(stdPars['BD+28d4211']['irafName']) + ',' + str(stdPars['Hiltner600']['irafName']) # Monochomatic mags in IRAF are absent beyond 8900 A, so use monochromatic # magnitudes from ESO instead. Place files in work directory. #'fl_inter':'yes','starname':stars,'caldir':'onedstds$spec50cal/', sensPars = { 'fl_inter':'yes','starname':stars,'caldir':'./', 'observatory':'Gemini-North','extinction':'./mk_extinct.txt', 'function':'chebyshev','order':6,'verbose':'no','logfile':'gsstdLog.txt' } #gsstandard('BD+28d4211,Hiltner600', 'std.txt', 'sens', **sensFlags) gsstandard(','.join(str(x) for x in stdPars.keys()), 'std.txt', 'sens', **sensPars) print ("=== Processing Science Targets ===") # Flata and an arc was obtained for each science exposure (more or less). # Small, relative y-offsets were imposed on the exposures in a sequence. # Capture the associations in a dict, though it would be better to do this # from a database table. # Reducing all the targets is overkill for this tutorial. Also, there are # exceptions to the otherwise simple associations between science files # and calibration files, which would make the flow a little messy. sciPars = {'d034.waa7_10': {'arc':'N20031028S0152', 'sci':['N20031028S0154','N20031028S0157'], 'flat':['N20031028S0153','N20031028S0158']}} sciPars['d085.waa5_16'] = {'arc':'N20031028S0163', 'sci':['N20031028S0165','N20031028S0167','N20031028S0170'], 'flat':['N20031028S0164','N20031028S0168','N20031028S0169']} sciPars['d115.wbb6_11'] = {'arc':'N20031028S0174', 'sci':['N20031028S0176','N20031028S0178','N20031028S0181'], 'flat':['N20031028S0175','N20031028S0179','N20031028S0180']} sciPars['e309.waa9_14'] = {'arc':'N20031123S0028', 'sci':['N20031123S0030','N20031123S0031','N20031123S0034'], 'flat':['N20031123S0029','N20031123S0032','N20031123S0033']} sciPars['f017.wdd9_10'] = {'arc':'N20031220S0192', 'sci':['N20031220S0194','N20031220S0196','N20031220S0199'], 'flat':['N20031220S0193','N20031220S0197','N20031220S0198']} sciPars['f213.wbb4_12'] = {'arc':'N20031219S0112', 'sci':['N20031219S0114','N20031219S0116'], 'flat':['N20031219S0113','N20031219S0117']} sciPars['f247.wbb8_10'] = {'arc':'N20040129S0114', 'sci':['N20040129S0116','N20040129S0117'], 'flat':['N20040129S0118','N20040129S0119']} sciPars['f509.wcc4_4'] = {'arc':'N20031224S0026', 'sci':['N20031224S0028','N20031224S0030'], 'flat':['N20031224S0027','N20031224S0031']} sciPars['f525.wcc4_12'] = {'arc':'N20031226S0015', 'sci':['N20031226S0017','N20031226S0019'], 'flat':['N20031226S0016','N20031226S0020']} sciFlatFlags = copy.deepcopy(stdFlatFlags) sciFlatFlags['fl_double'] = 'yes' sciFlatFlags['nshuffle'] = 1536 print (" -- Basic Processing --") # Set the task parameters. gmos.gsreduce.unlearn() sciFlags = { 'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_gscrrej':'no', 'fl_dark':'yes','fl_flat':'no','fl_gmosaic':'no','fl_fixpix':'no', 'fl_gsappwave':'no','fl_cut':'no','fl_oversize':'no', 'fl_vardq':'yes','fl_fulldq':'no','ovs_flinter':'no', 'rawpath':'./raw','fl_inter':'no','logfile':'gnsreduceLog.txt','verbose':'no' } epochs = ['Oct','NovDec','Jan'] for e in epochs: bias = 'MCbias'+e qs = qd[e] sciFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qs), qs) gmos.gsreduce(','.join(str(x) for x in sciFiles), bias=bias, dark='MCdark', **sciFlags) # Regrettably, gnscombine reads image offsets from a file, rather than from the image # header. The y-value can be derived from the DTAX keywords in an exposure sequence, # where the value for the first exposure must be subtracted from that of each exposure. qd['Sci'] = {'use_me':1, 'Instrument':'GMOS-N','CcdBin':'2 1','RoI':'Full', 'Disperser':'R400+_%','CentWave':740.0,'AperMask':'NS0.75arcsec', 'Object':'%.w%', 'DateObs':'*', } sciFlags.update({'fl_over':'no','fl_trim':'no','fl_bias':'no','fl_dark':'no', 'fl_flat':'yes','rawpath':'./','outpref':'f'}) sciFlags2 = { 'fl_over':'no','fl_trim':'no','fl_bias':'no','fl_gscrrej':'yes', 'fl_dark':'no','fl_flat':'no','fl_gmosaic':'yes','fl_fixpix':'yes', 'fl_gsappwave':'yes','fl_cut':'yes','fl_oversize':'no', 'fl_vardq':'yes','fl_fulldq':'no','fl_inter':'no','outpref':'m', 'rawpath':'./','logfile':'gnsreduceLog.txt','verbose':'no' } gmos.gnscombine.unlearn() gnscombPars = { 'fl_vardq':'yes','logfile':'gemcombineLog.txt','verbose':'no' } gmos.gstransform.unlearn() transPars = { 'fl_vardq':'yes','interptype':'linear','fl_flux':'yes', 'logfile':'gstransLog.txt' } # Loop through each target qs = qd['Sci'] for targ,pars in sciPars.iteritems(): qs['Object'] = targ + '%' sciFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qs), qs) shifts = fs.offsetQuery(dbFile, fs.createQuery('offset',qs),qs) if len(sciFiles) > 0: # Now we have a sequence of exposures for a single target. # Compute and write to a file the Y-offsets relative to the first exposure. shft0 = shifts[0][1] shiftFile = targ+'_shft.txt' with open (shiftFile, 'w') as o: for s in shifts: line = '0 %d %s\n' % (int(np.rint(s[1]-shft0)), s[0]) o.write(line) # Apply the matched (by index) flat-field. for f in sciFiles: flat = pars['flat'][pars['sci'].index(f)] gmos.gsreduce(f, flat=flat, **sciFlags) targFile = targ + '_comb.fits' # Combine, mosaic the exposures, and wavelength calibrate. gmos.gnscombine(','.join('fgs'+str(x) for x in sciFiles), shiftFile, targFile, **gnscombPars) gmos.gsreduce(targFile, **sciFlags2) gmos.gstransform('m'+targFile, wavtraname='gs'+pars['arc'], **transPars) # Perform extraction twice to extract both "positive" and "negative" spectra, # then combine and calibrate gmos.gsextract.unlearn() extrFlags = { 'apwidth':4.,'find':'yes','trace':'no', 'tfunction':'chebyshev','torder':2,'tnsum':100,'tstep':100, 'background':'none', 'fl_inter':'yes','fl_vardq':'no','logfile':'gsextrLog.txt' } gemtools.gemarith.unlearn() gemarithPars = { 'fl_vardq':'no','verbose':'no','logfile':'gemarithLog.txt' } gmos.gscalibrate.unlearn() calibPars = { 'extinction':'./mk_extinct.txt','fl_ext':'yes','fl_scale':'no', 'sfunc':'sens','fl_vardq':'no','logfile':'gscalibLog.txt' } for targ in sciPars.keys(): targFile = targ + '_comb' outPos = 'eptm'+ targFile outNeg = 'entm'+ targFile # Extract interactively gmos.gsextract('tm'+targFile, outimages=outPos, **extrPars) gmos.gsextract('tm'+targFile, outimages=outNeg, **extrPars) gemtools.gemarith(outPos, '-', outNeg, 'etm'+targFile, **gemarithPars) gmos.gscalibrate('etm'+targFile, **calibPars) print ("=== Finished Calibration Processing ===") if __name__ == "__main__": gmos_lsns_proc()