Gain

The gain of a camera is the number of electrons represented by one data number, and it has to be measured for each tap separately.

An \(^{55}\text{Fe}\) source provides a ruler. Its nucleus captures one of its own electrons and becomes \(^{55}\text{Mn}\), which emits an X-ray at 5.895 keV (the K-\(\alpha\) line) or, about one time in ten, at 6.490 keV (the K-\(\beta\) line). An X-ray of either energy releases a known number of electrons in silicon, so the position of the corresponding peak in a histogram of event charges measures the gain directly.

This report measures the gain of the ESIS channel 3 camera from four exposures gathered by the MSFC sounding rocket team on 2017-07-12, and then measures the gain of the channel 1 camera a second way, from the noise in a pair of uniformly illuminated images.

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
import astropy.visualization
import named_arrays as na
import msfc_ccd

astropy.visualization.quantity_support();

The energies of the two lines, and the charge each one releases, are given by msfc_ccd.Fe55.

fe55 = msfc_ccd.Fe55()

print(f"K-alpha  {fe55.energy_k_alpha:0.4f}  ->  {fe55.charge_k_alpha:0.1f}")
print(f"K-beta   {fe55.energy_k_beta:0.4f}  ->  {fe55.charge_k_beta:0.1f}")

K-alpha  5.8950 keV  ->  1602.3 electron
K-beta   6.4905 keV  ->  1764.2 electron

That charge is not a fixed number. It is the energy of the X-ray divided by the energy needed to create one electron-hole pair, and the latter falls as silicon warms, so it depends on how cold the sensor was. msfc_ccd.Fe55 defaults to minus 55 degrees Celsius, the temperature the ESIS cameras are designed to run at. At room temperature the same X-ray would release nearly a percent more charge, which would bias the gain by the same amount.

warm = msfc_ccd.Fe55(temperature=300 * u.K)

print(f"{fe55.charge_k_alpha:0.1f} at {fe55.temperature}")
print(f"{warm.charge_k_alpha:0.1f} at {warm.temperature}")

1602.3 electron at -55.0 deg_C
1616.9 electron at 300.0 K

Load the four sample images as a sequence and split them into taps.

axis_time = "time"

images = msfc_ccd.fits.open(
    path=na.ScalarArray(
        ndarray=np.array(msfc_ccd.samples.paths_fe55_esis3),
        axes=axis_time,
    ),
)

taps = images.taps

axis_x = taps.axis_x
axis_y = taps.axis_y
axis_tap_x = taps.axis_tap_x
axis_tap_y = taps.axis_tap_y

Finding the events

An X-ray releases its charge in one pixel and the pixels immediately around it, so an event is a pixel far above the dark level whose eight neighbors are not. msfc_ccd.TapData.hits() finds them and measures the charge of each one as the sum of the surrounding three by three region, which recovers the charge that spilled into the neighbors.

The result has the shape of the images, with the charge of each event in the pixel where it landed and numpy.nan everywhere else.

hits = taps.hits()

num = np.sum(np.isfinite(hits.outputs), axis=(axis_time, axis_x, axis_y))

print(f"events found in each tap, over {taps.shape[axis_time]} images:")
print(num.ndarray)

events found in each tap, over 4 images:
[[182 261]
 [126 205]]

The histogram

Pooling the events along the time axis and the two detector axes gives one histogram of event charges for each tap.

axis_charge = "charge"

hist = na.histogram(
    a=hits.outputs,
    bins={axis_charge: 81},
    axis=(axis_time, axis_x, axis_y),
    min=0 * u.DN,
    max=800 * u.DN,
)

edges = hist.inputs
centers = (edges[{axis_charge: slice(None, -1)}] + edges[{axis_charge: slice(1, None)}]) / 2

The K-\(\alpha\) peak is the tall one near 630 DN, with the K-\(\beta\) peak about 10% further right and about 10% as tall. The broad population below them is the events which lost part of their charge, either to a neighboring pixel or to the surface of the sensor.

fig, ax = na.plt.subplots(
    axis_rows=axis_tap_y,
    axis_cols=axis_tap_x,
    nrows=taps.shape[axis_tap_y],
    ncols=taps.shape[axis_tap_x],
    sharex=True,
    sharey=True,
    constrained_layout=True,
)
na.plt.stairs(
    edges,
    hist.outputs,
    axis=axis_charge,
    ax=ax,
    baseline=None,
)
na.plt.set_ylabel("number of events", ax[{axis_tap_x: 0}])
na.plt.set_xlabel("charge (DN)", ax=ax[{axis_tap_y: 0}])
na.plt.text(
    x=0.05,
    y=0.95,
    s=taps.label,
    ax=ax,
    transform=na.plt.transAxes(ax),
    ha="left",
    va="top",
);

An event from that peak looks like this. Almost all of the charge is in a single pixel, which is what makes the peak sharp; summing the three by three region around it collects the small remainder that reached the neighbors. Events which lost more than that to their neighbors, or to the surface of the sensor, are the ones spread out below the peak in the histogram.

tap = {axis_tap_x: 0, axis_tap_y: 0}
peak = centers[np.argmax(hist.outputs[tap], axis=axis_charge)]

# The event in the first image whose charge is closest to the peak
charge = hits.outputs[tap][{axis_time: 0}]
index = np.nanargmin(np.abs(charge - peak))
index = {ax: index[ax].ndarray for ax in index}

half = 4
region = {
    axis_x: slice(index[axis_x] - half, index[axis_x] + half + 1),
    axis_y: slice(index[axis_y] - half, index[axis_y] + half + 1),
}

fig, ax = plt.subplots(figsize=(4, 3.5), constrained_layout=True)
image = na.plt.pcolormesh(
    C=taps.unbiased.active.outputs[tap][{axis_time: 0}][region].value,
    axis_rgb=None,
    ax=ax,
)
ax.set_xlabel("detector $x$ (pix)")
ax.set_ylabel("detector $y$ (pix)")
plt.colorbar(image.ndarray.item(), ax=ax, label="signal (DN)");

The fit

msfc_ccd.gain.fe55() fits the two lines together. Their separation and their relative height are fixed by the line energies and the emission probabilities, so the only free parameters are the gain, the width shared by both lines, and the size of a flat background which absorbs the events that lost charge.

The fit is an unbinned maximum likelihood, so the histogram above is only for looking at: the bin width plays no part in the result. This matters because a single image gives only a few tens of events per tap.

gain = msfc_ccd.gain.fe55(taps).outputs

gain.ndarray

$[[2.5525477,~2.5781669],~ [2.5228956,~2.533233]] \; \mathrm{\frac{e^{-}}{DN}}$

Overlaying the fitted position of each line on the histogram shows where the fit placed them.

charge_k_alpha = (fe55.charge_k_alpha / gain).to(u.DN)
charge_k_beta = (fe55.charge_k_beta / gain).to(u.DN)

fig, ax = na.plt.subplots(
    axis_rows=axis_tap_y,
    axis_cols=axis_tap_x,
    nrows=taps.shape[axis_tap_y],
    ncols=taps.shape[axis_tap_x],
    sharex=True,
    sharey=True,
    constrained_layout=True,
)
na.plt.stairs(
    edges,
    hist.outputs,
    axis=axis_charge,
    ax=ax,
    baseline=None,
)
na.plt.axvline(
    x=charge_k_alpha,
    ax=ax,
    color="tab:orange",
    linestyle="--",
    label=r"fitted K-$\alpha$",
)
na.plt.axvline(
    x=charge_k_beta,
    ax=ax,
    color="tab:green",
    linestyle="--",
    label=r"fitted K-$\beta$",
)
na.plt.set_ylabel("number of events", ax[{axis_tap_x: 0}])
na.plt.set_xlabel("charge (DN)", ax=ax[{axis_tap_y: 0}])
handles, labels = ax.ndarray.flat[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="outside upper center", ncols=2);

The result

MSFC measured this camera, serial number 9, with an independent analysis of the same test, and published a gain for each of its four quadrants. Their quadrant numbers map onto our tap indices through the layout of the amplifiers across the readout frame, which msfc_ccd.TapData.amplifier records.

published = na.ScalarArray(
    ndarray=np.array([[2.57, 2.59], [2.53, 2.52]]) * u.electron / u.DN,
    axes=(axis_tap_y, axis_tap_x),
)

for i in range(taps.shape[axis_tap_y]):
    for j in range(taps.shape[axis_tap_x]):
        index = {axis_tap_y: i, axis_tap_x: j}
        print(
            f"  {taps.amplifier[index].ndarray}"
            f"   ours {gain[index].ndarray:0.3f}"
            f"   MSFC {published[index].ndarray:0.2f}"
            f"   {100 * (gain[index] / published[index] - 1).ndarray:+5.1f}%"
        )

  1/E   ours 2.553 electron / DN   MSFC 2.57 electron / DN    -0.7%
  4/F   ours 2.578 electron / DN   MSFC 2.59 electron / DN    -0.5%
  2/H   ours 2.523 electron / DN   MSFC 2.53 electron / DN    -0.3%
  3/G   ours 2.533 electron / DN   MSFC 2.52 electron / DN    +0.5%

Every tap agrees to better than one percent.

Four images are enough for this. Repeating the fit on forty images of the same run moves every tap by less than half a percent, and the requirement these cameras were built to, 2.2 to 5.0 electron / DN, is met by a wide margin either way.

The photon transfer curve

The second method needs no X-rays, only a pair of uniformly illuminated images, or flats. The charge in each pixel of a flat is a count of photons, so its noise is Poisson noise, whose variance in electrons equals the signal. In data numbers, the variance is

\[\sigma^2 = \frac{S}{g} + \sigma_r^2,\]

where \(S\) is the signal, \(g\) is the gain, and \(\sigma_r\) is the readout noise, so the ratio of the signal to the variance measures the gain.

The variance has to come from the noise alone. A single flat also varies from pixel to pixel because the illumination is uneven and because each pixel responds a little differently, and both are far larger than the noise. They are the same in two flats taken one after the other with the same illumination, so subtracting one from the other removes them, and half the variance of the difference is \(\sigma^2\).

The flats included with this package are from the channel 1 camera, gathered two seconds apart during the linearity test on 2017-07-12, along with a pair of darks gathered two minutes later.

flats = msfc_ccd.fits.open(
    path=na.ScalarArray(
        ndarray=np.array([
            msfc_ccd.samples.path_led_esis1,
            msfc_ccd.samples.path_led_esis1_next,
        ]),
        axes=axis_time,
    ),
).taps

darks = msfc_ccd.fits.open(
    path=na.ScalarArray(
        ndarray=np.array([
            msfc_ccd.samples.path_led_dark_esis1,
            msfc_ccd.samples.path_led_dark_esis1_next,
        ]),
        axes=axis_time,
    ),
).taps

msfc_ccd.noise.photon_transfer() computes the signal and the variance of each tap from each pair of images. The darks give the point at zero signal, where the variance is the square of the readout noise.

ptc_flats = msfc_ccd.noise.photon_transfer(flats, axis_time)
ptc_darks = msfc_ccd.noise.photon_transfer(darks, axis_time)

for i in range(flats.shape[axis_tap_y]):
    for j in range(flats.shape[axis_tap_x]):
        index = {axis_tap_y: i, axis_tap_x: j, axis_time: 0}
        print(
            f"  {flats.amplifier[index].ndarray}"
            f"   flats: signal {ptc_flats.inputs[index].ndarray:7.0f},"
            f" variance {ptc_flats.outputs[index].ndarray:6.0f}"
            f"   darks: variance {ptc_darks.outputs[index].ndarray:4.1f}"
        )

  1/E   flats: signal   15641 DN, variance   6200 DN2   darks: variance 16.2 DN2
  4/F   flats: signal   11732 DN, variance   4693 DN2   darks: variance 15.0 DN2
  2/H   flats: signal   24188 DN, variance   9579 DN2   darks: variance 17.4 DN2
  3/G   flats: signal   19076 DN, variance   7612 DN2   darks: variance 18.0 DN2

One pair of flats gives one point for each tap. The LED lights each tap unevenly, though, so small regions of the same pair sample a range of signals, which traces out a section of the curve. Split each tap into square tiles of 32 by 32 pixels, leaving out the masked rows at the start of each tap, and compute the signal and the variance of each tile.

axis_tile = "tile"
axis_pixel = "pixel"


def tiles(a: na.AbstractScalar, n: int = 32) -> na.ScalarArray:
    """Split each tap into square tiles of n by n pixels."""
    x = a.ndarray_aligned((axis_tap_y, axis_tap_x, axis_y, axis_x))
    num_y = x.shape[-2] // n
    num_x = x.shape[-1] // n
    x = x[..., : num_y * n, : num_x * n]
    x = x.reshape(x.shape[:2] + (num_y, n, num_x, n))
    x = np.moveaxis(x, 4, 3)
    x = x.reshape(x.shape[:2] + (num_y * num_x, n * n))
    return na.ScalarArray(x, axes=(axis_tap_y, axis_tap_x, axis_tile, axis_pixel))


rows = {axis_y: slice(flats.camera.sensor.num_masked, None)}
active = flats.unbiased.active.outputs[rows]
first = active[{axis_time: 0}]
second = active[{axis_time: 1}]

signal_tiles = tiles((first + second) / 2).mean(axis_pixel)
variance_tiles = np.var(tiles(second - first), axis=axis_pixel) / 2

msfc_ccd.gain.photon_transfer() turns the signal and the variance of the whole of each tap into a gain, measuring the readout noise from the blank columns of the same pair.

gain_ptc = msfc_ccd.gain.photon_transfer(flats, axis_time).outputs

gain_ptc.ndarray

$[[2.5292256,~2.5076198],~ [2.5297409,~2.5115486]] \; \mathrm{\frac{e^{-}}{DN}}$

The tiles fall along the line that this gain and the readout noise predict.

signal_line = na.linspace(0, 30000, axis="signal", num=2) * u.DN
variance_line = signal_line * u.electron / gain_ptc
variance_line = variance_line + ptc_darks.outputs[{axis_time: 0}]

fig, ax = na.plt.subplots(
    axis_rows=axis_tap_y,
    axis_cols=axis_tap_x,
    nrows=flats.shape[axis_tap_y],
    ncols=flats.shape[axis_tap_x],
    sharex=True,
    sharey=True,
    constrained_layout=True,
)
na.plt.plot(
    signal_tiles,
    variance_tiles,
    ax=ax,
    axis=axis_tile,
    linestyle="none",
    marker=".",
    markersize=2,
    alpha=0.3,
    label="tiles",
)
na.plt.plot(
    ptc_darks.inputs,
    ptc_darks.outputs,
    ax=ax,
    axis=axis_time,
    linestyle="none",
    marker="o",
    color="black",
    label="darks",
)
na.plt.plot(
    signal_line,
    variance_line,
    ax=ax,
    axis="signal",
    color="tab:red",
    label="fitted gain",
)
na.plt.set_ylabel("variance (DN$^2$)", ax[{axis_tap_x: 0}])
na.plt.set_xlabel("signal (DN)", ax=ax[{axis_tap_y: 0}])
na.plt.text(
    x=0.05,
    y=0.95,
    s=flats.label,
    ax=ax,
    transform=na.plt.transAxes(ax),
    ha="left",
    va="top",
)
handles, labels = ax.ndarray.flat[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="outside upper center", ncols=3);

These flats are bright, between about 12,000 and 24,000 DN. The full linearity test stepped the LED through 37 levels, and repeating the measurement at each of them shows that the ratio of the signal to the variance is constant up to about 5,000 DN and then rises, because the variance of the ESIS cameras grows more slowly than the signal as the pixels fill. For this camera, in electron / DN,

signal (DN)	1/E	4/F	2/H	3/G
400 to 900	2.49	2.48	2.45	2.45
1,700 to 3,700	2.49	2.48	2.45	2.45
4,000 to 8,800	2.49	2.48	2.47	2.46
11,000 to 24,000	2.52	2.50	2.53	2.51
21,000 to 46,000	2.60	2.54	2.69	2.62

where the 2/H and 3/G taps, which the LED lit more brightly, are at the upper end of each range. The gain of this camera is the value at low signal, 2.49, 2.48, 2.45 and 2.45 electron / DN, and the sample flats overestimate it by 1 to 3 percent.

Comparing the two methods

The two methods can only be compared on the same camera. On 2017-07-12 the \(^{55}\text{Fe}\) source was on the channel 3 and 4 cameras, and the LED lit all four. On the channel 4 camera, with 800 images of the \(^{55}\text{Fe}\) source and the photon transfer curve at low signal, in electron / DN,

tap	Fe 55	photon transfer	difference
1/E	2.587	2.559	+1.1%
4/F	2.570	2.542	+1.1%
2/H	2.594	2.563	+1.2%
3/G	2.529	2.506	+0.9%

The \(^{55}\text{Fe}\) gain is higher by about 1 percent on every tap. On the channel 3 camera of the samples above the difference is 2 to 3.5 percent. That camera has a grade 5 CCD with many columns that trap charge, which lowers the charge of the events in those columns; leaving out the worst of them lowers its \(^{55}\text{Fe}\) gain by 0.25 to 0.55 percent, and milder traps may account for more.

Neither measurement has an obvious flaw at that level. Pairs of flats ten seconds apart, and the variance of each pixel over a sequence of flats, give the same gain to within 0.2 percent, so the variance of the flats holds nothing but noise. Fitting a fuller model of the shape of the \(^{55}\text{Fe}\) lines, or fitting the K-\(\alpha\) line alone, moves the \(^{55}\text{Fe}\) gain by less than 0.2 percent, and the ring of pixels just outside the 3 by 3 region around each event holds no charge. A blue LED occasionally releases a second electron from a photon, which raises the variance of the flats, but by only 0.1 to 0.4 percent. The energy needed to create an electron-hole pair, which sets the charge of each \(^{55}\text{Fe}\) line, is itself uncertain at about the percent level, but a measurement with a CCD at 185 K that took its gain from flats, Kotov et al. [2018], points the other way.

So the gain of each tap is known to about 1 percent, the size of the difference between the two methods.