Algorithms

Algorithms to be implemented for the IRIS imager and Integral Field Spectrograph. Once the actual classes are implemented in liger_iris_pipeline, we will just link their implementations.

Generate master dark

The master dark is generated by the median of 5-10 individual dark frames taken at the same exposure. This removes any of the frame-to-frame noise variation. Individual dark frames are zero illumination calibration frames with a shutter or blocking filter to not allow any light to hit the detector.

# 4096x4096 at 100 seconds exposures, N-frames
# dark1, dark2, dark3, …, darkN
darks = np.array([dark1, dark2, dark3, darkN])
master_dark = np.median(darks, axis=0)

Dark subtraction

Each science and calibration frame will have dark current, residual electric current that is flowing through the array when there is no light. Dark current increases with time, so all science and calibration frames need equivalent exposure dark frames. To remove the dark current, we subtract a single dark in the real-time or a master dark in the final pipeline, at the same exposure time as the science and calibration frames.

dark_subtracted_science = science - dark
dark_subtracted_flat = flat - dark

Remove detector artifacts

The detector will have two types of artifacts; permanent/semi-permanent and transient artifacts. Dead pixels, hot pixels or “frozen” pixel fall with the permanent/semi-permanent artifacts while cosmic rays (CR) are within the transient artifacts. Dead pixel are pixels that no longer function. Hot pixels are sensitive pixels that can have a non-linear response to incoming photons. They can also turn on and off and can even respond normally until a certain flux is achieved in which they become non-linear. “Frozen” pixels are pixels that have a low response rate (opposite of the hot pixels), with similar types of problems. CRs are high energy photons from the sky that can hit the detector randomly and leave bright artifacts.

For all types of detector artifacts, they are generally difficult to remove or remove completely with flat-fielding. To deal with detector artifacts, they need to be either masked out or subtracted off. For permanent artifacts, a bad pixel mask is used to mask the frame. Bad pixel masks need to be throughout the lifetime of the detector. Over time more pixels may become dead, hot or frozen. In some cases, hot or frozen pixels might be able to be recovered, depending on their severity. Hot, dead or frozen pixels, can be found by taking various length N number of dark exposures and median combining them, if features are 10-20 sigma above and below the noise level, they will be clipped added to the masked. Since the pixel-to-pixel response will change, we may apply a set of flat-fielding before clipping. CRs on the other hand need to be removed from the frame. There are several methods for removing CRs such as L.A. Cosmic (Dokkem et al. 2001) which takes the Laplacian to find artifacts with steep slopes in individual frames. In the most severe cases it is possible find cosmic rays by taking a median of several science and calibration frames. In order to properly mitigate CRs in the median case, one needs 3 or more frames in multiples of odd numbers (i.e. 3, 5, 7…).

Flat Fielding

The response of individual pixels may vary significantly pixel-to-pixel. In order to correct for this in the imager and the IFS slicer, it is necessary illuminate the detector with a uniform light source. This can be accomplished by an internal flat field, twilight flats or science fields with a sparse number of sources. Internal flat fields are usually provided by the calibration unit, in this case the NSCU. Twilight flats are flats in which the telescope points to an area of the sky, while tracking, and dithering around while taking exposures. Science fields with sparse number of sources can be also used for flat fields, as long as only few bright sources need to be masked out and the there is a large enough dither such that every pixel will be on blank sky at some point. The exposures for all cases need be before the detectors become non-linear but with enough signal to illuminate each pixel, this is typically 10000-20000 counts for a Hawaii-2RG. The flats are processed by the normalize flats algorithm . Science frames are then divided by the normalized flats. Sources detected in sky flats or science flats will be masked out.

Scaled sky-subtraction

Imager: Science images can as be used sky in sparse or medium density fields (density compared to the dither size). In high density fields or a source that fills the detectors, a separate sky will be needed. The sky background will be estimated by measuring the median of each detector. For the real-time pipeline the median of the sky background will be subtracted off. The F-DRP will include an advanced algorithm that iteratively determines the sky background and subtracts it off following the procedure outlined in Clement et al. 2012;

  1. Subtract median from each image (detector)

    • Account for gain and sky, apply 5 sigma rejection (i.e. MAD)

    • Keep track of the values subtracted

  2. For each individual median subtracted image:

    • Run SExtractor to detect all objects in frame and generate an object mask.

      • Save “-object” check image

  3. Compute the median of each pixel using the mask computed previously but not including the pixel from the image you are working on

  4. Add background value from (3) to the median computed in (1), this will be the new sky background

  5. Subtract the new background from original image

The iterative process would include the mosaicking algorithm which would be the following

  1. Stack final images

  2. Redetect sources

  3. Mask sources

This process would repeat until no more additional sources are detected for masking.

IFS: The sky-subtraction algorithm will scale the sky frame to match each of the individual science frames, utilizing the Davies et al. 2007 methodology. Various OH lines arise from families of vibrational transitions. While sky lines can vary randomly throughout the night, these families fluctuate together. Using brighter sky lines, comparing the science and sky data cubes it is possible to determine the ratio between OH lines for each transition family. These ratios can be applied to the sky data cube in order to minimize the residuals in the subtracted cube. The scaling ratios are applied to the entire sky data cube, rather than to an extracted spectrum, such that any spatial or wavelength variations in the sky lines across the cube will still be accurately matched and cancelled out in the sky subtraction.

Flux calibration

Imager: To convert from DN to flux (erg/s/cm^2/Hz) or AB magnitude a standard star needs to be observed in the same instrument configuration, airmass, and close in time as the science observations.
For science fields and filters that overlap with SDSS (citation), Pan-STARRS (citation) or UKIDDS (citation) will be able to use the stars within the science frame for as standards (2MASS may be too bright and low resolution to use). For science fields outside these surveys or for more precise photometry will require observing a standard star from a standard field. Apertures of increasing radii will be used to determine the curve of growth and the appropriate aperture to use with the PSF and seeing, maximizing S/N. Once an aperture size is determined, the flux is integrated the flux for a given band to produce the flux of the star in DN. Aperture corrections will be applied based PSF and the seeing. For relative photometry, \(m_1-m_2=-2.5log_{10}\left(i \dfrac{f_1}{f_2}\right)\), where m1 and m2 are magnitudes of the sources and f1 and f1 are fluxes of the sources. This can be performed with a single source or the entire field with known sources to scale image. The zeropoints of the image can be determined from the known sources integrated flux and magnitude, (i.e. \(m=-2.5log_{10}\left( \dfrac{DN}{exptime}\right)+zeropoint\)). On sky tests will be required to determine the extinction corrected instrumental zeropoints. IFS: To convert from DN to flux units (erg/s/cm^2/Ang) a standard star needs to be observed in the same instrument configuration, airmass,and close in time as the science observations. In the near-IR the standard star at minimum needs to have zJHK photometry or ideally a spectrophotometric standard (in which a calibrated spectrum already exists). For a standard star with zJHK photometry, the photometry will be fit with a Planck law (or Rayleigh-Jeans approximation \(1/\lambda^4\)). Apertures of increasing radii will be used to determine the curve of growth and the appropriate aperture to use with the PSF and seeing. Once an aperture size is determined, the flux is integrated for a given wavelength to produce the spectrum of the standard star in DN. Aperture corrections based on the with growth curve and imager data. The science data cube and standard data cube are normalized by the exposure time such that they are each DN/s (count rate).
For the standard, we take the ratio of the flux (ergs/s/cm^2/Ang) over the count rate (DN/s). Each spaxel in the science data cube is multiplied by the ratio (flux/count rate) from the standard \(F_{sci}=\dfrac{F_{std}}{R_{std}}*R_{sci}\) , where F is flux (erg/s/cm^2) and R is count rate (DN/s)

Mosaic/Combine SCI

Imager: Mosaicking in the imager will be based on the dither pattern selected, and integer and non-integer pixel shifts will be supported. The dithers will be stored in the FITS header keywords and there will be support for an external file with the offsets. For integer pixel shits, frames will be combined using the median or mean, with sigma clipping to clip out deviant pixels. The clipping options will include using the standard deviation or median absolute deviation (MAD). For non-integer pixel shifts, there are widely used efficient software packages that handle drizzling and resampling, such as SWarp and DrizzlePac (previously known as AstroDrizzle).[j]

IFS: Mosaicking in the IFS will be relative to a source or the dither keywords in the FITS headers at a fixed PA. There will also be an option to stack the images based on an external offsets file. Currently, only integer pixel shifts will be supported. Frames will be combined using the median or mean, with sigma clipping to clip out deviant pixels. The clipping options will include using the standard deviation or median absolute deviation (MAD).

Imager Algorithms

Field distortion correction

The field distortion correction will correct the distortion in the imager field due to the optics of the system. These distortions can be chromatic and may need to be corrected per band. A calibration file with the distortion solution will be used using the distortion solution algorithm . The final image will need to be rectified and resampled based on the distortion solution. Software already exists to perform this task such as, SWarp, for rectifying and resampling the image based on the new distortion solution.

IFS Algorithms

Spectral extraction

Lenslet: Flux from an individual lenslet will be spread out into neighboring lenslets. Depending on the spacing between the lenslets, will determine how much flux falls into a neighboring lenslet. In order to recover the flux for an individual lenslet, it will be necessary to perform a deconvolution on the entire lenslet array, assigning flux to the appropriate lenslet. OSIRIS uses the Gauss-Seidel method to iteratively assign flux to individual spatial pixels (spaxels; Krabbe et al. 2004). The biggest assumption of the method is the knowledge of the PSF. In order to mitigate this problem, the PSF needs to be mapped in 2D and the structure of each lenslets PSF needs to be known precisely. Thus, the spectral extraction requires additional calibration files, rectification matrix (rectmat), which contains information about each lenslets PSF as a function of wavelength. Additional methods may be needed during INT. Slicer: Spectral extraction of the slicer will be similar to MOS (multi-object spectroscopy). The trace of each spectrum will be performed, typically fitting a low order polynomial. An aperture will be used over the spectrum, optimizing signal-to-noise (Horne 1986?). The extraction will be highly dependent on the extraction region and sky-subtraction algorithms .

Wavelength calibration

Wavelength calibration is performed using on arc lamps taken during daytime calibrations, typically Ar, Kr, and Xe. The arc lamps provide better velocity resolution and stability over OH skylines. A global wavelength solution is found for all of the spectra by fitting a low order polynomial. Legendre polynomials are preferred as they can be inverted (i.e. wavelength(pixel) → pixel(wavelength)) without significant errors in the coefficients. Using the global solution, a solution is found for each spaxel (spatial pixel). The solutions will be resampled to a common linear wavelength scale. These solutions are found be fairly stable in OSIRIS and we expect them to be similar. We anticipate checking the solution monthly for any changes. The solutions will be static based on the input lamp spectra and date they were taken.

Cube assembly

The spectral data cubes are assembled in this algorithm. The algorithm takes each extracted spectrum from spectral extraction routine and maps them to an x, y position on the sky (spatial rectification) based on the WCS information, and their z positions are shifted based on their individual wavelength solutions. The data cube format is (x, y, wavelength), which is common among data cubes with wavelength and frequency (i.e. VLT/SINFONI, ALMA and VLA).

Residual ADC

If necessary, implement residual ADC module (TBD). The ADC corrects the for the refraction caused by the atmosphere, at varying airmasses (or elevation). If the residuals from ADC correction are significant (like 4th order), it may be necessary to implement a module. To calibrate it, on-sky tests are required. One such test is to use a star to map the dispersion through the system at varying airmasses. Once the system is calibrated, temperature and pressure from the local weather, dome, telescope and instrument can be incorporated into the correction of the residuals per wavelength of light. With temperature/pressure lookup table, the DRS will have the correct spectral trace for the extraction. See instrument dispersion for how this is dealt with internally.

Telluric correction

Telluric absorption is caused by the Earths atmosphere, in which all spectra are attenuated by it. In order to correct for it, typically a featureless star is used to measure the attenuation carefully and apply a correction to the science spectra. Telluric correction as outlined by Vacca et al. 2003:

  1. Normalization of the observed A-type main sequence star spectrum (e.g. O, B, and A should be fine with “featureless” spectra, as well as white dwarfs) in the vicinity of a suitable absorption feature (as defined below);

  2. Determination of the radial velocity shift of the A-type star;

  3. Shifting the Vega model spectrum to the radial velocity of the A-type star;

  4. Scaling and reddening the Vega model spectrum to match the observed magnitudes of the A-type star;

  5. Construction of a convolution kernel from a small region around an absorption feature in the normalized observed A-type and model Vega spectra;

  6. Convolution of the kernel with the shifted, scaled, and reddened model of Vega;

  7. Scaling the equivalent widths of the various H lines to match those of the observed A-type star.

Finally, the convolved model is divided by the observed A-type spectrum and the resulting telluric correction spectrum is multiplied by the observed target spectrum.

Advanced Algorithms

Optimizing readouts

All of the algorithms used with ROP-DRS, including the various sampling techniques (i.e UTR, MCDS), will be available offline for an end user that wants extra control of optimizing the readouts of their science. For example, a user may want to include readouts with a specific seeing constraint (i.e. removing poor seeing frames).

PSF-reconstruction

Knowledge of the PSF is essential in the reduction of AO data. However, this is challenging because of changing conditions (seeing) and the rate at which they change as well as the structure of the PSF. In order to reconstruct the PSF for a given observation, a simulated PSF from the NFIRAOS PSF simulator will be used to do the deconvolution on the imager and IFS. Laurent Jolissaint et al. 2011 (AO4ELT 2011)

Calibration algorithms

Rectmats

IFS lenslets: Rectmat (or rectification matrix) is a calibration file used for the spectral extraction of the IFS lenslet data for a specific scale and filter. An individual rectmat contains information about each lenslets PSF as a function of wavelength. The rectmats are constructed from spectral white light scans, which scan each individual lenslet to determine their PSF and contribution to neighboring lensiets. This information can also be used to remove the variation from lenslet-to-lenslet, similar to a flat field.

Distortion solution

The distortion solution algorithm will determine the distortion of the image on the imager. It will be constructed by using a static uniform grid pinholes (pinhole mask) and on sky calibration using dense stellar field (i.e. globular cluster). The distortion solution will be determined by fitting some type of nth order 2D polynomial (surface) to the position of the pinholes. Software already exists to perform these tasks such as; (1) SExtractor, for detecting the sources and (2) SCAMP, for determining the distortion.

Super sky

Super sky frames are median combined sky frames. The purpose of combining them is increase the signal-to-noise of the sky. The super sky frames are used for scaled sky subtraction of the imager (in the case where the source fills the imager) and the IFS slicer.

Super dark

See generate master dark

Instrumental dispersion

The optics of IRIS (including from NFIRAOS) can cause spectral curvature, or instrumental chromatic dispersion. A white light fiber can be used to map the dispersion (x and y position of the spectra) in the system. In OSIRIS, most of the instrumental dispersion was caused by the dichroic used in the AO system.

Normalize flats

Imager: The normalize flats algorithm takes the imaging flats and generates a normalized flat (values 1 or close to one) which are used to correct the pixel-to-pixel variation. To normalize the imaging flats, N number of flats are median combined, subtracted by a dark (real-time) or master dark (F-DRP) and then divided by either the median or mode (depending on the distribution of pixel values on the detector) of the combined flats.

normalize_flat = (flat - dark)/np.median(flat - dark)

IFS slicer: The normalize flats algorithm takes the spectral flats and generates a normalized flat (values 1 or close to one) which are used to correct the pixel-to-pixel variation. The spectral flats median combined and subtracted by a dark (real-time) or master dark (F-DRP). To normalize the spectral flats, the spectral response is fit with a polynomial and subtracted off each flat, and then divided by either their median or mode (depending on the distribution of pixel values on the detector).