utilize dark current in variance calcs

docs/stcal/ramp_fitting/description.rst

@@ -195,7 +195,7 @@ time in seconds (from the keyword TGROUP).
 The variance of the slope in a segment due to Poisson noise is:
 
 .. math::
-   var^P_{s} = \frac{ slope_{est} }{  tgroup \times gain\ (ngroups_{s} -1)}  \,,
+   var^P_{s} = \frac{ slope_{est} + darkcurrent}{  tgroup \times gain\ (ngroups_{s} -1)}  \,,
 
 where :math:`gain` is the gain for the pixel (from the GAIN reference file),
 in e/DN. The :math:`slope_{est}` is an overall estimated slope of the pixel,

src/stcal/ramp_fitting/ols_fit.py

@@ -1138,7 +1138,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
             # 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.average_dark_current[rlo:rhi, :]
 
             # 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
@@ -1175,10 +1176,11 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans)
         # Huge variances correspond to non-existing segments, so are reset to 0
         #  to nullify their contribution.
         var_p3[var_p3 > utils.LARGE_VARIANCE_THRESHOLD] = 0.0
-        var_p3[:, med_rates <= 0.0] = 0.0
+        med_rate_mask = med_rates <= 0.0
+        var_p3[:, med_rate_mask] = 0.0
         warnings.resetwarnings()
 
-        var_p4[num_int, :, med_rates <= 0.0] = 0.0
+        var_p4[num_int, :, med_rate_mask] = ramp_data.average_dark_current[med_rate_mask][..., np.newaxis]
         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

@@ -54,9 +54,10 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False):
     ramp_data = ramp_fit_class.RampData()
 
     if isinstance(model.data, u.Quantity):
-        ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, model.pixeldq)
+        ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq,
+                             model.pixeldq, model.average_dark_current)
     else:
-        ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq)
+        ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq, model.average_dark_current)
 
     # Attribute may not be supported by all pipelines.  Default is NoneType.
     drop_frames1 = model.meta.exposure.drop_frames1 if hasattr(model, "drop_frames1") else None

src/stcal/ramp_fitting/ramp_fit_class.py

@@ -6,6 +6,7 @@ def __init__(self):
         self.err = None
         self.groupdq = None
         self.pixeldq = None
+        self.average_dark_current = None
 
         # Meta information
         self.instrument_name = None
@@ -41,7 +42,7 @@ def __init__(self):
 
         self.current_integ = -1
 
-    def set_arrays(self, data, err, groupdq, pixeldq):
+    def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current):
         """
         Set the arrays needed for ramp fitting.
 
@@ -68,6 +69,7 @@ def set_arrays(self, data, err, groupdq, pixeldq):
         self.err = err
         self.groupdq = groupdq
         self.pixeldq = pixeldq
+        self.average_dark_current = average_dark_current
 
     def set_meta(self, name, frame_time, group_time, groupgap, nframes, drop_frames1=None):
         """