JP-3463: Dark current noise floor

CHANGES.rst

@@ -1,7 +1,10 @@
 1.5.3 (unreleased)
 ==================
 
-- 
+ramp_fitting
+------------
+- Added the avg_dark_current parameter that is used to inlcude the
+  dark current in the calculation of the Poisson noise [#236]
 
 1.5.2 (2023-12-13)
 ==================

src/stcal/jump/jump.py

@@ -1,11 +1,13 @@
 import logging
 import multiprocessing
 import time
-
-import cv2 as cv
 import numpy as np
-from astropy import stats
-from astropy.convolution import Ring2DKernel, convolve
+import cv2 as cv
+import astropy.stats as stats
+from astropy.io import fits
+from astropy.convolution import Ring2DKernel
+from astropy.convolution import convolve
+
 
 from . import constants
 from . import twopoint_difference as twopt
@@ -834,6 +836,9 @@
     Total number of showers detected.
 
     """
+    log.info('Flagging Showers')
+    fits.writeto("shower_input_data.fits", indata, overwrite=True)
+    fits.writeto("shower_input_gdq.fits", gdq, overwrite=True)
     read_noise_2 = readnoise_2d**2
     data = indata.copy()
     data[gdq == sat_flag] = np.nan
@@ -842,6 +847,7 @@
     all_ellipses = []
     first_diffs = np.diff(data, axis=1)
     first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs))
+    fits.writeto("first_diffs.fits", first_diffs, overwrite=True)
     nints = data.shape[0]
     if nints > minimum_sigclip_groups:
         mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0)
@@ -881,6 +887,10 @@
             extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8)
             exty, extx = np.where(masked_smoothed_ratio > snr_threshold)
             extended_emission[exty, extx] = 1
+            if grp == 12:
+                fits.writeto("masked_ratio.fits", masked_ratio.filled(np.nan), overwrite=True)
+                fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True)
+                fits.writeto("extended_emission.fits", extended_emission, overwrite=True)
             #  find the contours of the extended emission
             contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
             #  get the contours that are above the minimum size

src/stcal/ramp_fitting/ols_fit.py

@@ -19,7 +19,9 @@
 log.setLevel(logging.DEBUG)
 
 
-def ols_ramp_fit_multi(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores):
+def ols_ramp_fit_multi(
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores,
+        avg_dark_current):
     """
     Setup the inputs to ols_ramp_fit with and without multiprocessing. The
     inputs will be sliced into the number of cores that are being used for
@@ -57,6 +59,9 @@
         'half', and 'all'. This is the fraction of cores to use for multi-proc. The
         total number of cores includes the SMT cores (Hyper Threading for Intel).
 
+    avg_dark_current: float (electrons)
+        The average dark current for this detector
+
     Returns
     -------
     image_info : tuple
@@ -96,15 +101,17 @@
     if number_slices == 1:
         # Single threaded computation
         image_info, integ_info, opt_info = ols_ramp_fit_single(
-            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting
-        )
+            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting,
+            avg_dark_current)
         if image_info is None or integ_info is None:
             return None, None, None
 
         return image_info, integ_info, opt_info
 
     # Call ramp fitting for multi-processor (multiple data slices) case
-    image_info, integ_info, opt_info = ols_ramp_fit_multiprocessing(
+
+    else:
+        image_info, integ_info, opt_info = ols_ramp_fit_multiprocessing(
         ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
     )
 
@@ -144,6 +151,8 @@
     number_slices: int
         The number of slices to partition the data into for multiprocessing.
 
+    avg_dark_current: float
+        The average dark current in electrons
     Return
     ------
     image_info: tuple
@@ -157,8 +166,8 @@
     """
     log.info("Number of processors used for multiprocessing: %s", number_slices)
     slices, rows_per_slice = compute_slices_for_starmap(
-        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
-    )
+        ramp_data, buffsize, save_opt,
+        readnoise_2d, gain_2d, weighting, avg_dark_current, number_slices)
 
     pool = Pool(processes=number_slices)
     pool_results = pool.starmap(ols_ramp_fit_single, slices)
@@ -432,8 +441,8 @@
 
 
 def compute_slices_for_starmap(
-    ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, number_slices
-):
+        ramp_data, buffsize, save_opt,
+        readnoise_2d, gain_2d, weighting, avg_dark_current, number_slices):
     """
     Creates the slices needed for each process for multiprocessing.  The slices
     for the arguments needed for ols_ramp_fit_single.
@@ -473,9 +482,11 @@
     start_row = 0
     for k in range(len(rslices)):
         ramp_slice = slice_ramp_data(ramp_data, start_row, rslices[k])
-        rnoise_slice = readnoise_2d[start_row : start_row + rslices[k], :].copy()
-        gain_slice = gain_2d[start_row : start_row + rslices[k], :].copy()
-        slices.insert(k, (ramp_slice, buffsize, save_opt, rnoise_slice, gain_slice, weighting))
+
+        rnoise_slice = readnoise_2d[start_row:start_row + rslices[k], :].copy()
+        gain_slice = gain_2d[start_row:start_row + rslices[k], :].copy()
+        slices.insert(
+            k, (ramp_slice, buffsize, save_opt, rnoise_slice, gain_slice, weighting, avg_dark_current))
         start_row = start_row + rslices[k]
 
     return slices, rslices
@@ -618,7 +629,8 @@
     ramp_data.one_groups_time = (ramp_data.nframes + 1) * ramp_data.frame_time / 2
 
 
-def ols_ramp_fit_single(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting):
+def ols_ramp_fit_single(
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, avg_dark_current):
     """
     Fit a ramp using ordinary least squares. Calculate the count rate for each
     pixel in all data cube sections and all integrations, equal to the weighted
@@ -645,6 +657,9 @@
     weighting : str
         'optimal' is the only valid value
 
+    avg_dark_current : float (electrons)
+        The average dark current for this detector
+
     Return
     ------
     image_info : tuple
@@ -697,7 +712,7 @@
     #   noise only, read noise only, and the combination of Poisson noise and
     #   read noise. The integration-specific variances are 3D arrays, and the
     #   segment-specific variances are 4D arrays.
-    variances_ans = ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
+    variances_ans = ramp_fit_compute_variances( ramp_data, gain_2d, readnoise_2d, fit_slopes_ans, avg_dark_current)
 
     # Now that the segment-specific and integration-specific variances have
     #   been calculated, the segment-specific, integration-specific, and
@@ -1017,7 +1032,8 @@
     )
 
 
-def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans):
+def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans,
+                               avg_dark_current):
     """
     In this 'Second Pass' over the data, loop over integrations and data
     sections to calculate the variances of the slope using the estimated
@@ -1049,6 +1065,8 @@
     fit_slopes_ans : tuple
         Contains intermediate values computed in the first pass over the data.
 
+    avg_dark_current : float (electrons)
+        The average dark current for this detector
     Return
     ------
     var_p3 : ndarray
@@ -1138,7 +1156,7 @@
             # Suppress harmless arithmetic warnings for now
             warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning)
             warnings.filterwarnings("ignore", ".*divide by zero.*", RuntimeWarning)
-            var_p4[num_int, :, rlo:rhi, :] = den_p3 * med_rates[rlo:rhi, :]
+            var_p4[num_int, :, rlo:rhi, :] = (den_p3 * med_rates[rlo:rhi, :]) + avg_dark_current
 
             # Find the segment variance due to read noise and convert back to DN
             var_r4[num_int, :, rlo:rhi, :] = num_r3 * den_r3 / gain_sect**2
@@ -1178,7 +1196,10 @@
         var_p3[:, med_rates <= 0.0] = 0.0
         warnings.resetwarnings()
 
-        var_p4[num_int, :, med_rates <= 0.0] = 0.0
+        # For pixels with zero or negative rates set the Poisson variance to the
+        # dark current rate.
+        var_p4[num_int, :, med_rates <= 0.0] = avg_dark_current
+
         var_both4[num_int, :, :, :] = var_r4[num_int, :, :, :] + var_p4[num_int, :, :, :]
         inv_var_both4[num_int, :, :, :] = 1.0 / var_both4[num_int, :, :, :]
 

src/stcal/ramp_fitting/ramp_fit.py

@@ -91,8 +91,9 @@ def ramp_fit(
     weighting,
     max_cores,
     dqflags,
+    avg_dark_current=0.0,
     suppress_one_group=False,
-):
+    ):
     """
     Calculate the count rate for each pixel in all data cube sections and all
     integrations, equal to the slope for all sections (intervals between
@@ -167,13 +168,12 @@ def ramp_fit(
     ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group)
 
     return ramp_fit_data(
-        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, algorithm, weighting, max_cores, dqflags
-    )
+        ramp_data, buffsize, save_opt, readnoise_2d, gain_2d,
+        algorithm, weighting, max_cores, dqflags, avg_dark_current)
 
 
-def ramp_fit_data(
-    ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, algorithm, weighting, max_cores, dqflags
-):
+def ramp_fit_data(ramp_data, buffsize, save_opt, readnoise_2d, gain_2d,
+                  algorithm, weighting, max_cores, dqflags, avg_dark_current):
     """
     This function begins the ramp fit computation after the creation of the
     RampData class.  It determines the proper path for computation to take
@@ -215,6 +215,8 @@ def ramp_fit_data(
         A dictionary with at least the following keywords:
         DO_NOT_USE, SATURATED, JUMP_DET, NO_GAIN_VALUE, UNRELIABLE_SLOPE
 
+    avg_dark_current : float (electrons)
+        The average dark current for this detector
     Returns
     -------
     image_info : tuple
@@ -247,7 +249,7 @@ def ramp_fit_data(
 
         # Compute ramp fitting using ordinary least squares.
         image_info, integ_info, opt_info = ols_fit.ols_ramp_fit_multi(
-            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores
+            ramp_data, buffsize, save_opt, readnoise_2d, gain_2d, weighting, max_cores, avg_dark_current
         )
         gls_opt_info = None
 

tests/test_jump.py

@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-
+from astropy.io import fits
 from stcal.jump.jump import (
     calc_num_slices,
     extend_saturation,
@@ -10,6 +10,7 @@
     point_inside_ellipse,
 )
 
+
 DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8}
 
 
@@ -351,4 +352,25 @@
     assert calc_num_slices(n_rows, "all", max_available_cores) == 10
     assert calc_num_slices(n_rows, "3/4", max_available_cores) == 1
     n_rows = 9
-    assert calc_num_slices(n_rows, "21", max_available_cores) == 9
+    assert (calc_num_slices(n_rows, '21', max_available_cores) == 9)
+
+def test_shower_slowmode():
+    data = fits.getdata("shower_input_data.fits")
+    gdq = fits.getdata("shower_input_gdq.fits")
+    nint = data.shape[0]
+    ngrps = data.shape[1]
+    ncols = data.shape[2]
+    nrows = data.shape[3]
+    gain = 5
+    readnoise = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_readnoise_0086.fits")
+#    readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain
+    rng = np.random.default_rng(12345)
+#    data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7
+#    data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise
+    gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100,
+                                           snr_threshold=1.3,
+                                           min_shower_area=40, inner=1,
+                                           outer=2, sat_flag=2, jump_flag=4,
+                                           ellipse_expand=1.1, num_grps_masked=3)
+    print("number of showers", num_showers)
+    fits.writeto("outgdq.fits", gdq, overwrite=True)