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 include 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.convolution import Ring2DKernel
+from astropy.convolution import convolve
+
 
 from . import constants
 from . import twopoint_difference as twopt
@@ -848,6 +850,7 @@ def find_faint_extended(
     Total number of showers detected.
 
     """
+    log.info('Flagging Showers')
     read_noise_2 = readnoise_2d**2
     data = indata.copy()
     data[gdq == sat_flag] = np.nan

src/stcal/ramp_fitting/ols_fit.py

@@ -19,7 +19,8 @@
 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):
     """
     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 +58,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 +100,16 @@
     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)
         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 +149,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 +164,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, number_slices)
 
     pool = Pool(processes=number_slices)
     pool_results = pool.starmap(ols_ramp_fit_single, slices)
@@ -432,8 +439,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 +480,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 +627,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):
     """
     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 +655,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
@@ -1049,6 +1062,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 +1153,8 @@
             # 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, :]) +
+                                              ramp_data.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 +1194,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] = ramp_data.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

@@ -30,7 +30,7 @@
 BUFSIZE = 1024 * 300000  # 300Mb cache size for data section
 
 
-def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False):
+def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False, avg_dark_current=0):
     """
     Create an internal ramp fit class from a data model.
 
@@ -77,6 +77,7 @@
     ramp_data.num_rows = ramp_data.data.shape[2]
 
     ramp_data.suppress_one_group_ramps = suppress_one_group
+    ramp_data.avg_dark_current = avg_dark_current
 
     return ramp_data
 
@@ -92,7 +93,8 @@
     max_cores,
     dqflags,
     suppress_one_group=False,
-):
+    avg_dark_current=0.0,
+    ):
     """
     Calculate the count rate for each pixel in all data cube sections and all
     integrations, equal to the slope for all sections (intervals between
@@ -141,6 +143,9 @@
         Find ramps with only one good group and treat it like it has zero good
         groups.
 
+    avg_dark_current : flt
+        The average dark current for this detector in units of electrons per second.
+
     Returns
     -------
     image_info : tuple
@@ -164,16 +169,15 @@
     # Create an instance of the internal ramp class, using only values needed
     # for ramp fitting from the to remove further ramp fitting dependence on
     # data models.
-    ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group)
+    ramp_data = create_ramp_fit_class(model, dqflags, suppress_one_group, avg_dark_current)
 
     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):
     """
     This function begins the ramp fit computation after the creation of the
     RampData class.  It determines the proper path for computation to take
@@ -211,10 +215,6 @@
         to use for multi-proc. The total number of cores includes the SMT cores
         (Hyper Threading for Intel).
 
-    dqflags : dict
-        A dictionary with at least the following keywords:
-        DO_NOT_USE, SATURATED, JUMP_DET, NO_GAIN_VALUE, UNRELIABLE_SLOPE
-
     Returns
     -------
     image_info : tuple
@@ -247,8 +247,7 @@
 
         # 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)
         gls_opt_info = None
 
     return image_info, integ_info, opt_info, gls_opt_info

src/stcal/ramp_fitting/ramp_fit_class.py

@@ -15,6 +15,7 @@ def __init__(self):
         self.groupgap = None
         self.nframes = None
         self.drop_frames1 = None
+        self.avg_dark_current = 0.0
 
         # Data quality flags
         self.flags_do_not_use = 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,
@@ -11,6 +11,7 @@
     detect_jumps,
 )
 
+
 DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8}
 
 
@@ -445,4 +446,26 @@
     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)
+
+@pytest.mark.skip(reason="local test only")
+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)