Edge glitch in template due to data resolution

[londumas]

The following plot shows the issue linked to applying the data resolution to the templates: spec.R.dot(mx). You can see that at each edge of the three spectrographs, we get a sudden drop of flux. We might be mixing zeros with data, or something of the kind.

Code to reproduce:

import os
import subprocess
from astropy.io import fits
import matplotlib.pyplot as plt
import scipy as sp

import desisim.templates
import redrock
from redrock.external import desi

def resolution_glitch():

    obj = "qso"
    simdir = os.getenv('SCRATCH') + '/desi/Resolution_glitch/'
    infile      = simdir + obj +'-input-spectra.fits'
    specoutfile = simdir + obj +'-observed-spectra.fits'

    ###
    qso_maker = desisim.templates.QSO()
    flux, wave, meta = qso_maker.make_templates(nmodel=1,seed=42)
    zTrue = meta['REDSHIFT'][0]

    ###
    os.makedirs(simdir, exist_ok=True)
    hdr = fits.Header()
    hdr['EXTNAME'] = 'WAVELENGTH'
    hdr['BUNIT'] = 'Angstrom'
    fits.writeto(infile, wave, header=hdr, overwrite=True)
    hdr['EXTNAME'] = 'FLUX'
    hdr['BUNIT'] = '10^-17 ergs/s/cm2/A'
    fits.append(infile, flux, header=hdr)

    ###
    cmd  = " quickspectra --input " + infile
    cmd += " --out-spectra        " + specoutfile
    cmd += " --seed 42"
    subprocess.call(cmd, shell=True)

    ### Data
    targets = desi.DistTargetsDESI(specoutfile)._my_data
    target = targets[0]
    specs_to_read = target.spectra

    #- Templates
    templates_path = redrock.templates.find_templates()
    templates = {}
    for el in templates_path:
        t = redrock.templates.Template(filename=el)
        templates[t.full_type] = t
    tp = templates['QSO']

    ###
    for i, spec in enumerate(specs_to_read):
        coef = sp.zeros(tp.flux.shape[0])
        coef[0] = 1
        mx = tp.eval(coef, spec.wave, zTrue) * (1.+zTrue)
        if i==0:
            plt.plot(spec.wave, mx, alpha=0.8,marker='o',markersize=2,color='black',label=r"$\mathrm{Template}$")
        else:
            plt.plot(spec.wave, mx, alpha=0.8,marker='o',markersize=2,color='black')
        model = spec.R.dot(mx)
        plt.plot(spec.wave, model, alpha=0.8,marker='o',markersize=2,label=r"$\mathrm{Template+Resolution}$")
    plt.xlabel(r"$\lambda_{\mathrm{Obs.}}$",fontsize=20)
    plt.ylabel(r"$Flux$",fontsize=20)
    plt.legend(fontsize=20, numpoints=1,ncol=1, loc=1)
    plt.grid()
    plt.show()

    return

resolution_glitch()

[londumas]

The bug is present in the fit and in the ploting. It might be the cause of the lines in https://github.com/desihub/redrock/issues/85

[ngbusca]

@londumas it looks like the resolution matrix is not normalized at the edges of the cameras, and there's some ambiguity on how to normalize there. A possible protection against this would be to check where spec.R.sum(axis=1) starts deviating from 1 and mask those bins.

[londumas]

@ngbusca, @sbailey, if we keep data only when spec.R.sum(axis=1)>0.99 we loose only the first two pixels and the last two. Does that seem reasonable to you ?

[sbailey]

@londumas yes, that seems reasonable. I also confirmed on "real" extractions of pixsim data that this would only lose 1-2 pixels on each end of each spectrum. Please make sure that those are not included in the NPIXELS output.