save the trace shift offsets in the psf file

py/desispec/qa/qa_quicklook.py

@@ -341,7 +341,7 @@ def run_qa(self,image,inputs):
         from desispec.trace_shifts import compute_dy_using_boxcar_extraction as compute_dy
         fibers=np.arange(500) #RS: setting nfibers to 500 for now
         ox,oy,odx,oex,of,ol=compute_dx(xcoef,ycoef,wavemin,wavemax,image,fibers=fibers)
-        x_for_dy,y_for_dy,ody,ey,fiber_for_dy,wave_for_dy=compute_dy(psftrace,image,fibers)
+        x_for_dy,y_for_dy,ody,ey,fiber_for_dy,wave_for_dy,dwave,dwave_err=compute_dy(psftrace,image,fibers)
 
         # return average shifts in x and y
         dx=np.mean(odx)

py/desispec/scripts/qproc.py

@@ -17,7 +17,6 @@
 from desispec.io.fluxcalibration import read_average_flux_calibration
 from desispec.io.xytraceset import read_xytraceset,write_xytraceset
 import desispec.scripts.trace_shifts as trace_shifts_script
-from desispec.trace_shifts import write_traces_in_psf
 from desispec.calibfinder import CalibFinder
 from desispec.qproc.qframe import QFrame
 from desispec.qproc.io import read_qframe,write_qframe

py/desispec/scripts/trace_shifts.py

@@ -94,7 +94,15 @@ def read_specter_psf(filename) :
 
 
 def fit_trace_shifts(image, args):
-
+    """
+    Perform the fitting of shifts of spectral traces 
+    This consists of two steps, one is internal, by 
+    cross-correlating spectra to themselves, and then
+    cross-correlating to external (ususally sky) spectrum
+
+    Return updated traceset and two dictionaies with offset information 
+    to be written in the PSF file
+    """
     global psfs
 
     log=get_logger()
@@ -168,6 +176,7 @@ def fit_trace_shifts(image, args):
         nfibers = args.nfibers # FOR DEBUGGING
 
     fibers=np.arange(nfibers)
+    internal_offset_info = None
 
     if lines is not None :
 
@@ -198,11 +207,14 @@ def fit_trace_shifts(image, args):
         x_for_dx,y_for_dx,dx,ex,fiber_for_dx,wave_for_dx = compute_dx_from_cross_dispersion_profiles(xcoef,ycoef,wavemin,wavemax, image=image, fibers=fibers, width=args.width, deg=args.degxy,image_rebin=args.ccd_rows_rebin)
         if internal_wavelength_calib :
             # measure y shifts
-            x_for_dy,y_for_dy,dy,ey,fiber_for_dy,wave_for_dy = compute_dy_using_boxcar_extraction(tset, image=image, fibers=fibers, width=args.width, continuum_subtract=continuum_subtract)
+            x_for_dy,y_for_dy,dy,ey,fiber_for_dy,wave_for_dy,dwave,dwave_err  = compute_dy_using_boxcar_extraction(tset, image=image, fibers=fibers, width=args.width, continuum_subtract=continuum_subtract)
             mdy = np.median(dy)
             log.info("Subtract median(dy)={}".format(mdy))
             dy -= mdy # remove median, because this is an internal calibration
-
+            internal_offset_info = dict(wave=wave_for_dy,
+                                        fiber=fiber_for_dy,
+                                        dwave=dwave,
+                                        dwave_err=dwave_err)
         else :
             # duplicate dx results with zero shift to avoid write special case code below
             x_for_dy = x_for_dx.copy()
@@ -211,7 +223,7 @@ def fit_trace_shifts(image, args):
             ey       = 1.e-6*np.ones(ex.shape)
             fiber_for_dy = fiber_for_dx.copy()
             wave_for_dy  = wave_for_dx.copy()
-
+            
     degxx=args.degxx
     degxy=args.degxy
     degyx=args.degyx
@@ -341,11 +353,16 @@ def fit_trace_shifts(image, args):
         # the psf is used only to convolve the input spectrum
         # the traceset of the psf is not used here
         psf = read_specter_psf(args.psf)
-        tset.y_vs_wave_traceset._coeff = shift_ycoef_using_external_spectrum(psf=psf,xytraceset=tset,
-                                                                             image=image,fibers=fibers,
+        (tset.y_vs_wave_traceset._coeff,
+         wave_external, dwave_external, dwave_err_external) = shift_ycoef_using_external_spectrum(psf=psf, xytraceset=tset,
+                                                                             image=image, fibers=fibers,
                                                                              spectrum_filename=spectrum_filename,
-                                                                             degyy=args.degyy,width=7)
-
+                                                                             degyy=args.degyy, width=7)
+        external_offset_info = dict(wave=wave_external,
+                                    dwave=dwave_external,
+                                    dwave_err=dwave_err_external)
+    else:
+        external_offset_info = None
     x = np.zeros(x0.shape)
     y = np.zeros(x0.shape)
     for s in range(tset.nspec) :
@@ -361,7 +378,7 @@ def fit_trace_shifts(image, args):
     tset.meta["MINDY"]=np.min(dy)
     tset.meta["MAXDY"]=np.max(dy)
 
-    return tset
+    return tset, internal_offset_info, external_offset_info
 
 def main(args=None) :
 
@@ -382,10 +399,11 @@ def main(args=None) :
         log.critical(f"Entire {os.path.basename(args.image)} image is masked; can't fit traceshifts")
         sys.exit(1)
 
-    tset = fit_trace_shifts(image=image, args=args)
+    tset, internal_offset_info, external_offset_info = fit_trace_shifts(image=image, args=args)
     tset.meta['IN_PSF'] = shorten_filename(args.psf)
     tset.meta['IN_IMAGE'] = shorten_filename(args.image)
 
     if args.outpsf is not None :
-        write_traces_in_psf(args.psf,args.outpsf,tset)
+        write_traces_in_psf(args.psf,args.outpsf,tset, internal_offset_info=internal_offset_info,
+                            external_offset_info=external_offset_info)
         log.info("wrote modified PSF in %s"%args.outpsf)

py/desispec/trace_shifts.py

@@ -25,14 +25,17 @@
 from desispec.interpolation import resample_flux
 from desispec.qproc.qextract import qproc_boxcar_extraction
 
-def write_traces_in_psf(input_psf_filename,output_psf_filename,xytraceset) :
+def write_traces_in_psf(input_psf_filename,output_psf_filename,xytraceset, internal_offset_info=None,
+                        external_offset_info=None) :
     """
     Writes traces in a PSF.
 
     Args:
         input_psf_filename : Path to input fits file which has to contain XTRACE and YTRACE HDUs
         output_psf_filename : Path to output fits file which has to contain XTRACE and YTRACE HDUs
         xytraceset : xytraceset
+        internal_offset_info: dictionary of internal offsets (i.e. fiber vs 'median  fiber') in wavelength
+        external_offset_info: dictionary of external offsets (i.e. 'median fiber' vs external spectrum) in wavelength
     """
 
     xcoef=xytraceset.x_vs_wave_traceset._coeff
@@ -104,7 +107,24 @@ def write_traces_in_psf(input_psf_filename,output_psf_filename,xytraceset) :
     if (xytraceset.meta is not None) and ("PSF" in psf_fits):
         for k in xytraceset.meta.keys() :
             psf_fits["PSF"].header[k] = xytraceset.meta[k]
-
+    if internal_offset_info is not None:
+        data = {}
+        dwave, dwave_err, fiber, wave=[internal_offset_info[_] for _ in ['dwave','dwave_err','fiber','wave']]
+        data  = np.rec.fromarrays((fiber, wave, dwave, dwave_err),
+                                  dtype=np.dtype([('FIBER','i4'),
+                                         ('WAVE','f4'),
+                                         ('DWAVE','f4'),
+                                         ('DWAVE_ERR','f4')]))
+        psf_fits.append(pyfits.BinTableHDU(data, name='INTOFF'))
+    if external_offset_info is not None:
+        data = {}
+        dwave,dwave_err,wave=[external_offset_info[_] for _ in ['dwave','dwave_err','wave']]
+        data  = np.rec.fromarrays((wave, dwave, dwave_err),
+                                  dtype=np.dtype([
+                                         ('WAVE','f4'),
+                                         ('DWAVE','f4'),
+                                         ('DWAVE_ERR','f4')]))
+        psf_fits.append( pyfits.BinTableHDU(data, name='EXTOFF'))
 
     tmpfile = get_tempfilename(output_psf_filename)
     psf_fits.writeto(tmpfile, overwrite=True)
@@ -336,6 +356,8 @@ def compute_dy_from_spectral_cross_correlations_of_frame(flux, ivar, wave , xcoe
         ey : 1D array of uncertainties on dy
         fiber : 1D array of fiber ID (first fiber = 0)
         wave  : 1D array of wavelength
+        dwave : 1D array of wavelength offsets
+        dwave_err: 1D array of wavelength offset uncertainties
 
     """
     log=get_logger()
@@ -346,7 +368,9 @@ def compute_dy_from_spectral_cross_correlations_of_frame(flux, ivar, wave , xcoe
     ey=np.array([])
     fiber_for_dy=np.array([])
     wave_for_dy=np.array([])
-
+    dwave_list = np.array([])
+    dwave_err_list = np.array([])
+
     nfibers = flux.shape[0]
 
     for fiber in range(nfibers) :
@@ -375,7 +399,8 @@ def compute_dy_from_spectral_cross_correlations_of_frame(flux, ivar, wave , xcoe
 
             if err > 1 :
                 continue
-
+            dwave_list = np.append(dwave_list, dwave)
+            dwave_err_list = np.append(dwave_err_list, err)
             rw = legx(block_wave,wavemin,wavemax)
             tx = legval(rw,xcoef[fiber])
             ty = legval(rw,ycoef[fiber])
@@ -392,7 +417,7 @@ def compute_dy_from_spectral_cross_correlations_of_frame(flux, ivar, wave , xcoe
             fiber_for_dy=np.append(fiber_for_dy,fiber)
             wave_for_dy=np.append(wave_for_dy,block_wave)
 
-    return x_for_dy,y_for_dy,dy,ey,fiber_for_dy,wave_for_dy
+    return x_for_dy,y_for_dy,dy,ey,fiber_for_dy,wave_for_dy, dwave_list, dwave_err_list
 
 def compute_dy_using_boxcar_extraction(xytraceset, image, fibers, width=7, degyy=2,
                                            continuum_subtract = False):
@@ -444,7 +469,7 @@ def compute_dy_using_boxcar_extraction(xytraceset, image, fibers, width=7, degyy
     wavemax = xytraceset.wavemax
     xcoef   = xytraceset.x_vs_wave_traceset._coeff
     ycoef   = xytraceset.y_vs_wave_traceset._coeff
-
+    
     return compute_dy_from_spectral_cross_correlations_of_frame(flux=flux, ivar=ivar, wave=wave, xcoef=xcoef, ycoef=ycoef, wavemin=wavemin, wavemax=wavemax, reference_flux = mflux , n_wavelength_bins = degyy+4)
 
 @numba.jit
@@ -800,6 +825,9 @@ def shift_ycoef_using_external_spectrum(psf, xytraceset, image, fibers,
 
     Returns:
         ycoef  : 2D np.array of same shape as input, with modified Legendre coefficients for each fiber to convert wavelength to YCCD
+        wave:  : 1D array of wavelengths for which offsets are computed
+        dwave: : 1D array of wavelength offsets
+        dwave_err: 1D array of wavelength offset uncertainties
 
     """
     log = get_logger()
@@ -834,6 +862,8 @@ def shift_ycoef_using_external_spectrum(psf, xytraceset, image, fibers,
     dy = np.array([])
     ey = np.array([])
     wave_for_dy = np.array([])
+    dwave_list = np.array([])
+    dwave_err_list = np.array([])
     fiber_for_psf_evaluation = flux.shape[0] //2
     wavelength_bins = np.linspace(wave[0], wave[-1], n_wavelength_bins+1)
     for b in range(n_wavelength_bins) :
@@ -857,6 +887,8 @@ def shift_ycoef_using_external_spectrum(psf, xytraceset, image, fibers,
             dy = np.append(dy, -dwave * dydw)
             ey = np.append(ey, err*dydw)
             wave_for_dy = np.append(wave_for_dy,bin_wave)
+            dwave_list = np.append(dwave_list, dwave)
+            dwave_err_list = np.append(dwave_err_list, err)
             y_for_dy=np.append(y_for_dy,y)
             log.info(f"wave = {bin_wave}A , y={y}, measured dwave = {dwave} +- {err} A")
 
@@ -895,7 +927,8 @@ def shift_ycoef_using_external_spectrum(psf, xytraceset, image, fibers,
     for fiber in range(ycoef.shape[0]) :
         ycoef[fiber] += dycoef
 
-    return ycoef
+    return ycoef, wave_for_dy, dwave_list, dwave_err_list
+
 
 
 # end of routines for cross-correlation method for trace shifts