RCAL-509: Update jump step docs (#1035)

Co-authored-by: Eddie Schlafly <[email protected]>

CHANGES.rst

@@ -1,3 +1,11 @@
+0.13.1 (unreleased)
+===================
+
+jump detection
+--------------
+
+- Added uneven ramp-jump detection docs. [#1035]
+
 0.13.0 (2023-11-28)
 ===================
 

docs/roman/jump/arguments.rst

@@ -46,3 +46,5 @@ The ``jump`` step has five optional arguments that can be set by the user:
   the flagging of neighbors of marginal detections. Any primary jump below this value will
   not have its neighbors flagged. The goal is to prevent flagging jumps that would be too
   small to significantly affect the slope determination.  The default value is 10.
+
+* ``--use_ramp_jump_detection``: See the description in :ref:`ramp fitting <rampfit-arguments>`.

docs/roman/jump/even_ramp.rst

@@ -0,0 +1,17 @@
+=============================
+Jump Detection for Even Ramps
+=============================
+
+.. warning::
+
+    This step can only be run on evenly-spaced ramps. Using turning this step
+    will turn off the jump detection algorithm integrated into ramp fitting.
+
+.. toctree::
+   :maxdepth: 2
+
+   description.rst
+   arguments.rst
+   reference_files.rst
+
+.. automodapi:: romancal.jump

docs/roman/jump/index.rst

@@ -5,10 +5,7 @@ Jump Detection
 ==============
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 3
 
-   description.rst
-   arguments.rst
-   reference_files.rst
-
-.. automodapi:: romancal.jump
+   uneven_ramp.rst
+   even_ramp.rst

docs/roman/jump/uneven_ramp.rst

@@ -0,0 +1,13 @@
+===============================
+Jump Detection for Uneven Ramps
+===============================
+
+.. note::
+
+    This version of jump detection will work on both unevenly-spaced and
+    evenly-spaced ramps. So there is no need to worry about which type of
+    of ramp you have.
+
+Description
+-----------
+This "step" is actually integrated into the :ref:`unevenly-spaced ramp fitting algorithm<rampfit-algorithm-ols22>`.

docs/roman/ramp_fitting/arguments.rst

@@ -1,3 +1,5 @@
+.. _rampfit-arguments:
+
 Arguments
 =========
 The ramp fitting step has the following optional argument that can be set by the user:
@@ -23,3 +25,9 @@ The following optional arguments are valid only if using the `ols` algorithm.
   six real cores and 6 virtual cores setting maximum_cores to 'half' results in a
   decrease of a factor of six in the clock time for the step to run. Depending on the system
   the clock time can also decrease even more with maximum_cores is set to 'all'.
+
+* ``--use_ramp_jump_detection``: A True/False value that specifies whether to use
+  the unevenly-spaced jump detection integrated into the ramp fitting algorithm.
+  If True, then the jump detection step will be skipped and then revisited alongside
+  the ramp fitting step. If False, then the jump detection step will be run. The
+  default is True.

docs/roman/ramp_fitting/description.rst

@@ -193,6 +193,129 @@ Upon successful completion of this step, the status attribute ramp_fit will be s
 to "COMPLETE".
 
 
+Jump Detection
+==============
+
+For most pixels, the ramp steadily accumulates flux from the sky as an integration
+proceeds.  However, in relatively rare cases, a cosmic ray can pass through the
+detector which instantaneously deposits a large amount of charge in a pixel.
+This leads the resulting ramp to have a discontinuous *jump* in a particular read,
+and accordingly to discontinuities in the resultants downlinked from the telescope.
+The jump detection algorithm attempts to identify uncontaminated segments of ramps
+for ramp fitting, so that the underlying astronomical signal can be extracted without
+contamination from these jumps.
+
+If the uneven-ramp jump detection algorithm is turned on (the default), the ramp fitting
+algorithm is then run iteratively on a "queue" (list) of ramps. The queue is initialized
+with the ramp(s).
+Then following the algorithm presented in Sharma et al (2023) (in preparation),
+the jump detection algorithm picks a ramp, say :math:`[R_0, \dots, R_M]`, out of the
+queue and runs the ramp fitting algorithm on it. It then checks the resulting ramp for jumps.
+If a jump is detected, then two sub-ramps are created from the passed in ramp which exclude
+the resultants predicted to be affected by the jump. These sub-ramps are then added to the
+queue. This process continues until the queue is empty.
+
+.. note::
+
+   It may not always be possible to create two sub-ramps around the resultants predicted to
+   be part of a jump. For example if these jump resultants include the first, second, second-to-last,
+   or last resultant of the ramp then it is not possible to create two meaningful sub-ramps, as one
+   cannot run the ramp fitting algorithm on a ramp with zero or only one resultant. Therefore, in
+   these cases, only one sub-ramp is created and added to the queue.
+
+The method use for determining if and where a jump occurs is divided into two parts. First,
+a *statistic*, :math:`S` and possible jump resultants are determined for the fitted ramp.
+Then the statistic is compared against a threshold function, :math:`S_{\text{threshold}}` to determine
+if a jump has occurred.
+
+
+Statistic and Possible Jump Resultants
+++++++++++++++++++++++++++++++++++++++
+
+The statistic used to determine if a jump has occurred in the ramp, :math:`[R_0, \dots, R_M]`,
+is computed from a list of statistics computed for each *single* and *double-difference* of
+the resultants in the ramp.  By single-difference we mean the difference between two adjacent
+resultants in the ramp, while double-difference refers to the difference between a resultant
+and a resultant two steps away (the resultant adjacent to a resultant adjacent to the resultant
+in question).
+
+To compute these statistics, the single-difference excess slope :math:`\delta_{i, i+1}` and
+the double-difference excess slope :math:`\delta_{i, i+2}` are computed as:
+
+.. math::
+
+   \delta_{i, i+1} &= \frac{R_{i+1} - R_i} {\bar t_{i+1} - \bar t_i} - \hat \alpha
+
+   \delta_{i, i+2} &= \frac{R_{i+2} - R_i} {\bar t_{i+2} - \bar t_i} - \hat \alpha
+
+where :math:`\hat \alpha` is the slope computed by the ramp fitting algorithm. The
+variance in the excess slope:
+
+.. math::
+
+   Var(\delta_{i, j}) &= \frac {Var(R_j - R_i)} {(\bar t_j - \bar t_i)^2} + f_{corr}(\hat \alpha)
+
+   Var(R_j - R_i) &= \sigma_{RN} \left( \frac{1}{N_j} + \frac{1}{N_i} \right) + \hat \alpha \left[\tau_j + \tau_i - \min(\bar t_j, \bar t_i) \right]
+
+   f_{corr}(\hat \alpha) &= - \frac{\hat \alpha}{t_{M - 1} - t_0}
+
+where :math:`\sigma_{RN}` is the read noise. The single-difference statistic, :math:`s_i^\prime`,
+and double-difference statistic, :math:`s_i^{\prime\prime}` are,
+
+.. math::
+
+   s_i^\prime &= \frac{\delta_{i, i+1}} {\sqrt{Var(\delta_{i, i+1})}}
+
+   s_i^{\prime\prime} &= \frac{\delta_{i, i+2}} {\sqrt{Var(\delta_{i, i+2})}}.
+
+The statistic :math:`s_i` for each resultants :math:`0 \leq i \leq M - 1` (no differences from the last
+resultant are possible) is then computed as:
+
+.. math::
+   :nowrap:
+
+   \[
+   s_i =
+   \begin{cases}
+   s_i^\prime & \text{if } i = M - 2\\
+   \max(s_i^\prime, s_i^{\prime\prime}) & \text{otherwise}
+   \end{cases}
+   \]
+
+
+Finally, :math:`S = \max(s_i)` is the statistic used to determine if a jump has occurred in the fitted
+ramp. The possible jump resultants for this ramp are the resultants :math:`R_i` and :math:`R_{i+1}`,
+where :math:`i = \arg\max(s_i)`.
+
+Two possible jump resultants are necessary, because the statistics cannot determine which of the two
+adjacent resultants is the one affected by the jump. This is because if the jump occurs near the last
+read making up :math:`R_i`, then it might appear that :math:`R_{i+1}` has the jump, this jump will be
+picked up the :math:`s_i^{\prime\prime}` statistic. Using just the :math:`s_i^\prime` statistic, the
+jump would be incorrectly identified in :math:`R_{i+1}`.
+
+
+Threshold Function
+++++++++++++++++++
+
+Similarly to the statistic, the threshold function relies on the slope computed by the ramp fitting
+algorithm, :math:`\hat \alpha`. The function itself was determined empirically by running simulations
+of ramps with jumps and ramps without jumps. The threshold function was determined to be:
+
+.. math::
+
+   S_{\text{threshold}}(\hat \alpha) = 5.5 - \frac{1}{3}\log_{10}(\hat \alpha)
+
+This corresponds to identifying jumps at 5.5 sigma when the count rate is 1 electron per second, and
+4.5 sigma when the count rate is 1000 electrons per second. The decision was made to have the threshold
+depend on the count rate because the pixels with lots of signal have larger uncertainties; meaning that
+lower amplitude cosmic rays get identified in these cases.
+
+
+A jump is determined to have occurred for a ramp fit with statistic, :math:`S`, with possible jump
+resultants :math:`R_i,\ R_{i+1}`, if :math:`S \geq S_{\text{threshold}}(\hat \alpha)`. This results
+in the ramp being split into two sub-ramps :math:`[R_0, \dots R_{i-1}]` and :math:`[R_{i+2}, \dots R_M]`,
+which are then added to the ramp queue.
+
 Error Propagation
 =================