JP-3470: Add ifu covariance scaling to extract1d (#8457)

CHANGES.rst

@@ -60,12 +60,13 @@ extract_1d
 - Correct the output slit name for non-primary slit extractions in the
   spec3 pipeline, for NIRSpec fixed slit mode. [#8470]
 
+- Add ``ifu_covar_scale`` parameter to help correct for IFU cube covariance. [#8457]
+
 extract_2d
 ----------
 
 - Added handling for NIRCam GRISM time series pointing offsets. [#8449]
 
-
 flat_field
 ----------
 

docs/jwst/extract_1d/arguments.rst

@@ -121,7 +121,13 @@ The ``extract_1d`` step has the following step-specific arguments.
    is a MIRI MRS only paramenter. Values accepted are between 0.5 to 3.0. The default extraction
    size is set to 2 * FWHM. Values below 2 will result in a smaller
    radius, a value of 2 results in no change to radius and a value above 2 results in a larger
-   extraction radius.  
+   extraction radius.
+
+``--ifu_covar_scale``
+   A float to be multiplied into the error arrays of the extracted spectra to account
+   for covariance between adjacent spaxels in the IFU data cube.  The default value is
+   1.0 (i.e., no correction) unless set by a user or a parameter reference file.  This
+   parameter only affects MIRI and NIRSpec IFU spectroscopy.
 
 ``--soss_atoca``
   This is a NIRISS-SOSS algorithm-specific parameter; if True, use the ATOCA

docs/jwst/extract_1d/description.rst

@@ -297,11 +297,12 @@ is used to set the resampling value. The default value is 10.
 
 For IFU cubes the error information is contained entirely in the ERR array, and is not broken out into the
 VAR_POISSON, VAR_RNOISE, and VAR_FLAT arrays.  As such, ``extract_1d`` only propagates this
-non-differentiated error term.  Note that while covariance is also extremely important for IFU data cubes
-(as the IFUs themselves are significantly undersampled) this term is not presently computed or taken
-into account in the ``extract_1d`` step.  As such, the error estimates should be taken as a rough
-approximation that will be characterized and improved as flight data become available. 
-
+non-differentiated error term.  Since covariance is also extremely important for undersampled IFU data
+(see discussion by Law et al. 2023; AJ, 166, 45) the optional parameter ``ifu_covar_scale``
+will multiply all ERR arrays in the extracted spectra by a constant prefactor to account
+for this covariance.  As discussed by Law et al. 2023, this prefactor provides
+a reasonable first-order correction for the vast majority of use cases.  Values for the prefactor
+are provided in the ``extract_1d`` parameter reference files for MIRI and NIRSpec.
 
 .. _MIRI-MRS-1D-residual-fringe:
 

jwst/extract_1d/extract.py

@@ -2549,6 +2549,7 @@ def run_extract1d(
         ifu_rfcorr: Union[bool, None],
         ifu_set_srctype: str,
         ifu_rscale: float,
+        ifu_covar_scale: float,
         was_source_model: bool = False,
 ) -> DataModel:
     """Extract 1-D spectra.
@@ -2619,6 +2620,10 @@ def run_extract1d(
         radius, a value of 2 results in no change to the radius and a value above 2 results in a larger
         extraction radius.
 
+    ifu_covar_scale : float
+        Scaling factor by which to multiply the ERR values in extracted spectra to account
+        for covariance between adjacent spaxels in the IFU data cube.
+
     was_source_model : bool
         True if and only if `input_model` is actually one SlitModel
         obtained by iterating over a SourceModelContainer.  The default
@@ -2669,6 +2674,7 @@ def run_extract1d(
         ifu_rfcorr,
         ifu_set_srctype,
         ifu_rscale,
+        ifu_covar_scale,
         was_source_model,
     )
 
@@ -2733,6 +2739,7 @@ def do_extract1d(
         ifu_rfcorr: Union[bool, None] = None,
         ifu_set_srctype: str = None,
         ifu_rscale: float = None,
+        ifu_covar_scale: float = 1.0,
         was_source_model: bool = False
 ) -> DataModel:
     """Extract 1-D spectra.
@@ -2813,6 +2820,10 @@ def do_extract1d(
         radius, a value of 2 results in no change to the radius and a value above 2 results in a larger
         extraction radius.
 
+    ifu_covar_scale : float
+        Scaling factor by which to multiply the ERR values in extracted spectra to account
+        for covariance between adjacent spaxels in the IFU data cube.
+
     was_source_model : bool
         True if and only if `input_model` is actually one SlitModel
         obtained by iterating over a SourceModelContainer.  The default
@@ -3028,7 +3039,7 @@ def do_extract1d(
                 log.info(f"Overriding source type and setting it to = {ifu_set_srctype}")
             output_model = ifu.ifu_extract1d(
                 input_model, extract_ref_dict, source_type, subtract_background,
-                bkg_sigma_clip, apcorr_ref_model, center_xy, ifu_autocen, ifu_rfcorr, ifu_rscale
+                bkg_sigma_clip, apcorr_ref_model, center_xy, ifu_autocen, ifu_rfcorr, ifu_rscale, ifu_covar_scale
             )
 
         else:

jwst/extract_1d/extract_1d_step.py

@@ -91,6 +91,10 @@ class Extract1dStep(Step):
         in a smaller radius, a value of 2 results in no change to the radius and a value above
         2 results in a larger extraction radius.
 
+    ifu_covar_scale : float
+        Scaling factor by which to multiply the ERR values in extracted spectra to account
+        for covariance between adjacent spaxels in the IFU data cube.
+
     soss_atoca : bool, default=False
         Switch to toggle extraction of SOSS data with the ATOCA algorithm.
         WARNING: ATOCA results not fully validated, and require the photom step
@@ -161,6 +165,7 @@ class Extract1dStep(Step):
     ifu_rfcorr = boolean(default=False) # Apply 1d residual fringe correction
     ifu_set_srctype = option("POINT", "EXTENDED", None, default=None) # user-supplied source type
     ifu_rscale = float(default=None, min=0.5, max=3) # Radius in terms of PSF FWHM to scale extraction radii
+    ifu_covar_scale = float(default=1.0) # Scaling factor to apply to errors to account for IFU cube covariance
     soss_atoca = boolean(default=True)  # use ATOCA algorithm
     soss_threshold = float(default=1e-2)  # TODO: threshold could be removed from inputs. Its use is too specific now.
     soss_n_os = integer(default=2)  # minimum oversampling factor of the underlying wavelength grid used when modeling trace.
@@ -296,6 +301,7 @@ def process(self, input):
                         self.ifu_rfcorr,
                         self.ifu_set_srctype,
                         self.ifu_rscale,
+                        self.ifu_covar_scale,
                         was_source_model=was_source_model
                     )
                     # Set the step flag to complete
@@ -334,6 +340,7 @@ def process(self, input):
                             self.ifu_rfcorr,
                             self.ifu_set_srctype,
                             self.ifu_rscale,
+                            self.ifu_covar_scale,
                             was_source_model=was_source_model,
                         )
                         # Set the step flag to complete in each MultiSpecModel
@@ -375,6 +382,7 @@ def process(self, input):
                     self.ifu_rfcorr,
                     self.ifu_set_srctype,
                     self.ifu_rscale,
+                    self.ifu_covar_scale,
                     was_source_model=was_source_model,
                 )
 
@@ -513,6 +521,7 @@ def process(self, input):
                     self.ifu_rfcorr,
                     self.ifu_set_srctype,
                     self.ifu_rscale,
+                    self.ifu_covar_scale,
                     was_source_model=False,
                 )
 

jwst/extract_1d/ifu.py

@@ -37,7 +37,7 @@
 
 def ifu_extract1d(input_model, ref_dict, source_type, subtract_background,
                   bkg_sigma_clip, apcorr_ref_model=None, center_xy=None,
-                  ifu_autocen=False, ifu_rfcorr=False, ifu_rscale=None):
+                  ifu_autocen=False, ifu_rfcorr=False, ifu_rscale=None, ifu_covar_scale=1.0):
     """Extract a 1-D spectrum from an IFU cube.
 
     Parameters
@@ -84,6 +84,10 @@ def ifu_extract1d(input_model, ref_dict, source_type, subtract_background,
         radius, a value of 2 results in no change to radius and a value above 2 results in a larger
         extraction radius.
 
+    ifu_covar_scale : float
+        Scaling factor by which to multiply the ERR values in extracted spectra to account
+        for covariance between adjacent spaxels in the IFU data cube.
+
     Returns
     -------
     output_model : MultiSpecModel
@@ -120,6 +124,7 @@ def ifu_extract1d(input_model, ref_dict, source_type, subtract_background,
     extract_params['ifu_autocen'] = ifu_autocen
     extract_params['ifu_rfcorr'] = ifu_rfcorr
     extract_params['ifu_rscale'] = ifu_rscale
+    extract_params['ifu_covar_scale'] = ifu_covar_scale
 
     # If the user supplied extraction center coords,
     # load them into extract_params for use later.
@@ -239,6 +244,22 @@ def ifu_extract1d(input_model, ref_dict, source_type, subtract_background,
         sb_var_poisson *= 0
         b_var_poisson *= 0
 
+    # Deal with covariance in the IFU cube by multiplying 1d spectra errors by a scaling factor
+    if (extract_params['ifu_covar_scale'] != 1.0):
+        log.info("Applying ERR covariance prefactor of %g", extract_params['ifu_covar_scale'])
+        error *= extract_params['ifu_covar_scale']
+        sb_error *= extract_params['ifu_covar_scale']
+        berror *= extract_params['ifu_covar_scale']
+        f_var_poisson *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        f_var_rnoise *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        f_var_flat *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        sb_var_poisson *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        sb_var_rnoise *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        sb_var_flat *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        b_var_poisson *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        b_var_rnoise *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+        b_var_flat *= extract_params['ifu_covar_scale'] * extract_params['ifu_covar_scale']
+
     otab = np.array(
         list(
             zip(wavelength, flux, error, f_var_poisson, f_var_rnoise, f_var_flat,