Also scan for large trace shifts in continuum/LED exposures

py/desispec/scripts/proc.py

@@ -558,7 +558,7 @@ def main(args=None, comm=None):
                     cmd += " --psf {}".format(inpsf)
                     cmd += " --degxx 2 --degxy 0"
                     if args.obstype in ['FLAT', 'TESTFLAT', 'TWILIGHT'] :
-                        cmd += " --continuum"
+                        cmd += " --continuum --no-large-shift-scan"
                     else :
                         cmd += " --degyx 2 --degyy 0"
                     if args.obstype in ['SCIENCE', 'SKY']:

py/desispec/scripts/trace_shifts.py

@@ -22,7 +22,7 @@
 from desispec.io import shorten_filename
 from desiutil.log import get_logger
 from desispec.util import parse_int_args
-from desispec.trace_shifts import write_traces_in_psf,compute_dx_from_cross_dispersion_profiles,compute_dy_from_spectral_cross_correlation,monomials,polynomial_fit,compute_dy_using_boxcar_extraction,compute_dx_dy_using_psf,shift_ycoef_using_external_spectrum,recompute_legendre_coefficients,recompute_legendre_coefficients_for_x,recompute_legendre_coefficients_for_y,list_of_expected_spot_positions
+from desispec.trace_shifts import write_traces_in_psf,compute_dx_from_cross_dispersion_profiles,compute_dy_from_spectral_cross_correlation,monomials,polynomial_fit,compute_dy_using_boxcar_extraction,compute_dx_dy_using_psf,shift_ycoef_using_external_spectrum,recompute_legendre_coefficients,recompute_legendre_coefficients_for_x,recompute_legendre_coefficients_for_y,list_of_expected_spot_positions,compute_x_offset_from_central_band_cross_dispersion_profile
 from desispec.large_trace_shifts import detect_spots_in_image,match_same_system
 
 def parse(options=None):
@@ -60,7 +60,7 @@ def parse(options=None):
     parser.add_argument('--degyy', type = int, default = 0, required=False,
                         help = 'polynomial degree for y shifts along y')
     parser.add_argument('--continuum', action='store_true',
-                        help = 'only fit shifts along x for continuum input image')
+                        help = 'only fit shifts along x for continuum or LED input image')
     parser.add_argument('--auto', action='store_true',
                         help = 'choose best method (sky,continuum or just internal calib) from the FLAVOR keyword in the input image header')
 
@@ -72,7 +72,8 @@ def parse(options=None):
                         help = "width of cross-dispersion profile")
     parser.add_argument('--ccd-rows-rebin', type = int, default = 4 , required=False,
                         help = "rebinning of CCD rows to run faster")
-
+    parser.add_argument('--no-large-shift-scan', action="store_true",
+                        help = "do not perform a large shift scan for arc lamp or continuum exposures")
     args = parser.parse_args(options)
 
     return args
@@ -459,15 +460,17 @@ def fit_trace_shifts(image, args):
             args.sky = True
         log.info("wavelength calib, internal={}, sky={} , arc_lamps={}".format(internal_wavelength_calib,args.sky,args.arc_lamps))
 
-    if args.arc_lamps :
+    cfinder = CalibFinder([image.meta])
+    fibers=np.arange(tset.nspec)
+    if cfinder.haskey("BROKENFIBERS") :
+        brokenfibers=parse_int_args(cfinder.value("BROKENFIBERS"))%500
+        log.debug(f"brokenfibers={brokenfibers}")
+        fibers=fibers[np.isin(fibers, brokenfibers, invert=True)]
+
+    if args.arc_lamps and (not args.no_large_shift_scan) :
 
         log.info("for arc lamps, find a first solution by comparing expected spots positions with detections over the whole CCD")
-        cfinder = CalibFinder([image.meta])
-        fibers=np.arange(tset.nspec)
-        if cfinder.haskey("BROKENFIBERS") :
-            brokenfibers=parse_int_args(cfinder.value("BROKENFIBERS"))%500
-            log.debug(f"brokenfibers={brokenfibers}")
-            fibers=fibers[np.isin(fibers, brokenfibers, invert=True)]
+
         xref,yref = list_of_expected_spot_positions(tset,fibers)
         xin,yin   = detect_spots_in_image(image)
 
@@ -510,6 +513,12 @@ def fit_trace_shifts(image, args):
         xcoef[:,0] += delta_xref
         ycoef[:,0] += delta_yref
 
+    if args.continuum  and (not args.no_large_shift_scan) :
+        log.info("for continuum or LED exposures, find a first solution on delta_X by comparing a wide cross-dispersion profile with expectations")
+        delta_xref = compute_x_offset_from_central_band_cross_dispersion_profile(tset, image, fibers=fibers)
+        log.info(f"apply best shift delta x = {delta_xref} to traceset")
+        xcoef[:,0] += delta_xref
+
     spectrum_filename = args.spectrum
     continuum_subtract = False
     if args.sky :
@@ -557,6 +566,7 @@ def fit_trace_shifts(image, args):
 
     else :
 
+
         # internal calibration method that does not use the psf
         # nor a prior set of lines. this method is much faster
 

py/desispec/trace_shifts.py

@@ -1507,3 +1507,42 @@ def list_of_expected_spot_positions(traceset,fibers=None,max_number_of_lines=50)
     xspot=np.hstack(xspot)
     yspot=np.hstack(yspot)
     return xspot,yspot
+
+
+def compute_x_offset_from_central_band_cross_dispersion_profile(tset, image, fibers=None,halfwidth=50) :
+    log=get_logger()
+    bands=7 # measure delta_x in 7 horizontal bands of 100 pixel wide each and return the median offset
+    bandwidth=100
+    hb=bands//2
+    n0=image.pix.shape[0]
+    n1=image.pix.shape[1]
+    ivar=image.ivar*(image.mask==0)
+    if fibers is None :
+        fibers = np.arange(tset.nspec)
+    delta_x_array=[]
+    for b in range(-hb,hb+1) :
+        begin=int(n0//2+(b-0.5)*bandwidth)
+        end=begin+bandwidth
+
+        sw=np.sum(ivar[begin:end,:],axis=0)
+        swf=np.sum(image.pix[begin:end,:]*ivar[begin:end,:])
+        prof=(swf/(sw+(sw==0)))
+        y=int(n0//2+b*bandwidth)
+        xc=np.array([tset.x_vs_y(fiber,y) for fiber in fibers])
+        oversample=4
+        ox=np.arange(n1*oversample)/oversample
+        xx=np.arange(n1)
+        sigma=1.5 # pixel, approximate
+        mprof=np.exp(-(ox[:,None]-xc[None,:])**2/2./sigma**2).sum(axis=1)
+        norm = 1./np.sqrt(np.sum(prof**2)*(np.sum(mprof**2)/oversample))
+        ddx=np.linspace(-halfwidth,halfwidth,4*halfwidth+1)
+        xcorr=np.zeros(ddx.size)
+        for i,dx in enumerate(ddx) :
+            xcorr[i]=np.sum(prof*np.interp(xx,ox+dx,mprof))
+        xcorr *= norm
+        imax=np.argmax(xcorr)
+        log.info("horizontal band {}:{} delta x = {:.1f} pixels".format(begin,end,ddx[imax]))
+        delta_x_array.append(ddx[imax])
+    delta_x = np.median(delta_x_array)
+    log.info("median delta x = {:.1f} pixels".format(delta_x))
+    return delta_x