L2 images

L2 images are constructed from L1 images, which are in turn constructed from idealized count images. This means that even when constructing L2 images, one must go through the process of simulating how counts get apportioned among the various reads. It is challenging to realistically model the statistics in the noise of L2 images without going through this process.

L2 images are constructed by doing “ramp fitting” on the level 1 images. Each pixel of a level 1 image is a series of “resultants”, giving the measured value of that pixel averaged over a series on non-destructive reads as the exposure is being observed. A simple model for a pixel is that its flux as a function of time is simply two numbers: a pedestal and a linear ramp representing the rate at which photons are detected by the pixel. Ramp fitting turns these one-dimensional series of resultants into a “slope” image that is of interest astronomically. Due to details of the H4RG detectors, the pedestals of the ramp vary widely from exposure to exposure, and so current fitting completely throws away as non-astronomical any information in the ramp that is sensitive to the pedestal.

The package contains two algorithms for ramp fitting. The first uses “optimal” weighting, considering the full covariance matrix between each of the resultants stemming from read & Poisson noise from the sources. The covariance is inverted and combined with the design matrix in the usual least-squares approach to solve for the optimal slope and pedestal measurements for each pixel. This approach is naively expensive, but because the covariance matrices for each pixel for a one-dimensional family depending only on the ratio of the flux in the pixel to the read variance, the relevant matrices can be precomputed. These are then interpolated between for each pixel and summed over to get the parameters and variances for each pixel.

This approach does not handle cosmic rays or saturated pixels well, though for modest sized sets of resultants introducing an additional series of fits for the roughly \(2 n_\mathrm{resultant}\) sub-ramps would be straightforward. That approach would also naturally handle saturated ramps. Even explicitly computing every possible \(n_\mathrm{resultant} (n_\mathrm{resultant} - 1) / 2\) subramp would likely still be quite inexpensive for modestly sized ramps.

The second approach follows Casertano+2022. In this approach, a diagonal set of weights is used in place of the full covariance matrix. The choice of weights depend on the particular pattern of reads assigned to each resultant and the amount of flux in the ramp, allowing them to interpolate from simply differencing the first and last resultants when the flux is very large to weighting the resultants by the number of reads when the flux is zero. This approach more efficiently handles dealing with ramps that have been split by cosmic rays, and obtaining uncertainties within a few percent of the “optimal” weighting approach. For these cases, we report final ramp slopes and variances derived from the inverse variance weighted subramp slopes and variances, using the read-noise derived variances.

This is a fairly faithful representation of how level two image construction works, so there are not many additional effects to add here. But mentioning some limitations:

  • We have a simplistic saturation treatment, just clipping resultants that exceed the saturation level to the saturation level and throwing a flag.

romanisim.ramp Module

Ramp fitting routines.

The simulator need not actually fit any ramps, but we would like to do a good job simulating the noise induced by ramp fitting. That requires computing the covariance matrix coming out of ramp fitting. But that’s actually a big part of the work of ramp fitting.

There are a few different proposed ramp fitting algorithms, differing in their weights. The final derived covariances are all somewhat similarly difficult to compute, however, since we ultimately end up needing to compute

\[(A^T C^{-1} A)^{-1}\]

for the “optimal” case, or

\[(A^T W^{-1} A)^{-1} A^T W^{-1} C W^{-1} A (A^T W^{-1} A)^{-1}\]

for some alternative weighting.

We start trying the “optimal” case below.

For the “optimal” case, a challenge is that we don’t want to compute \(C^{-1}\) for every pixel individually. Fortunately, we only need \((A^T C^{-1} A)^{-1}\) (which is only a 2x2 matrix) for variances, and only \((A^T C^{-1} A)^{-1} A^T C^{-1}\) for ramp fitting, which is 2xn. Both of these matrices are effectively single parameter families, depending after rescaling by the read noise only on the ratio of the read noise and flux.

So the routines in these packages construct these different matrices, store them, and interpolate between them for different different fluxes and ratios.

Functions

construct_covar(read_noise, flux, read_pattern)

Constructs covariance matrix for first finite differences of unevenly sampled resultants.

construct_ki_and_variances(atcinva, atcinv, ...)

Construct the \(k_i\) weights and variances for ramp fitting.

construct_ramp_fitting_matrices(covar, ...)

Construct \(A^T C^{-1} A\) and \(A^T C^{-1}\), the matrices needed to fit ramps from resultants.

fit_ramps_casertano(resultants, dq, ...)

Fit ramps following Casertano+2022, including averaging partial ramps.

fit_ramps_casertano_no_dq(resultants, ...)

Fit ramps following Casertano+2022, only using full ramps.

ki_and_variance_grid(read_pattern, ...)

Construct a grid of \(k\) and covariances for the values of flux_on_readvar.

read_pattern_to_tau(read_pattern)

Construct the tau for each resultant from a read_pattern.

read_pattern_to_tbar(read_pattern)

Construct the mean times for each resultant from a read_pattern.

resultants_to_differences(resultants)

Convert resultants to their finite differences.

simulate_many_ramps([ntrial, flux, ...])

Simulate many ramps with a particular flux, read noise, and read_pattern.

Classes

RampFitInterpolator(read_pattern[, ...])

Ramp fitting tool aiding efficient fitting of large number of ramps.

Class Inheritance Diagram

Inheritance diagram of romanisim.ramp.RampFitInterpolator