Up-the-Ramp Readout Simulation — Liger Detector

from liger_iris_sim.utr.create_ramp import create_ramp

import numpy as np
import matplotlib.pyplot as plt

Detector and readout parameters

np.random.seed(42)

shape = (2048, 2048)
n_reads   = 32
readtime  = 2.7      # seconds between reads
n_channels = 32      # parallel readout channels

gain       = np.ones(shape, dtype=np.float32)
flat       = np.ones(shape, dtype=np.float32)
dark       = np.zeros(shape, dtype=np.float32)   # e-/s
bias       = np.full(shape, 1000.0, dtype=np.float32)  # e-
rn_map     = np.full(shape, 9.0, dtype=np.float32)     # e- RMS
kTC_noise  = 50.0    # e- RMS

# Uniform source: 1000 e-/s everywhere
electron_rate = np.full(shape, 1000.0, dtype=np.float32)

Run the UTR simulation

ramp_data = create_ramp(
    electron_rate,
    readtime=readtime,
    n_reads=n_reads,
    nonlin_coeffs=[1E-5, 1, 0],
    gain=gain,
    flat=flat,
    dark=dark,
    bias=bias,
    kTC_noise=kTC_noise,
    poisson_noise=True,
    read_noise=rn_map,
    convert_to_uint16=True,
    clip_ramps=True,
    max_cores=1,
    n_channels=n_channels,
)

data = ramp_data['data']   # (n_reads, 2048, 2048), uint16
print(f"Ramp shape: {data.shape}, dtype: {data.dtype}")

Single-pixel ramp and fitted slope

The simple slope estimate below is wrong due to saturated pixels and nonlinearity.

cy, cx = shape[0] // 2, shape[1] // 2
t = np.arange(n_reads) * readtime
y = data[:, cy, cx].astype(np.float32)
pfit = np.polyfit(t, y, deg=1)

print(f"Injected rate : 1000.0 e-/s")
print(f"Recovered rate: {pfit[0]:.1f} e-/s  (centre pixel)")

plt.figure(figsize=(7, 4))
plt.plot(t, y, 'o', ms=4, label='reads')
plt.plot(t, np.polyval(pfit, t), '-', lw=1.5, label=f'fit: {pfit[0]:.0f} e/s')
plt.xlabel('Time (s)')
plt.ylabel('Signal (ADU)')
plt.title('Up-the-ramp — centre pixel')
plt.legend()
plt.tight_layout()
plt.show()

Slope map from linear fits

# Vectorised least-squares slope over the ramp axis
t_norm = t - t.mean()
slopes = np.tensordot(t_norm, data.astype(np.float32), axes=([0], [0])) / np.dot(t_norm, t_norm)

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

im0 = axes[0].imshow(slopes, origin='lower', cmap='viridis')
plt.colorbar(im0, ax=axes[0], label='Slope (ADU/s)')
axes[0].set_title('Slope map')

axes[1].hist(slopes.ravel(), bins=80, color='C0', edgecolor='none')
axes[1].axvline(1000.0, color='r', lw=1.5, ls='--', label='injected')
axes[1].set_xlabel('Slope (ADU/s)')
axes[1].set_ylabel('Pixels')
axes[1].set_title('Slope distribution')
axes[1].legend()

plt.tight_layout()
plt.show()

print(f"Slope  mean : {slopes.mean():.2f} ADU/s")
print(f"Slope  std  : {slopes.std():.2f} ADU/s")