From 2c32e5557ce596db4c1a27fb9bb570c046983202 Mon Sep 17 00:00:00 2001
From: David Law <dlaw@stsci.edu>
Date: Fri, 5 Jul 2024 15:03:45 -0400
Subject: [PATCH] JP-3678: Fix C poisson variance calculation (#269)

* Fix poisson variance calculation when median rate below zero

* Bugfix to C code poisson variance calculation

* Add change log

---------

Co-authored-by: Nadia Dencheva <nadia.dencheva@gmail.com>
---
 CHANGES.rst                               |  7 +++++++
 src/stcal/ramp_fitting/src/slope_fitter.c | 23 ++++++++++++++---------
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 09f78790..00b16791 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -14,6 +14,13 @@ Changes to API
 Bug Fixes
 ---------
 
+ramp_fitting
+~~~~~~~~~~~~
+
+- Fix bugs in the C algorithm Poisson variance calculation when provided with
+  an average dark current. [#269]
+
+=======
 
 ramp_fitting
 ~~~~~~~~~~~~
diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c
index 1fbd5462..b27a2f35 100644
--- a/src/stcal/ramp_fitting/src/slope_fitter.c
+++ b/src/stcal/ramp_fitting/src/slope_fitter.c
@@ -2365,7 +2365,7 @@ ramp_fit_pixel(
         pr->rate.dq |= rd->sat;
     }
 
-    if ((pr->median_rate > 0.) && (pr->rate.var_poisson > 0.)) {
+    if (pr->rate.var_poisson > 0.) {
         pr->rate.var_poisson = 1. / pr->rate.var_poisson;
     }
     if ((pr->rate.var_poisson >= LARGE_VARIANCE_THRESHOLD)
@@ -2479,7 +2479,7 @@ ramp_fit_pixel_integration_fit_slope(
         // DBG_DEFAULT_SEG; /* XXX */
 
         invvar_r += (1. / current->var_r);
-        if (pr->median_rate > 0.) {
+	if (current->var_p > 0.) {
             invvar_p += (1. / current->var_p);
         }
 
@@ -2488,7 +2488,7 @@ ramp_fit_pixel_integration_fit_slope(
     }
 
     /* Get rateints computations */
-    if (pr->median_rate > 0.) {
+    if (invvar_p > 0.) {
         pr->rateints[integ].var_poisson = 1. / invvar_p;
     }
 
@@ -2520,7 +2520,7 @@ ramp_fit_pixel_integration_fit_slope(
     // DBG_INTEG_INFO;  /* XXX */
 
     /* Get rate pre-computations */
-    if (pr->median_rate > 0.) {
+    if (invvar_p > 0.) {
         pr->rate.var_poisson += invvar_p;
     }
     pr->rate.var_rnoise += invvar_r;
@@ -2617,8 +2617,13 @@ ramp_fit_pixel_integration_fit_slope_seg_len1(
 
     seg->slope = pr->data[idx] / timing;
 
+    /* Segment Poisson variance */
     pden = (timing * pr->gain);
-    seg->var_p = (pr->median_rate + pr->dcurrent) / pden;
+    if (pr->median_rate > 0.) {
+        seg->var_p = (pr->median_rate + pr->dcurrent) / pden;
+    } else {
+        seg->var_p = (pr->dcurrent)/pden;
+    }
 
     /* Segment read noise variance */
     rnum = pr->rnoise / timing;
@@ -2664,11 +2669,11 @@ ramp_fit_pixel_integration_fit_slope_seg_len2(
     seg->slope = data_diff / rd->group_time;
 
     /* Segment Poisson variance */
+    pden = (rd->group_time * pr->gain);
     if (pr->median_rate > 0.) {
-        pden = (rd->group_time * pr->gain);
         seg->var_p = (pr->median_rate + pr->dcurrent) / pden;
     } else {
-        seg->var_p = pr->dcurrent;
+        seg->var_p = (pr->dcurrent)/pden;
     }
 
     /* Segment read noise variance */
@@ -2825,11 +2830,11 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg(
     seglen = (float)seg->length;
 
     /* Segment Poisson variance */
+    pden = (rd->group_time * pr->gain * (seglen - 1.));
     if (pr->median_rate > 0.) {
-        pden = (rd->group_time * pr->gain * (seglen - 1.));
         seg->var_p = (pr->median_rate + pr->dcurrent) / pden;
     } else {
-        seg->var_p = pr->dcurrent;
+        seg->var_p = (pr->dcurrent) / pden;
     }
 
     /* Segment read noise variance */