Source code for kcwidrp.primitives.MakeInvsens

from keckdrpframework.primitives.base_primitive import BasePrimitive
from kcwidrp.core.bokeh_plotting import bokeh_plot
from kcwidrp.core.kcwi_correct_extin import kcwi_correct_extin
from kcwidrp.core.kcwi_get_std import kcwi_get_std
from kcwidrp.core.kcwi_plotting import set_plot_lims, save_plot
from kcwidrp.primitives.kcwi_file_primitives import kcwi_fits_writer, \
    kcwi_fits_reader, strip_fname

from bokeh.plotting import figure, ColumnDataSource
from scipy.signal import find_peaks
from scipy.interpolate import interp1d
from astropy.io import fits as pf
from astropy.nddata import CCDData
from astropy import units as u

import time
import os
import numpy as np


[docs]class MakeInvsens(BasePrimitive): """ Generate inverse sensitivity curve from a standard star observation. Uses object name to determine if the observation is a standard star. Then checks that the image has been processed through DAR correction. Reads in the reference spectrum and compares it with the observed spectrum to generate the inverse sensitivity curve that can be used to flux calibrate science observations. Outputs FITS inverse sensitivity spectrum along with diagnostic plots that include residuals between reference spectrum and calibrated observed spectrum and effective area and efficiency plots as a function of wavelength. """ def __init__(self, action, context): BasePrimitive.__init__(self, action, context) self.logger = context.pipeline_logger def _pre_condition(self): """Checks if object is a standard star""" self.logger.info("Checking precondition for MakeInvsens") stdfile = None stdname = None if 'object' in self.action.args.imtype.lower(): self.logger.info("Checking OBJECT keyword") stdfile, stdname = kcwi_get_std( self.action.args.ccddata.header['OBJECT'], self.logger) if not stdfile: self.logger.info("Checking TARGNAME keyword") stdfile, stdname = kcwi_get_std( self.action.args.ccddata.header['TARGNAME'], self.logger) else: self.logger.warning("Not object type: %s" % self.action.args.imtype) self.action.args.stdfile = stdfile self.action.args.stdname = stdname # have we been processed correctly? if 'DARCOR' in self.action.args.ccddata.header: if self.action.args.ccddata.header['DARCOR']: # does the standard file exist? if self.action.args.stdfile is not None: # have we been stacked? if 'NSTACK' in self.action.args.ccddata.header: nstack = self.action.args.ccddata.header['NSTACK'] else: nstack = 0 # does file already exist? ofn = self.action.args.name msname = strip_fname(ofn) + '_invsens.fits' rdir = self.config.instrument.output_directory invsensf = os.path.join(self.config.instrument.cwd, rdir, msname) if os.path.exists(invsensf): if nstack <= 0: self.logger.warning("Master cal already exists: %s" % invsensf) return False else: os.unlink(invsensf) self.logger.info("Master cal will re-generated " "from stacked image") return True else: self.logger.info("Master cal will be generated.") return True else: self.logger.warning("Not a KCWI standard star observation.") return False else: self.logger.warning("DAR not corrected, cannot generate " "inverse sensitivity") else: self.logger.warning("Not processed enough to generate " "inverse sensitivity") return False def _perform(self): self.logger.info("Making inverse sensitivity curve") suffix = 'invsens' stdname = self.action.args.stdname do_plots = self.config.instrument.plot_level >= 3 # get size sz = self.action.args.ccddata.data.shape # default pixel ranges z = np.arange(sz[0]) z0 = 175 z1 = sz[0] - 175 # get exposure time expt = self.action.args.ccddata.header['XPOSURE'] if expt == 0.: self.logger.warning("No exposure time found, setting to 1s") expt = 1. else: self.logger.info("Using exposure time of %.1f" % expt) # get wavelength scale w0 = self.action.args.ccddata.header['CRVAL3'] dw = self.action.args.ccddata.header['CD3_3'] crpixw = self.action.args.ccddata.header['CRPIX3'] # get all good wavelength range wgoo0 = self.action.args.ccddata.header['WAVGOOD0'] if wgoo0 < 3650: wgoo0 = 3650. wgoo1 = self.action.args.ccddata.header['WAVGOOD1'] # get all-inclusive wavelength range wall0 = self.action.args.ccddata.header['WAVALL0'] wall1 = self.action.args.ccddata.header['WAVALL1'] # get DAR padding in y pad_y = self.action.args.ccddata.header['DARPADY'] # get sky subtraction status if 'SKYCOR' in self.action.args.ccddata.header: skycor = self.action.args.ccddata.header['SKYCOR'] else: skycor = False # get telescope and atm. correction if 'TELESCOP' in self.action.args.ccddata.header: tel = self.action.args.ccddata.header['TELESCOP'] else: tel = 'KeckI' if 'Keck' in tel: area = 760000.0 else: area = -1.0 # compute good y pixel ranges if w0 > 0. and dw > 0. and wgoo0 > 0. and wgoo1 > 0.: z0 = int((wgoo0 - w0) / dw) + 10 z1 = int((wgoo1 - w0) / dw) - 10 # wavelength scale w = w0 + z * dw # good spatial range gy0 = pad_y if pad_y > 1 else 1 gy1 = sz[1] - (pad_y if pad_y > 2 else 2) # log results self.logger.info("Invsen. Pars: Y0, Y1, Z0, Z1, Wav0, Wav1: " "%d, %d, %d, %d, %.3f, %.3f" % (gy0, gy1, z0, z1, w[z0], w[z1])) # central wavelength cwv = self.action.args.cwave # find standard in slices # sum over wavelength tot = np.sum(self.action.args.ccddata.data[z0:z1, gy0:gy1, :], 0) # yy = np.arange(gy1-gy0) + gy0 mxsl = -1. mxsg = 0. # for each slice for i in range(sz[2]): tstd = float(np.nanstd(tot[:, i])) if tstd > mxsg: mxsg = tstd mxsl = i # plot slice data, if requested if do_plots: py = tot[:, i] px = np.arange(len(py)) p = figure( title=self.action.args.stdlabel + ' Std Slice %d (DAR padded)' % i, x_axis_label="Position along slice", y_axis_label="Flux summed over WLs", plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) p.scatter(px, py, marker='x') bokeh_plot(p, self.context.bokeh_session) print("sl, std: %02d, %.3f" % (i, tstd)) qstr = input("Next? <cr>, q - quit: ") if 'Q' in qstr.upper(): do_plots = False # relevant slices if self.action.args.ifunum == 1: # Large slicer slset = 3 elif self.action.args.ifunum == 2: # Medium slicer slset = 5 elif self.action.args.ifunum == 3: # Small slicer slset = 12 else: slset = 3 self.logger.error("Assuming Large slicer: id undefined! - %d" % self.action.args.ifunum) # Get set of slices needed for calculation sl0 = (mxsl - slset) if mxsl >= slset else 0 sl1 = (mxsl + slset) if (mxsl + slset) <= sz[2]-1 else sz[2]-1 # get y position of std cy, _ = find_peaks(tot[:, mxsl], height=np.nanmean(tot[:, mxsl])) cy = int(cy[0]) + gy0 # log results self.logger.info("Std slices (DAR padded): " "max, sl0, sl1, spatial cntrd: " "%d, %d, %d, %.2f" % (mxsl, sl0, sl1, cy)) # get dwave spectrum ofn = self.action.args.name delfn = strip_fname(ofn) + '_dcubed.fits' full_path = os.path.join( self.config.instrument.cwd, self.config.instrument.output_directory, delfn) if os.path.exists(full_path): dew = kcwi_fits_reader(full_path)[0] dwspec = dew.data[:, cy, mxsl] zeros = np.where(dwspec == 0) if len(zeros) > 0: dwspec[zeros] = dw else: dwspec = np.zeros(sz[0]) + dw # copy of input cube scub = self.action.args.ccddata.data.copy() # sky window width in pixels skywin = int(self.config.instrument.psfwid / self.action.args.xbinsize) # do sky subtraction, if needed if not skycor: self.logger.warning("Sky should have been subtraced already") # apply extinction correction ucub = scub.copy() # uncorrected cube kcwi_correct_extin(scub, self.action.args.ccddata.header, logger=self.logger) # get slice spectra y limits sy0 = (cy - skywin) if (cy - skywin) > 0 else 0 sy1 = (cy + skywin) if (cy + skywin) < (sz[1]-1) else (sz[1]-1) # sum over y range slspec = np.sum(scub[:, sy0:sy1, :], 1) ulspec = np.sum(ucub[:, sy0:sy1, :], 1) # sum over slices obsspec = np.sum(slspec[:, sl0:sl1], 1) ubsspec = np.sum(ulspec[:, sl0:sl1], 1) # convert to e-/second obsspec /= expt ubsspec /= expt # check for zeros zeros = np.where(obsspec == 0.) if len(zeros) > 0: obsmean = np.nanmean(obsspec) obsspec[zeros] = obsmean # read in standard star spectrum stdfile = self.action.args.stdfile hdul = pf.open(stdfile) swl = hdul[1].data['WAVELENGTH'] sflx = hdul[1].data['FLUX'] sfw = hdul[1].data['FWHM'] hdul.close() # get region of interest sroi = [i for i, v in enumerate(swl) if w[0] <= v <= w[-1]] nsroi = len(sroi) if nsroi <= 0: self.logger.error("No standard wavelengths in common") return self.action.args # expand range after checking for edges if sroi[0] > 0: sroi.insert(0, sroi[0]-1) nsroi += 1 if sroi[-1] < (len(swl)-1): sroi.append(sroi[-1]+1) nsroi += 1 # how many points? self.logger.info("Number of standard points = %d" % nsroi) # very sparsely sampled w.r.t. object if nsroi <= 5: self.logger.error("Not enough standard points") return self.action.args swl = swl[sroi] sflx = sflx[sroi] sfw = sfw[sroi] fwhm = np.max(sfw) self.logger.info("Reference spectrum FWHM used = %.1f (A)" % fwhm) # resample standard onto our wavelength grid rsint = interp1d(swl, sflx, kind='cubic', fill_value='extrapolate') rsflx = rsint(w) # get effective inverse sensitivity invsen = rsflx / obsspec # convert to photons/s/cm^2/(wl bin = dw) rspho = 5.03411250e7 * rsflx * w * dw # get effective area earea = ubsspec / rspho # correct to native bins earea *= dw/dwspec # Line masks lmasks = [] # Header for line mask file lmhdr = [] # default values (for BM) ford = 9 # fit order if 'BL' in self.action.args.grating: ford = 7 elif 'BH' in self.action.args.grating: ford = 9 # Adjust for dichroic fraction try: dichroic_fraction = self.action.args.ccddata.header['DICHFRAC'] except KeyError: dichroic_fraction = 1. ford = int(ford * dichroic_fraction) if ford < 3: ford = 3 self.logger.info("Fitting Invsens with polynomial order %d" % ford) # fit inverse sensitivity and effective area # get initial region of interest wl_good = [i for i, v in enumerate(w) if wgoo0 <= v <= wgoo1] nwl_good = len(wl_good) if nwl_good <= 0: self.logger.error("No good wavelengths to fit") return self.action.args wlm0 = wgoo0 wlm1 = wgoo1 # interactively set wavelength limits if self.config.instrument.plot_level >= 1: print("CHECKING WAVELENGTH LIMITS") print("Current WL limits: %.1f - %.1f Angstroms " "(blue vertical lines)" % (wlm0, wlm1)) print("A <cr> will accept current values, " "or enter new values to avoid edge problems (if present).") print("Hover the cursor over spectrum to find wavelengths.") # yran = [np.min(obsspec), np.max(obsspec)] yran = [np.min(obsspec[wl_good]), np.max(obsspec[wl_good])] source = ColumnDataSource(data=dict(x=w, y=obsspec)) done = False while not done: p = figure( tooltips=[("x", "@x{0,0.0}"), ("y", "@y{0,0.0}")], title=self.action.args.stdlabel + ' Obs Spec', x_axis_label='Wave (A)', y_axis_label='Intensity (e-)', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) p.line('x', 'y', line_color='black', source=source) p.line([wgoo0, wgoo0], yran, line_color='green', legend_label='WAVGOOD') p.line([wgoo1, wgoo1], yran, line_color='green') p.line([wlm0, wlm0], yran, line_color='blue', legend_label='LIMITS') p.line([wlm1, wlm1], yran, line_color='blue') p.line([cwv, cwv], yran, line_color='red', legend_label='CWAV') set_plot_lims(p, xlim=[wall0, wall1], ylim=yran) bokeh_plot(p, self.context.bokeh_session) qstr = input("New values? <float> <float> (or <cr> - done): ") if len(qstr) <= 0: done = True else: try: wlm0 = float(qstr.split()[0]) wlm1 = float(qstr.split()[1]) if wlm1 < wlm0 or wlm0 < wall0 or wlm1 > wall1: wlm0 = wgoo0 wlm1 = wgoo1 print("range/order error, try again") except (IndexError, ValueError): wlm0 = wgoo0 wlm1 = wgoo1 print("format error, try again") # update region of interest wl_good = [i for i, v in enumerate(w) if wlm0 <= v <= wlm1] nwl_good = len(wl_good) # END: interactively set wavelength limits # Look for mask data file # first in local directory local_lmfile = os.path.basename(stdfile).split('.fits')[0] + '.lmsk' if not os.path.exists(local_lmfile): # then in stds directory lmfile = stdfile.split('.fits')[0] + '.lmsk' if not os.path.exists(lmfile): lmfile = None else: lmfile = local_lmfile # if we found a mask file, read it in if lmfile is None: self.logger.info("No line mask file found") else: self.logger.info("Using line mask file %s" % lmfile) with open(lmfile) as lmf: lmasks_str = lmf.readlines() # parse file into mask list for lmws in lmasks_str: lmws = lmws.strip() # Collect header lines if lmws.startswith('#'): lmhdr.append(lmws) continue # Skip blank lines if len(lmws.strip()) < 1: continue try: data = lmws.split('#')[0] # Parse comment on line if '#' in lmws: comm = lmws.split('#')[1].lstrip() else: comm = '' # Parse line mask range lm = [float(v) for v in data.split()] except ValueError: print("bad line: %s" % lmws) continue # Only collect good lines if len(lm) == 2 and lm[0] < lm[1]: lmdict = {'w0': lm[0], 'w1': lm[1], 'com': comm} lmasks.append(lmdict) else: print("bad line: %s" % lmws) # Now interactively identify lines if requested if self.config.instrument.plot_level >= 1: yran = [np.min(obsspec[wl_good]), np.max(obsspec[wl_good])] # source = ColumnDataSource(data=dict(x=w, y=obsspec)) print("MASKING SHARP FEATURES: ABSORPTION LINES/COSMIC RAYS") print("To mask, enter starting and stopping wavelengths " "followed by an optional comment (string) for each line " "you want to mask.") print("Current masks are shown as vertical dashed yellow lines.") print("Hover the cursor over spectrum to find wavelengths.") print("Mask values and comments are saved to this file: %s." % local_lmfile) print("A <cr> with no line limits will accept current masks.") done = False while not done: p = figure( tooltips=[("x", "@x{0.0}"), ("y", "@y{0.0}")], title=self.action.args.stdlabel + ' Obs Spec', x_axis_label='Wave (A)', y_axis_label='Intensity (e-)', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) p.line(w, obsspec, line_color='black') p.line([wgoo0, wgoo0], yran, line_color='green', legend_label='WAVGOOD') p.line([wgoo1, wgoo1], yran, line_color='green') p.line([wlm0, wlm0], yran, line_color='blue', legend_label='LIMITS') p.line([wlm1, wlm1], yran, line_color='blue') p.line([cwv, cwv], yran, line_color='red', legend_label='CWAV') for ml in lmasks: if wall0 < ml['w0'] < wall1 and wall0 < ml['w1'] < wall1: p.line([ml['w0'], ml['w0']], yran, line_color='orange', line_dash='dashed') p.line([ml['w1'], ml['w1']], yran, line_color='orange', line_dash='dashed') set_plot_lims(p, xlim=[wall0, wall1], ylim=yran) bokeh_plot(p, self.context.bokeh_session) qstr = input("Mask line: wavelength start stop (Ang) comment " "(str) <float> <float> <str> (<cr> - done): ") if len(qstr) <= 0: done = True else: try: newl = {'w0': float(qstr.split()[0]), 'w1': float(qstr.split()[1]), 'com': " ".join(qstr.split()[2:]).strip()} # newl = [float(val) for val in qstr.split()] except ValueError: print("bad line: %s" % qstr) continue if wlm0 < newl['w0'] < wlm1 and wlm0 < newl['w1'] < wlm1: lmasks.append(newl) else: print("line mask outside wl range: %s" % qstr) # END: interactively identify lines # Write out a local copy of line masks so user can edit if len(lmasks) > 0: local_lmfile = os.path.basename( stdfile).split('.fits')[0] + '.lmsk' with open(local_lmfile, 'w') as lmf: for lmh in lmhdr: lmf.write(lmh) for lm in lmasks: if len(lm['com']) == 0: lmf.write("%.2f %.2f\n" % (lm['w0'], lm['w1'])) else: lmf.write("%.2f %.2f # %s\n" % (lm['w0'], lm['w1'], lm['com'])) # set up fitting vectors, flux, waves, measure errors sf = invsen[wl_good] # dependent variable af = earea[wl_good] # effective area wf = w[wl_good] # independent variable bf = obsspec[wl_good] # input electrons rf = rsflx[wl_good] # reference flux mw = np.ones(nwl_good, dtype=float) # weights use = np.ones(nwl_good, dtype=int) # toggles for usage # loop over line masks and apply them for ml in lmasks: roi = [i for i, v in enumerate(wf) if ml['w0'] <= v <= ml['w1']] nroi = len(roi) if nroi > 0: use[roi] = 0 self.logger.info("Masking line at %.2f to %.2f (A)" % (ml['w0'], ml['w1'])) # ignore bad points by setting large errors mf = [] ef = [] ww = [] used = [] not_used = [] for i in range(len(use)): if use[i] == 1: used.append(i) else: mw[i] = 1.e-9 mf.append(sf[i]) ef.append(100.*af[i]/area) ww.append(wf[i]) not_used.append(i) # initial polynomial fit of inverse sensitivity wf0 = np.min(wf) res = np.polyfit(wf-wf0, sf, deg=ford, w=mw) finvsen = np.polyval(res, w-wf0) sinvsen = np.polyval(res, wf-wf0) calspec = obsspec * finvsen scalspec = bf * sinvsen # initial polynomial fit of effective area res = np.polyfit(wf-wf0, af, ford, w=mw) fearea = np.polyval(res, w-wf0) # calculate residuals resid = 100.0 * (scalspec - rf) / rf if len(not_used) > 0: rbad = resid[not_used] else: rbad = None rsd_mean = float(np.nanmean(resid[used])) rsd_stdv = float(np.nanstd(resid[used])) self.logger.info("Calibration residuals = %f +- %f %%" % (rsd_mean, rsd_stdv)) # plots peff = None pivs = None pcal = None prsd = None # interactively adjust fit if self.config.instrument.plot_level >= 1: done = False while not done: yran = [np.min(100.*af/area), np.max(100.*af/area)] effmax = np.nanmax(100.*fearea/area) effmean = np.nanmean(100.*fearea/area) peff = figure( title=self.action.args.stdlabel + ' Efficiency', x_axis_label='Wave (A)', y_axis_label='Effective Efficiency (%)', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) peff.line(wf, 100.*af/area, line_color='black', legend_label='Data') peff.line(w, 100.*fearea/area, line_color='red', legend_label='Fit') peff.scatter(ww, ef, marker='x', legend_label='Rej') peff.line([wlm0, wlm0], yran, line_color='green', legend_label='WL lims') peff.line([wlm1, wlm1], yran, line_color='green') peff.line([wall0, wall1], [effmax, effmax], line_color='black', line_dash='dashed') peff.line([wall0, wall1], [effmean, effmean], line_color='black', line_dash='dashdot') set_plot_lims(peff, xlim=[wall0, wall1], ylim=yran) bokeh_plot(peff, self.context.bokeh_session) if self.config.instrument.plot_level >= 1: input("Next? <cr>: ") else: time.sleep(2. * self.config.instrument.plot_pause) yran = [np.min(sf), np.max(sf)] pivs = figure( title=self.action.args.stdlabel + ' Inverse sensitivity', x_axis_label='Wave (A)', y_axis_label='Invserse Sensitivity (Flux/e-/s)', y_axis_type='log', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) pivs.line(wf, sf, line_color='black', legend_label='Data') pivs.line(w, finvsen, line_color='red', legend_label='Fit') pivs.scatter(ww, mf, marker='x', legend_label='Rej') pivs.line([wlm0, wlm0], yran, line_color='green', legend_label='WL lims') pivs.line([wlm1, wlm1], yran, line_color='green') set_plot_lims(pivs, xlim=[wall0, wall1], ylim=yran) bokeh_plot(pivs, self.context.bokeh_session) if self.config.instrument.plot_level >= 1: input("Next? <cr>: ") else: time.sleep(2. * self.config.instrument.plot_pause) yran = [np.min(calspec[wl_good]), np.max(calspec[wl_good])] pcal = figure( title=self.action.args.stdlabel + ' Calibrated', x_axis_label='Wave (A)', y_axis_label='Flux (ergs/s/cm^2/A)', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) pcal.line(w, calspec, line_color='black', legend_label='Obs') pcal.line(w, rsflx, line_color='red', legend_label='Ref') pcal.line([wlm0, wlm0], yran, line_color='green', legend_label='WL lims') pcal.line([wlm1, wlm1], yran, line_color='green') set_plot_lims(pcal, xlim=[wall0, wall1], ylim=yran) bokeh_plot(pcal, self.context.bokeh_session) if self.config.instrument.plot_level >= 1: input("Next? <cr>: ") else: time.sleep(2. * self.config.instrument.plot_pause) yran = [np.min(resid), np.max(resid)] prsd = figure( title=self.action.args.stdlabel + ' Residuals = %.1f +- %.1f (%%)' % (rsd_mean, rsd_stdv), x_axis_label='Wave (A)', y_axis_label='Obs - Ref / Ref (%)', plot_width=self.config.instrument.plot_width, plot_height=self.config.instrument.plot_height) prsd.line(wf, resid, line_color='black', legend_label='Obs - Ref (%)') if len(not_used) > 0: prsd.scatter(ww, rbad, marker='x', legend_label='Rej') prsd.line([wlm0, wlm0], yran, line_color='green', legend_label='WL lims') prsd.line([wlm1, wlm1], yran, line_color='green') prsd.line([wall0, wall1], [rsd_mean, rsd_mean], line_color='red') prsd.line([wall0, wall1], [rsd_mean+rsd_stdv, rsd_mean+rsd_stdv], line_color='black', line_dash='dashed') prsd.line([wall0, wall1], [rsd_mean-rsd_stdv, rsd_mean-rsd_stdv], line_color='black', line_dash='dashed') set_plot_lims(prsd, xlim=[wall0, wall1], ylim=yran) bokeh_plot(prsd, self.context.bokeh_session) if self.config.instrument.plot_level >= 1: qstr = input("Current fit order = %d, " "New fit order? <int>, <cr> - done: " % ford) if len(qstr) <= 0: done = True else: try: ford = int(qstr) # update fit of inverse sensitivity res = np.polyfit(wf - wf0, sf, deg=ford, w=mw) finvsen = np.polyval(res, w - wf0) sinvsen = np.polyval(res, wf - wf0) calspec = obsspec * finvsen scalspec = bf * sinvsen # update polynomial fit of effective area res = np.polyfit(wf - wf0, af, ford, w=mw) fearea = np.polyval(res, w - wf0) # re-calculate residuals resid = 100.0 * (scalspec - rf) / rf if len(not_used) > 0: rbad = resid[not_used] else: rbad = None rsd_mean = float(np.nanmean(resid[used])) rsd_stdv = float(np.nanstd(resid[used])) self.logger.info( "Calibration residuals = %f +- %f %%" % (rsd_mean, rsd_stdv)) except ValueError: print("Bad fit order, try again") else: done = True time.sleep(2. * self.config.instrument.plot_pause) # log results effmax = float(np.nanmax(100. * fearea / area)) effmean = float(np.nanmean(100. * fearea / area)) self.logger.info("Peak, mean efficiency (%%): %.1f, %.1f" % (effmax, effmean)) self.logger.info("Fit order = %d" % ford) # output plots pfname = "std_%05d_%s_%s_%s_%d" % ( self.action.args.ccddata.header['FRAMENO'], stdname, self.action.args.grating.strip(), self.action.args.ifuname.strip(), int(self.action.args.cwave)) # Save plots save_plot(peff, filename=pfname + '_eff.png') # Efficiency save_plot(pivs, filename=pfname + '_invsens.png') # Inv. Sens. save_plot(prsd, filename=pfname + '_resid.png') # Residuals save_plot(pcal, filename=pfname + '_cal.png') # Calibrated log_string = MakeInvsens.__module__ self.action.args.ccddata.header['HISTORY'] = log_string self.logger.info(log_string) # write out effective inverse sensitivity # update inverse sensitivity header hdr = self.action.args.ccddata.header.copy() hdr['HISTORY'] = log_string hdr['INVSENS'] = (True, 'effective inv. sens. spectrum?') hdr['INVSW0'] = (w[z0], 'low wavelength for eff. inv. sens.') hdr['INVSW1'] = (w[z1], 'high wavelength for eff. inv. sens.') hdr['INVSZ0'] = (z0, 'low wave pixel for eff. inv. sens.') hdr['INVSZ1'] = (z1, 'high wave pixel for eff. inv. sens.') hdr['INVSY0'] = (gy0, 'low spatial pixel for eff. inv. sens.') hdr['INVSY1'] = (gy1, 'high spatial pixel for eff. inv. sens.') hdr['INVSLMX'] = (mxsl, 'brightest std star slice') hdr['INVSL0'] = (sl0, 'lowest std star slice summed') hdr['INVSL1'] = (sl1, 'highest std star slice summed') hdr['INVSLY'] = (cy, 'spatial pixel position of std within slice') hdr['INVFW0'] = (wlm0, 'low wavelength for fits') hdr['INVFW1'] = (wlm1, 'high wavelength for fits') hdr['INVFORD'] = (ford, 'fit order') hdr['EXPTIME'] = (1., 'effective exposure time (seconds)') hdr['XPOSURE'] = (1., 'effective exposure time (seconds)') # remove old WCS del hdr['RADESYS'] del hdr['EQUINOX'] del hdr['LONPOLE'] del hdr['LATPOLE'] del hdr['NAXIS2'] if 'NAXIS3' in hdr: del hdr['NAXIS3'] del hdr['CTYPE1'] del hdr['CTYPE2'] del hdr['CTYPE3'] del hdr['CUNIT1'] del hdr['CUNIT2'] del hdr['CUNIT3'] del hdr['CNAME1'] del hdr['CNAME2'] del hdr['CNAME3'] del hdr['CRVAL1'] del hdr['CRVAL2'] del hdr['CRVAL3'] del hdr['CRPIX1'] del hdr['CRPIX2'] del hdr['CRPIX3'] del hdr['CD1_1'] del hdr['CD1_2'] del hdr['CD2_1'] del hdr['CD2_2'] del hdr['CD3_3'] # set wavelength axis WCS values hdr['WCSDIM'] = 1 hdr['CTYPE1'] = ('AWAV', 'Air Wavelengths') hdr['CUNIT1'] = ('Angstrom', 'Wavelength units') hdr['CNAME1'] = ('KCWI INVSENS Wavelength', 'Wavelength name') hdr['CRVAL1'] = (w0, 'Wavelength zeropoint') hdr['CRPIX1'] = (crpixw, 'Wavelength reference pixel') hdr['CDELT1'] = (dw, 'Wavelength Angstroms per pixel') ofn = self.action.args.name invsname = strip_fname(ofn) + '_' + suffix + '.fits' eaname = strip_fname(ofn) + '_ea.fits' # set units invsens_u = u.erg / (u.angstrom * u.cm ** 2 * u.s * u.electron) # output inverse sensitivity out_invsens = CCDData(np.asarray([invsen, finvsen, obsspec]), meta=hdr, unit=invsens_u) kcwi_fits_writer(out_invsens, output_file=invsname, output_dir=self.config.instrument.output_directory) self.context.proctab.update_proctab(frame=out_invsens, suffix=suffix, newtype='INVSENS', filename=self.action.args.name) # output effective area ea_u = u.cm ** 2 / u.angstrom out_ea = CCDData(np.asarray([earea, fearea]), meta=hdr, unit=ea_u) kcwi_fits_writer(out_ea, output_file=eaname, output_dir=self.config.instrument.output_directory) self.context.proctab.update_proctab(frame=out_ea, suffix='ea', newtype='EAREA', filename=self.action.args.name) return self.action.args
# END: class MakeInvsens()