# GMOS Data Reduction Cookbook companion script to the chapter: # "Reduction of IFU Spectral Images with IRAF" # # IRAF CL command script to: # Process IFU exposures for J2240-0927, in programs GS-2012B-Q-26 and GN-2012B-Q-226. # 2016-Apr-06 shaw@noao.edu # # 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: # Define tasks: # task lacos_spec=/path/to/lacos_spec.cl # task gscrspec=/path/to/gscrspec.cl # Then: # cd /path/to/science_files print ("### Begin Processing GMOS/IFU Spectra ###") print ("###") cd ./raw #### Select files for processing # print ("===Exosure Selection===") print (" --Selecting Bias Exposures--") string obsSelect string ccdSelect # All bias & flat files match a single CCD RoI and binning. ccdSelect = "&& detrO1ys>1024 && ccdsum?='1 1' " print (" --Selecting Bias Exposures--") # # Bias exposures with a common instrument, observation type & class: obsSelect = "instrume?='GMOS-N' && obstype?='BIAS' && obsclass?='dayCal' " # Select bias exposures within ~10 days of the target observations. s2 = " && @'date-obs' > '2012-10-09' && @'date-obs' < '2012-10-26'" s1 = obsSelect // ccdSelect // s2 # Select using information from both PHDU and HDU-1 hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "bias_tmp.txt") # Now remove the IRAF index and kernel parameters from the file names string omit = "exten,index,kernel" gemextn ("@bias_tmp.txt", omit=omit, outfile="biasFilesN.txt") # Now for GMOS-S # Taken w/science targets obsSelect = "instrume?='GMOS-S' && obstype?='BIAS' && obsclass?='dayCal' " s2 = " && @'date-obs' > '2012-08-20' && @'date-obs' < '2012-09-01'" s1 = obsSelect // ccdSelect // s2 hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "bias_tmp.txt") gemextn ("@bias_tmp.txt", omit=omit, outfile="biasFilesS1.txt") # Taken w/std program at a much later date. s2 = " && @'date-obs' > '2012-10-26' && @'date-obs' < '2012-11-03'" s1 = obsSelect // ccdSelect // s2 hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "bias_tmp.txt") gemextn ("@bias_tmp.txt", omit=omit, outfile="biasFilesS2.txt") print (" --Selecting Flat-field Expsoures--") string bandSelect # Start with GMOS-N flats # GCAL B600/499: obsSelect = "instrume?='GMOS-N' && obstype?='FLAT' " bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5307' && grwlen=499." s1 = obsSelect // ccdSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "flt_tmp.txt") gemextn ("@flt_tmp.txt", omit=omit, outfile="flatGcal_B6-499.txt") # R831/853: bandSelect = "&& maskname?='IFU-2' && grating?='R831+_G5302' && grwlen=853." s1 = obsSelect // ccdSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "flt_tmp.txt") gemextn ("@flt_tmp.txt", omit=omit, outfile="flatGcal_R8-853.txt") # GMOS-S flats # GCAL B600/625: bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5323' && grwlen=625." obsSelect = "instrume?='GMOS-S' && obstype?='FLAT' " s1 = obsSelect // ccdSelect // bandSelect hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "flt_tmp.txt") gemextn ("@flt_tmp.txt", omit=omit, outfile="flatGcal_B6-625.txt") dele ("*_tmp.txt") rename ("bias*.txt,flat*.txt", "../") ### Comparison arcs print (" --Selecting comparison arcs--") ## Arcs # Only one arc is useful for each configuration, so choose filenames from the log print ("N20121005S0902.fits", > "arc_B6-499.txt") print ("S20120828S0005.fits", > "arc_B6-625.txt") print ("N20121018S0197.fits", > "arc_R8-853.txt") ## Science exposures print (" --Selecting science exposures--") # Select via object name, observation class, mask, grating & cenWave. # GMOS-N: obsSelect = "instrume?='GMOS-N' && i_title?='J2240-0927' && obsclass?='science' " # R831-853 bandSelect = "&& maskname?='IFU-2' && grating?='R831+_G5302' && grwlen=853." s1 = obsSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "sci_tmp.txt") gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_R8-853.txt") # B600-499 bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5307' && grwlen=499." s1 = obsSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "sci_tmp.txt") gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_B6-499.txt") # GMOS-S: obsSelect = "instrume?='GMOS-S' && i_title?='J2240-0927' && obsclass?='science' " # B600-625 bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5323' && grwlen=625." s1 = obsSelect // bandSelect hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "sci_tmp.txt") gemextn ("@sci_tmp.txt", omit=omit, outfile="sci_B6-625.txt") ## Standard stars print (" --Selecting standard star exposures--") # Select via observation type & class, mask, grating & cenWave. ## GMOS-N, program Cal: obsSelect = "instrume?='GMOS-N' && obstype?='OBJECT' && obsclass?='progCal' " # B600-499 bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5307' && grwlen=499." s1 = obsSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "std_tmp.txt" gemextn ("@std_tmp.txt", omit=omit, outfile="std_B6-499.txt") # R831-853 bandSelect = "&& maskname?='IFU-2' && grating?='R831+_G5302' && grwlen=853." s1 = obsSelect // bandSelect hselect ("N*.fits[1,inherit=yes]", "$I", s1, > "std_tmp.txt" gemextn ("@std_tmp.txt", omit=omit, outfile="std_R8-853.txt") ## GMOS-S: # B600-625 bandSelect = "&& maskname?='IFU-2' && grating?='B600+_G5323' && grwlen=625." # EG131 obsSelect = "instrume?='GMOS-S' && obstype?='OBJECT' && i_title?='EG131' " s1 = obsSelect // bandSelect hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "std_tmp.txt" gemextn ("@std_tmp.txt", omit=omit, outfile="EG131_B6-625.txt") # LTT9239 obsSelect = "instrume?='GMOS-S' && obstype?='OBJECT' && i_title?='LTT9239' " s1 = obsSelect // bandSelect hselect ("S*.fits[1,inherit=yes]", "$I", s1, > "std_tmp.txt" gemextn ("@std_tmp.txt", omit=omit, outfile="LTT9239_B6-625.txt") ## Move the file lists to the parent directory, and clean up. dele ("*_tmp.txt") rename ("arc*.txt,std*.txt,sci*.txt", "../") rename ("EG131_B6-625.txt,LTT9239_B6-625.txt", "../") print ("===Creating Bias MasterCals===") cd .. ## Create Bias Residual and Flat-field MasterCals, including VAR and DQ arrays. print (" --Creating Bias Residual MasterCals--") unlearn gbias unlearn gemextn # Disarm a bug in gbias gbias.logfile="gbiasLog.txt" gbias.rawpath="./raw" gbias.verbose=no print (" -GMOS-N Bias-") gbias ("@biasfilesN.txt", "MCbiasN.fits", fl_vardq+) print (" -GMOS-S Bias 1-") gbias ("@biasfilesS1.txt", "MCbiasS1.fits", fl_vardq+) print (" -GMOS-S Bias 2-") gbias ("@biasfilesS2.txt", "MCbiasS2.fits", fl_vardq+) # Clean up imdele gS*.fits,gN*.fits print (" --Creating Trace MasterCals--") # # Use primarily the default task parameters, which include: # fl_addmdf+ fl_bias+ fl_extract+ fl_fluxcal+ fl_fixgaps- fl_fulldq- # fl_gnsskysub- fl_gsappwave+ fl_gscrrej+ fl_nodsuffle- fl_novlap+ # fl_over+ fl_skysub+ fl_trim+ fl_wavtran+ fl_vardq- # Use primarily the default task parameters. unlearn gfreduce unlearn gireduce gfreduce.logfile="gfreduceLog.txt" # Note the trailing slash, which is for some reason required for gfreduce: gfreduce.rawpath="./raw/" gfreduce.verbose=no # Construct a list of output trace files from contemporaneous flat names # Start with basic reductions for GCAL flat-fields. gfreduce ("@flatGcal_B6-625-1.txt", bias="MCbiasS1", fl_extract+, \ fl_gscrrej-, fl_gsappwave-, fl_wavtran-, fl_skysub-, fl_fluxcal-, \ fl_vardq+, fl_inter-) gfreduce ("@flatGcal_B6-625-2.txt", bias="MCbiasS2", fl_extract+, \ fl_gscrrej-, fl_gsappwave-, fl_wavtran-, fl_skysub-, fl_fluxcal-, \ fl_vardq+, fl_inter-) # Insert the static BPM in-place. # The addmdf task seems not to work, so insert the BPM by hand. string bpm = "bpm_gmos-s_EEV_v2_1x1_spec_MEF.fits" string extn = "[DQ," sections ("rg//@flatGcal_B6-625.txt", option="root", >"rgFlat_B6-625.txt") list = "rgFlat_B6-625.txt" while (fscan (list, s1) != EOF) { for (i=1; i<=3; i+=1) { s2=bpm//extn//i//"]" imcopy (s2, s1//extn//i//",overwrite+]") } } ## Interpolate over bad columns unlearn gemfix gemfix.logfile="gemfixLog.txt" gemfix ("@rgFlat_B6-625.txt", "p@rgFlat_B6-625.txt", method="fit1d", bitmask=1, order=32) ## Extract flats from interpolated exposures gfreduce.rawpath="./" gfreduce ("p@rgFlat_B6-625.txt", fl_addmdf-, fl_over-, fl_trim-, fl_bias-, \ fl_extract+, fl_gscrrej-, fl_gsappwave+, fl_wavtran-, fl_skysub-, \ fl_fluxcal-, fl_vardq+, fl_inter-) # Clean up imdele ("gS2012*.fits,gN2012*.fits") ## Subtracting scattered light, then re-extracting flats, skipped. ## Apply QE correction here ## Normalize flats to create the response functions. unlearn gfresponse gfresponse.verbose=no gfresponse.logfile="gfresponseLog.txt" sections ("ep//@rgFlat_B6-625.txt", option="root", >"eprgFlats_B6-625.txt") list = "eprgFlats_B6-625.txt" while (fscan (list, s1) != EOF) { gfresponse (s1, s1//'_flat', skyimage="", function='spline3', order=47, sample='1:1,30:2855,2879:2879') } ### Process the Arcs. print (" --Creating Wavelength Calibration--") # Only need one trace file for this. gfreduce.rawpath="./raw/" gfreduce ("@arc_B6-625.txt", reference="eprgS20120829S0062.fits", \ fl_bias-, fl_extract+, trace-, recen-, fl_gsappwave+, fl_wavtran-, \ fl_gscrrej-, fl_skysub-, fl_fluxcal-, fl_vardq-, fl_inter-) # Perform wavelength calibration unlearn gswavelength gswavelength.logfile="gswaveLog.txt" gswavelength.verbose=no # Hi-res linelist is problematic because of the weakly exposed dayCals. Use default. gswavelength ("erg@arc_B6-625.txt", fl_inter-, fwidth=8, minsep=2.5) print ("=== General Basic Processing ===") # # Create a list of all Std and Sci files to process concatenate ("EG131_B6-625.txt,s*_B6-625.txt", "allSfiles.txt", out_type="in_type") concatenate ("s*_B6-499.txt,s*_R8-853.txt", "allNfiles.txt", out_type="in_type") gfreduce.rawpath="./raw/" gfreduce.verbose=no #gfreduce ("@allNfiles.txt", bias="MCbiasN", reference="eprgS20120829S00XX", \ #fl_extract-, fl_gscrrej-, fl_gsappwave-, fl_wavtran-, fl_skysub-, \ #fl_fluxcal-, fl_vardq+, fl_inter-) gfreduce ("@allSfiles.txt", bias="MCbiasS1", reference="eprgS20120829S0062", \ fl_extract-, fl_gscrrej-, fl_gsappwave-, fl_wavtran-, fl_skysub-, \ fl_fluxcal-, fl_vardq+, fl_inter-) gfreduce ("@LTT9239_B6-625.txt", bias="MCbiasS2", reference="eprgS20121031S0025.fits", \ fl_extract-, fl_gscrrej-, fl_gsappwave-, fl_wavtran-, fl_skysub-, \ fl_fluxcal-, fl_vardq+, fl_inter-) # Insert the updated GMOS-S BPM in-place, by hand. concatenate ("allSfiles.txt,LTT9239_B6-625.txt", "allFiles.txt", out_type="in_type") sections ("rg//@allFiles.txt", option="root", >"rgAllfiles.txt") list = "rgAllfiles.txt" while (fscan (list, s1) != EOF) { for (i=1; i<=3; i+=1) { s2=bpm//extn//i//"]" imcopy (s2, s1//extn//i//",overwrite+]") } } # Clean up imdele ("gS2012*.fits,gN2012*.fits") ### Artifact rejection ### This will take a bloody long time... print (" -Removing CRs-") # Clean cosmic rays. unlearn gemcrspec gemcrspec.logfile="gemcrspecLog.txt" gemcrspec ("@rgAllfiles.txt", "x@rgAllfiles.txt", sigfrac=0.32, niter=3, fl_vardq+) # Replace bad pixels and flagged CRs with interpolated values. unlearn gemfix gemfix.logfile="gemfixLog.txt" gemfix.verbose=yes gemfix ("x@rgAllfiles.txt", "px@rgAllfiles.txt", method="fixpix") ## Perform QE correction. # [TBD] imdele ("rg*.fits") # Resume Std processing to apply flat-field & wavelength calibrations, and subtract sky. # print ("=== Standard Star Processing ===") gfreduce.rawpath="./" sections ("pxrg@EG131_B6-625.txt", >"cleanEG131_B6-625.txt") sections ("pxrg@LTT9239_B6-625.txt", >"cleanLTT9239_B6-625.txt") gfreduce ("@cleanEG131_B6-625.txt", fl_addmdf-, fl_over-, fl_trim-, fl_bias-, \ fl_extract+, reference="eprgS20120829S0062", trace-, recenter-, \ response="eprgS20120829S0062_flat", \ fl_gsappwave+, fl_wavtran+, wavtraname="ergS20120828S0005", \ w1=5618., w2=INDEF, dw=0.4622, nw=2822, \ fl_skysub+, sepslits+, fl_fluxcal-, fl_gscrrej-, fl_vardq+, fl_inter-) gfreduce ("@cleanLTT9239_B6-625.txt", fl_addmdf-, fl_over-, fl_trim-, fl_bias-, \ fl_extract+, reference="eprgS20121031S0025", trace-, recenter-, \ response="eprgS20121031S0025_flat", \ fl_gsappwave+, fl_wavtran+, wavtraname="ergS20120828S0005", \ w1=5618., w2=INDEF, dw=0.4622, nw=2822, \ fl_skysub+, sepslits+, fl_fluxcal-, fl_gscrrej-, fl_vardq+, fl_inter-) # Combine all the fibers from the object field into a 1-D spectrum. unlearn gfapsum gfapsum.logfile="gfapsumLog.txt" sections ("ste@cleanEG131_B6-625.txt", option="root", >"steEG131_B6-625.txt") list = "steEG131_B6-625.txt" while (fscan (list, s1) != EOF) { gfapsum (s1, lthreshold=0., fl_inter-) } copy astepxrgS20120829S0061.fits EG131sum_B6-625.fits sections ("ste@cleanLTT9239_B6-625.txt", option="root", >"steLTT9239_B6-625.txt") list = "steLTT9239_B6-625.txt" while (fscan (list, s1) != EOF) { gfapsum (s1, lthreshold=0., fl_inter-) } print (" -Combining spectra for LTT9239-") unlearn gemcombine gemcombine.logfile="gemcombineLog.txt" sections ("a@steLTT9239_B6-625.txt", option="root", >"asteLTT9239_B6-625.txt") gemcombine ("@asteLTT9239_B6-625.txt", "LTT9239sum_B6-625", reject="none", scale="exposure") print (" -Deriving the sensitivity function-") # Derive the sensitivity function. # Use the Mauna Kea extinction function for GMOS-N. copy onedstds$ctionewcal/l9239.dat ./ copy gmos$calib/eg131.dat ./ unlearn gsstandard gsstandard.logfile="gsstandardLog.txt" gsstandard.observatory="Gemini-South" gsstandard ("EG131sum_B6-625.fits,LTT9239sum_B6-625.fits", "std_B6-625", \ "Sens_B6-625", starname="eg131,l9239", caldir="./", \ extinction="onedstds$ctioextinct.dat", order=11, fl_inter+) #### Process science targets # print ("=== Science Target Processing ===") print (" --Extract spectra & Apply Calibrations-- ") gfreduce.rawpath="./" sections ("pxrg@sci_B6-625.txt", >"cleanSci_B6-625.txt") gfreduce ("@cleanSci_B6-625.txt", fl_addmdf-, fl_over-, fl_trim-, fl_bias-, fl_extract+, reference="eprgS20120827S0069", trace-, recenter-, \ response="eprgS20120827S0069_flat", \ fl_gsappwave+, fl_wavtran+, wavtraname="ergS20120828S0005", \ w1=5618., w2=INDEF, dw=0.4622, nw=2822, \ fl_skysub+, sepslits+, fl_fluxcal+, sfunction="Sens_B6-625", \ extinction="onedstds$ctioextinct.dat", \ fl_gscrrej-, fl_vardq+, fl_inter-) # ...and similarly for other grating settings. print (" --Wavelength Zeropoint Correction-- ") unlearn rvidlines rv.observatory="gemini-south" rvidlines.coordlist="skylines.txt" rvidlines.ftype="emission" rvidlines.logfile="rvLog.txt" rvidlines ("stepxrgS20120827S0066.fits[SKY]", threshold=7., \ nsum=1, maxfeatures=10, fwidth=10., cradius=10., threshold=10., minsep=5.) # Results borrowed from JT measurements: sections ("cste//@cleanSci_B6-625.txt//.fits[SCI]", >"calibSci_B6-625.txt") unlearn hedit hedit.update=yes hedit.verify=no hedit ("cstepxrgS20120827S0066.fits[sci]", "CRVAL1", 5618.164) hedit ("cstepxrgS20120827S0067.fits[sci]", "CRVAL1", 5618.268) hedit ("cstepxrgS20120827S0068.fits[sci]", "CRVAL1", 5618.387) hedit ("cstepxrgS20120827S0070.fits[sci]", "CRVAL1", 5618.560) hedit ("cstepxrgS20120827S0071.fits[sci]", "CRVAL1", 5618.667) hedit ("cstepxrgS20120827S0072.fits[sci]", "CRVAL1", 5618.787) print (" --Advanced Product Generation-- ") # Create (x,y,lambda) datacube. unlearn gfcube gfcube.logfile="gfcubeLog.txt" list = "calibSci_B6-625.txt" while (fscan (list, s1) != EOF) { gfcube (s1, fl_atmdisp+, fl_flux+, fl_var+, fl_dq+) } # Adjust the spatial WCS reference values. # The offsets were derived by James Turner using his pyfalign software. hedit ("dcstepxrgS20120827S0070.fits[sci]", "CRVAL1", 65.48042, update+, verify-) hedit ("dcstepxrgS20120827S0071.fits[sci]", "CRVAL1", 65.47042, update+, verify-) hedit ("dcstepxrgS20120827S0071.fits[sci]", "CRVAL2", 0.045000, update+, verify-) hedit ("dcstepxrgS20120827S0072.fits[sci]", "CRVAL1", 65.50042, update+, verify-) hedit ("dcstepxrgS20120827S0072.fits[sci]", "CRVAL2", 0.035000, update+, verify-) # Simple combine of datacubes imcombine ("dcstepxrgS20120827S00*.fits[SCI]", "j2240-097.fits") ## FIN