Nonlinearity correction + datacube spec1d extractions + datacube fluxing improvements

deprecated/datacube.py

@@ -443,3 +443,71 @@ def make_whitelight_frompixels(all_ra, all_dec, all_wave, all_sci, all_wghts, al
             whitelight_ivar[:, :, ff] = ivar_img.copy()
     return whitelight_Imgs, whitelight_ivar, whitelightWCS
 
+
+def make_sensfunc(ss_file, senspar, blaze_wave=None, blaze_spline=None, grating_corr=False):
+    """
+    Generate the sensitivity function from a standard star DataCube.
+
+    Args:
+        ss_file (:obj:`str`):
+            The relative path and filename of the standard star datacube. It
+            should be fits format, and for full functionality, should ideally of
+            the form :class:`~pypeit.coadd3d.DataCube`.
+        senspar (:class:`~pypeit.par.pypeitpar.SensFuncPar`):
+            The parameters required for the sensitivity function computation.
+        blaze_wave (`numpy.ndarray`_, optional):
+            Wavelength array used to construct blaze_spline
+        blaze_spline (`scipy.interpolate.interp1d`_, optional):
+            Spline representation of the reference blaze function (based on the illumflat).
+        grating_corr (:obj:`bool`, optional):
+            If a grating correction should be performed, set this variable to True.
+
+    Returns:
+        `numpy.ndarray`_: A mask of the good sky pixels (True = good)
+    """
+    # TODO :: This routine has not been updated to the new spec1d plan of passing in a sensfunc object
+    #      :: Probably, this routine should be removed and the functionality moved to the sensfunc object
+    msgs.error("coding error - make_sensfunc is not currently supported.  Please contact the developers")
+    # Check if the standard star datacube exists
+    if not os.path.exists(ss_file):
+        msgs.error("Standard cube does not exist:" + msgs.newline() + ss_file)
+    msgs.info(f"Loading standard star cube: {ss_file:s}")
+    # Load the standard star cube and retrieve its RA + DEC
+    stdcube = fits.open(ss_file)
+    star_ra, star_dec = stdcube[1].header['CRVAL1'], stdcube[1].header['CRVAL2']
+
+    # Extract a spectrum of the standard star
+    wave, Nlam_star, Nlam_ivar_star, gpm_star = extract_standard_spec(stdcube)
+
+    # Extract the information about the blaze
+    if grating_corr:
+        blaze_wave_curr, blaze_spec_curr = stdcube['BLAZE_WAVE'].data, stdcube['BLAZE_SPEC'].data
+        blaze_spline_curr = interp1d(blaze_wave_curr, blaze_spec_curr,
+                                     kind='linear', bounds_error=False, fill_value="extrapolate")
+        # Perform a grating correction
+        grat_corr = correct_grating_shift(wave, blaze_wave_curr, blaze_spline_curr, blaze_wave, blaze_spline)
+        # Apply the grating correction to the standard star spectrum
+        Nlam_star /= grat_corr
+        Nlam_ivar_star *= grat_corr ** 2
+
+    # Read in some information above the standard star
+    std_dict = flux_calib.get_standard_spectrum(star_type=senspar['star_type'],
+                                                star_mag=senspar['star_mag'],
+                                                ra=star_ra, dec=star_dec)
+    # Calculate the sensitivity curve
+    # TODO :: This needs to be addressed... unify flux calibration into the main PypeIt routines.
+    msgs.warn("Datacubes are currently flux-calibrated using the UVIS algorithm... this will be deprecated soon")
+    zeropoint_data, zeropoint_data_gpm, zeropoint_fit, zeropoint_fit_gpm = \
+        flux_calib.fit_zeropoint(wave, Nlam_star, Nlam_ivar_star, gpm_star, std_dict,
+                                 mask_hydrogen_lines=senspar['mask_hydrogen_lines'],
+                                 mask_helium_lines=senspar['mask_helium_lines'],
+                                 hydrogen_mask_wid=senspar['hydrogen_mask_wid'],
+                                 nresln=senspar['UVIS']['nresln'],
+                                 resolution=senspar['UVIS']['resolution'],
+                                 trans_thresh=senspar['UVIS']['trans_thresh'],
+                                 polyorder=senspar['polyorder'],
+                                 polycorrect=senspar['UVIS']['polycorrect'],
+                                 polyfunc=senspar['UVIS']['polyfunc'])
+    wgd = np.where(zeropoint_fit_gpm)
+    sens = np.power(10.0, -0.4 * (zeropoint_fit[wgd] - flux_calib.ZP_UNIT_CONST)) / np.square(wave[wgd])
+    return interp1d(wave[wgd], sens, kind='linear', bounds_error=False, fill_value="extrapolate")

doc/calibrations/calibrations.rst

@@ -127,6 +127,7 @@ The primary calibration procedures are, in the order they're performed:
    flexure
    wave_calib
    Slit Alignment (IFU only) <alignment>
+   non-linear correction <nonlinear>
    flat_fielding
    scattlight
 

doc/calibrations/nonlinear.rst

@@ -0,0 +1,93 @@
+
+=====================
+Non-Linear correction
+=====================
+
+The non-linear correction is a simple way to correct the non-linear behavior of the
+sensor. Imagine a perfect light source that emits photons with a constant count rate.
+An ideal detector would be able to measure this count rate independently of the exposure
+time. However, in practice, detectors have a non-linear response to the incoming photons.
+This means that the count rate measured by the detector is not proportional to the count
+rate of the light source. The non-linear correction is a simple way to correct this
+non-linear behavior. Most modern CCD cameras have a linear response to about 1% over
+most of the dynamic range, so this correction is not necessary for most applications.
+The correction parameters can be measured by the user, but care should be taken with
+how you collect the calibration data, since the internal continuum light source used
+for most flat-field calibrations is not a perfect light source.
+
+At present, the non-linear correction is only implemented for KCWI, which is already very
+close to linear over most of the dynamic range of the CCD. If you have collected non-linear
+calibration data with another instrument, please contact the developers to see if we can
+implement this correction for your instrument.
+
+Calibrations required
+---------------------
+
+The non-linear correction requires a set of dedicated calibration data that can be used to measure
+the non-linear response of the detector. This calibration data should be collected with
+a continuum light source (usually internal to the instrument). We currently recommend that
+you turn on the internal lamp, and leave it for 30 minutes to warm up and settle. The lamp
+count rate will likely decay over time, so we recommend that you collect a series of exposures
+that are bracketed by `reference` exposures of constant exposure time (also known as the
+`bracketed repeat exposure` technique). The reference exposures are used to monitor the
+temporal decay of the lamp count rate. Here is a screenshot to show the lamp decay as a
+function of time. The points are the `reference exposures` and are colour-coded by the
+wavelength of the light, and the solid lines are the best-fit exponential + 2D polynomial
+decay curves. The bottom panel shows a residual of the fit. Thus, both the relative shape
+and the intensity of the lamp decay can be monitored as a function of time.
+
+.. image:: ../figures/nonlinear_lamp_decay.png
+
+Between a set of reference exposures, you should
+acquire a series of exposures with increasing exposure time. The exposure time should be
+increased by a factor of 2 between each exposure. The exposure time should be chosen so that
+the count rate is not too high, but also not too low. The counts should extend to at least 10%
+of the maximum counts of the detector. The `reference` exposures should have an exposure time
+that is close to the middle of the range of exposure times used for the other exposures. It is
+often good practice to collect an even number of exposures at each exposure time; some shutters
+move in opposite directions for adjacent exposures, and the true exposure time may depend on
+the direction that the shutter is moving.
+
+To determine the non-linearity, we assume that the correction is quadratic, such that the
+quadratic term is a small perturbation from linear response. The functional form of the
+true count rate is then:
+
+.. math::
+
+    C_T \cdot t = C_M \cdot t (1 + b \cdot C_M \cdot t)
+
+where :math:`C_T` is the true count rate, :math:`C_M` is the measured count rate, :math:`t`
+is the exposure time, and :math:`b` is the non-linearity coefficient.
+The non-linearity coefficient can be determined by fitting
+the measured count rate as a function of the exposure time. That that as the true count rate
+tends to zero, the measured count rate tends to zero. The non-linearity coefficient can be
+measured for each pixel on the detector, and it is a good idea to consider the non-linearity
+coefficient for different amplifiers separately. The above equation can be written as a matrix
+equation:
+
+.. math::
+
+    M \cdot x = e
+    M = \begin{bmatrix}
+    C_{M,1}t_1 & (C_{M,1}t_1)^2 \\
+    C_{M,2}t_2 & (C_{M,2}t_2)^2 \\
+    \vdots & \vdots \\
+    C_{M,n}t_n & (C_{M,n}t_n)^2 \\
+    \end{bmatrix}
+    x = \begin{bmatrix}
+    1/C_T \\
+    b/C_T \\
+    \end{bmatrix}
+    e = \begin{bmatrix}
+    t_1 \\
+    t_2 \\
+    \vdots \\
+    t_n \\
+    \end{bmatrix}
+
+where :math:`M` is a matrix of measured count rates, :math:`x` is a vector of the non-linearity
+coefficient and the true count rate, and :math:`e` is a vector of exposure times. The non-linearity
+coefficient can be determined by inverting the matrix :math:`M`, and solving for :math:`x`. The
+non-linearity coefficient can be determined for each pixel on the detector. The central value of
+the distribution of non-linearity coefficients can be used as a fixed non-linearity coefficient
+for each amplifier of the detector.