desihub / desispec

DESI spectral pipeline
BSD 3-Clause "New" or "Revised" License
36 stars 24 forks source link

Resolution matrix not normalized #959

Open dkirkby opened 4 years ago

dkirkby commented 4 years ago

The FRAME resolution matrix documented here does not appear to be normalized in the SV0 reductions I have checked so far.

Specifically, I expected the sum of the 11 resolution function bins tabulated for each (fiber, wlen) to be close to 1, except possibly for some leakage beyond the 11-bin window.

The actual sums have a mean close to 1 but a lot of dispersion that depends mostly on the wavelength (rather than the fiber number). There might also be some correlation of R > 1 with cosmics. The spread in sums is largest in z and smallest in b.

To demonstrate the issue, I used:

def frame_res_test(night='20200315', expid='00055589', spectro=9, camera='z', window=(0.8, 1.2)):
    path = f'/global/cfs/cdirs/desi/spectro/redux/daily/exposures/{night}/{expid}/frame-{camera}{spectro}-{expid}.fits'
    resdata = fitsio.read(path, ext='RESOLUTION')
    Rnorm = resdata.sum(axis=1)
    fig, ax = plt.subplots(1, 2, figsize=(12, 5))
    img = ax[0].imshow(Rnorm, origin='lower', aspect='auto', vmin=window[0], vmax=window[1])
    plt.colorbar(img, ax=ax[0]).set_label('R sum')
    ax[0].set_xlabel('Wavelength bin')
    ax[0].set_ylabel('Fiber number')
    ax[1].hist(Rnorm.reshape(-1), bins=np.linspace(*window, 101))
    ax[1].set_xlim(*window)
    ax[1].set_xlabel('R sum')
    ax[1].set_yticks([])
    plt.tight_layout()
    plt.savefig(f'res_test_{night}_{expid}_{spectro}_{camera}.png')

which produces: res_test_20200315_00055589_9_z

sbailey commented 4 years ago

The non-normalization of R is a "feature" of the spectroperfectionism math from Bolton & Schelgel 2009, roughly analogous to a weighted- vs. unweighted- fit and how they treat low-weight outlier points. I think you will find the largest Rsum outliers to be near noisy sky lines, which is why it is a bigger effect in Z than B.

The normalization comes from whether you sum over rows of columns in Bolton & Schelgel 2009 equations 11 and 12. If you swap the normalization axis (ex2d.py line 532) you can get something that does conserve flux in R with the somewhat surprising (but perhaps better) side effect of oscillating noise: neighboring bins of the extracted spectra can have radically different noise levels in a ringing-like oscillatory manner. But at least R is a real convolution matrix then. We explored that a bit years ago (see section 4.2.2 of DESI-0890) but didn't converge on a detailed study and specter currently implements the normalization as defined in the original paper.

All that being said, the effect you show here is larger than I expected, and perhaps more work could be done on

  1. Ringing noise could be better (or less surprising?) than R not really being a convolution matrix. I don't think we've come to grips with the consequences in either direction.
  2. The size of this might be a side effect of some poorly constrained eigenvectors in the matrix algebra, and could be better constrained.