Repository: spacetelescope/jwst

Issue:
Updated NIRCam WFSS photom algorithm to account for varying dispersion 
_Issue [JP-3238](https://jira.stsci.edu/browse/JP-3238) was created on JIRA by [Alicia Canipe](https://jira.stsci.edu/secure/ViewProfile.jspa?name=acanipe):_

NIRCam is working on addressing an issue with the interaction of the flight WFSS/TSGRISM flux cal files with the current design of the photom step, which is causing the flux calibration to be off by a factor of 10. The near-term fix will be to divide the PHOTMJSR values by 10. However, this triggered a discussion that dividing by 10 is not quite right – the dispersion varies across the detector, and is also wavelength dependent, so we can't just create a 2048x2048 array of dispersions and extract 2D boxes. Instead, the pipeline needs to know the wavelength at each pixel in order to find the dispersion. This is something that needs to be done individually for each source in the FOV. Forcing a pixel size/sr/bin size to the sensitivity will always be an approximation and introduces an error in the flux calibration. Because it's a small variation, it will probably not be noticed until the contamination and background subtraction are resolved, but we should begin investigating this issue to improve the quality of WFSS data products. 

The photom step already creates the 2d array of wavelength values for each source. A 2d array of dispersion values would have to be constructed on the fly for each source and the pipeline would have to know the function for calculating dispersion.

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/8207?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`dac66be`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/dac66be8645b8830690109e9726bdd833b0a398f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 75.25% compared to head [(`92dac8b`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/8207?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 77.36%.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #8207      +/-   ##
==========================================
+ Coverage   75.25%   77.36%   +2.11%     
==========================================
  Files         470      451      -19     
  Lines       38435    35285    -3150     
==========================================
- Hits        28923    27300    -1623     
+ Misses       9512     7985    -1527     
```

| [Flag](https://app.codecov.io/gh/spacetelescope/jwst/pull/8207/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | | *Carryforward flag |
|---|---|---|---|
| [nightly](https://app.codecov.io/gh/spacetelescope/jwst/pull/8207/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | `77.36% <ø> (-0.01%)` | :arrow_down: | Carriedforward from [649b28f](https://app.codecov.io/gh/spacetelescope/jwst/commit/649b28f33560796614b472ee805f42b8ee4d0266?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

*This pull request uses carry forward flags. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to find out more.


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/8207?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

hbushouse: Is there a regression test run from earlier that has these changes applied? If so, can you provide the link here?
stscirij: Link to developer regression tests:

https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/1139/

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 19bd53764..bf5a0f4d3 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -46,6 +46,12 @@ photom
   ``MultiSlitModel`` data models, forcing information on units to only come
   from individual slit meta data. [#8189]
 
+- Updated photom step to include spectral dispersion when applying NIRCam WFSS and TSGRISM
+  flux calibration, because the dispersion varies with location in the field and wavelength.
+  The PHOTOM reference files have the dispersion factored out, requiring that
+  the pipeline put the (variable, calculated per pixel) dispersion back in.  Assumes that
+  the dispersion needs to be in Angstroms/pixel to match the required factor of ~10. [#8207]
+
 refpix
 ------
 
diff --git a/jwst/photom/photom.py b/jwst/photom/photom.py
index 345db4a4d..6e60a43d5 100644
--- a/jwst/photom/photom.py
+++ b/jwst/photom/photom.py
@@ -9,6 +9,7 @@ from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import dqflags
 
 from .. lib.wcs_utils import get_wavelengths
+from .. lib.dispaxis import get_dispersion_direction
 from . import miri_mrs
 from . import miri_imager
 
@@ -469,7 +470,7 @@ class DataSet():
                 row = find_row(ftab.phot_table, fields_to_match)
                 if row is None:
                     return
-                
+
             # Check to see if the reference file contains the coefficients for the
             # time-dependent correction of the PHOTOM value
             try:
@@ -853,9 +854,11 @@ class DataSet():
                     # First, compute 2D array of photom correction values using
                     # uncorrected wavelengths, which is appropriate for a uniform source
                     conversion_2d_uniform, no_cal = self.create_2d_conversion(slit,
-                                                                   self.exptype, conversion_uniform,
-                                                                   waves, relresps, order,
-                                                                   use_wavecorr=False)
+                                                                              self.exptype,
+                                                                              conversion_uniform,
+                                                                              waves, relresps,
+                                                                              order,
+                                                                              use_wavecorr=False)
                     slit.photom_uniform = conversion_2d_uniform  # store the result
 
                     # Now repeat the process using corrected wavelength values,
@@ -867,6 +870,13 @@ class DataSet():
                                                                    use_wavecorr=True)
                     slit.photom_point = conversion  # store the result
 
+                elif self.exptype in ["NRC_WFSS", "NRC_TSGRISM"]:
+                    log.info("Including spectral dispersion in 2-d flux calibration")
+                    conversion, no_cal = self.create_2d_conversion(slit,
+                                                                   self.exptype, conversion,
+                                                                   waves, relresps, order,
+                                                                   include_dispersion=True)
+
                 else:
                     conversion, no_cal = self.create_2d_conversion(slit,
                                                                    self.exptype, conversion,
@@ -879,8 +889,17 @@ class DataSet():
                 conversion, no_cal = self.create_1d_conversion(self.input.spec[self.specnum],
                                                                conversion, waves, relresps)
             else:
-                conversion, no_cal = self.create_2d_conversion(self.input, self.exptype, conversion,
-                                                               waves, relresps, order)
+                # NRC_TSGRISM data produces a SpecModel, which is handled here
+                if self.exptype in ["NRC_WFSS", "NRC_TSGRISM"]:
+                    log.info("Including spectral dispersion in 2-d flux calibration")
+                    conversion, no_cal = self.create_2d_conversion(self.input,
+                                                                   self.exptype, conversion,
+                                                                   waves, relresps, order,
+                                                                   include_dispersion=True)
+
+                else:
+                    conversion, no_cal = self.create_2d_conversion(self.input, self.exptype, conversion,
+                                                                   waves, relresps, order)
 
         # Apply the conversion to the data and all uncertainty arrays
         if isinstance(self.input, datamodels.MultiSlitModel):
@@ -974,7 +993,7 @@ class DataSet():
         return
 
     def create_2d_conversion(self, model, exptype, conversion, waves, relresps,
-                             order, use_wavecorr=None):
+                             order, use_wavecorr=None, include_dispersion=False):
         """
         Create a 2D array of photometric conversion values based on
         wavelengths per pixel and response as a function of wavelength.
@@ -997,6 +1016,9 @@ class DataSet():
         use_wavecorr : bool or None
             Flag indicating whether or not to use corrected wavelengths.
             Typically only used for NIRSpec fixed-slit data.
+        include_dispersion : bool or None
+            Flag indicating whether the dispersion needs to be incorporated
+            into the 2-d conversion factors.
 
         Returns
         -------
@@ -1013,8 +1035,18 @@ class DataSet():
 
         # Interpolate the photometric response values onto the
         # 2D wavelength grid
+        # waves is in microns, so wl_array must be in microns too
         conv_2d = np.interp(wl_array, waves, relresps, left=np.nan, right=np.nan)
 
+        if include_dispersion:
+            dispaxis = get_dispersion_direction(self.exptype, self.grating, self.filter, self.pupil)
+            if dispaxis is not None:
+                dispersion_array = self.get_dispersion_array(wl_array, dispaxis)
+                # Convert dispersion from micron/pixel to angstrom/pixel
+                dispersion_array *= 1.0e4
+                conv_2d /= np.abs(dispersion_array)
+            else:
+                log.warning("Unable to get dispersion direction, so cannot calculate dispersion array")
         # Combine the scalar and 2D conversion factors
         conversion = conversion * conv_2d
         no_cal = np.isnan(conv_2d)
@@ -1022,6 +1054,34 @@ class DataSet():
 
         return conversion, no_cal
 
+    def get_dispersion_array(self, wavelength_array, dispaxis):
+        """Create an array of dispersion values from the wavelength array
+
+        Parameters
+        ----------
+        wavelength_array : float
+            2-d array of wavelength values, assumed to be in microns
+        dispaxis : integer
+            Direction along which light is dispersed: 1 = along rows, 2 = along columns
+
+        Returns
+        -------
+        dispersion_array : float
+            2-d array of dispersion values, in microns/pixel
+
+        """
+        nrows, ncols = wavelength_array.shape
+        dispersion_array = np.zeros(wavelength_array.shape)
+        if dispaxis == 1:
+            for row in range(nrows):
+                dispersion_array[row] = np.gradient(wavelength_array[row])
+        elif dispaxis == 2:
+            for column in range(ncols):
+                dispersion_array[:, column] = np.gradient(wavelength_array[:, column])
+        else:
+            log.warning(f"Can't process data with DISPAXIS={dispaxis}")
+        return dispersion_array
+
     def create_1d_conversion(self, model, conversion, waves, relresps):
         """Create a 1D array of photometric conversion values based on
         wavelength array of input spectrum and response as a function of wavelength.


Test patch:
diff --git a/jwst/photom/tests/test_photom.py b/jwst/photom/tests/test_photom.py
index f10105cd8..fa987f5ce 100644
--- a/jwst/photom/tests/test_photom.py
+++ b/jwst/photom/tests/test_photom.py
@@ -9,6 +9,7 @@ from stdatamodels.jwst import datamodels
 from stdatamodels.jwst.datamodels import SpecModel, MultiSpecModel
 
 from jwst.photom import photom
+from jwst.lib.dispaxis import get_dispersion_direction
 
 MJSR_TO_UJA2 = (u.megajansky / u.steradian).to(u.microjansky / (u.arcsecond**2))
 
@@ -1488,9 +1489,14 @@ def test_nircam_spec():
         ix = shape[1] // 2
         iy = shape[0] // 2
         wl = slit.wavelength[iy, ix]
+        # Incude the dispersion in the correction, as per JP-3238
+        dispaxis = get_dispersion_direction(ds.exptype, ds.grating, ds.filter, ds.pupil)
+        dispersion_array = ds.get_dispersion_array(slit.wavelength, dispaxis)
+        # Convert dispersion in micron/pixel to Angstrom/pixel
+        disp = dispersion_array[iy, ix] * 10000.0
         rel_resp = np.interp(wl, wavelength, relresponse,
                              left=np.nan, right=np.nan)
-        compare = photmjsr * rel_resp
+        compare = photmjsr * rel_resp / disp
         # Compare the values at the center pixel.
         ratio = output[iy, ix] / input[iy, ix]
         assert np.allclose(ratio, compare, rtol=1.e-7)