Add option to preserve reference pixels after IRS2 processing (#8255)

Co-authored-by: Howard Bushouse <[email protected]>

CHANGES.rst

@@ -61,6 +61,10 @@ extract_1d
   pixels in the 2nd-order spectrum are flagged and would cause the step
   to fail. [#8265]
 
+- Fixed ifu auto-centroiding to only use wavelengths shortward of 26 microns
+  to avoid failures for moderate-brightness sources due to extremely low
+  throughput at the long wavelength end of MRS band 4C. [#8199]
+
 - Replaced instances of deprecated interp2d with 
   RectBivariateSpline in ``apply_apcorr``. [#8291]
 
@@ -120,9 +124,8 @@ refpix
   all groups in each integration and robustly replace bad values from their
   nearest neighbors. [#8197, #8214]
 
-- Fixed ifu auto-centroiding to only use wavelengths shortward of 26 microns
-  to avoid failures for moderate-brightness sources due to extremely low
-  throughput at the long wavelength end of MRS band 4C. [#8199]
+- Add option for NIRSpec IRS2 to preserve interleaved reference pixels in the
+  output file, for calibration and diagnostic purposes. [#8255]
 
 resample
 --------

docs/jwst/refpix/arguments.rst

@@ -1,7 +1,7 @@
 Step Arguments
 ==============
 
-The reference pixel correction step has five step-specific arguments:
+The reference pixel correction step has seven step-specific arguments:
 
 *  ``--odd_even_columns``
 
@@ -44,3 +44,11 @@ This is a factor to avoid overcorrection of intermittently bad reference
 pixels in the IRS2 algorithm. This factor is the number of sigmas away
 from the mean. The default value is 3.0, and this argument applies
 only to NIRSpec data taken with IRS2 mode.
+
+*  ``--preserve_irs2_refpix``
+
+If the ``preserve_irs2_refpix`` argument is set, interleaved reference pixels
+in IRS2 mode will be processed along with the normal pixels and preserved
+in the output.  This option is intended for calibration or diagnostic reductions
+only. For normal science operation, this argument should always be False,
+so that interleaved pixels are stripped before continuing processing.

jwst/refpix/irs2_subtract_reference.py

@@ -10,8 +10,8 @@
 log.setLevel(logging.DEBUG)
 
 
-def correct_model(input_model, irs2_model,
-                  scipix_n_default=16, refpix_r_default=4, pad=8):
+def correct_model(input_model, irs2_model, scipix_n_default=16, refpix_r_default=4,
+                  pad=8, preserve_refpix=False):
     """Correct an input NIRSpec IRS2 datamodel using reference pixels.
 
     Parameters
@@ -35,6 +35,11 @@ def correct_model(input_model, irs2_model,
         of each row (new-row overhead).  The padding is needed to preserve
         the phase of temporally periodic signals.
 
+    preserve_refpix: bool
+        If True, reference pixels will be preserved in the output.
+        This is not used in the science pipeline, but is necessary to
+        create new bias files for IRS2 mode.
+
     Returns
     -------
     output_model: ramp model
@@ -163,12 +168,19 @@ def correct_model(input_model, irs2_model,
         # Y axis.  This is the reason for the slice `nx-ny:` that is used
         # below.  The last axis of output_model.data should be 2048.
         data0 = data[integ, :, :, :]
-        data0 = subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad)
-        data[integ, :, :, nx - ny:] = data0
+        data0 = subtract_reference(data0, alpha, beta, irs2_mask, scipix_n,
+                                   refpix_r, pad, preserve_refpix=preserve_refpix)
+        if not preserve_refpix:
+            data[integ, :, :, nx - ny:] = data0
+        else:
+            data[integ, :, :, :] = data0
 
     # Convert corrected data back to sky orientation
     output_model = input_model.copy()
-    temp_data = data[:, :, :, nx - ny:]
+    if not preserve_refpix:
+        temp_data = data[:, :, :, nx - ny:]
+    else:
+        temp_data = data
     if detector == "NRS1":
         output_model.data = np.swapaxes(temp_data, 2, 3)
     elif detector == "NRS2":
@@ -177,7 +189,8 @@ def correct_model(input_model, irs2_model,
         output_model.data = temp_data
 
     # Strip interleaved ref pixels from the PIXELDQ, GROUPDQ, and ERR extensions.
-    strip_ref_pixels(output_model, irs2_mask)
+    if not preserve_refpix:
+        strip_ref_pixels(output_model, irs2_mask)
 
     return output_model
 
@@ -642,7 +655,8 @@ def flag_bad_refpix(datamodel, n_sigma=3.0, flag_only=False, replace_only=False)
                               | dqflags.pixel['DO_NOT_USE'])
 
 
-def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
+def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n,
+                       refpix_r, pad, preserve_refpix=False):
     """Subtract reference output and pixels for the current integration.
 
     Parameters
@@ -682,6 +696,11 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
         The effective number of pixels sampled during the pause at the end
         of each row (new-row overhead).
 
+    preserve_refpix: bool
+        If True, reference pixels will be preserved in the output.
+        This is not used in the science pipeline, but is necessary to
+        create new bias files for IRS2 mode.
+
     Returns
     -------
     data0: ramp data
@@ -722,6 +741,8 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
     hnorm1 = ind_n + (refpix_r + 2) * ((ind_n + scipix_n // 2) // scipix_n)
     href1 = ind_ref + (scipix_n + 2) * (ind_ref // refpix_r) + scipix_n // 2 + 1
 
+    unpad = np.sort(np.hstack([hnorm1, href1]))
+
     # Subtract the average over the ramp for each pixel.
     # b_offset is saved so that it can be added back in at the end.
     b_offset = data0.sum(axis=0, dtype=np.float64) / float(ngroups)
@@ -852,7 +873,9 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
         # FFT.  This really shouldn't matter.
         normalization = float(shape_d[2] * shape_d[3])
 
-        if beta is not None:
+        # Set up refout if alpha was provided
+        refout0 = None
+        if alpha is not None:
             # IDL:  refout0 = reform(data0[*,*,*,0], sd[1] * sd[2], sd[3])
             refout0 = data0[0, :, :, :].reshape((shape_d[1], shape_d[2] * shape_d[3]))
 
@@ -872,7 +895,7 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
 
         # IDL:  for k=0,3 do oBridge[k]->Execute,
         #           "for i=0, s3-1 do r0[*,i] += beta * refout0[*,i]"
-        if beta is not None:
+        if alpha is not None:
             r0k_fft += (alpha[k - 1] * refout0)
         del refout0
 
@@ -888,20 +911,31 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
         # IDL:  r0 = reform(r0, sd[1], sd[2], sd[3], 5, /over)
         r0k = r0k.reshape(shape_d[1], shape_d[2], shape_d[3])
         r0k = r0k.real
-        r0k = r0k[:, :, hnorm1]
+        if not preserve_refpix:
+            r0k = r0k[:, :, hnorm1]
+        else:
+            r0k = r0k[:, :, unpad]
 
         # Subtract the correction from the data in this sector
-        data0[k, :, :, hnorm1] -= np.transpose(r0k, (2, 0, 1))
+        if not preserve_refpix:
+            data0[k, :, :, hnorm1] -= np.transpose(r0k, (2, 0, 1))
+        else:
+            data0[k, :, :, unpad] -= np.transpose(r0k, (2, 0, 1))
         del r0k
 
     # End of loop over 4 sectors
 
     # Original data0 array has shape (5, ngroups, 2048, 712). Now that
     # correction has been applied, remove the interleaved reference pixels.
     # This leaves data0 with shape (5, ngroups, 2048, 512).
-    data0 = data0[:, :, :, hnorm1]
+    if not preserve_refpix:
+        data0 = data0[:, :, :, hnorm1]
+    else:
+        data0 = data0[:, :, :, unpad]
 
     # Unflip the data in the sectors that have opposite readout direction
+    if preserve_refpix:
+        data0[0, :, :, :] = data0[0, :, :, ::-1]
     data0[2, :, :, :] = data0[2, :, :, ::-1]
     data0[4, :, :, :] = data0[4, :, :, ::-1]
 
@@ -920,13 +954,23 @@ def subtract_reference(data0, alpha, beta, irs2_mask, scipix_n, refpix_r, pad):
     # the interleaved reference pixels stripped out.
     # IDL:  data0 = reform(data0[*, 1:*, *, *], s[2], s[2], s[3], /over)
     # Note:  ny x ny, not ny x nx.
-    data0 = data0[:, :, 1:, :].reshape((ngroups, ny, ny))
+    if not preserve_refpix:
+        data0 = data0[:, :, 1:, :].reshape((ngroups, ny, ny))
+    else:
+        data0 = data0.reshape((ngroups, ny, nx))
 
     # b_offset is the average over the ramp that we subtracted near the
     # beginning; add it back in.
     # Shape of b_offset is (2048, 3200), but data0 is (ngroups, 2048, 2048),
     # so a mask is applied to b_offset to remove the reference pix locations.
-    data0 += b_offset[..., irs2_mask]
+    if not preserve_refpix:
+        data0 += b_offset[..., irs2_mask]
+    else:
+        # add in only data value -
+        # reference mean should be subtracted if not stripped,
+        # except in reference sector
+        data0[..., irs2_mask] += b_offset[..., irs2_mask]
+        data0[..., :nx // 5] += b_offset[..., :nx // 5]
 
     return data0
 

jwst/refpix/refpix_step.py

@@ -23,6 +23,7 @@ class RefPixStep(Step):
         side_gain = float(default=1.0) # Multiplicative factor for side reference signal before subtracting from rows
         odd_even_rows = boolean(default=True) # Compute reference signal separately for even- and odd-numbered rows
         ovr_corr_mitigation_ftr = float(default=3.0) # Factor to avoid overcorrection of bad reference pixels for IRS2
+        preserve_irs2_refpix = boolean(default=False) # Preserve reference pixels in output
     """
 
     reference_file_types = ['refpix']
@@ -54,7 +55,8 @@ def process(self, input):
                 irs2_model = datamodels.IRS2Model(self.irs2_name)
 
                 # Apply the IRS2 correction scheme
-                result = irs2_subtract_reference.correct_model(datamodel, irs2_model)
+                result = irs2_subtract_reference.correct_model(
+                    datamodel, irs2_model, preserve_refpix=self.preserve_irs2_refpix)
 
                 if result.meta.cal_step.refpix != 'SKIPPED':
                     result.meta.cal_step.refpix = 'COMPLETE'