desihub / desispec

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

Unflagged cosmics in desi_preproc #762

Open londumas opened 5 years ago

londumas commented 5 years ago

The following data are using SP2, but the results are similar for SP0. I preproc using the following for the series of 10 zeros (199,200,201,202,203,204,205,207,208,209), i.e. not using the 206:

desi_preproc -i /exposures/desi/sps/20190227/00000208/sp2-00000*.fits.fz
--outdir /n/home/datasystems/users/hdumasde/20190227/00000208/
--cameras b2,r2,z2

Here is the histograms of data. You can clearly see that there are a lot of unflagged cosmics, spetially in r2 and z2.

histo_b2_zero

histo_r2_zero

histo_z2_zero

Looking at the images, most of them seem to be associated with flagged cosmics. Some are isolated though. In white good data with value<40, in violet masked data, in color the log of the value of the unflagged cosmics.

r2_202_cosmic_3

r2_202_cosmic_2

r2_202_cosmic

r2_199_cosmic_isolated

With the following code:

import fitsio
import scipy as sp
import matplotlib.pyplot as plt
import subprocess

from desispec.interpolation import resample_flux
from desispec.preproc import _parse_sec_keyword

def look_cosmics():

    ### image
    for c in ['b2','r2','z2']:
        for exp in range(199,210):
            if exp==206: continue
            path = '/n/home/datasystems/users/hdumasde/20190227/00000{}/preproc-{}-00000{}.fits'.format(exp,c,exp)
            h = fitsio.FITS(path)

            d = h['IMAGE'].read()
            w = h['MASK'].read()==0.
            ww = (h['MASK'].read()==0.) & (d<40.)
            www = (h['MASK'].read()==0.) & (d>=40.)
            if (w.sum()==0) or (www.sum()==0): continue

            i = sp.ones(d.shape[0])[:,None]*sp.arange(d.shape[1])
            j = sp.arange(d.shape[0])[:,None]*sp.ones(d.shape[1])
            print(list(zip(i[www],j[www])))

            td = d.copy()
            td[~w] = 1.
            td[ww] = sp.nan

            td[~sp.isnan(td)] = sp.log10(td[~sp.isnan(td)])

            plt.imshow(td)
            plt.colorbar()
            plt.title(c+' '+str(exp))
            plt.grid()
            plt.show()

    ### histo
    for c in ['b2','r2','z2']:
        for exp in range(199,210):
            if exp==206: continue
            path = '/n/home/datasystems/users/hdumasde/20190227/00000{}/preproc-{}-00000{}.fits'.format(exp,c,exp)
            h = fitsio.FITS(path)

            d = h['IMAGE'].read()
            w = h['MASK'].read()==0.
            if(w.sum()==0): continue

            plt.hist(d[w],bins=sp.arange(d[w].min(),d[w].max(),1.),histtype='step')
        plt.xlabel('IMAGE')
        plt.ylabel('#')
        plt.title(c)
        plt.yscale('log')
        plt.grid()
        plt.show()

    return

[EDIT] Adding the code to reproduce results.

londumas commented 5 years ago

Updated results using 50 zeros: [1159,1208]:

b0-50-bias

b2-50-bias

r0-50-bias

r2-50-bias

z0-50-bias

z2-50-bias

londumas commented 5 years ago

From @julienguy, maybe the answer is in updating the hardcodded parameters of the gradient: https://github.com/desihub/desispec/blob/36e5b82d40f06841cccd6a7a293d61c1c89fee77/py/desispec/cosmics.py#L364

julienguy commented 5 years ago

I worked a bit on this for the purpose of studying cosmics in MOSAIC CCDs. This is in branch improve_cosmics.

The following plot shows from top to bottom the fraction of remaining unmasked cosmic ray pixels, the rejection efficiency, and the signal masking rate that we want to keep as low as possible.

efficiency

The following shows the example for the NIR CCD image which combines exposure #1284 which is a 20min dark, with exposure #1450 and #1451 which are 10s and 1s arc lamp exposures. The first image is the data, the second the data with masked pixels in gray, and the last image the mask. One can see that most cosmics are masked while the data (the spectral lines) are not.

Screenshot from 2019-04-19 17-04-24