sdss / lvmdrp

Local Volume Mapper (LVM) Data Reduction Pipeline
BSD 3-Clause "New" or "Revised" License
2 stars 0 forks source link

lvmSFrame sky and sky_error extensions #81

Closed ymamay closed 2 months ago

ymamay commented 2 months ago

The current quick_sky_subtraction routine is saving the SkyE flux and error as the sky and sky_error extensions for the lvmSFrame. These should be the sky that was subtracted from each fiber in the lvmCFrame (skysub) and associated error (np.sqrt(error^2+skysub_error^2)).

In master branch, the current code is (note the TODOs):

# TODO - deal with error in sky-subtracted data ; replaced "error_c"
# TODO - deal properly with the error, sky, and sky_error
# this is incorrect
error_c = cframe._error
sky_c = cframe._sky_east
sky_error = cframe._sky_east_error

print("writing lvmSFrame")
sframe = lvmSFrame(data=skydata, error=error_c, mask=cframe._mask, sky=sky_c, sky_error=sky_error,
                   wave=cframe._wave, lsf=cframe._lsf, header=cframe._header, slitmap=cframe._slitmap)
kslong commented 2 months ago

I think we have several issues to sort out with regard to the code here. Here is more of the region that Amy cited in SkyMethod.

 # create spectrum for sky subtraction
    scisky = create_skysub_spectrum(sky_hdu, tel='sci', method=skymethod)
    skyesky = create_skysub_spectrum(sky_hdu, tel="skye", method=skymethod)
    skywsky = create_skysub_spectrum(sky_hdu, tel="skyw", method=skymethod)
   # select correct fibers
    data = cframe._data
    slitmap = Table(cframe._slitmap)
    scifibers = slitmap[slitmap["telescope"] == "Sci"]["fiberid"] - 1
    skyefibers = slitmap[slitmap["telescope"] == "SkyE"]["fiberid"] - 1
    skywfibers = slitmap[slitmap["telescope"] == "SkyW"]["fiberid"] - 1

    # sky subtract each spectrum
    skydata = data.copy()
    skydata[scifibers] = data[scifibers] - scisky
    skydata[skyefibers] = data[skyefibers] - skyesky
    skydata[skywfibers] = data[skywfibers] - skywsky

    # write out sky table to ancillary file
    mjd = cframe._header['MJD']
    expnum = cframe._header['EXPOSURE']
    tileid = cframe._header['TILE_ID']
    skytable = path.full('lvm_anc', mjd=mjd, tileid=tileid, drpver=drpver,
                         kind='sky', camera='brz', imagetype='table', expnum=expnum)

    sky_hdu.writeto(skytable, overwrite=True)

    # TODO - deal with error in sky-subtracted data ; replaced "error_c"
    # TODO - deal properly with the error, sky, and sky_error
    # this is incorrect
    error_c = cframe._error
    sky_c = cframe._sky_east
    sky_error = cframe._sky_east_error

    print("writing lvmSFrame")
    sframe = lvmSFrame(data=skydata, error=error_c, mask=cframe._mask, sky=sky_c, sky_error=sky_error,
                       wave=cframe._wave, lsf=cframe._lsf, header=cframe._header, slitmap=cframe._slitmap)
    sframe.writeFitsData(out_sframe)
    # TODO: check on expnum=7632 for halpha emission in sky fibers