JP-3830: MRS residual cosmic ray shower straylight correction

changes/9126.straylight.rst

@@ -0,0 +1 @@
+Add an additional optional step to the MIRI MRS straylight correction code to measure and remove residual cosmic ray showers using information from the inter-slice pixels.

docs/jwst/references_general/references_general.rst

@@ -169,6 +169,8 @@ documentation on each reference file.
 |                                                       | :ref:`ABVEGAOFFSET <abvegaoffset_reffile>`       |
 +-------------------------------------------------------+--------------------------------------------------+
 | :ref:`straylight <straylight_step>`                   | :ref:`MRSXARTCORR <mrsxartcorr_reffile>`         |
++                                                       +--------------------------------------------------+
+|                                                       | :ref:`REGIONS <regions_reffile>`                 |
 +-------------------------------------------------------+--------------------------------------------------+
 | :ref:`spectral_leak <spectral_leak_step>`             | :ref:`MRSPTCORR <mrsptcorr_reffile>`             |
 +-------------------------------------------------------+--------------------------------------------------+
@@ -273,6 +275,8 @@ documentation on each reference file.
 | :ref:`REGIONS <regions_reffile>`                 | :ref:`assign_wcs <assign_wcs_step>`                   |
 +                                                  +-------------------------------------------------------+
 |                                                  | :ref:`residual_fringe <residual_fringe_step>`         |
++                                                  +-------------------------------------------------------+
+|                                                  | :ref:`straylight <straylight_step>`                   |
 +--------------------------------------------------+-------------------------------------------------------+
 | :ref:`RESET <reset_reffile>`                     | :ref:`reset <reset_step>`                             |
 +--------------------------------------------------+-------------------------------------------------------+

docs/jwst/references_general/regions_reffile.rst

@@ -28,6 +28,8 @@ The IFU takes a region reference file that defines the pixels on the detector th
 IFU slice they correspond to.  The reference file is a three-dimensional array
 measuring 1032 x 1024 x 9.  Each of the nine planes within the array represents
 a two-dimensional detector image for which each pixel has a value indicating the slice to which it belongs.
+Negative values represent slice pixels that receive light from the sky but for data quality reasons
+are not mapped to the sky using the world coordinate solution.
 The nine planes each correspond to different throughput thresholds ranging from 10% - 90%; these differ slightly
 in the effective along-slice size.
 

docs/jwst/straylight/arguments.rst

@@ -1,3 +1,32 @@
 Step Arguments
 ==============
-The ``straylight`` step has no step-specific arguments.
+The ``straylight`` step has the following optional arguments.
+
+``--clean_showers`` (boolean, default=False)
+  A flag to indicate whether to remove straylight due to residual cosmic
+  ray showers using pixels between the IFU slices.
+
+``--shower_plane`` (int, default=3)
+  Identifies the throughput plane to use from the MRS regions reference
+  files to identify between-slice pixels.  Lower values identify pixels
+  with less science signal, but result in fewer (or no) pixels between
+  slices in some bands.
+
+``--shower_x_stddev`` (float, default=18)
+   Standard deviation to use for gaussian kernel convolution in the X
+   (between-slice) direction.  Must be large enough to reach the middle
+   of the IFU slices.
+
+``--shower_y_stddev`` (float, default=5)
+   Standard deviation to use for gaussian kernel convolution in the Y
+   (along-slice) direction.
+
+``--shower_low_reject`` (float, default=0.1)
+   Percentile of low-valued pixels to reject when performing gaussian
+   kernel convolution.  Important to ensure unmasked bad pixels do
+   not bias the kernel convolution.
+
+``--shower_high_reject`` (float, default=99.9)
+   Percentile of high-valued pixels to reject when performing gaussian
+   kernel convolution.  Important to ensure unmasked bad pixels do
+   not bias the kernel convolution.

docs/jwst/straylight/main.rst

@@ -10,18 +10,19 @@ The ``straylight`` correction is only valid for MIRI MRS data.
 
 Overview
 --------
-This routine removes contamination of MIRI MRS spectral data
-by the MIRI cross artifact feature produced by internal reflections
-within the detector arrays.  As discussed in depth for the MIRI Imager
-by A. Gáspár et al. 2021 (PASP, 133, 4504), the cross artifact manifests
+This routine removes contamination of MIRI MRS spectral data due to a variety of stray light
+sources including the MIRI cross artifact, variable detector dark current, and residual cosmic ray showers.
+
+As discussed in depth for the MIRI Imager
+by A. Gáspár et al. 2021 (PASP, 133, 4504), the MIRI cross artifact feature is produced by internal reflections
+within the detector arrays.  For the MIRI MRS, the cross artifact manifests
 as a signal extending up to hundreds of pixels along the detector column and row directions from
 bright sources.  This signal has both smooth and structured components whose
 profiles vary as a function of wavelength.
 Although the peak intensity of the cross artifact is at
 most 1% of the source intensity in Channel 1 (decreasing toward longer wavelengths),
 the total integrated light in this feature can be of order 20% of the total light from a given source.
 
-
 In the MIRI MRS, such a signal extending along detector rows is more disruptive
 than for the MIRI imager.
 Since the individual IFU slices are interleaved on the detector
@@ -32,11 +33,18 @@ The purpose of this routine is thus to model the cross artifact feature in a giv
 and subtract it at the detector level prior to reformatting
 the data into three-dimensional cubes.
 
-At the same time, this step also ensures that the median count rate (in DN/s) in regions of the detector that
-see no direct light from the sky is zero for consistency with the applied flux calibration vectors.
+This step also ensures that the median count rate (in DN/s) in regions of the detector that
+see no direct light from the sky is zero for consistency with the applied flux calibration vectors.  This
+correction is necessary because the MRS darks are strongly time-variable, and significant signal can remain
+after subtraction of the reference dark in the ``dark_current`` step.
 
-Algorithm
----------
+Finally, this step optionally applies a correction for residual cosmic ray showers that can produce
+blotchy straylight features across the detector.  While the ``jump`` detection step masks the cores of the largest
+such artifacts, the extended halos of showers can cover hundreds of pixels and are best corrected as
+coherent extended stray light artifacts.
+
+Algorithm: Cross-Artifact Correction
+------------------------------------
 The basic idea of the cross artifact correction is to convolve a given science detector image with a
 kernel function that has been pre-calibrated based on observations
 of isolated sources and subtract the corresponding convolved image.
@@ -85,3 +93,25 @@ as the magnitude of the correction is typically too small to be visible from poi
 channels use only a Lorentzian kernel with the Gaussian amplitudes set to zero as such structured components are less
 obvious at these longer wavelengths.  In Channel 4 no correction appears to be necessary,
 and the amplitudes of all model components are set equal to zero.
+
+Algorithm: Residual Shower Correction
+-------------------------------------
+The residual cosmic ray shower correction relies upon the fact that cosmic rays shower artifacts tend to
+have extremely large, blobby halos that extend across one or more IFU slices and the inter-slice pixels between
+them.  These inter-slice pixels receive little to no direct light from the sky, and can be used to construct
+a model of residual showers which are coherent in detector space.
+
+Algorithmically, this routine first uses the :ref:`REGIONS <regions_reffile>` reference file to locate and
+mask all pixels on the detector that receive direct light from the sky.  Pixels that are flagged as
+DO_NOT_USE or REFERENCE_PIXEL are similarly masked, as are some highest and lowest percentile of pixel values
+to help ensure robustness against unmasked bad pixels.
+
+The unmasked remaining pixels are then used to construct a model of the detector rates in DN/s across
+both the masked and unmasked regions using a 2D Gaussian kernel.  This model of the residual shower artifacts
+is then subtracted from the data.
+
+Since the inter-slice pixels never reach truly zero throughput from the sky (particularly at short wavelengths)
+this correction is most reliable in cases for which the observations contain only faint sources and
+have had a pixel-based dedicated background subtraction applied.  While the correction can be applied to brighter
+sources, this may result in artifacts as the subtracted model may be dominated more by spilled source flux
+than by residual cosmic ray showers.

docs/jwst/straylight/reference_files.rst

@@ -6,4 +6,9 @@ for each MRS band.  In Channel 1 these vectors include power in a broad Lorentzi
 plus a pair of double-Gaussian profiles.  In Channels 2 and 3 these vectors include only
 power in the broad Lorentzian, while in Channel 4 there is no correction.
 
+Optionally, the ``straylight`` step can also use the :ref:`REGIONS <regions_reffile>`
+reference files during correction of residual cosmic ray showers.  These reference files
+describe the footprint of the dispersed IFU traces on the detector and are used to
+identify inter-slice pixels that do not see direct light from the sky.
+
 .. include:: ../references_general/mrsxartcorr_reffile.inc

jwst/straylight/straylight.py

@@ -11,6 +11,7 @@
 
 from stdatamodels.jwst.datamodels import dqflags
 from astropy.stats import sigma_clipped_stats as scs
+from astropy.convolution import convolve_fft, Gaussian2DKernel
 from .calc_xart import xart_wrapper  # c extension
 
 log = logging.getLogger(__name__)
@@ -99,7 +100,7 @@
     nrows, ncols = input_model.data.shape
 
     # Create a copy of the input data array that will be modified
-    # for use in the straylight calculations
+    # for use in the cross artifact calculations
     usedata = input_model.data.copy()
 
     # mask is same size as image - set = 1 everywhere to start
@@ -199,27 +200,25 @@
         right_model[:, :] = 0
         log.info("No parameters for right detector half, not applying Cross-Artifact correction.")
 
-    model = left_model + right_model
+    xartifact_model = left_model + right_model
     # remove the straylight correction for the reference pixels
-    model[:, 1028:1032] = 0.0
-    model[:, 0:4] = 0.0
+    xartifact_model[:, 1028:1032] = 0.0
+    xartifact_model[:, 0:4] = 0.0
 
-    # Create output as a copy of the real data prior to replacement of NaNs with zeros
-    output = input_model.copy()
-    # Subtract the model from the data
-    output.data = output.data - model
-    usedata = usedata - model
+    # Subtract the xartifact model from the original data
+    input_model.data = input_model.data - xartifact_model
+    usedata = usedata - xartifact_model
 
     # Now measure and remove the pedestal dark count rate measured between the channels
     # Embed in a try/except block to catch unusual failures
     try:
         _, themed, therms = scs(usedata[yd1:yd2, xd1:xd2])
-        pedestal = np.zeros_like(output.data) + themed
+        pedestal = np.zeros_like(input_model.data) + themed
         # remove the pedestal correction for the reference pixels
         pedestal[:, 1028:1032] = 0.0
         pedestal[:, 0:4] = 0.0
 
-        output.data = output.data - pedestal
+        input_model.data = input_model.data - pedestal
         log.info("Derived pedestal correction " + str(themed) + " DN/s")
     except Exception:
         log.info("Straylight pedestal correction failed.")
@@ -228,4 +227,78 @@
     del usedata
 
     log.info("Cross-artifact model complete.")
-    return output
+    return input_model
+
+def clean_showers(input_model, allregions, shower_plane=3, shower_x_stddev=18.0, shower_y_stddev=5.0,
+                  shower_low_reject=0.1, shower_high_reject=99.9):
+    """
+    Corrects the MIRI MRS data for straylight produced by residual cosmic ray showers.
+
+    Parameters
+    ----------
+    input_model : `~jwst.datamodels.IFUImageModel`
+        Science data to be corrected.
+
+    allregions : numpy array
+        Holds the regions information mapping MRS pixels to slices (3-D, planes for different throughput)
+
+    shower_plane : integer, optional
+        Throughput plane for identifying inter-slice regions
+
+    shower_x_stddev : float, optional
+        X standard deviation for shower model
+
+    shower_y_stddev : float, optional
+        Y standard deviation for shower model
+
+    shower_low_reject : float, optional
+        Low percentile of pixels to reject
+
+    shower_high_reject : float, optional
+        High percentile of pixels to reject
+
+    Returns
+    -------
+    output : `~jwst.datamodels.IFUImageModel`
+        Straylight-subtracted science data.
+
+    """
+
+    log.info("Applying correction for residual cosmic ray showers.")
+
+    # Create a copy of the input data array that will be modified
+    # for use in the shower calculations
+    usedata = input_model.data.copy()
+    mask_dq = input_model.dq
+
+    # Which throughput plane of the slice map should be used?
+    regions = allregions[shower_plane,:,:]
+
+    # NaN-out the science pixels by using the slice footprint regions
+    usedata[regions != 0] = np.nan
+
+    # NaN-out pixels that should not be used for computation
+    all_flags = (dqflags.pixel['DO_NOT_USE'] | dqflags.pixel['REFERENCE_PIXEL'])
+    # where are pixels set to any one of the all_flags cases
+    testflags = mask_dq & all_flags
+    # where are testflags ne 0 and mask == 1
+    bad_flags = (testflags != 0)
+    usedata[bad_flags] = np.nan
+
+    # Apply a thresholding analysis and mask out any pixels that do not pass it
+    lowcut = np.nanpercentile(usedata, shower_low_reject)
+    hicut = np.nanpercentile(usedata, shower_high_reject)
+    badpix = (usedata < lowcut) | (usedata > hicut)
+    usedata[badpix] = np.nan
+
+    # Construct a 2d gaussian convolution kernel with specified parameters
+    gauss = Gaussian2DKernel(x_stddev=shower_x_stddev, y_stddev=shower_y_stddev)
+    shower_model = convolve_fft(usedata, gauss)
+
+    # Subtract the shower model from the original data
+    input_model.data = input_model.data - shower_model
+
+    # Delete our temporary working copy of the data
+    del usedata
+
+    return input_model

jwst/straylight/straylight_step.py

@@ -14,41 +14,59 @@
 
     class_alias = "straylight"
 
-    reference_file_types = ['mrsxartcorr']
+    spec = """
+        clean_showers = boolean(default=False)  # Clean up straylight from residual cosmic ray showers
+        shower_plane = integer(default=3)  # Throughput plane for identifying inter-slice regions
+        shower_x_stddev = float(default=18) # X standard deviation for shower model
+        shower_y_stddev = float(default=5) # Y standard deviation for shower model
+        shower_low_reject = float(default=0.1) # Low percentile of pixels to reject
+        shower_high_reject = float(default=99.9) # High percentile of pixels to reject
+    """  # noqa: E501
+
+    reference_file_types = ['mrsxartcorr', 'regions']
 
     def process(self, input):
 
         with datamodels.open(input) as input_model:
+            # Set up the output result
+            result = input_model.copy()
 
             # check the data is an IFUImageModel (not TSO)
-
             if isinstance(input_model, (datamodels.ImageModel, datamodels.IFUImageModel)):
                 # Check for a valid mrsxartcorr reference file
                 self.straylight_name = self.get_reference_file(input_model, 'mrsxartcorr')
 
                 if self.straylight_name == 'N/A':
                     self.log.warning('No MRSXARTCORR reference file found')
                     self.log.warning('Straylight step will be skipped')
-                    result = input_model.copy()
                     result.meta.cal_step.straylight = 'SKIPPED'
                     return result
 
                 self.log.info('Using mrsxartcorr reference file %s', self.straylight_name)
 
                 modelpars = datamodels.MirMrsXArtCorrModel(self.straylight_name)
 
-                # Apply the correction
-                result = straylight.correct_xartifact(input_model, modelpars)
+                # Apply the cross artifact correction
+                result = straylight.correct_xartifact(result, modelpars)
 
                 modelpars.close()
+
+                # Apply the cosmic ray droplets correction if desired
+                if self.clean_showers:
+                    self.regions_name = self.get_reference_file(input_model, 'regions')
+                    with datamodels.RegionsModel(self.regions_name) as f:
+                        allregions = f.regions
+                        result = straylight.clean_showers(result, allregions, self.shower_plane,
+                                                          self.shower_x_stddev, self.shower_y_stddev,
+                                                          self.shower_low_reject, self.shower_high_reject)
+
                 result.meta.cal_step.straylight = 'COMPLETE'
 
             else:
                 if isinstance(input_model, (datamodels.ImageModel, datamodels.IFUImageModel)) is False:
                     self.log.warning('Straylight correction not defined for datatype %s',
                                      input_model)
                 self.log.warning('Straylight step will be skipped')
-                result = input_model.copy()
                 result.meta.cal_step.straylight = 'SKIPPED'
 
         return result