Fixing 1 km/s radial velocity systematic

[segasai]

Hi,

One of the long-standing issues we notice in MWS when looking at RVs wrt external data is the offset of ~ 1km/s that has been present for long time. (see i.e. figure

from
MWS EDR VAC paper K+2024).

We have already done some improvements in RV accuracy with trace_shifts #2106 (but mostly for backup data), but I thought I'd look at it again. And I've looked at the MEANDY from trace_shifts across the whole kibo, and this is how it looks like

There is a clear bias of 0.03 pixels. This corresponds to about 1km/s.
The sign is also in the right direction for +1km/s offst see figures in the #2106 for the plots of meandy vs velocity offset.

From what I can see the most likely cause is the sky spectrum sky-spec.dat, rather than anything algorithmic (that's my hypothesis though, not the fact). So my proposal/plan is to try to either

1. substitute sky-spec.dat by similar synthesized spectrum bazed on Hanuschik+2003 atlas
2. update the algorithm compute_dy_from_spectral_cross_correlation() to something a bit more sophisticated where a sky spectrum is not compared to a fixed template, but for each offset we on the fly solve the line amplitudes (given the line-list from Hanuschik and LSF).

In #2106 I briefly looked at Hanuschik, but then decided to deal with that later as it was tangential to the issue with backup data and large MEANDY offsets.

I don't currently have a patch, but have written a code somewhere for 2nd method. I'll aim to try it on some subset of data to see if that gets rid of the offset and/or also improves the meandy accuracy.

This is just an issue to keep track of this problem, avoid forgetting it, as well as solicit any help/feedback.

[segasai]

A quick attempt of rerunning trace_shifts with new Hanuschik synthesized sky template.
Using only lines with >.1 flux from Hanuschik and their fluxes and synthesized with sigma=0.15 A on the same grid as the existing sky-spec.dat
To my disappointment that does not seem to be improving things much, just maybe a bit smaller shift...

One possible theory is that the line blends from weaker lines are producing this, or that
line ratios needs to be changed or alternatively it's not the problem with the sky template at all.

(this was run on a random subset of exposures from kibo)

(UPDATE to avoid spamming)

These are dwave vs wavelength as returned by shift_ycoef_using_external_spectrum
The left one is the one with the new spectrum, the right is original one.

The x-axis is the wavelength of the wavelength bin

desispec/py/desispec/trace_shifts.py

Line 757 in 4eae355

for b in range(n_wavelength_bins) :

and the y-axis is the measurement by compute_dy_from_spectral_cross_correlation

desispec/py/desispec/trace_shifts.py

Line 767 in 4eae355

dwave,err = compute_dy_from_spectral_cross_correlation(mflux[ok],

(focusing on the >8000 region)
We can clearly see persistent systematics here that different wavelength bins are shifted in the same direction.

(UPDATE2)
Previously I've been told that the provenance of the spec-sky.dat is not quite clear, but apparently it was based on Hanuschik+2003 convolved to desi resolution (according to Guy+2023 paper), hence it's not surprising we are getting the same results.

This makes me think that the root cause is not the sky-spectrum in fact.

(UPDATE3)
Now investigating whether the PSF assymmetries could be the cause here.
This shows the Average(y) for PSF - center of the PSF (in pix) as a function wavelength for one frame
(all 10 petals and for set of fibers) (the code is below)
The conclusion from this is that there are considerable assymmetries of roughly correct amplitude and possibly showing similar trends as in one of the plots above. So it is possible this is one of the reasons for the offsets

import numpy as np
import matplotlib.pyplot as plt
import desispec.scripts.trace_shifts

plt.clf()
for petal in range(10):
    plt.subplot(3, 4, petal + 1)
    PP = desispec.scripts.trace_shifts.read_specter_psf(
        f'/global/cfs/cdirs/desi/spectro/redux/kibo/exposures/20231008/00200048/psf-z{petal}-00200048.fits'
    )
    ws = np.linspace(7520, 9824, 20)
    fibers = np.arange(1, 499, 20)
    colors = (fibers - fibers.min()) / fibers.ptp()
    for i, fiber in enumerate(fibers):
        sk = []
        for w in ws:
            x1 = PP.x(fiber, w)
            y1 = PP.y(fiber, w)
            yg = np.linspace(y1 - 40, y1 + 40, 800)
            x, y = np.meshgrid(np.linspace(x1 - 20, x1 + 20, 400), yg)
            q = (PP._value(x, y, fiber, w).sum(axis=1))
            mean = (yg * q).sum()
            '''std = (q*(yg-mean)**2).sum()**.5'''
            '''sk.append((q*(yg-mean)**3).sum()/std**3 )'''
            sk.append(mean - y1)
        plt.plot(ws, sk, color=plt.cm.turbo(colors[i]))
    plt.xlabel('Wavelength')

[julienguy]

@segasai , thank for digging into this.

• Yes the sky spectrum is actually using the Hanuschik catalog.
• I also looked at that and came to the same conclusion that offsets from line to line are probably caused by blends.
• We can have genuine wavelength calibration systematics because of tree rings on the CCDs, it would be interesting to look at the mean dY per spectrograph.
• One complication: we also have a PCA correction for the sky lines (that partially fixes the tree rings issue). This could be part of the problem. One had to look at the values in the HDU DWAVE_MEAN of the files

ls $DESI_SPECTRO_CALIB/spec/sm*/skycorr-pca-sm*.fits

[segasai]

Thanks for the suggestions/comments @julienguy
I was not aware of the PCA correction, I'll look into it and what it does.
Regarding the dependence on the spectrograph, as far as I can see the meandy is basically the same across petals.

Also I still cannot get my head around why the dy offsets show roughly constant pattern as a f-n of wavelength (one of the plots above). My first thought was it was related to psf asymmetries/dependence on wavelength, but it does not seem to be the case, since the improvement in the #2386 do not fix that.

[segasai]

Just looked at the DWAVE_MEAN in skycorr-pca for different spectrographs (see below) and to first order I do not see the same patterns observed in DWAVE from trace_shifts, so this skycorr-pca is probably not the root-cause of the offsets.

[segasai]

I will reopen this issue, as while #2386 should have improved things, but it does not fully fix this. The root cause for 1km/s offset is still unclear.
As a next step I'd like to have #2411 implemented so we can properly monitor the offsets and then hopefully fix them.