Making L1 images¶
An L1 (level 1) image is a “raw” image received from the detectors. The actual measurements made on the spacecraft consist of a number of non-destructive reads of the pixels of the H4RG detectors. These reads have independent read noise but because the pixels count the total number of photons having entered each pixel, the Poisson noise in different reads of the same pixel is correlated.
Because the telescope has limited bandwidth, every read is not transferred to ground stations. Instead, reads are averaged into “resultants” according to a specification called a MultiAccum table, and these resultants are transferred, archived, and analyzed. These resultants make up an L1 image, which romanisim simulates.
L1 images are created using an idealized counts image described here, which contains the number of photons each pixel of the detector would receive absent any instrumental systematics. To transform this into an L1 image, these counts must be apportioned into reads and averaged into resultants, and instrumental effects must be added.
This process proceeds by simulating each read, drawing the appropriate number of photons from the total number of photons for each read following a binomial distribution. These photons are added to a running sum that is then averaged into a resultant according to the MultiAccum table specification. This process requires drawing random numbers from the binomial distribution for every read of every pixel, and so can take on the order of a minute, but it allows detailed simulation of the statistics of the noise in each resultant together with their correlations. It also makes it straightforward to add various instrumental effects into the simulation accurately, since these usually apply to individual reads rather than to resultants (e.g., cosmic rays affect individual reads, and their affect on a resultant depends on the read in the resultant to which they apply).
After apportioning counts to resultants, systematic effects are added to the resultants. Presently only read noise is added. The read noise is averaged down like \(1/\sqrt{N}\), where \(N\) is the number of reads contributing to the resultant.
Nonlinearity¶
Non-linearity is considered when L1 images are constructed and a non-linearity model is provided (e.g., from CRDS). We treat non-linearity as a difference between the electrons captured in the detector and the amount of signal read out. This function is modeled as a high order polynomial, and the coefficients of this polynomial and its inverse are stored in CRDS (linearity, inverselinearity reference files for each detector). When assigning counts to each read, these are transformed through the inverselinearity polynomial for each pixel and then added to the resultant buffer to account for this effect. The linearity polynomial then corrects for this effects as part of calibrating an L1 file and eventually producing an L2 file.
The linearity polynomials are occasionally poorly determined or cannot be computed. When marked as problematic in the reference file, we use trivial polynomials (i.e., the identity), and mark the affected pixels with a DQ bit indicating a problematic linearity correction.
Interpixel Capacitance¶
Interpixel capacitance (IPC) is added following non-linearity and before read-out. Read noise remains independent among different pixels but the Poisson noise is correlated between pixels by the IPC. We simply convolve the resultants by a 3x3 kernel after apportioning counts to resultants and applying non-linearity but before adding read noise.
This is slightly different than including IPC in the PSF kernel because including IPC in the PSF kernel leaves the Poisson noise uncorrelated.
Persistence¶
Persistence is implemented in the simulator following Sanchez+2023. This follows the Fermi description of persistence implemented in GalSim, where the flux in electrons per second recorded in a pixel is parameterized in terms of the total number of counts recorded in an earlier frame.
Here \(P(x, t)\) is the rate in electrons per second that the pixel records \(t\) seconds following receiving a total number of electrons \(x\). The parameters \(A\), \(x_0\), \(\delta x\), \(\alpha\), \(\gamma\) may vary from pixel to pixel, though are presently fixed to global constants. This equation for the rate only applies to pixels which were illuminated more than to fill more than their half-well. We follow GalSim and linearly increase the persistence from 0 to the half-well value for illuminations between 0 and half-well.
This persistence rate is sampled with a Poisson distribution and added to each pixel read-by-read and incorporated into the resultants in the L1 images.
Persistence-affected pixels are expected to be rare, and are tracked sparsely via a list of the indices of affected pixels, the amount of the illumination, and the times of their illumination. Pixels are dropped from persistence tracking when their persistence rate is less than one electron per 100 seconds. If the same pixel is receives large fluxes multiple times, these are treated as two independent events and the resulting persistence flux is handled by summing the persistence rates given above over each event.
Cosmic rays¶
Cosmic rays are added to the simulation read-by-read. The cosmic ray parameters follow Wu et al. (2023). The locations of cosmic rays are chosen at random to sample the focal plane uniformly. Lengths are chosen according to a power law distribution \(p(l) \\sim l^{-4.33}\), with lengths between 10 and 10,000 microns. Charge deposition rates per micron are selected from a Moyal distribution located at 120 electrons per micron with a width of 50 electrons per micron. An idealized charge is computed for each pixel in a read according to the product of the deposition rate per micron and the length of the cosmic ray’s path within that pixel. This idealized charge is Poisson sampled and added to the relevant pixels in a read.
romanisim.l1 Module¶
Convert images into L1 images, with ramps.
We imagine starting with an image that gives the total number of counts from all Poisson processes (at least: sky, sources, dark current). We then need to redistribute these counts over the resultants of an L1 image.
The easiest thing to do, and probably where I should start, is to sample the image read-by-read with a binomial draw from the total counts weighted by the chance that each count landed in this particular time window. Then gather those for each resultant and average, as done on the spacecraft.
It’s tempting to go straight to making the appropriate resultants. Following Casertano (2022?), the variance in each resultant is:
where f is the count rate, N is the number of reads in the resultant, and \(\tau\) is the ‘variance-based resultant time’
where the t_k is the time of the kth read in the resultant.
For uniformly spaced reads,
where t_0 is the time of the first read in the resultant and d is the spacing of the reads.
So that gives the variance from Poisson counts in resultant. But how should we draw random numbers to get that variance and the right mean? I can separately control the mean and variance by scaling the Poisson distribution, but I’m not sure that’s doing the right thing with the higher order moments.
It probably isn’t that expensive to just work out all of the reads, individually, which will also allow more natural incorporation of cosmic rays down the road. So let’s take that approach instead for the moment.
How do we want to specify an L1 image? An L1 image is defined by a total count image and a list of lists \(t_{i, j}\), where \(t_{i, j}\) is the time at which the jth read in the ith resultant is made. We demand \(t_{i, j} > t_{k, l}\) whenever i > k or whenever i = k and j > l.
Things this doesn’t allow neatly:
jitter in telescope pointing: the rate image is the same for each read/resultant
weird non-linear systematics in darks?
Some systematics need to be applied to the individual reads, rather than to the final image. Currently linearity, persistenc, and CRs are implemented at individual read time. I need to think about when in the chain things like IPC, etc., come in. But it still seems correct to first generate the total number of counts that an ideal detector would measure from sources, and then apply these effects read-by-read as we build up the resultants—i.e., I expect the current framework will be able to handle this without major issue.
This approach is not super fast. For a high latitude set of resultants, generating all of the random numbers to determine the apportionment takes 43 s on the machine I’m currently using; this will scale linearly with the number of reads. That’s longer than the actual image production for the dummy scene I’m using (though only ~2x longer).
I don’t have a good way to speed this up. Explicitly doing the Poisson noise on each read from a rate image (rather than the apportionment of an image that already has Poisson noise) is 2x slower—generating billions of random numbers just takes a little while.
Plausibly I could figure out how to draw numbers directly from what a resultant is rather than looking at each read individually. That would likely bring a ~10x speed-up. The read noise there is easy. The poisson noise is a sum of scaled Poisson variables:
where \(c_i\) is a Poisson-distributed variable. The sum of Poisson-distributed variables is Poisson-distributed, but I wasn’t immediately able to find anything about the sum of scaled Poisson-distributed variables. The result is clearly not Poisson-distributed, but maybe there’s some special way to sample from that directly.
If we did sample from that directly, we’d still also need to get the sum of the counts in the reads comprising the resultant. So I think you’d need a separate draw for that, conditional on the number you got for the resultant. Or, reversing that, you might want to draw the total number of counts first, e.g., via the binomial distribution, and then you’d want to draw a number for what the average number of counts was among the reads comprising the resultant, conditional on the total number of counts. Then
is some kind of statistic of the multinomial distribution. That sounds a little more tractable?
We want to draw from \(\sum (N-i) c_i\). I think the probabilities are always \(1/N\), with the possible small but important exception of ‘skipped’ or ‘dropped’ reads, in which case the first read would be more like \(2/(N+1)\) and all the others \(1/(N+1)\). If the probabilities were always 1/N, this looks vaguely like it could have a nice analytic solution. Otherwise, I don’t immediately see a route forward. So I am not going to pursue this avenue further.
Functions¶
|
Add IPC to resultants. |
|
Adds read noise to resultants. |
|
Apportion counts to resultants given read times. |
|
Package and optionally write out an L1 frame. |
|
Make an L1 image from a counts image. |
|
Get the times of each read going into resultants for a read_pattern. |
|
Convert a set of times tij to corresponding probabilities for sampling. |
|
Verify that a set of times tij for a valid resultant. |
romanisim.nonlinearity Module¶
Routines to handle non-linearity in simulating ramps.
The approach taken here is straightforward. The detector is accumulating photons, but the capacitance of the pixel varies with flux level and so the mapping between accumulated photons and read-out digital numbers changes with flux level. The CRDS linearity and inverse-linearity reference files describe the mapping between linear DN and observed DN. This module implements that mapping. When simulating an image, the photons entering each pixel are simulated, and then before being “read out” into a buffer, are transformed with this mapping into observed counts. These are then averaged and emitted as resultants.
Functions¶
|
Correct the observed counts for non-linearity. |
|
Fix cases of zeros and NaNs in non-linearity coefficients. |
Classes¶
|
Keep track of non-linearity and inverse non-linearity coefficients. |
Class Inheritance Diagram¶
romanisim.persistence Module¶
Persistence module.
This module implements a persistence simulation following Sanchez+2023.
Functions¶
|
The Fermi model for persistence: A * (x/x0)**alpha * (t/1000.)**(-gamma) / (exp(-(x-x0)/dx) + 1) For influence level below the half well, the persistence is linear in x. |
Classes¶
|
Track persistence information. |
Class Inheritance Diagram¶
romanisim.cr Module¶
Functions¶
|
A function for performing inverse transform sampling. |
|
Return unnormalized Moyal distribution, which approximates a Landau distribution and is used to describe the energy loss probability distribution of a charged particle through a detector. |
|
Return unnormalized power-law distribution parameterized by a log-log slope, used to describe the cosmic ray path lengths. |
|
Generates cosmic ray parameters randomly sampled from distribution. |
|
Adds CRs to an existing image. |
|
Given a starting and ending pixel, returns a list of pixel coordinates (ii, jj) and their traversed path lengths. |