Coincidences and timing jitter (delta-t)

A daplis-style pipeline for coincidence / timing-jitter analysis of the Kelpie ORT (timestamp) program, built on three layers:

1. unpack decodes each .bin file into the per-frame, per-pixel TDC code (time_series, the MATLAB time_mat).
2. calc_diff.calculate_differences turns structured per-pixel timestamps into per-window differences for pixel pairs.
3. delta_t drives those two over a folder of files, pools the differences and saves them to a .feather table, then histograms, fits and plots the coincidence peak.

Timing model

Each acquisition window is ~9 µs, but the ring-oscillator TDC measures only within one period of the reference clock. The recorded time_series value is the oscillator code on that range, and the timestamp is

t = code * time_unit_s        with  time_unit_s ≈ 71.6 ps  (core.timing.CODE_S)

on average. The ring sets that code width, so it does not change with the clock — only the number of codes does, as 1/f:

clock	full scale	codes
50 MHz	20 ns	~279
10 MHz	100 ns	~1397
5 MHz	200 ns	~2794
2 MHz	500 ns	~6985

That full scale is what every entry point needs, as open_shutter_time_s (one clock period) or full_scale_s. It bounds a difference, so it is also the accidental triangle's half-base — see the triangle background.

The exact per-code bin widths come from a density-test LUT — pass it as time_lut once you have one; see TDC calibration. A LUT is specific to its clock; the loaders check that for you and refuse a mismatch, because the silent version of that mistake costs 96 % of the coincidences at 2 MHz.

Units

Every time parameter in the package is in seconds, with an _s suffix: bin_width_s, plot_window_s, fit_window_s, support_s, time_unit_s, delta_window_s, open_shutter_time_s. Histogram bin centres, edges and plot axes are in picoseconds — that is the scale the TDC lives on — and the conversion happens on entry.

There is one timestamp per photon, per window, per macropixel.

Fitting the peak: the background is a triangle

This is the part that most affects the number you report.

The accidental-coincidence background under an ORT delta-t peak is triangular, not flat, because the free-running oscillator has no cycle counter. Fitting a Gaussian on a flat background over a wide window therefore inflates the fitted sigma — you measure the background's slope as if it were jitter.

Two models are provided:

model	fit	when to use
gaussian / fit_gaussian_peak	Gaussian on flat background	valid only on a narrow window around the peak
gaussian_plus_triangle / fit_gaussian_on_triangle	Gaussian on triangular background	the correct model for ORT delta-t over a wide window

The full derivation, with the demonstration that the triangle is intrinsic rather than an artefact, is in the triangular background note.

Before you analyse: check the data

Run data_quality.plot_time_code_histogram on a pixel or two first. When the TDC never stops, unpack still returns a full array of plausible integers and this whole pipeline still produces a peak-shaped histogram — it is simply meaningless. See Data quality.

Large runs: the differences are streamed, not accumulated

calculate_and_save_timestamp_differences writes as it goes. Whenever the buffered differences would exceed the 64 MB part cap they are flushed to a part file, and the parts are concatenated into the final feather one Arrow record batch at a time:

processed/
  <dataset>_delta_t_parts/
    <dataset>_delta_t_0.feather      <- flushed during the run
    <dataset>_delta_t_1.feather
    ...
  <dataset>_delta_t.feather          <- the combined result

This matters at scale. Holding a whole 10 000-file run in memory and writing it once fails twice over: the buffers exhaust RAM, and the single write then asks pyarrow for a multi-GB contiguous allocation on top of what is already held. Streaming makes the loop's memory ceiling one part plus one unpacked file, whatever the file count.

Two consequences worth knowing:

• A crashed run is recoverable. Every part flushed before the crash is a complete, readable feather. Assemble them with combine_delta_t_parts, or plot a single part directly via feather_path.
• The parts are always kept. They cost disk rather than memory, they are the crash insurance, and — see below — the read side streams them just as happily as the combined file, so there is nothing to gain from deleting them automatically. Remove the folder by hand once you have checked the combined feather.

Why 64 MB and not 10: unpack already holds a (32, 32, nframes) float64 array — about 82 MB at 10 000 frames — so below that the part buffer is not the binding constraint, and a smaller cap only multiplies the number of part files. The cap is the module constant _PART_MB, not a parameter: sharding is a memory strategy and cannot change the numbers, so no call site needs to tune it.

Reading it back: the histogram is streamed, not loaded

The write side streams, and so must the read side. A 10 000-file run's combined feather reaches ~135 GB; loading it into a DataFrame and pooling the columns needs it resident several times over, which simply does not happen.

compute_and_save_delta_histogram reads it one Arrow record batch at a time instead. A histogram over fixed edges is additive, so the counts can be summed batch by batch with one batch resident — and the result is bin-for-bin what histogramming the whole pooled array at once would have given. This is the exact baseline: no sampling, no re-binning, no rounding.

Nothing downstream ever needed the raw differences. Both fit models take (bin_centers, counts), and so does the plot; the pooled array was a pure intermediate. So the counts are saved as a stage-1 artifact:

processed/
  <dataset>_delta_t.feather        <- ~135 GB, read once
  <dataset>_delta_hist.npy         <- the counts, a few MB
  <dataset>_delta_hist.meta.json   <- binning, source fingerprint, statistics

The sidecar records bin_width_s, plot_window_s, time_unit_s, pairs and a fingerprint of the source files. A saved histogram is reused only when all of it matches, so changing the binning or rewriting the feather re-streams automatically; pass reuse=False (or reuse_histogram=False to the plot driver) to force the full pass. Re-fitting, switching between the flat and triangular background, or redrawing therefore costs a small load rather than another read of the feather.

collect_and_plot_timestamp_differences goes through this, and returns (counts, centers, fit) — enough to re-fit by hand without touching disk:

counts, centers, fit = dt.collect_and_plot_timestamp_differences(
    path,
    background="triangle",
    open_shutter_time_s=100e-9,
    plot_window_s=100e-9,
)
narrow = np.abs(centers) < 1000
dt.fit_gaussian_peak(centers[narrow], counts[narrow])

The .meta.json also reports padding_fraction: the share of the cells read that were padding rather than data. The parts are written wide and padded to the longest pair column, so an unbalanced blob pair inflates the feather — that number says how much of the 135 GB is real.

Plotting without the combined file

feather_path accepts a *_delta_t_parts folder as well as a .feather, and a run whose combined feather is missing falls back to its parts on its own. The parts are complete feathers, so the combine step is optional for plotting:

dt.collect_and_plot_timestamp_differences(
    feather_path=".../processed/SPDC_delta_t_parts"
)

Parts or combined file is not a numerical choice — the histogram is a sum over record batches either way, and a sum does not care which file a batch came out of, so the counts are identical. It is only an I/O choice, and the parts win it: histogramming them skips the combine pass, which reads and rewrites every byte of the run to produce a second copy of it. Working from the parts is the reason they are always kept. Keep the combined feather when you want one portable file per run, or when the data set is small enough to read whole into a DataFrame.

The other stage 1: histogram straight off the .bin

Everything above keeps every difference and histograms it later. That is the baseline — nothing is thrown away — but the artifact grows with the run, and the binning is still chosen at plot time, so every change of mind costs another pass over hundreds of GB.

calculate_and_save_delta_counts runs the same loop and bins the differences onto a fixed native grid as it finds them, keeping only the counts. It writes no feather at all:

processed/
  <dataset>_delta_counts.npy         <- (n_pairs, n_cells) int64
  <dataset>_delta_counts.meta.json   <- grid, pair labels, totals, progress

The size is set by the grid, not by the data. A 256-pair run on the default grid is 106 MB whether it came from 100 files or 10 000 — measured against a 3.22 GB feather for the same run, and the same 106 MB for one that reaches 135 GB. That is what makes every pair plottable at once instead of whichever columns you could afford to stream.

Choosing the grid

Two arguments, grid_s (cell width) and support_s (half-range):

data	default cell	why
raw codes (apply_TDC_calibration=False)	one code, ~71.6 ps	delta_code is an integer, so counts-per-code is a re-encoding, not a summary — exactly lossless
calibrated (unit='ps')	1/10 code, ~7.2 ps	the per-pixel LUT smears the code lattice; 7.2 ps is ~5× finer than the jitter

support_s defaults to twice the TDC full scale — read off the LUT when there is one, else from open_shutter_time_s. A difference cannot exceed the full scale, so that covers the whole physical support at any clock; a fixed 200 ns was 10× too wide at 50 MHz and 2.5× too narrow at 2 MHz. Anything landing outside is reported as n_outside rather than dropped silently; if that is not ~0, widen it.

Re-binning is free, and exact

plot_window_s becomes a slice and bin_width_s a reduce, both in memory:

dt.collect_and_plot_delta_counts(
    path,
    pairs=["6,10-21,22"],
    bin_width_s=71e-12,
    plot_window_s=10e-9,
    background="triangle",
    full_scale_s=500e-9,        # or leave it to open_shutter_time_s
)

By default the binning is snapped onto the grid, and the snap is not cosmetic:

• the width goes to a whole number of cells, or bins hold 9 then 10 then 9 cells and a smooth distribution picks up an ~11% sawtooth;
• that number is made odd. Bins are centred on zero, so their edges sit at half-integer multiples of the width — odd puts those edges on cell boundaries, even puts them through cell centres, splitting a cell at every boundary. Measured on a real 3.22 GB run: exactly zero difference from the feather path at 9, 11 and 13 cells per bin, ~1% at 8, 10 and 12;
• the window goes to a whole number of bins, so zero is a bin centre.

So with snapping on, this is not an approximation of the feather histogram — it is the same histogram. Pass snap_to_grid=False to take your numbers literally, and you will be warned about the sawtooth you are asking for.

The one honest disagreement is the outermost bin at each end. The feather path cuts data at abs(delta) <= plot_window_s but draws edges half a bin beyond it, so its first and last bins can only ever be part filled; the counts path fills them. That, not the method, is why fitted sigma differs by ~1%.

What you give up

No per-event data survives, so there is no re-calibrating with a different LUT, no re-cutting delta_window, and nothing that needs individual differences — time ordering, splitting by acquisition, higher-order correlations. It also does not make stage 1 meaningfully faster: unpacking dominates either way. The win is artifact size and replot latency. Keep the feather until you are sure you want none of that back.

Crash insurance without a part folder

The grid is rewritten every flush_every files (default 50), so a run that dies leaves a valid artifact covering everything up to its last checkpoint, and the next call resumes from it. There is only ever one file, so two runs cannot interleave — the failure mode combine_delta_t_parts exists to guard against simply does not arise here.

The sidecar records files_done and total_in_grid. The array is moved into place before the sidecar, so a crash between the two moves could pair a newer array with an older sidecar; resuming on that would re-count files and inflate the total. The resume path checks the array against total_in_grid and starts over rather than continue from something inconsistent.

Measuring the background instead of fitting it

subtract_background=True histograms the frame-shifted pixel pairs alongside the real ones, so the accidental triangle can be subtracted bin-for-bin rather than fitted. Both stage-1 paths take it, and both give the same residual — see Removing the accidental background by shifting frames.

# counts path: the lags cost one small grid plane each
dt.calculate_and_save_delta_counts(path, [signal, idler], nframes=10_000,
                                   subtract_background=True)
counts, centers, fit = dt.collect_and_plot_delta_counts(
    path, subtract_background=True, peak_window_s=2e-9)

# feather path: same call, but the lags cost one column set each - so ~9x the
# feather at the default eight lags. Use fewer, or use the counts path.
dt.calculate_and_save_timestamp_differences(path, [signal, idler], nframes=10_000,
                                            subtract_background=True,
                                            background_lags=(1, 2, 3, 4))
counts, centers, fit = dt.collect_and_plot_timestamp_differences(
    path, subtract_background=True, peak_window_s=2e-9)

The feather's lag columns are named <pair>@lag<k>; lag 0 keeps the bare pair name, and pooling "every pair column" never picks up the shifted ones.

Overwriting: the five-second countdown

Every rewrite=True in this module routes through store.confirm_rewrite. On a fresh run it does nothing at all. When something really is about to be destroyed it names each target with its size and timestamp, then counts down five seconds before proceeding:

!!  rewrite=True - the following data WILL BE OVERWRITTEN:
!!    .../processed/SPDC_delta_t.feather  (1.4 GB, written 2026-07-30 02:11)
!!    .../processed/SPDC_delta_t_parts/  (23 file(s), 1.4 GB)
!!  Press Ctrl-C within 5 s to abort and copy it somewhere safe.

The pause is a plain time.sleep, so Ctrl-C raises straight out of the analysis with the data on disk untouched. This is aimed at the common accident: a rewrite=True left in the script from the previous acquisition, silently eating an overnight run. Set store.REWRITE_DELAY_S = 0 in an unattended pipeline that has its own guard.

Note the part folder is named in the warning as well — a rewrite clears the previous run's parts, otherwise they would be pooled into the new output.

Pooling runs

combine_delta_t_feathers pools the per-run delta-t feathers under a parent folder (many acquisitions of the same measurement) into one combined .feather in a combined sub-folder of processed/. It never folds a previous combined output back into itself, and it streams batch-by-batch just as the part combiner does, so pooling a hundred runs costs one batch of memory rather than the whole pooled data set.

Part files are not picked up by the search: they end in _<i>.feather, which the *_delta_t.feather glob does not match, so nothing is ever counted twice.

API

delta_t

Coincidence / timing-jitter (delta-t) pipeline for the ORT program.

Drives 'unpack' and 'calc_diff' over a folder, pools the per-pixel-pair timestamp differences into a '.feather' table, then histograms, fits and plots the coincidence peak.

Both halves stream. Stage 1 flushes size-capped part files as it goes, and the read side histograms the feather one Arrow record batch at a time rather than loading it - a 10 000-file run's combined feather reaches ~135 GB, which no amount of RAM makes loadable. Since a histogram over fixed edges is additive, the streamed counts are bin-for-bin what loading the whole table would have given; they are saved beside the feather so re-fitting costs a small load.

NOTE on fitting: the ORT accidental background is TRIANGULAR, not flat, because the free-running oscillator has no cycle counter. Use 'fit_gaussian_on_triangle' over a wide window; 'fit_gaussian_peak' (flat background) is valid only on a narrow window around the peak, and a flat fit over a wide one inflates the fitted sigma. Or do not model it at all: subtract_background=True measures it from frame-shifted pairs and subtracts it, leaving the peak on a zero baseline ('dapkel.functions.background_subtraction', docs/guide/background_subtraction.md).

The triangle's half-base is the TDC full scale, which is one period of the reference clock - 20 ns at 50 MHz, 100 ns at 10 MHz, 500 ns at 2 MHz - so it is not a constant and neither is the number of TDC codes. Every entry point here takes open_shutter_time_s (or full_scale_s) rather than assuming 100 ns, and the calibration LUT is checked against it: feeding a 10 MHz LUT to a 2 MHz run discards every code past 100 ns, which is 79 % of the singles and 96 % of the coincidences, and the surviving peak still looks plausible.

UNITS: every time parameter in this module is in seconds and named with an _s suffix, per 'dapkel.core.timing'. Histogram edges, bin centres and plot axes are in picoseconds - that is the scale the TDC lives on - and the conversion happens on entry with 'timing.to_ps'.

Check a fresh acquisition with 'dapkel.functions.data_quality' first. See docs/guide/coincidences.md and docs/ort_triangle_background.md.

calculate_and_save_timestamp_differences

calculate_and_save_timestamp_differences(path: str, pixels: Sequence[Sequence[Pixel]], rewrite: bool = False, *, nframes: int, open_shutter_time_s: float, tag: str = 'ORT', mode: str = 'all_pairs', delta_window_s: float | None = None, apply_TDC_calibration: bool = True, daughterboard_number: str | None = None, motherboard_number: str | None = None, spad: int | str = 'average', max_files: int | None = None, subtract_background: bool = False, background_lags: int | Sequence[int] = BACKGROUND_LAGS) -> str

Unpack ORT data, compute blob-vs-blob delta-t, and save to feather.

Writes path/processed/<name>_delta_t.feather with one column per pixel pair (named "ra,ca-rb,cb"), padded to a common length; every requested pair gets a column even if it saw no coincidences.

The differences are streamed to disk as they are found, not held until the end: the running buffer is flushed to processed/<name>_delta_t_parts/<name>_delta_t_<run>_<i>.feather whenever it would exceed '_PART_MB', and the parts are concatenated into the final feather one record batch at a time. Memory therefore stays flat in the number of '.bin' files, and a run that dies at file 9 000 of 10 000 leaves every completed part on disk - recover them with 'combine_delta_t_parts'.

The parts are kept. They are complete feathers in their own right and every read-side entry point takes a part folder as readily as a single file ('compute_and_save_delta_histogram' streams either one a record batch at a time), so they are the cheaper of the two things on disk to work from as well as the crash insurance. Delete the folder by hand once the combined feather has been checked.

Each run tags its parts with its own <run> id and records them in a manifest, and a run refuses to start while another run's parts are still in the folder. Part numbering restarts at 0 every run, so without both of those a re-run would overwrite a dead run's first parts, leave its later ones, and the combine step would pool the two - counting some differences twice and inflating the feather.

With apply_TDC_calibration=True (default) the board's density-test LUT converts timestamps to picoseconds before differencing, saving delta_ps; otherwise raw delta_code is saved. Which one it is goes into the feather's schema metadata as delta_unit.

Parameters:

Name	Type	Description	Default
path	str	Path to the folder with the ORT '.bin' data files.	required
pixels	Sequence[Sequence[tuple[int, int]]]	Two groups [group_a, group_b] — the signal and idler blobs read off the hitmap, each a list of (row, col) pixels.	required
rewrite	bool	Overwrite an existing '.feather' (and part folder) for this data set. A guard against accidental overwrites: the default is False, which raises instead. When True and something really is about to be destroyed, 'store.confirm_rewrite' names it and counts down five seconds first, so a stale rewrite=True can still be caught with Ctrl-C.	False
nframes	int	Number of frames stored in each '.bin' file.	required
open_shutter_time_s	float	Open shutter time, i.e. one reference-clock period, in seconds (e.g. 100e-9 at 10 MHz, 500e-9 at 2 MHz). Selects and validates the calibration LUT, and converts delta_window_s into raw code units when no LUT is applied. There is no default: the TDC's full scale is this window, so nothing downstream is right without it.	required
tag	str	Filename fragment selecting the files. The default is 'ORT'.	'ORT'
mode	str	'all_pairs' (default) or '1v1'; see 'calc_diff.calculate_differences'.	'all_pairs'
delta_window_s	float | None	Keep only differences with abs(delta) <= delta_window_s while accumulating (bounds memory), in seconds. Converted to whichever unit the timestamps are in - ps with a LUT, raw codes without, using 'core.timing.CODE_S'. The default is None (keep all).	None
apply_TDC_calibration	bool	Apply the density-test TDC calibration (per-pixel code -> ps LUT) before differencing, so the saved delta_ps is already calibrated. The default is True, in which case the packaged board LUT is loaded via daughterboard_number / motherboard_number (a ValueError is raised if either is missing). Set False to save raw delta_code and calibrate later.	True
daughterboard_number	str | None	Daughterboard id (e.g. 'D0') used to pick the packaged board LUT via 'tdc_calibration.load_board_lut'. The default is None.	None
motherboard_number	str | None	Motherboard id (e.g. 'M0') used with daughterboard_number to pick the packaged board LUT. The default is None.	None
spad	int | str	Which SPAD's LUT to use: 'average' (default) for the mean of the board's four SPAD LUTs, an int 0..3 for a specific SPAD, or a tag fragment like 'SPAD0_S0'. ORT cannot tell which micropixel fired, so 'average' is the safe default; choose a SPAD only when the optics fix the micropixel. A missing specific SPAD falls back to the average (with a warning).	'average'
max_files	int | None	Process at most this many files. The default is None (all files).	None
subtract_background	bool	Also difference the frame-shifted pairs, so the accidental background can be measured rather than fitted. The default is False, which writes exactly the feather it always has. When True each pair gains one column per lag, named "<pair>@lag<k>", while lag 0 keeps the bare pair name - so the feather is a strict superset and every existing reader still sees the same columns. It also multiplies the feather by roughly 1 + len(background_lags): prefer the counts path ('calculate_and_save_delta_counts'), where the same lags cost one small grid plane each. See 'dapkel.functions.background_subtraction' and docs/guide/background_subtraction.md.	False
background_lags	int | Sequence[int]	Frame offsets used for the accidental estimate when subtract_background is True. The default is (1, ..., 8) - which is 9x the feather, so consider fewer here. A bare int is taken as a single lag. Lag 0 is always included and must not appear here. Ignored when subtract_background is False.	BACKGROUND_LAGS

Returns:

Type	Description
str	Path to the saved '.feather' file.

Raises:

Type	Description
FileExistsError	Raised when the '.feather' already exists and rewrite is False, or when an earlier run's parts are still in the part folder.
ValueError	Raised on a calibration request with no board, or a 0 / negative / repeated entry in background_lags.

Source code in src/dapkel/functions/delta_t.py

def calculate_and_save_timestamp_differences(
    path: str,
    pixels: Sequence[Sequence[Pixel]],
    rewrite: bool = False,
    *,
    nframes: int,
    open_shutter_time_s: float,
    tag: str = "ORT",
    mode: str = "all_pairs",
    delta_window_s: float | None = None,
    apply_TDC_calibration: bool = True,
    daughterboard_number: str | None = None,
    motherboard_number: str | None = None,
    spad: int | str = "average",
    max_files: int | None = None,
    subtract_background: bool = False,
    background_lags: int | Sequence[int] = bs.BACKGROUND_LAGS,
) -> str:
    """Unpack ORT data, compute blob-vs-blob delta-t, and save to feather.

    Writes ``path/processed/<name>_delta_t.feather`` with **one column per
    pixel pair** (named ``"ra,ca-rb,cb"``), padded to a common length; every
    requested pair gets a column even if it saw no coincidences.

    The differences are **streamed to disk as they are found**, not held until
    the end: the running buffer is flushed to
    ``processed/<name>_delta_t_parts/<name>_delta_t_<run>_<i>.feather``
    whenever it would exceed '_PART_MB', and the parts are concatenated
    into the final feather one record batch at a time. Memory therefore stays
    flat in the number of '.bin' files, and a run that dies at file 9 000 of
    10 000 leaves every completed part on disk - recover them with
    'combine_delta_t_parts'.

    The parts are **kept**. They are complete feathers in their own right and
    every read-side entry point takes a part folder as readily as a single
    file ('compute_and_save_delta_histogram' streams either one a record batch
    at a time), so they are the cheaper of the two things on disk to work
    from as well as the crash insurance. Delete the folder by hand once the
    combined feather has been checked.

    Each run tags its parts with its own ``<run>`` id and records them in a
    manifest, and a run refuses to start while another run's parts are still
    in the folder. Part numbering restarts at 0 every run, so without both of
    those a re-run would overwrite a dead run's first parts, leave its later
    ones, and the combine step would pool the two - counting some differences
    twice and inflating the feather.

    With ``apply_TDC_calibration=True`` (default) the board's density-test LUT
    converts timestamps to picoseconds *before* differencing, saving
    ``delta_ps``; otherwise raw ``delta_code`` is saved. Which one it is goes
    into the feather's schema metadata as ``delta_unit``.

    Parameters
    ----------
    path : str
        Path to the folder with the ORT '.bin' data files.
    pixels : Sequence[Sequence[tuple[int, int]]]
        Two groups ``[group_a, group_b]`` — the signal and idler blobs read
        off the hitmap, each a list of ``(row, col)`` pixels.
    rewrite : bool, optional
        Overwrite an existing '.feather' (and part folder) for this data set.
        A guard against accidental overwrites: the default is False, which
        raises instead. When True and something really is about to be
        destroyed, 'store.confirm_rewrite' names it and counts down five
        seconds first, so a stale ``rewrite=True`` can still be caught with
        Ctrl-C.
    nframes : int
        Number of frames stored in each '.bin' file.
    open_shutter_time_s : float
        Open shutter time, i.e. one reference-clock period, in seconds (e.g.
        ``100e-9`` at 10 MHz, ``500e-9`` at 2 MHz). Selects and validates the
        calibration LUT, and converts ``delta_window_s`` into raw code units
        when no LUT is applied. There is no default: the TDC's full scale is
        this window, so nothing downstream is right without it.
    tag : str, optional
        Filename fragment selecting the files. The default is ``'ORT'``.
    mode : str, optional
        ``'all_pairs'`` (default) or ``'1v1'``; see
        'calc_diff.calculate_differences'.
    delta_window_s : float | None, optional
        Keep only differences with ``abs(delta) <= delta_window_s`` while
        accumulating (bounds memory), in seconds. Converted to whichever unit
        the timestamps are in - ps with a LUT, raw codes without, using
        'core.timing.CODE_S'. The default is None (keep all).
    apply_TDC_calibration : bool, optional
        Apply the density-test TDC calibration (per-pixel code -> ps LUT)
        before differencing, so the saved ``delta_ps`` is already calibrated.
        The default is True, in which case the packaged board LUT is loaded
        via ``daughterboard_number`` / ``motherboard_number`` (a ValueError
        is raised if either is missing). Set False to save raw ``delta_code``
        and calibrate later.
    daughterboard_number : str | None, optional
        Daughterboard id (e.g. ``'D0'``) used to pick the packaged board LUT
        via 'tdc_calibration.load_board_lut'. The default is None.
    motherboard_number : str | None, optional
        Motherboard id (e.g. ``'M0'``) used with ``daughterboard_number`` to
        pick the packaged board LUT. The default is None.
    spad : int | str, optional
        Which SPAD's LUT to use: ``'average'`` (default) for the mean of the
        board's four SPAD LUTs, an int ``0..3`` for a specific SPAD, or a tag
        fragment like ``'SPAD0_S0'``. ORT cannot tell which micropixel fired,
        so ``'average'`` is the safe default; choose a SPAD only when the
        optics fix the micropixel. A missing specific SPAD falls back to the
        average (with a warning).
    max_files : int | None, optional
        Process at most this many files. The default is None (all files).
    subtract_background : bool, optional
        Also difference the *frame-shifted* pairs, so the accidental
        background can be measured rather than fitted. The default is False,
        which writes exactly the feather it always has. When True each pair
        gains one column per lag, named ``"<pair>@lag<k>"``, while lag 0 keeps
        the bare pair name - so the feather is a strict superset and every
        existing reader still sees the same columns. **It also multiplies the
        feather by roughly ``1 + len(background_lags)``**: prefer the counts
        path ('calculate_and_save_delta_counts'), where the same lags cost one
        small grid plane each. See
        'dapkel.functions.background_subtraction' and
        ``docs/guide/background_subtraction.md``.
    background_lags : int | Sequence[int], optional
        Frame offsets used for the accidental estimate when
        ``subtract_background`` is True. The default is ``(1, ..., 8)`` - which
        is 9x the feather, so consider fewer here. A bare ``int`` is taken as a
        single lag. Lag 0 is always included and must not appear here. Ignored
        when ``subtract_background`` is False.

    Returns
    -------
    str
        Path to the saved '.feather' file.

    Raises
    ------
    FileExistsError
        Raised when the '.feather' already exists and ``rewrite`` is False, or
        when an earlier run's parts are still in the part folder.
    ValueError
        Raised on a calibration request with no board, or a 0 / negative /
        repeated entry in ``background_lags``.
    """
    lags = _background_lags(background_lags) if subtract_background else [0]

    files = io.find_bin_files(path, tag)
    if max_files is not None:
        files = files[:max_files]

    name = os.path.basename(os.path.normpath(path))
    stem = f"{name}_delta_t"
    processed = store.processed_dir(path)
    out_path = os.path.join(processed, f"{stem}.feather")
    parts_dir = os.path.join(processed, f"{stem}{_PARTS_SUFFIX}")

    if os.path.isfile(out_path) and not rewrite:
        raise FileExistsError(
            f"{out_path} already exists. Pass rewrite=True to overwrite."
        )
    if rewrite:
        # Name the part folder too: a previous run's parts are also data, and
        # stale ones left in place would be folded into this run's output.
        store.confirm_rewrite([out_path, parts_dir])
        _clear_parts(parts_dir)
    else:
        # A run that died leaves its parts behind, and part numbering restarts
        # at 0, so this run would overwrite the dead one's first N parts and
        # leave the rest - a folder holding two runs at once. Stop instead of
        # producing a feather that counts some differences twice. Nothing is
        # deleted here; the parts are the crash insurance and they are the
        # caller's to keep, combine or discard.
        stale = _all_part_files(parts_dir)
        if stale:
            raise FileExistsError(
                f"{len(stale)} part file(s) from an earlier run are still in\n"
                f"  {parts_dir}\n"
                "Starting now would interleave both runs' parts. Either\n"
                "  - recover them:  combine_delta_t_parts(parts_dir)\n"
                "  - or discard them by re-running with rewrite=True."
            )

    # Resolve the density-test LUT to apply (per-pixel, code -> ps). None
    # leaves the raw codes in place (saved as delta_code).
    time_lut = None
    if apply_TDC_calibration:
        if daughterboard_number is not None and motherboard_number is not None:
            time_lut = tc.load_board_lut(
                daughterboard_number,
                motherboard_number,
                spad,
                open_shutter_time_s=open_shutter_time_s,
            )
        else:
            raise ValueError(
                "apply_TDC_calibration=True but no board was given. Pass "
                "daughterboard_number and motherboard_number, or set "
                "apply_TDC_calibration=False to save raw codes."
            )

    unit = "ps" if time_lut is not None else "code"
    # The window is a time; the timestamps are ps with a LUT and raw codes
    # without, so convert it into whichever unit the differences will be in.
    delta_window_native = _native_window(delta_window_s, unit)

    # Build the stable, ordered list of pair labels up front (matching the
    # keys 'calculate_differences' produces) so every requested pair gets a
    # column, in a deterministic order, even if it never fires. Shared with
    # the counts path through 'core.pairs' - the two stage-1 paths are only
    # comparable if they agree on which pairs exist and in what order.
    labels, all_pixels = pixel_pairs.pair_labels(pixels, mode)
    # Lag 0 first and under its bare name, then one block of columns per
    # shifted lag: a reader that knows nothing about lags finds precisely the
    # feather it has always found, at the front.
    columns = [_lag_column(lbl, k) for k in lags for lbl in labels]

    print(
        f"\n> > > Collecting delta-t (mode='{mode}', tag='{tag}') from "
        f"{len(files)} file(s) into {len(columns)} pixel-pair column(s)"
        + (f" at lags {lags}" if subtract_background else "")
        + f", streaming {_PART_MB:.0f} MB parts to "
        f"{os.path.basename(parts_dir)}{os.sep} < < <\n"
    )

    writer = _PartWriter(
        parts_dir,
        stem,
        columns,
        unit,
        _PART_MB,
        extra_meta=_lag_meta(lags, nframes) if subtract_background else None,
    )
    for fp in tqdm(files, desc=tag or "delta_t"):
        ts, _ = unpack(fp, nframes, compute_time_series=True)
        pixel_ts = _structure_pixel_timestamps(ts, all_pixels, time_lut)
        if subtract_background:
            deltas = {}
            for k in lags:
                lagged = bs.compute_lagged_differences(
                    pixel_ts, pixels, delta_window_native, mode, frame_lag=k
                )
                deltas.update(
                    {_lag_column(lbl, k): arr for lbl, arr in lagged.items()}
                )
        else:
            # Left as it was: lag 0 is tested to be identical to this, but the
            # plain path has no reason to route through the lag machinery.
            deltas = cd.calculate_differences(
                pixel_ts,
                pixels,
                delta_window_native=delta_window_native,
                mode=mode,
            )
        writer.add(deltas)
    parts = writer.close()

    if parts:
        print(
            f"\n> > > Combining {len(parts)} part file(s) into "
            f"{os.path.basename(out_path)} < < <"
        )
        total, _, _ = _combine_delta_feathers(parts, out_path)
    else:
        # Not one coincidence anywhere. Still write the empty table so the
        # requested pair columns and the unit are on record.
        total = _write_delta_feather(
            out_path,
            {lbl: np.empty(0, dtype=np.float64) for lbl in columns},
            columns,
            unit,
            _lag_meta(lags, nframes) if subtract_background else None,
        )

    print(
        f"\n> > > {total} differences across {len(columns)} pixel-pair "
        f"column(s) saved to {out_path} (unit '{unit}') < < <"
    )
    return out_path

combine_delta_t_parts

combine_delta_t_parts(parts_dir: str, out_path: str | None = None, *, rewrite: bool = False, run: str | None = None) -> str

Assemble one run's delta-t part files into its combined feather.

'calculate_and_save_timestamp_differences' does this itself at the end of a successful run. Call it by hand to recover a run that died partway through: every part flushed before the crash is complete and usable, so the differences up to the last flush are not lost.

Parameters:

Name	Type	Description	Default
parts_dir	str	The processed/<dataset>_delta_t_parts folder to assemble.	required
out_path	str | None	Where to write the combined '.feather'. The default is None, meaning processed/<dataset>_delta_t.feather - the same path the run would have written, beside the part folder.	None
rewrite	bool	Overwrite an existing combined feather. The default is False, which raises instead. When True, 'store.confirm_rewrite' counts down first.	False
run	str | None	Which run's parts to assemble, when the folder holds more than one. The default is None, which requires the folder to describe exactly one run and raises otherwise rather than pooling them.	None

Returns:

Type	Description
str	Path to the saved combined '.feather'.

Raises:

Type	Description
FileNotFoundError	Raised when the folder holds no part files.
FileExistsError	Raised when the output exists and rewrite is False.
ValueError	Raised when the folder holds several runs and run is None.

Source code in src/dapkel/functions/delta_t.py

def combine_delta_t_parts(
    parts_dir: str,
    out_path: str | None = None,
    *,
    rewrite: bool = False,
    run: str | None = None,
) -> str:
    """Assemble one run's delta-t part files into its combined feather.

    'calculate_and_save_timestamp_differences' does this itself at the end of
    a successful run. Call it by hand to recover a run that died partway
    through: every part flushed before the crash is complete and usable, so
    the differences up to the last flush are not lost.

    Parameters
    ----------
    parts_dir : str
        The ``processed/<dataset>_delta_t_parts`` folder to assemble.
    out_path : str | None, optional
        Where to write the combined '.feather'. The default is None, meaning
        ``processed/<dataset>_delta_t.feather`` - the same path the run would
        have written, beside the part folder.
    rewrite : bool, optional
        Overwrite an existing combined feather. The default is False, which
        raises instead. When True, 'store.confirm_rewrite' counts down first.
    run : str | None, optional
        Which run's parts to assemble, when the folder holds more than one.
        The default is None, which requires the folder to describe exactly
        one run and raises otherwise rather than pooling them.

    Returns
    -------
    str
        Path to the saved combined '.feather'.

    Raises
    ------
    FileNotFoundError
        Raised when the folder holds no part files.
    FileExistsError
        Raised when the output exists and ``rewrite`` is False.
    ValueError
        Raised when the folder holds several runs and ``run`` is None.
    """
    parts = _find_parts(parts_dir, run)
    if not parts:
        raise FileNotFoundError(
            f"No *_delta_t_<i>.feather part files in:\n  {parts_dir}"
        )

    if out_path is None:
        stem = os.path.basename(os.path.normpath(parts_dir)).removesuffix(
            _PARTS_SUFFIX
        )
        out_path = os.path.join(
            os.path.dirname(os.path.normpath(parts_dir)), f"{stem}.feather"
        )

    if os.path.isfile(out_path):
        if not rewrite:
            raise FileExistsError(
                f"{out_path} already exists. Pass rewrite=True to overwrite."
            )
        store.confirm_rewrite(out_path)

    print(
        f"\n> > > Assembling {len(parts)} part file(s) from "
        f"{os.path.basename(os.path.normpath(parts_dir))} into "
        f"{os.path.basename(out_path)} < < <"
    )
    total, labels, unit = _combine_delta_feathers(parts, out_path)
    print(
        f"\n> > > {total} differences across {len(labels)} pixel-pair "
        f"column(s) -> {out_path} (unit '{unit}') < < <"
    )
    return out_path

combine_delta_t_feathers

combine_delta_t_feathers(path: str, *, out_name: str | None = None, rewrite: bool = False, feathers: Sequence[str] | None = None, combined_dirname: str = 'combined') -> str

Pool the delta-t feathers under a folder into one combined feather.

Finds every *_delta_t.feather beneath path (each written by 'calculate_and_save_timestamp_differences' for one acquisition), pools them per pixel-pair column — concatenating the differences of matching pairs and taking the union of pairs across runs — and writes a single combined '.feather' into path/<combined_dirname>. The combined table keeps the one-column-per-pair, padded layout and the delta_unit metadata, so it plots exactly like a single-run feather (pass it to 'collect_and_plot_timestamp_differences' via feather_path).

The pooling streams a record batch at a time (see '_combine_delta_feathers'), so combining a hundred runs costs one batch of memory rather than the whole pooled data set. Per-run part folders are not picked up: their files end in _<i>.feather, so the *_delta_t.feather search skips them and nothing is counted twice.

Parameters:

Name	Type	Description	Default
path	str	Parent folder holding the per-run subfolders (e.g. .../SPDC_external). Searched recursively for *_delta_t.feather; anything already inside the combined subfolder is skipped.	required
out_name	str | None	Base name for the combined file (<out_name>_delta_t.feather). The default is None (<folder name>_combined).	None
rewrite	bool	Overwrite an existing combined feather. The default is False, which raises instead. When True, 'store.confirm_rewrite' counts down first.	False
feathers	Sequence[str] | None	Explicit list of feather paths to combine, bypassing the recursive search under path. The default is None (auto-discover).	None
combined_dirname	str	Name of the output subfolder created under path. The default is 'combined'.	'combined'

Returns:

Type	Description
str	Path to the saved combined '.feather'.

Raises:

Type	Description
FileNotFoundError	Raised when no *_delta_t.feather files are found.
FileExistsError	Raised when the combined feather exists and rewrite is False.
ValueError	Raised when the inputs mix delta_unit ('code' and 'ps').

Source code in src/dapkel/functions/delta_t.py

def combine_delta_t_feathers(
    path: str,
    *,
    out_name: str | None = None,
    rewrite: bool = False,
    feathers: Sequence[str] | None = None,
    combined_dirname: str = "combined",
) -> str:
    """Pool the delta-t feathers under a folder into one combined feather.

    Finds every ``*_delta_t.feather`` beneath ``path`` (each written by
    'calculate_and_save_timestamp_differences' for one acquisition), pools
    them per pixel-pair column — concatenating the differences of matching
    pairs and taking the union of pairs across runs — and writes a single
    combined '.feather' into ``path/<combined_dirname>``. The combined table
    keeps the one-column-per-pair, padded layout and the ``delta_unit``
    metadata, so it plots exactly like a single-run feather (pass it to
    'collect_and_plot_timestamp_differences' via ``feather_path``).

    The pooling streams a record batch at a time (see
    '_combine_delta_feathers'), so combining a hundred runs costs one batch of
    memory rather than the whole pooled data set. Per-run part folders are not
    picked up: their files end in ``_<i>.feather``, so the ``*_delta_t.feather``
    search skips them and nothing is counted twice.

    Parameters
    ----------
    path : str
        Parent folder holding the per-run subfolders (e.g.
        ``.../SPDC_external``). Searched recursively for
        ``*_delta_t.feather``; anything already inside the ``combined``
        subfolder is skipped.
    out_name : str | None, optional
        Base name for the combined file (``<out_name>_delta_t.feather``). The
        default is None (``<folder name>_combined``).
    rewrite : bool, optional
        Overwrite an existing combined feather. The default is False, which
        raises instead. When True, 'store.confirm_rewrite' counts down first.
    feathers : Sequence[str] | None, optional
        Explicit list of feather paths to combine, bypassing the recursive
        search under ``path``. The default is None (auto-discover).
    combined_dirname : str, optional
        Name of the output subfolder created under ``path``. The default is
        ``'combined'``.

    Returns
    -------
    str
        Path to the saved combined '.feather'.

    Raises
    ------
    FileNotFoundError
        Raised when no ``*_delta_t.feather`` files are found.
    FileExistsError
        Raised when the combined feather exists and ``rewrite`` is False.
    ValueError
        Raised when the inputs mix ``delta_unit`` ('code' and 'ps').
    """
    combined_dir = os.path.join(store.processed_dir(path), combined_dirname)
    combined_abs = os.path.abspath(combined_dir)

    if feathers is None:
        found = sorted(
            glob.glob(os.path.join(path, "**", "*_delta_t.feather"),
                      recursive=True)
        )
        # Never fold a previous combined output back into itself.
        feathers = [
            f for f in found
            if not os.path.abspath(f).startswith(combined_abs + os.sep)
        ]
    if not feathers:
        raise FileNotFoundError(
            f"No *_delta_t.feather files found under:\n  {path}"
        )

    name = out_name or f"{os.path.basename(os.path.normpath(path))}_combined"
    out_path = os.path.join(combined_dir, f"{name}_delta_t.feather")
    if os.path.isfile(out_path):
        if not rewrite:
            raise FileExistsError(
                f"{out_path} already exists. Pass rewrite=True to overwrite."
            )
        store.confirm_rewrite(out_path)

    print(
        f"\n> > > Combining {len(feathers)} delta-t feather(s) under "
        f"{path} into {os.path.basename(out_path)} < < <\n"
    )

    total, labels, unit = _combine_delta_feathers(feathers, out_path)

    print(
        f"\n> > > Combined {len(feathers)} feather(s): {total} differences "
        f"across {len(labels)} pixel-pair column(s) -> {out_path} "
        f"(unit '{unit}') < < <"
    )
    return out_path

compute_and_save_delta_histogram

compute_and_save_delta_histogram(path: str | None = None, *, feather_path: str | None = None, time_unit_s: float = CODE_S, bin_width_s: float = CODE_S, plot_window_s: float = 3e-09, pairs: Sequence[str] | None = None, subtract_background: bool = False, reuse: bool = True, quiet: bool = False) -> tuple[ndarray, ndarray, dict]

Stream a delta-t feather down to a saved histogram.

This is the expensive pass - it reads every record batch of every source - and it is the only thing standing between a hundreds-of-GB feather and a plot. The counts are saved to processed/<name>_delta_hist.npy with a '.meta.json' sidecar, so replotting, re-fitting or changing the background model costs a ~5 MB load rather than another full read.

The saved histogram is lossless for everything downstream: both fit models take (bin_centers, counts), and so does the plot. Binning at bin_width_s=timing.CODE_S (one TDC code) is the native resolution of the data, so nothing is thrown away at the default settings.

The sidecar records the binning and a fingerprint of the sources; a saved histogram is only reused when all of it matches, so changing bin_width_s, plot_window_s, time_unit_s or pairs - or rewriting the feather - re-streams automatically.

Parameters:

Name	Type	Description	Default
path	str | None	Path to the raw-data folder. Its processed/ holds the feather and receives the histogram. May be None when feather_path is given.	None
feather_path	str | None	Explicit delta-t '.feather' to read, or a *_delta_t_parts folder to read the parts directly (which skips needing the combined file at all). When given, path is ignored. The default is None.	None
time_unit_s	float	Seconds per raw TDC code, applied to a 'code' feather. The default is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent.	CODE_S
bin_width_s	float	Histogram bin width in seconds. The default is one TDC code.	CODE_S
plot_window_s	float	Histogram half-range in seconds around zero. The default is 3e-9 (3 ns).	3e-09
pairs	Sequence[str] | None	Restrict pooling to these pixel-pair column names. The default is None (pool every pair column). Frame-shifted columns are never pooled into the signal histogram.	None
subtract_background	bool	Histogram every frame lag separately instead of pooling one signal histogram, giving an (n_lags, n_bins) array with lag 0 - the ordinary histogram, bin-for-bin - as row 0. The default is False. Requires a feather written with calculate_and_save_timestamp_differences(..., subtract_background=True).	False
reuse	bool	Load a matching saved histogram instead of re-streaming. The default is True. Pass False to force the full pass.	True
quiet	bool	Suppress the progress bar and the summary. The default is False.	False

Returns:

Name	Type	Description
counts	ndarray	Counts per bin, or (n_lags, n_bins) with subtract_background.
edges	ndarray	The bin edges the counts belong to.
info	dict	The stats from '_stream_delta_histogram' plus name, root, sources, hist_path and reused. With subtract_background also lags and frames_used.

Raises:

Type	Description
ValueError	Raised when subtract_background is asked for on a feather with no frame-shifted columns.

Source code in src/dapkel/functions/delta_t.py

def compute_and_save_delta_histogram(
    path: str | None = None,
    *,
    feather_path: str | None = None,
    time_unit_s: float = timing.CODE_S,
    bin_width_s: float = timing.CODE_S,
    plot_window_s: float = 3e-9,
    pairs: Sequence[str] | None = None,
    subtract_background: bool = False,
    reuse: bool = True,
    quiet: bool = False,
) -> tuple[np.ndarray, np.ndarray, dict]:
    """Stream a delta-t feather down to a saved histogram.

    This is the expensive pass - it reads every record batch of every source -
    and it is the only thing standing between a hundreds-of-GB feather and a
    plot. The counts are saved to ``processed/<name>_delta_hist.npy`` with a
    '.meta.json' sidecar, so replotting, re-fitting or changing the background
    model costs a ~5 MB load rather than another full read.

    The saved histogram is *lossless for everything downstream*: both fit
    models take ``(bin_centers, counts)``, and so does the plot. Binning at
    ``bin_width_s=timing.CODE_S`` (one TDC code) is the native resolution of
    the data, so nothing is thrown away at the default settings.

    The sidecar records the binning and a fingerprint of the sources; a saved
    histogram is only reused when all of it matches, so changing
    ``bin_width_s``, ``plot_window_s``, ``time_unit_s`` or ``pairs`` - or
    rewriting the feather - re-streams automatically.

    Parameters
    ----------
    path : str | None, optional
        Path to the raw-data folder. Its ``processed/`` holds the feather and
        receives the histogram. May be None when ``feather_path`` is given.
    feather_path : str | None, optional
        Explicit delta-t '.feather' to read, or a ``*_delta_t_parts`` folder
        to read the parts directly (which skips needing the combined file at
        all). When given, ``path`` is ignored. The default is None.
    time_unit_s : float, optional
        Seconds per raw TDC code, applied to a ``'code'`` feather. The default
        is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent.
    bin_width_s : float, optional
        Histogram bin width in seconds. The default is one TDC code.
    plot_window_s : float, optional
        Histogram half-range in seconds around zero. The default is ``3e-9``
        (3 ns).
    pairs : Sequence[str] | None, optional
        Restrict pooling to these pixel-pair column names. The default is
        None (pool every pair column). Frame-shifted columns are never pooled
        into the signal histogram.
    subtract_background : bool, optional
        Histogram every frame lag separately instead of pooling one signal
        histogram, giving an ``(n_lags, n_bins)`` array with lag 0 - the
        ordinary histogram, bin-for-bin - as row 0. The default is False.
        Requires a feather written with
        ``calculate_and_save_timestamp_differences(..., subtract_background=True)``.
    reuse : bool, optional
        Load a matching saved histogram instead of re-streaming. The default
        is True. Pass False to force the full pass.
    quiet : bool, optional
        Suppress the progress bar and the summary. The default is False.

    Returns
    -------
    counts : np.ndarray
        Counts per bin, or ``(n_lags, n_bins)`` with ``subtract_background``.
    edges : np.ndarray
        The bin edges the counts belong to.
    info : dict
        The stats from '_stream_delta_histogram' plus ``name``, ``root``,
        ``sources``, ``hist_path`` and ``reused``. With
        ``subtract_background`` also ``lags`` and ``frames_used``.

    Raises
    ------
    ValueError
        Raised when ``subtract_background`` is asked for on a feather with no
        frame-shifted columns.
    """
    sources, name, root = _locate_delta_target(path, feather_path)
    lag_meta = _read_lag_meta(sources) if subtract_background else None
    if subtract_background and lag_meta is None:
        raise ValueError(
            "subtract_background=True needs a feather written with "
            "frame-shifted columns. Recompute stage 1 with "
            "'calculate_and_save_timestamp_differences(..., "
            "subtract_background=True)', or use the counts path."
        )
    lags = lag_meta["lags"] if lag_meta else None

    edges = _delta_hist_edges(bin_width_s, plot_window_s)
    signature = _hist_signature(
        sources,
        time_unit_s=time_unit_s,
        bin_width_s=bin_width_s,
        plot_window_s=plot_window_s,
        pairs=pairs,
        lags=lags,
    )
    hist_path = os.path.join(
        store.processed_dir(root, create=False), f"{name}{_HIST_SUFFIX}"
    )

    def _lag_info(info: dict) -> dict:
        """Attach the lag bookkeeping the subtraction needs."""
        if lag_meta is not None:
            info["lags"] = lags
            info["frames_used"] = [
                int(v) for v in bs.frames_per_lag(lag_meta["nframes"], lags)
            ]
        return info

    if reuse and os.path.isfile(hist_path):
        counts, meta = store.load_map(npy_path=hist_path)
        if (
            meta.get("signature") == signature
            and counts.shape[-1] == len(edges) - 1
        ):
            if not quiet:
                # 'n' lives under 'stats', not at the top level - reading it
                # from the wrong place made every cache hit report "0".
                n_saved = meta.get("stats", {}).get("n", 0)
                print(
                    f"\n> > > Reusing the saved delta-t histogram "
                    f"({n_saved} differences in the window) from "
                    f"{os.path.basename(hist_path)} - pass reuse=False to "
                    "re-stream < < <"
                )
            info = dict(meta.get("stats", {}))
            info.update(
                name=name,
                root=root,
                sources=sources,
                hist_path=hist_path,
                reused=True,
            )
            return counts, edges, _lag_info(info)

    if not quiet:
        print(
            f"\n> > > Streaming {len(sources)} delta-t feather(s) into a "
            f"{len(edges) - 1}-bin histogram "
            f"({bin_width_s * 1e12:.0f} ps bins, "
            f"+/-{plot_window_s * 1e9:.2f} ns)"
            + (f" at lags {lags}" if lags else "")
            + " < < <\n"
        )

    counts, stats = _stream_delta_histogram(
        sources,
        edges=edges,
        time_unit_s=time_unit_s,
        plot_window_s=plot_window_s,
        pairs=pairs,
        lags=lags,
        quiet=quiet,
    )

    hist_path = store.save_map(
        counts,
        root,
        kind=_KIND,
        file_name=f"{name}{_HIST_SUFFIX}",
        meta={
            "signature": signature,
            "stats": stats,
            "sources": [os.path.basename(s) for s in sources],
        },
        quiet=True,
    )

    if not quiet:
        print(
            f"\n> > > {stats['total']} difference(s) across "
            f"{stats['n_columns']} pair column(s); {stats['n']} inside "
            f"+/-{plot_window_s * 1e9:.2f} ns. Padding was "
            f"{100 * stats['padding_fraction']:.1f}% of the cells read. "
            f"Histogram saved to {hist_path} < < <"
        )
    return (
        counts,
        edges,
        _lag_info(
            {
                **stats,
                "name": name,
                "root": root,
                "sources": sources,
                "hist_path": hist_path,
                "reused": False,
            }
        ),
    )

calculate_and_save_delta_counts

calculate_and_save_delta_counts(path: str, pixels: Sequence[Sequence[Pixel]], rewrite: bool = False, *, nframes: int, open_shutter_time_s: float, tag: str = 'ORT', mode: str = 'all_pairs', delta_window_s: float | None = None, apply_TDC_calibration: bool = True, daughterboard_number: str | None = None, motherboard_number: str | None = None, spad: int | str = 'average', max_files: int | None = None, grid_s: float | None = None, support_s: float | None = None, subtract_background: bool = False, background_lags: int | Sequence[int] = BACKGROUND_LAGS, flush_every: int = _FLUSH_EVERY, resume: bool = True) -> str

Unpack ORT data and histogram delta-t directly - no feather.

The alternative to 'calculate_and_save_timestamp_differences'. Same loop (unpack -> per-pixel timestamps -> 'calc_diff.calculate_differences'), but the differences are binned onto a fixed native grid as they are found and only the counts are kept. Writes processed/<name>_delta_counts.npy: one row per pixel pair, one column per grid cell - and, with subtract_background, one plane per frame lag on top of that.

The size is set by the grid, not by the data - a 10 000-file run and a 100-file run produce the same ~100 MB artifact, against 3-135 GB of feather. That is what makes every pair plottable at once, and it makes re-binning free: 'rebin_delta_counts' turns the grid into any coarser histogram in memory, so bin_width_s and plot_window_s stop being compute-time decisions.

What it gives up. No per-event data survives, so there is no re-calibrating with a different LUT, no re-cutting delta_window_s, and no analysis that needs individual differences (time ordering, splitting by acquisition, higher-order correlations). Keep the feather until you have decided you never want those.

How exact it is. For raw codes (apply_TDC_calibration=False) the default grid is one code and the result is lossless - integers binned onto their own lattice. For calibrated data the LUT smears the lattice, so each difference is displaced by up to grid_s / 2 (~3.6 ps at the default), which is far below the detector jitter but is not bit-exact against the feather.

Crash insurance is one file rather than a part folder: the grid is rewritten every flush_every files, so a run that dies leaves a valid artifact covering everything up to its last checkpoint, and resume picks it up. Two runs cannot interleave, because there is only ever one file to write.

Parameters:

Name	Type	Description	Default
path	str	Path to the folder with the ORT '.bin' data files.	required
pixels	Sequence[Sequence[tuple[int, int]]]	Two groups [group_a, group_b] - the signal and idler blobs.	required
rewrite	bool	Overwrite a complete counts artifact for this data set. The default is False, which raises instead. An incomplete one is resumed rather than overwritten (see resume).	False
nframes	int	Number of frames stored in each '.bin' file. Part of the resume signature - it sets how much of each file is read, so a checkpoint only resumes against the same value.	required
open_shutter_time_s	float	Open shutter time, i.e. one reference-clock period, in seconds (e.g. 100e-9 at 10 MHz, 500e-9 at 2 MHz). Selects and validates the calibration LUT, and sets the default grid support - the TDC full scale is this window, so the accidental triangle reaches exactly this far.	required
tag	str	Filename fragment selecting the files. The default is 'ORT'.	'ORT'
mode	str	'all_pairs' (default) or '1v1'.	'all_pairs'
delta_window_s	float | None	Keep only differences with abs(delta) <= delta_window_s, in seconds. The default is None (keep all). Note this is frozen into the artifact - unlike the feather, it cannot be narrowed afterwards, so prefer leaving it None and cutting at plot time.	None
apply_TDC_calibration	bool	Apply the per-pixel code -> ps LUT before differencing. The default is True, which requires daughterboard_number and motherboard_number.	True
daughterboard_number	str | None	Daughterboard id (e.g. 'D0'). The default is None.	None
motherboard_number	str | None	Motherboard id (e.g. 'M0'). The default is None.	None
spad	int | str	Which SPAD's LUT to use. The default is 'average'.	'average'
max_files	int | None	Process at most this many files. The default is None (all files).	None
grid_s	float | None	Width of one grid cell, in seconds. The default is None, meaning one TDC code ('core.timing.CODE_S', ~71.6 ps) for raw-code data - where that is exact - and one tenth of a code (~7.2 ps) for calibrated data.	None
support_s	float | None	Half-range of the grid, in seconds. The default is None, which is twice the TDC full scale - read off the LUT when there is one, else from open_shutter_time_s. A difference cannot exceed the full scale, so that covers the whole physical support at any clock; differences beyond it are counted in n_outside.	None
subtract_background	bool	Also histogram the differences between frame-shifted pixel pairs, so the accidental background can be measured rather than fitted. The default is False, which writes exactly the 2D artifact it always has. When True the artifact gains a leading lag axis with lag 0 - the ordinary within-frame histogram, bin-for-bin identical to the 2D one - as plane 0; see 'dapkel.functions.background_subtraction' and docs/guide/background_subtraction.md. Costs one plane of grid per lag and one differencing pass per lag; the '.bin' decoding, which dominates the loop, is paid once either way.	False
background_lags	int | Sequence[int]	Frame offsets used for the accidental estimate when subtract_background is True. The default is (1, ..., 8). A bare int is taken as a single lag. Lag 0 is always included and must not appear here. Ignored when subtract_background is False.	BACKGROUND_LAGS
flush_every	int	Files between checkpoint writes. The default is 50.	_FLUSH_EVERY
resume	bool	Continue an incomplete artifact for the same inputs instead of starting over. The default is True.	True

Returns:

Type	Description
str	Path to the saved '.npy'.

Raises:

Type	Description
FileExistsError	Raised when a complete artifact exists and rewrite is False.
ValueError	Raised on a bad mode, unequal '1v1' groups, a calibration request with no board, a non-positive grid, or a 0 / negative / repeated entry in background_lags.

Source code in src/dapkel/functions/delta_t.py

def calculate_and_save_delta_counts(
    path: str,
    pixels: Sequence[Sequence[Pixel]],
    rewrite: bool = False,
    *,
    nframes: int,
    open_shutter_time_s: float,
    tag: str = "ORT",
    mode: str = "all_pairs",
    delta_window_s: float | None = None,
    apply_TDC_calibration: bool = True,
    daughterboard_number: str | None = None,
    motherboard_number: str | None = None,
    spad: int | str = "average",
    max_files: int | None = None,
    grid_s: float | None = None,
    support_s: float | None = None,
    subtract_background: bool = False,
    background_lags: int | Sequence[int] = bs.BACKGROUND_LAGS,
    flush_every: int = _FLUSH_EVERY,
    resume: bool = True,
) -> str:
    """Unpack ORT data and histogram delta-t directly - no feather.

    The alternative to 'calculate_and_save_timestamp_differences'. Same loop
    (unpack -> per-pixel timestamps -> 'calc_diff.calculate_differences'), but
    the differences are binned onto a fixed native grid as they are found and
    only the counts are kept. Writes
    ``processed/<name>_delta_counts.npy``: one row per pixel pair, one column
    per grid cell - and, with ``subtract_background``, one *plane* per frame
    lag on top of that.

    The size is set by the *grid*, not by the data - a 10 000-file run and a
    100-file run produce the same ~100 MB artifact, against 3-135 GB of
    feather. That is what makes every pair plottable at once, and it makes
    re-binning free: 'rebin_delta_counts' turns the grid into any coarser
    histogram in memory, so ``bin_width_s`` and ``plot_window_s`` stop being
    compute-time decisions.

    **What it gives up.** No per-event data survives, so there is no
    re-calibrating with a different LUT, no re-cutting ``delta_window_s``, and
    no analysis that needs individual differences (time ordering, splitting by
    acquisition, higher-order correlations). Keep the feather until you have
    decided you never want those.

    **How exact it is.** For raw codes (``apply_TDC_calibration=False``) the
    default grid is one code and the result is *lossless* - integers binned
    onto their own lattice. For calibrated data the LUT smears the lattice, so
    each difference is displaced by up to ``grid_s / 2`` (~3.6 ps at the
    default), which is far below the detector jitter but is not bit-exact
    against the feather.

    Crash insurance is one file rather than a part folder: the grid is
    rewritten every ``flush_every`` files, so a run that dies leaves a valid
    artifact covering everything up to its last checkpoint, and ``resume``
    picks it up. Two runs cannot interleave, because there is only ever one
    file to write.

    Parameters
    ----------
    path : str
        Path to the folder with the ORT '.bin' data files.
    pixels : Sequence[Sequence[tuple[int, int]]]
        Two groups ``[group_a, group_b]`` - the signal and idler blobs.
    rewrite : bool, optional
        Overwrite a *complete* counts artifact for this data set. The default
        is False, which raises instead. An incomplete one is resumed rather
        than overwritten (see ``resume``).
    nframes : int
        Number of frames stored in each '.bin' file. Part of the resume
        signature - it sets how much of each file is read, so a checkpoint
        only resumes against the same value.
    open_shutter_time_s : float
        Open shutter time, i.e. one reference-clock period, in seconds (e.g.
        ``100e-9`` at 10 MHz, ``500e-9`` at 2 MHz). Selects and validates the
        calibration LUT, and sets the default grid support - the TDC full scale
        is this window, so the accidental triangle reaches exactly this far.
    tag : str, optional
        Filename fragment selecting the files. The default is ``'ORT'``.
    mode : str, optional
        ``'all_pairs'`` (default) or ``'1v1'``.
    delta_window_s : float | None, optional
        Keep only differences with ``abs(delta) <= delta_window_s``, in seconds.
        The default is None (keep all). Note this is frozen into the artifact -
        unlike the feather, it cannot be narrowed afterwards, so prefer leaving
        it None and cutting at plot time.
    apply_TDC_calibration : bool, optional
        Apply the per-pixel code -> ps LUT before differencing. The default is
        True, which requires ``daughterboard_number`` and
        ``motherboard_number``.
    daughterboard_number : str | None, optional
        Daughterboard id (e.g. ``'D0'``). The default is None.
    motherboard_number : str | None, optional
        Motherboard id (e.g. ``'M0'``). The default is None.
    spad : int | str, optional
        Which SPAD's LUT to use. The default is ``'average'``.
    max_files : int | None, optional
        Process at most this many files. The default is None (all files).
    grid_s : float | None, optional
        Width of one grid cell, in seconds. The default is None, meaning one
        TDC code ('core.timing.CODE_S', ~71.6 ps) for raw-code data - where
        that is exact - and one tenth of a code (~7.2 ps) for calibrated data.
    support_s : float | None, optional
        Half-range of the grid, in seconds. The default is None, which is twice
        the TDC full scale - read off the LUT when there is one, else from
        ``open_shutter_time_s``. A difference cannot exceed the full scale, so
        that covers the whole physical support at any clock; differences beyond
        it are counted in ``n_outside``.
    subtract_background : bool, optional
        Also histogram the differences between *frame-shifted* pixel pairs, so
        the accidental background can be measured rather than fitted. The
        default is False, which writes exactly the 2D artifact it always has.
        When True the artifact gains a leading lag axis with lag 0 - the
        ordinary within-frame histogram, bin-for-bin identical to the 2D one -
        as plane 0; see 'dapkel.functions.background_subtraction' and
        ``docs/guide/background_subtraction.md``. Costs one plane of grid per
        lag and one differencing pass per lag; the '.bin' decoding, which
        dominates the loop, is paid once either way.
    background_lags : int | Sequence[int], optional
        Frame offsets used for the accidental estimate when
        ``subtract_background`` is True. The default is ``(1, ..., 8)``. A bare
        ``int`` is taken as a single lag. Lag 0 is always included and must not
        appear here. Ignored when ``subtract_background`` is False.
    flush_every : int, optional
        Files between checkpoint writes. The default is 50.
    resume : bool, optional
        Continue an incomplete artifact for the same inputs instead of
        starting over. The default is True.

    Returns
    -------
    str
        Path to the saved '.npy'.

    Raises
    ------
    FileExistsError
        Raised when a complete artifact exists and ``rewrite`` is False.
    ValueError
        Raised on a bad ``mode``, unequal ``'1v1'`` groups, a calibration
        request with no board, a non-positive grid, or a 0 / negative / repeated
        entry in ``background_lags``.
    """
    lags = _background_lags(background_lags) if subtract_background else [0]

    files = io.find_bin_files(path, tag)
    if max_files is not None:
        files = files[:max_files]

    labels, all_pixels = pixel_pairs.pair_labels(pixels, mode)
    time_lut = _resolve_time_lut(
        apply_TDC_calibration,
        daughterboard_number,
        motherboard_number,
        spad,
        open_shutter_time_s,
    )
    unit = "ps" if time_lut is not None else "code"

    if grid_s is None:
        grid_s = (
            timing.CODE_S / _GRID_SUBDIV if unit == "ps" else timing.CODE_S
        )
    if support_s is None:
        # The support is a *bound*, not a model parameter, so the shutter goes in
        # as a stated full scale rather than a guess at one: over-wide costs a
        # few empty cells, and 'resolve_tdc_full_scale_s' still prefers the LUT
        # when there is one. The half-base that actually gets fitted is resolved
        # at plot time, where being wrong matters.
        full_scale_s, _ = timing.resolve_tdc_full_scale_s(
            time_lut=time_lut, full_scale_s=open_shutter_time_s
        )
        support_s = _SUPPORT_FULL_SCALES * full_scale_s

    # 'n_bins' and 'zero_index' depend only on the support/grid ratio, so they
    # are unit-free; the cell width handed to '_accumulate_delta_counts' is not,
    # because it divides the raw difference values.
    n_bins, zero_index, _ = _delta_grid(grid_s, support_s)
    grid_native = (
        timing.to_ps(grid_s) if unit == "ps" else grid_s / timing.CODE_S
    )
    delta_window_native = _native_window(delta_window_s, unit)

    name = os.path.basename(os.path.normpath(path))
    counts_path = os.path.join(
        store.processed_dir(path), f"{name}{_COUNTS_SUFFIX}"
    )
    signature = _counts_signature(
        files,
        labels=labels,
        mode=mode,
        nframes=nframes,
        delta_window_s=delta_window_s,
        unit=unit,
        grid_s=grid_s,
        support_s=support_s,
        lags=lags,
    )

    # 2D without the lags, exactly as before, so an artifact written by an
    # older run still resumes and every downstream reader keeps working.
    counts = np.zeros(
        (len(labels), n_bins)
        if len(lags) == 1
        else (len(lags), len(labels), n_bins),
        dtype=np.int64,
    )
    rows = {label: i for i, label in enumerate(labels)}
    start = 0
    total = 0
    outside = 0

    existing = store.read_meta(counts_path) if os.path.isfile(counts_path) else {}
    if existing.get("signature") == signature:
        if existing.get("complete"):
            if not rewrite:
                raise FileExistsError(
                    f"{counts_path} already covers all {len(files)} file(s). "
                    "Pass rewrite=True to recompute."
                )
            store.confirm_rewrite([counts_path])
        elif resume:
            saved, _ = store.load_map(npy_path=counts_path)
            # The sidecar is moved into place after the array, so a crash
            # between the two moves can pair a newer array with an older
            # sidecar. Resuming on that would re-count files, and a silently
            # inflated coincidence total is worse than starting again.
            if (
                saved.shape == counts.shape
                and int(saved.sum()) == existing.get("total_in_grid")
            ):
                counts = saved
                start = int(existing.get("files_done", 0))
                total = int(existing.get("total", 0))
                outside = int(existing.get("n_outside", 0))
                print(
                    f"\n> > > Resuming {os.path.basename(counts_path)} at "
                    f"file {start} of {len(files)} < < <"
                )
            else:
                print(
                    f"\n  {os.path.basename(counts_path)} is inconsistent "
                    "with its sidecar (interrupted mid-write); starting over."
                )
    elif os.path.isfile(counts_path) and not rewrite:
        raise FileExistsError(
            f"{counts_path} exists but was computed from different inputs "
            "(files, pairs, grid or calibration). Pass rewrite=True to "
            "replace it."
        )
    elif os.path.isfile(counts_path):
        store.confirm_rewrite([counts_path])

    def _meta(files_done: int, complete: bool) -> dict:
        meta = {
            "signature": signature,
            "labels": labels,
            "unit": unit,
            "grid_s": float(grid_s),
            # The cell width in the *stored* unit - codes when uncalibrated.
            # Keeping both means the ps-per-code assumption is not baked in:
            # a code grid can be re-scaled at load time, exactly as the
            # feather path's 'time_unit_s' does.
            "grid_native": float(grid_native),
            "support_s": float(support_s),
            "open_shutter_time_s": float(open_shutter_time_s),
            "n_bins": int(n_bins),
            "zero_index": int(zero_index),
            "tag": tag,
            "files_done": int(files_done),
            "n_files": len(files),
            "complete": bool(complete),
            "total": int(total),
            "n_outside": int(outside),
            "total_in_grid": int(counts.sum()),
        }
        if subtract_background:
            # 'frames_used' is what 'background_subtraction' needs to rescale
            # the shifted lags. It is deterministic - lag k loses the last k
            # frames of every file - so it is derived rather than tracked, and
            # stays right on a resume.
            meta["lags"] = lags
            meta["frames_used"] = [
                int(v) for v in bs.frames_per_lag(nframes, lags, files_done)
            ]
        return meta

    print(
        f"\n> > > Histogramming delta-t straight from {len(files) - start} "
        f"'.bin' file(s) (mode='{mode}', tag='{tag}') onto a "
        f"{'x'.join(str(d) for d in counts.shape)} grid of "
        f"{grid_s * 1e12:.2f} ps cells spanning "
        f"+/-{support_s * 1e9:.1f} ns - "
        f"{counts.nbytes / 1024**2:.0f} MB, no feather"
        + (f"; lags {lags} < < <\n" if subtract_background else " < < <\n")
    )

    for i, fp in enumerate(
        tqdm(files[start:], desc=tag or "delta_counts"), start=start + 1
    ):
        ts, _ = unpack(fp, nframes, compute_time_series=True)
        pixel_ts = _structure_pixel_timestamps(ts, all_pixels, time_lut)
        if subtract_background:
            for li, lag in enumerate(lags):
                deltas = bs.compute_lagged_differences(
                    pixel_ts,
                    pixels,
                    delta_window_native,
                    mode,
                    frame_lag=lag,
                )
                seen, out = _accumulate_delta_counts(
                    counts[li],
                    deltas,
                    rows,
                    grid_native=grid_native,
                    zero_index=zero_index,
                )
                total += seen
                outside += out
        else:
            # Left exactly as it was: 'compute_lagged_differences' at lag 0 is
            # tested to be identical to this, but the plain path has no reason
            # to route through it.
            deltas = cd.calculate_differences(
                pixel_ts,
                pixels,
                delta_window_native=delta_window_native,
                mode=mode,
            )
            seen, out = _accumulate_delta_counts(
                counts,
                deltas,
                rows,
                grid_native=grid_native,
                zero_index=zero_index,
            )
            total += seen
            outside += out
        if i % flush_every == 0:
            _write_counts_checkpoint(counts_path, counts, _meta(i, False))

    _write_counts_checkpoint(counts_path, counts, _meta(len(files), True))

    print(
        f"\n> > > {total} difference(s) across {len(labels)} pixel-pair "
        f"row(s)"
        + (f" and {len(lags)} lag plane(s)" if subtract_background else "")
        + f" binned onto {n_bins} cells; {outside} outside "
        f"+/-{support_s * 1e9:.1f} ns. Saved to {counts_path} "
        f"(unit '{unit}') < < <"
    )
    if subtract_background:
        per_lag = counts.sum(axis=(1, 2))
        print(
            f"  lag 0: {per_lag[0]}; shifted lags: {per_lag[1:].mean():.0f} on "
            f"average -> a raw excess of {per_lag[0] - per_lag[1:].mean():.0f} "
            "over the whole grid."
        )
    if outside:
        print(
            f"  NOTE: {100 * outside / max(total, 1):.2f}% of differences "
            "fell outside the grid and are counted but not binned. Raise "
            "support_s to keep them."
        )
    return counts_path

load_delta_counts

load_delta_counts(path: str | None = None, *, counts_path: str | None = None, pairs: Sequence[str] | None = None, pool: bool = True, time_unit_s: float = CODE_S) -> tuple[ndarray, ndarray, dict]

Load a saved counts grid and its cell centres.

Parameters:

Name	Type	Description	Default
path	str | None	Data folder whose processed/ holds the artifact. May be None when counts_path is given.	None
counts_path	str | None	Explicit '.npy' to load, ignoring path. The default is None.	None
pairs	Sequence[str] | None	Restrict to these pixel-pair labels. The default is None (all).	None
pool	bool	Sum the selected pairs into one histogram. The default is True; pass False to keep the pair axis.	True
time_unit_s	float	Seconds per raw TDC code, used to place the cell centres of an uncalibrated ('code') grid. The default is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent. Ignored for a calibrated grid, whose cells already carry a time.	CODE_S

Returns:

Name	Type	Description
counts	ndarray	Pooled counts, or per-pair counts when pool=False. An artifact written with subtract_background keeps its leading lag axis, so the shape is (n_bins) / (n_pairs, n_bins) without lags and (n_lags, n_bins) / (n_lags, n_pairs, n_bins) with them - lag 0 first.
centers	ndarray	Grid cell centres, in picoseconds whatever the stored unit.
info	dict	The sidecar, plus name, root, counts_path, the selected labels, the effective cell_s and any requested pairs that were missing. A lagged artifact's sidecar also carries lags and frames_used, which is what 'background_subtraction.subtract_background' needs.

Raises:

Type	Description
ValueError	Raised when neither path nor counts_path is given.
FileNotFoundError	Raised when the artifact does not exist.

Source code in src/dapkel/functions/delta_t.py

def load_delta_counts(
    path: str | None = None,
    *,
    counts_path: str | None = None,
    pairs: Sequence[str] | None = None,
    pool: bool = True,
    time_unit_s: float = timing.CODE_S,
) -> tuple[np.ndarray, np.ndarray, dict]:
    """Load a saved counts grid and its cell centres.

    Parameters
    ----------
    path : str | None, optional
        Data folder whose ``processed/`` holds the artifact. May be None when
        ``counts_path`` is given.
    counts_path : str | None, optional
        Explicit '.npy' to load, ignoring ``path``. The default is None.
    pairs : Sequence[str] | None, optional
        Restrict to these pixel-pair labels. The default is None (all).
    pool : bool, optional
        Sum the selected pairs into one histogram. The default is True; pass
        False to keep the pair axis.
    time_unit_s : float, optional
        Seconds per raw TDC code, used to place the cell centres of an
        *uncalibrated* (``'code'``) grid. The default is 'core.timing.CODE_S'
        (~71.6 ps), which is clock-independent. Ignored for a calibrated grid,
        whose cells already carry a time.

    Returns
    -------
    counts : np.ndarray
        Pooled counts, or per-pair counts when ``pool=False``. An artifact
        written with ``subtract_background`` keeps its leading lag axis, so
        the shape is ``(n_bins)`` / ``(n_pairs, n_bins)`` without lags and
        ``(n_lags, n_bins)`` / ``(n_lags, n_pairs, n_bins)`` with them - lag 0
        first.
    centers : np.ndarray
        Grid cell centres, in **picoseconds** whatever the stored unit.
    info : dict
        The sidecar, plus ``name``, ``root``, ``counts_path``, the selected
        ``labels``, the effective ``cell_s`` and any requested pairs that
        were ``missing``. A lagged artifact's sidecar also carries ``lags``
        and ``frames_used``, which is what
        'background_subtraction.subtract_background' needs.

    Raises
    ------
    ValueError
        Raised when neither ``path`` nor ``counts_path`` is given.
    FileNotFoundError
        Raised when the artifact does not exist.
    """
    if counts_path is None:
        if path is None:
            raise ValueError("Provide either path or counts_path.")
        name = os.path.basename(os.path.normpath(path))
        counts_path = os.path.join(
            store.processed_dir(path, create=False), f"{name}{_COUNTS_SUFFIX}"
        )
    counts_path = os.path.abspath(counts_path)
    counts, meta = store.load_map(npy_path=counts_path)

    name = os.path.basename(counts_path).removesuffix(_COUNTS_SUFFIX)
    root = os.path.dirname(os.path.dirname(counts_path))

    # The pair axis is the last but one, whether or not a lag axis precedes it.
    pair_axis = counts.ndim - 2

    labels = list(meta.get("labels", []))
    missing: list[str] = []
    if pairs is not None:
        wanted = [p for p in pairs if p in labels]
        missing = [p for p in pairs if p not in labels]
        counts = np.take(
            counts, [labels.index(p) for p in wanted], axis=pair_axis
        )
        labels = wanted

    zero_index = int(meta.get("zero_index", (counts.shape[-1] - 1) // 2))
    grid_native = float(meta.get("grid_native", meta.get("grid_ps", 1.0)))
    # A code grid's cells are codes; scale them to ps here rather than at save
    # time, so a better ps-per-code can be applied without recomputing.
    cell_ps = (
        grid_native
        if meta.get("unit", "ps") == "ps"
        else grid_native * timing.to_ps(time_unit_s)
    )
    centers = (np.arange(counts.shape[-1]) - zero_index) * cell_ps

    if pool:
        counts = counts.sum(axis=pair_axis)

    info = {
        **meta,
        "name": name,
        "root": root,
        "counts_path": counts_path,
        "labels": labels,
        "missing": missing,
        "cell_s": timing.to_s(cell_ps),
    }
    if not meta.get("complete", True):
        print(
            f"  NOTE: {os.path.basename(counts_path)} is incomplete - "
            f"{meta.get('files_done')} of {meta.get('n_files')} file(s)."
        )
    return counts, centers, info

rebin_delta_counts

rebin_delta_counts(counts: ndarray, centers: ndarray, *, bin_width_s: float, plot_window_s: float) -> tuple[ndarray, ndarray]

Reduce a native counts grid onto a coarser histogram, in memory.

Each grid cell is assigned whole to the target bin its centre falls in. That is exact whenever bin_width_s is a whole multiple of the grid cell and the two lattices line up; otherwise a cell straddling a target edge goes entirely one way, which is the rebinning, not the data. Bin on a multiple of the grid when comparing against the feather path.

Parameters:

Name	Type	Description	Default
counts	ndarray	1D pooled counts on the native grid, from 'load_delta_counts'.	required
centers	ndarray	The grid cell centres in ps, from 'load_delta_counts'.	required
bin_width_s	float	Target bin width, in seconds.	required
plot_window_s	float	Target half-range, in seconds.	required

Returns:

Name	Type	Description
counts	ndarray	Counts per target bin.
edges	ndarray	The target bin edges in ps - the same '_delta_hist_edges' the feather path uses, so the two histograms are directly comparable.

Source code in src/dapkel/functions/delta_t.py

def rebin_delta_counts(
    counts: np.ndarray,
    centers: np.ndarray,
    *,
    bin_width_s: float,
    plot_window_s: float,
) -> tuple[np.ndarray, np.ndarray]:
    """Reduce a native counts grid onto a coarser histogram, in memory.

    Each grid cell is assigned whole to the target bin its *centre* falls in.
    That is exact whenever ``bin_width_s`` is a whole multiple of the grid
    cell and the two lattices line up; otherwise a cell straddling a target
    edge goes entirely one way, which is the rebinning, not the data. Bin on a
    multiple of the grid when comparing against the feather path.

    Parameters
    ----------
    counts : np.ndarray
        1D pooled counts on the native grid, from 'load_delta_counts'.
    centers : np.ndarray
        The grid cell centres in ps, from 'load_delta_counts'.
    bin_width_s : float
        Target bin width, in seconds.
    plot_window_s : float
        Target half-range, in seconds.

    Returns
    -------
    counts : np.ndarray
        Counts per target bin.
    edges : np.ndarray
        The target bin edges in ps - the same '_delta_hist_edges' the feather
        path uses, so the two histograms are directly comparable.
    """
    edges = _delta_hist_edges(bin_width_s, plot_window_s)
    counts = np.asarray(counts).ravel()
    idx = np.searchsorted(edges, centers, side="right") - 1
    inside = (idx >= 0) & (idx < len(edges) - 1)
    return (
        np.bincount(
            idx[inside], weights=counts[inside], minlength=len(edges) - 1
        ).astype(np.int64),
        edges,
    )

gaussian

gaussian(x: ndarray, amp: float, mu: float, sigma: float, bkg: float) -> ndarray

Gaussian peak on a flat background.

Parameters:

Name	Type	Description	Default
x	ndarray	Input positions.	required
amp	float	Peak amplitude above the background.	required
mu	float	Peak centre.	required
sigma	float	Standard deviation of the Gaussian.	required
bkg	float	Flat background level.	required

Returns:

Type	Description
ndarray	Model values amp * exp(-(x - mu)**2 / (2 sigma**2)) + bkg.

Source code in src/dapkel/functions/delta_t.py

def gaussian(
    x: np.ndarray, amp: float, mu: float, sigma: float, bkg: float
) -> np.ndarray:
    """Gaussian peak on a flat background.

    Parameters
    ----------
    x : np.ndarray
        Input positions.
    amp : float
        Peak amplitude above the background.
    mu : float
        Peak centre.
    sigma : float
        Standard deviation of the Gaussian.
    bkg : float
        Flat background level.

    Returns
    -------
    np.ndarray
        Model values ``amp * exp(-(x - mu)**2 / (2 sigma**2)) + bkg``.
    """
    return amp * np.exp(-((x - mu) ** 2) / (2.0 * sigma**2)) + bkg

gaussian_plus_triangle

gaussian_plus_triangle(x: ndarray, amp: float, mu: float, sigma: float, bkg: float, half_base: float) -> ndarray

Gaussian peak on a symmetric triangular background.

The correct accidental model for ORT: uncorrelated free-running phases difference into a triangle (zero at +/- half_base), with correlated pairs adding a Gaussian at mu ~ 0. Why: docs/ort_triangle_background.md.

Parameters:

Name	Type	Description	Default
x	ndarray	Input positions (timestamp differences, same units as half_base).	required
amp	float	Gaussian amplitude above the triangle.	required
mu	float	Gaussian centre.	required
sigma	float	Gaussian standard deviation.	required
bkg	float	Triangle height at its apex (x = 0).	required
half_base	float	Half-width of the triangle base; the background is zero for abs(x) >= half_base (physically ~one oscillator period, ~100 ns).	required

Returns:

Type	Description
ndarray	Model values amp * exp(-(x-mu)^2 / 2 sigma^2) + bkg * tri(x) where tri(x) = max(0, 1 - abs(x) / half_base).

Source code in src/dapkel/functions/delta_t.py

def gaussian_plus_triangle(
    x: np.ndarray,
    amp: float,
    mu: float,
    sigma: float,
    bkg: float,
    half_base: float,
) -> np.ndarray:
    """Gaussian peak on a symmetric triangular background.

    The correct accidental model for ORT: uncorrelated free-running phases
    difference into a triangle (zero at ``+/- half_base``), with correlated
    pairs adding a Gaussian at ``mu ~ 0``. Why:
    ``docs/ort_triangle_background.md``.

    Parameters
    ----------
    x : np.ndarray
        Input positions (timestamp differences, same units as ``half_base``).
    amp : float
        Gaussian amplitude above the triangle.
    mu : float
        Gaussian centre.
    sigma : float
        Gaussian standard deviation.
    bkg : float
        Triangle height at its apex (``x = 0``).
    half_base : float
        Half-width of the triangle base; the background is zero for
        ``abs(x) >= half_base`` (physically ~one oscillator period, ~100 ns).

    Returns
    -------
    np.ndarray
        Model values ``amp * exp(-(x-mu)^2 / 2 sigma^2) + bkg * tri(x)`` where
        ``tri(x) = max(0, 1 - abs(x) / half_base)``.
    """
    tri = np.clip(1.0 - np.abs(x) / half_base, 0.0, None)
    return amp * np.exp(-((x - mu) ** 2) / (2.0 * sigma**2)) + bkg * tri

two_gaussians_on_triangle

two_gaussians_on_triangle(x: ndarray, amp_n: float, mu: float, sigma_n: float, amp_b: float, dmu: float, sigma_ratio: float, bkg: float, half_base: float) -> ndarray

Narrow + broad Gaussian on a triangular background.

A single Gaussian fitted to the ORT coincidence peak reports the area-weighted mixture of two physically distinct timing populations, which is why it lands near 500-700 ps while the designer's focused-laser number is tens of ps:

• narrow - photons absorbed inside the SPAD depletion region, timed by the avalanche build-up alone (tens of ps, at or below this TDC's ~77 ps LSB, so it is resolution-limited here);
• broad - photons absorbed in the field-free epi/substrate below the junction, which reach the multiplication region by diffusion. The true shape of that component is one-sided exponential, not Gaussian; a Gaussian is the cheap two-parameter stand-in that separates the scales without committing to the tail shape. If the residuals stay structured in the flanks, that is the signal to move to an exponentially modified Gaussian instead.

The broad width is parameterised as a ratio rather than an absolute sigma so the two components cannot swap roles mid-fit: sigma_ratio >= 1 makes sigma_n the narrow one by construction. sigma_ratio = 1 collapses the model to a single Gaussian, so a data set that does not need two components degenerates gracefully instead of fitting noise.

Parameters:

Name	Type	Description	Default
x	ndarray	Input positions (timestamp differences, same units as half_base).	required
amp_n	float	Amplitude of the narrow component above the triangle.	required
mu	float	Centre of the narrow component.	required
sigma_n	float	Standard deviation of the narrow component.	required
amp_b	float	Amplitude of the broad component above the triangle.	required
dmu	float	Centre of the broad component relative to mu. Free because a diffusion tail is one-sided: with unequal signal/idler wavelengths the two arms' tails differ and the broad component sits off-centre.	required
sigma_ratio	float	sigma_broad / sigma_n, constrained to >= 1 by the fitter.	required
bkg	float	Triangle height at its apex (x = 0).	required
half_base	float	Half-width of the triangle base; the background is zero for abs(x) >= half_base (~one oscillator period, ~100 ns).	required

Returns:

Type	Description
ndarray	Model values: narrow Gaussian + broad Gaussian + bkg * tri(x), with tri(x) = max(0, 1 - abs(x) / half_base).

Source code in src/dapkel/functions/delta_t.py

def two_gaussians_on_triangle(
    x: np.ndarray,
    amp_n: float,
    mu: float,
    sigma_n: float,
    amp_b: float,
    dmu: float,
    sigma_ratio: float,
    bkg: float,
    half_base: float,
) -> np.ndarray:
    """Narrow + broad Gaussian on a triangular background.

    A single Gaussian fitted to the ORT coincidence peak reports the
    area-weighted mixture of two physically distinct timing populations, which
    is why it lands near 500-700 ps while the designer's focused-laser number
    is tens of ps:

    * **narrow** - photons absorbed *inside* the SPAD depletion region, timed
      by the avalanche build-up alone (tens of ps, at or below this TDC's
      ~77 ps LSB, so it is resolution-limited here);
    * **broad** - photons absorbed in the field-free epi/substrate *below* the
      junction, which reach the multiplication region by diffusion. The true
      shape of that component is one-sided exponential, not Gaussian; a
      Gaussian is the cheap two-parameter stand-in that separates the *scales*
      without committing to the tail shape. If the residuals stay structured
      in the flanks, that is the signal to move to an exponentially modified
      Gaussian instead.

    The broad width is parameterised as a **ratio** rather than an absolute
    sigma so the two components cannot swap roles mid-fit: ``sigma_ratio >= 1``
    makes ``sigma_n`` the narrow one by construction. ``sigma_ratio = 1``
    collapses the model to a single Gaussian, so a data set that does not need
    two components degenerates gracefully instead of fitting noise.

    Parameters
    ----------
    x : np.ndarray
        Input positions (timestamp differences, same units as ``half_base``).
    amp_n : float
        Amplitude of the narrow component above the triangle.
    mu : float
        Centre of the narrow component.
    sigma_n : float
        Standard deviation of the narrow component.
    amp_b : float
        Amplitude of the broad component above the triangle.
    dmu : float
        Centre of the broad component *relative to* ``mu``. Free because a
        diffusion tail is one-sided: with unequal signal/idler wavelengths the
        two arms' tails differ and the broad component sits off-centre.
    sigma_ratio : float
        ``sigma_broad / sigma_n``, constrained to ``>= 1`` by the fitter.
    bkg : float
        Triangle height at its apex (``x = 0``).
    half_base : float
        Half-width of the triangle base; the background is zero for
        ``abs(x) >= half_base`` (~one oscillator period, ~100 ns).

    Returns
    -------
    np.ndarray
        Model values: narrow Gaussian + broad Gaussian + ``bkg * tri(x)``,
        with ``tri(x) = max(0, 1 - abs(x) / half_base)``.
    """
    sigma_b = sigma_n * sigma_ratio
    tri = np.clip(1.0 - np.abs(x) / half_base, 0.0, None)
    narrow = amp_n * np.exp(-((x - mu) ** 2) / (2.0 * sigma_n**2))
    broad = amp_b * np.exp(-((x - mu - dmu) ** 2) / (2.0 * sigma_b**2))
    return narrow + broad + bkg * tri

fit_gaussian_on_triangle

fit_gaussian_on_triangle(bin_centers: ndarray, counts: ndarray, half_base: float, fit_half_base: bool = True) -> dict

Fit a Gaussian peak on a triangular background (Kelpie ORT model).

Use this rather than 'fit_gaussian_peak' for ORT delta-t over a wide window: the background is triangular, and a flat-background fit lets the peak swallow the triangle and reports an inflated sigma.

Parameters:

Name	Type	Description	Default
bin_centers	ndarray	Histogram bin centres (should span roughly the full +/- half_base).	required
counts	ndarray	Histogram counts.	required
half_base	float	Initial triangle half-base (~one oscillator period; ~100 ns in ps, or ~1300 in code units). Refined when fit_half_base is True.	required
fit_half_base	bool	Let half_base float in the fit. The default is True.	True

Returns:

Type	Description
dict	Fit results with keys amp, mu, sigma, sigma_err, bkg (triangle apex), half_base, and ok.

Raises:

Type	Description
ValueError	Raised when bin_centers and counts differ in length.

Source code in src/dapkel/functions/delta_t.py

def fit_gaussian_on_triangle(
    bin_centers: np.ndarray,
    counts: np.ndarray,
    half_base: float,
    fit_half_base: bool = True,
) -> dict:
    """Fit a Gaussian peak on a triangular background (Kelpie ORT model).

    Use this rather than 'fit_gaussian_peak' for ORT delta-t over a wide
    window: the background is triangular, and a flat-background fit lets the
    peak swallow the triangle and reports an inflated sigma.

    Parameters
    ----------
    bin_centers : np.ndarray
        Histogram bin centres (should span roughly the full +/- half_base).
    counts : np.ndarray
        Histogram counts.
    half_base : float
        Initial triangle half-base (~one oscillator period; ~100 ns in ps, or
        ~1300 in code units). Refined when ``fit_half_base`` is True.
    fit_half_base : bool, optional
        Let ``half_base`` float in the fit. The default is True.

    Returns
    -------
    dict
        Fit results with keys ``amp``, ``mu``, ``sigma``, ``sigma_err``,
        ``bkg`` (triangle apex), ``half_base``, and ``ok``.

    Raises
    ------
    ValueError
        Raised when ``bin_centers`` and ``counts`` differ in length.
    """
    x = np.asarray(bin_centers, dtype=np.float64)
    y = np.asarray(counts, dtype=np.float64)
    _check_hist_lengths(x, y)

    result = {
        "amp": np.nan,
        "mu": 0.0,
        "sigma": np.nan,
        "sigma_err": np.nan,
        "bkg": np.nan,
        "half_base": half_base,
        "ok": False,
    }
    if len(x) < 6:
        return result

    # Triangle apex seed: median just off zero (skip the peak-contaminated
    # central few bins), scaled back up to x=0 along the triangle.
    dx = float(np.median(np.diff(x))) if len(x) > 1 else 1.0
    near = (np.abs(x) > 3 * dx) & (np.abs(x) < 0.15 * half_base)
    apex0 = (
        float(
            np.median(y[near]) / (1.0 - np.median(np.abs(x[near])) / half_base)
        )
        if np.any(near)
        else float(np.median(y))
    )
    tri_seed = apex0 * np.clip(1.0 - np.abs(x) / half_base, 0.0, None)
    excess = np.clip(y - tri_seed, 0.0, None)
    amp0 = float(np.max(excess)) if excess.size else float(np.max(y))
    mu0 = float(x[np.argmax(excess)]) if excess.size else 0.0
    # Seed sigma from the excess *near the peak only*: far from zero the
    # apex-seed mismatch leaves a broad positive residual that would otherwise
    # inflate the moment to ~half_base and push p0 outside the bounds.
    core = np.abs(x - mu0) < 0.05 * half_base
    w = np.where(core, excess, 0.0)
    sigma0 = (
        float(np.sqrt(np.sum(w * (x - mu0) ** 2) / np.sum(w)))
        if w.sum() > 0
        else 5 * dx
    )
    if not sigma0 > 0:
        sigma0 = 5 * dx

    lo = [0.0, x.min(), 0.0, 0.0, 0.5 * half_base]
    hi = [
        np.inf,
        x.max(),
        0.5 * half_base,
        np.inf,
        2.0 * half_base if fit_half_base else half_base * 1.0001,
    ]
    if not fit_half_base:
        lo[4] = half_base * 0.9999

    # Clamp every seed strictly inside its bounds (curve_fit rejects p0 == or
    # outside a bound).
    p0 = [amp0, mu0, sigma0, apex0, half_base]
    p0 = [
        min(
            max(v, lo_i + 1e-9 * (abs(lo_i) + 1)),
            hi_i - 1e-9 * (abs(hi_i) + 1),
        )
        for v, lo_i, hi_i in zip(p0, lo, hi, strict=True)
    ]

    try:
        popt, pcov = curve_fit(
            gaussian_plus_triangle, x, y, p0=p0, bounds=(lo, hi), maxfev=20000
        )
    except (RuntimeError, ValueError):
        return result

    amp, mu, sigma, bkg, hb = popt
    if amp <= 0 or not np.isfinite(sigma) or sigma <= 0:
        return result
    sigma_err = (
        float(np.sqrt(pcov[2, 2])) if np.all(np.isfinite(pcov)) else np.nan
    )
    result.update(
        amp=float(amp),
        mu=float(mu),
        sigma=float(abs(sigma)),
        sigma_err=sigma_err,
        bkg=float(bkg),
        half_base=float(hb),
        ok=True,
    )
    return result

fit_two_gaussians_on_triangle

fit_two_gaussians_on_triangle(bin_centers: ndarray, counts: ndarray, half_base: float, fit_half_base: bool = True, poisson_weights: bool = True) -> dict

Fit a narrow + broad Gaussian on a triangular background.

Separates the two timing populations a single Gaussian averages together (see 'two_gaussians_on_triangle'), so the fast avalanche core can be compared against a focused-laser jitter number and the slow diffusion component quantified instead of inflating one sigma.

Seeded from 'fit_gaussian_on_triangle': its sigma is the area-weighted mixture of both components, which makes it a sound seed for the broad one (it carries most of the area) and the narrow component is then seeded from what is left over in the central few bins. Fitting eight free parameters from cold seeds does not converge reliably; from these it does.

Parameters:

Name	Type	Description	Default
bin_centers	ndarray	Histogram bin centres. Should span well past the peak so the triangle is constrained - a window of a few ns cannot separate the components.	required
counts	ndarray	Histogram counts.	required
half_base	float	Initial triangle half-base (~one oscillator period; ~100 ns in ps, or ~1300 in code units). Refined when fit_half_base is True.	required
fit_half_base	bool	Let half_base float in the fit. The default is True.	True
poisson_weights	bool	Weight each bin by 1 / sqrt(max(counts, 1)). The default is True, and it matters more here than for the single-Gaussian fits: the narrow core is a small fraction of the total counts, so unweighted least squares is dominated by the triangle and can leave the core badly fitted while still reporting convergence. Set False to reproduce the unweighted behaviour of the other fitters.	True

Returns:

Type	Description
dict	amp_n, mu, sigma_n, sigma_n_err, amp_b, dmu, sigma_b, sigma_b_err, sigma_ratio, bkg, half_base, frac_narrow (narrow share of the total peak area), fwhm_n, fwhm_b and ok. sigma_ratio at its lower bound of 1.0 means the data did not support two components.

Raises:

Type	Description
ValueError	Raised when bin_centers and counts differ in length.

Source code in src/dapkel/functions/delta_t.py

def fit_two_gaussians_on_triangle(
    bin_centers: np.ndarray,
    counts: np.ndarray,
    half_base: float,
    fit_half_base: bool = True,
    poisson_weights: bool = True,
) -> dict:
    """Fit a narrow + broad Gaussian on a triangular background.

    Separates the two timing populations a single Gaussian averages together
    (see 'two_gaussians_on_triangle'), so the fast avalanche core can be
    compared against a focused-laser jitter number and the slow diffusion
    component quantified instead of inflating one sigma.

    Seeded from 'fit_gaussian_on_triangle': its sigma is the area-weighted
    mixture of both components, which makes it a sound seed for the broad one
    (it carries most of the area) and the narrow component is then seeded from
    what is left over in the central few bins. Fitting eight free parameters
    from cold seeds does not converge reliably; from these it does.

    Parameters
    ----------
    bin_centers : np.ndarray
        Histogram bin centres. Should span well past the peak so the triangle
        is constrained - a window of a few ns cannot separate the components.
    counts : np.ndarray
        Histogram counts.
    half_base : float
        Initial triangle half-base (~one oscillator period; ~100 ns in ps, or
        ~1300 in code units). Refined when ``fit_half_base`` is True.
    fit_half_base : bool, optional
        Let ``half_base`` float in the fit. The default is True.
    poisson_weights : bool, optional
        Weight each bin by ``1 / sqrt(max(counts, 1))``. The default is True,
        and it matters more here than for the single-Gaussian fits: the narrow
        core is a small fraction of the total counts, so unweighted least
        squares is dominated by the triangle and can leave the core badly
        fitted while still reporting convergence. Set False to reproduce the
        unweighted behaviour of the other fitters.

    Returns
    -------
    dict
        ``amp_n``, ``mu``, ``sigma_n``, ``sigma_n_err``, ``amp_b``, ``dmu``,
        ``sigma_b``, ``sigma_b_err``, ``sigma_ratio``, ``bkg``, ``half_base``,
        ``frac_narrow`` (narrow share of the total peak area), ``fwhm_n``,
        ``fwhm_b`` and ``ok``. ``sigma_ratio`` at its lower bound of 1.0 means
        the data did not support two components.

    Raises
    ------
    ValueError
        Raised when ``bin_centers`` and ``counts`` differ in length.
    """
    x = np.asarray(bin_centers, dtype=np.float64)
    y = np.asarray(counts, dtype=np.float64)
    _check_hist_lengths(x, y)

    result = {
        "amp_n": np.nan,
        "mu": 0.0,
        "sigma_n": np.nan,
        "sigma_n_err": np.nan,
        "amp_b": np.nan,
        "dmu": 0.0,
        "sigma_b": np.nan,
        "sigma_b_err": np.nan,
        "sigma_ratio": np.nan,
        "bkg": np.nan,
        "half_base": half_base,
        "frac_narrow": np.nan,
        "fwhm_n": np.nan,
        "fwhm_b": np.nan,
        "ok": False,
    }
    if len(x) < 10:
        return result

    dx = float(np.median(np.diff(x))) if len(x) > 1 else 1.0

    seed = fit_gaussian_on_triangle(
        x, y, half_base, fit_half_base=fit_half_base
    )
    if seed["ok"]:
        mu0 = seed["mu"]
        sigma_b0 = seed["sigma"]
        apex0 = seed["bkg"]
        hb0 = seed["half_base"]
        amp_b0 = seed["amp"]
    else:
        mu0 = float(x[np.argmax(y)])
        sigma_b0 = 10.0 * dx
        apex0 = float(np.median(y))
        hb0 = half_base
        amp_b0 = max(float(np.max(y)) - apex0, 1.0)

    # Narrow seed: whatever the one-Gaussian fit failed to account for in the
    # central few bins. That under-shoot at x=0 is the whole reason for this
    # model, so it is exactly the right thing to seed the second component on.
    resid = y - gaussian_plus_triangle(x, amp_b0, mu0, sigma_b0, apex0, hb0)
    core = np.abs(x - mu0) <= 3.0 * dx
    amp_n0 = float(np.max(resid[core])) if np.any(core) else 0.0
    amp_n0 = max(amp_n0, 0.05 * amp_b0)
    # The core is expected at or below the TDC LSB, so seed it a couple of bins
    # wide rather than anywhere near the mixture sigma.
    sigma_n0 = max(1.5 * dx, 1e-3 * sigma_b0)
    ratio0 = float(np.clip(sigma_b0 / sigma_n0, 1.05, 400.0))

    lo = [
        0.0,
        x.min(),
        dx / 10.0,
        0.0,
        -0.05 * half_base,
        1.0,
        0.0,
        0.5 * half_base,
    ]
    hi = [
        np.inf,
        x.max(),
        0.5 * half_base,
        np.inf,
        0.05 * half_base,
        500.0,
        np.inf,
        2.0 * half_base if fit_half_base else half_base * 1.0001,
    ]
    if not fit_half_base:
        lo[7] = half_base * 0.9999

    p0 = [amp_n0, mu0, sigma_n0, amp_b0, 0.0, ratio0, apex0, hb0]
    p0 = [
        min(
            max(v, lo_i + 1e-9 * (abs(lo_i) + 1)),
            hi_i - 1e-9 * (abs(hi_i) + 1),
        )
        for v, lo_i, hi_i in zip(p0, lo, hi, strict=True)
    ]

    kwargs = {}
    if poisson_weights:
        kwargs["sigma"] = np.sqrt(np.maximum(y, 1.0))
        kwargs["absolute_sigma"] = False

    try:
        popt, pcov = curve_fit(
            two_gaussians_on_triangle,
            x,
            y,
            p0=p0,
            bounds=(lo, hi),
            maxfev=40000,
            **kwargs,
        )
    except (RuntimeError, ValueError):
        return result

    amp_n, mu, sigma_n, amp_b, dmu, ratio, bkg, hb = popt
    if not np.isfinite(sigma_n) or sigma_n <= 0 or (amp_n <= 0 and amp_b <= 0):
        return result

    sigma_n = abs(sigma_n)
    sigma_b = sigma_n * ratio

    finite_cov = np.all(np.isfinite(pcov))
    sigma_n_err = float(np.sqrt(pcov[2, 2])) if finite_cov else np.nan
    if finite_cov:
        # sigma_b = sigma_n * ratio, so both parameters and their covariance
        # enter; taking sqrt(pcov[5,5]) alone would understate it.
        var_b = (
            ratio**2 * pcov[2, 2]
            + sigma_n**2 * pcov[5, 5]
            + 2.0 * ratio * sigma_n * pcov[2, 5]
        )
        sigma_b_err = float(np.sqrt(var_b)) if var_b > 0 else np.nan
    else:
        sigma_b_err = np.nan

    area_n = max(amp_n, 0.0) * sigma_n
    area_b = max(amp_b, 0.0) * sigma_b
    frac_narrow = (
        float(area_n / (area_n + area_b)) if (area_n + area_b) > 0 else np.nan
    )

    fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0))
    result.update(
        amp_n=float(amp_n),
        mu=float(mu),
        sigma_n=float(sigma_n),
        sigma_n_err=sigma_n_err,
        amp_b=float(amp_b),
        dmu=float(dmu),
        sigma_b=float(sigma_b),
        sigma_b_err=sigma_b_err,
        sigma_ratio=float(ratio),
        bkg=float(bkg),
        half_base=float(hb),
        frac_narrow=frac_narrow,
        fwhm_n=float(fwhm * sigma_n),
        fwhm_b=float(fwhm * sigma_b),
        ok=True,
    )
    return result

fit_gaussian_peak

fit_gaussian_peak(bin_centers: ndarray, counts: ndarray, fit_window: float | None = None) -> dict

Fit a Gaussian-plus-background peak with scipy.curve_fit.

The flat background is seeded from the histogram tails and the peak centre/width from the background-subtracted moments, then all four parameters (amp, mu, sigma, bkg) are refined by a least-squares fit of 'gaussian'.

Parameters:

Name	Type	Description	Default
bin_centers	ndarray	Histogram bin centres.	required
counts	ndarray	Histogram counts.	required
fit_window	float | None	If given, restrict the fit to abs(x - x_peak) <= fit_window around the tallest bin. The default is None (use every bin).	None

Returns:

Type	Description
dict	Fit results with keys amp, mu, sigma, sigma_err, bkg, and ok (False when the fit did not converge to a real peak).

Source code in src/dapkel/functions/delta_t.py

def fit_gaussian_peak(
    bin_centers: np.ndarray,
    counts: np.ndarray,
    fit_window: float | None = None,
) -> dict:
    """Fit a Gaussian-plus-background peak with ``scipy.curve_fit``.

    The flat background is seeded from the histogram tails and the peak
    centre/width from the background-subtracted moments, then all four
    parameters (``amp``, ``mu``, ``sigma``, ``bkg``) are refined by a
    least-squares fit of 'gaussian'.

    Parameters
    ----------
    bin_centers : np.ndarray
        Histogram bin centres.
    counts : np.ndarray
        Histogram counts.
    fit_window : float | None, optional
        If given, restrict the fit to ``abs(x - x_peak) <= fit_window``
        around the tallest bin. The default is None (use every bin).

    Returns
    -------
    dict
        Fit results with keys ``amp``, ``mu``, ``sigma``, ``sigma_err``,
        ``bkg``, and ``ok`` (False when the fit did not converge to a real
        peak).
    """
    x = np.asarray(bin_centers, dtype=np.float64)
    y = np.asarray(counts, dtype=np.float64)

    n = len(x)
    tail = np.r_[y[: n // 3], y[-n // 3 :]]
    bkg = float(np.median(tail)) if tail.size else 0.0

    x_peak = float(x[np.argmax(y)])
    if fit_window is not None:
        sel = np.abs(x - x_peak) <= fit_window
        x, y = x[sel], y[sel]

    result = {
        "amp": float(np.max(y) - bkg),
        "mu": x_peak,
        "sigma": np.nan,
        "sigma_err": np.nan,
        "bkg": bkg,
        "ok": False,
    }
    if len(x) < 4:
        return result

    w = np.clip(y - bkg, 0, None)
    if w.sum() <= 0:
        return result
    mu0 = float(np.sum(w * x) / np.sum(w))
    sigma0 = float(np.sqrt(np.sum(w * (x - mu0) ** 2) / np.sum(w)))
    if not sigma0 > 0:
        sigma0 = (x[1] - x[0]) if len(x) > 1 else 1.0
    p0 = [float(np.max(y) - bkg), mu0, sigma0, bkg]

    try:
        popt, pcov = curve_fit(
            gaussian,
            x,
            y,
            p0=p0,
            bounds=(
                [0.0, x.min(), 0.0, 0.0],
                [np.inf, x.max(), np.inf, np.inf],
            ),
            maxfev=10000,
        )
    except (RuntimeError, ValueError):
        return result

    amp, mu, sigma, bkg_fit = popt
    if amp <= 0 or not np.isfinite(sigma) or sigma <= 0:
        return result
    sigma_err = (
        float(np.sqrt(pcov[2, 2])) if np.all(np.isfinite(pcov)) else np.nan
    )
    result.update(
        amp=float(amp),
        mu=float(mu),
        sigma=float(abs(sigma)),
        sigma_err=sigma_err,
        bkg=float(bkg_fit),
        ok=True,
    )
    return result

collect_and_plot_timestamp_differences

collect_and_plot_timestamp_differences(path: str | None = None, *, open_shutter_time_s: float | None = None, time_unit_s: float = CODE_S, bin_width_s: float = CODE_S, plot_window_s: float = 3e-09, fit_window_s: float | None = 1e-09, background: str = 'flat', full_scale_s: float | None = None, feather_path: str | None = None, pairs: Sequence[str] | None = None, label: str = 'spdc', cmap_title: str | None = None, reuse_histogram: bool = True, subtract_background: bool = False, peak_window_s: float = 1e-09, wide_rebin: int = 25, support_s: float | None = None) -> tuple[ndarray, ndarray, dict]

Histogram the saved delta-t feather, fit the peak, and plot it.

Pools the selected pair columns of a saved feather, converts to picoseconds (per the delta_unit metadata), histograms, fits the coincidence peak and saves the figure. Old single-column feathers (delta_code / delta_ps) are still read.

The feather is streamed through 'compute_and_save_delta_histogram', one Arrow record batch at a time, and never held in memory - a full run's combined feather is far too large to load. The counts are identical to histogramming the whole pooled array at once; only the memory profile differs. The histogram is cached in processed/, so re-fitting or changing the background model does not re-read the feather.

Locate the feather either from path (its processed/) or directly via feather_path, which ignores path entirely - that is how a combined feather is plotted. Either may point at a *_delta_t_parts folder to plot straight from the part files.

With subtract_background the histogram that gets fitted is the lag-0 minus shifted-frames residual: the accidental pedestal is measured from frame-shifted pairs and subtracted, rather than modelled. All the lags are histogrammed in the same single streaming pass. Needs a feather written with calculate_and_save_timestamp_differences(..., subtract_background=True).

Parameters:

Name	Type	Description	Default
path	str | None	Path to the raw-data folder (its processed holds the feather, results receives the plot). May be None when feather_path is given. Required otherwise.	None
open_shutter_time_s	float | None	Open shutter time, i.e. one reference-clock period, in seconds. Used as the fallback for the TDC full scale, hence for the triangle half-base and the default streaming support. The default is None, which is fine only for background='flat' without subtract_background; the triangle models need a full scale from somewhere.	None
time_unit_s	float	Seconds per raw TDC code, applied to a delta_code feather. The default is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent.	CODE_S
bin_width_s	float	Histogram bin width in seconds. The default is one TDC code.	CODE_S
plot_window_s	float	Histogram/plot half-range in seconds around zero. The default is 3e-9 (3 ns).	3e-09
fit_window_s	float | None	Fit half-range in seconds around the tallest bin. Only used for background='flat'. The default is 1e-9 (1 ns).	1e-09
background	str	Peak-and-background model. 'flat' (default) fits Gaussian + constant over fit_window_s — the legacy daplis model. 'triangle' fits Gaussian + triangular background via 'fit_gaussian_on_triangle' over the whole histogram, which is the physically correct model for the Kelpie ORT free-running phase and does not inflate sigma. 'two_gaussians' adds a second, broader Gaussian on the same triangle via 'fit_two_gaussians_on_triangle' — after TDC calibration the peak shows a narrow core on wide tails that one Gaussian averages into a single inflated sigma, and this separates the two. The plot then also draws both components and the bare triangle. For 'triangle' and 'two_gaussians' set plot_window_s wide (ideally the whole full scale) so the triangle is well constrained.	'flat'
full_scale_s	float | None	TDC full scale in seconds, which is the triangle half-base: the accidental background peaks at 0 and reaches zero at +/- full_scale_s. The default is None, which resolves it from open_shutter_time_s via 'core.timing.resolve_tdc_full_scale_s'. Only used when background is 'triangle' or 'two_gaussians'; the fit refines it. Pass it explicitly when the shutter is longer than one clock period, because then the timestamps fold and the triangle is narrower than the shutter.	None
feather_path	str | None	Explicit feather to read and plot directly, or a *_delta_t_parts folder to read the parts; when given, path is ignored and the raw-data folder is not needed (the plot is saved beside the feather). The default is None (derive from path).	None
pairs	Sequence[str] | None	Restrict pooling to these pixel-pair column names (e.g. ["16,16-20,20"]). The default is None (pool every pair column). Ignored for old single-column feathers.	None
label	str	Short label for the title and output file ('spdc' / 'hbt'). The default is 'spdc'.	'spdc'
cmap_title	str | None	Optional override for the plot title. The default is None.	None
reuse_histogram	bool	Reuse a saved histogram whose binning and sources match, instead of re-streaming the feather. The default is True. Pass False to force the full pass.	True
subtract_background	bool	Subtract the frame-shifted accidentals before fitting, and write the three-panel diagnostic beside the fit. The default is False. Requires a feather written with subtract_background=True, and only makes sense with background='flat' - there is no triangle left to fit. Note the streamed histogram then spans support_s rather than plot_window_s, since the residual out at the edge of the full scale is the evidence that the pedestal was accidental.	False
peak_window_s	float	Half-range the model-free coincidence count is summed over, in seconds. subtract_background only. The default is 1e-9 (1 ns).	1e-09
wide_rebin	int	Histogram bins per bin of the diagnostic's two full-support panels. subtract_background only. The default is 25.	25
support_s	float | None	Half-range the lag histograms are streamed over, in seconds. subtract_background only. The default is None, i.e. one TDC full scale - the whole physical support of a difference.	None

Returns:

Name	Type	Description
counts	ndarray	Counts per bin - the histogram that was fitted and plotted, or the residual when subtract_background is set.
centers	ndarray	Bin centres in picoseconds, to re-fit without re-reading the feather.
fit	dict	The fit results plus fwhm (2.355 sigma), jitter_per_detector (sigma / sqrt(2)), car (peak-to-background ratio) and n. With subtract_background also model_free, subtraction and subtraction_png.

Raises:

Type	Description
FileNotFoundError	Raised when the delta-t '.feather' cannot be found.
ValueError	Raised when subtract_background is asked for on a feather with no frame-shifted columns, or alongside a non-flat background model.

Source code in src/dapkel/functions/delta_t.py

def collect_and_plot_timestamp_differences(
    path: str | None = None,
    *,
    open_shutter_time_s: float | None = None,
    time_unit_s: float = timing.CODE_S,
    bin_width_s: float = timing.CODE_S,
    plot_window_s: float = 3e-9,
    fit_window_s: float | None = 1e-9,
    background: str = "flat",
    full_scale_s: float | None = None,
    feather_path: str | None = None,
    pairs: Sequence[str] | None = None,
    label: str = "spdc",
    cmap_title: str | None = None,
    reuse_histogram: bool = True,
    subtract_background: bool = False,
    peak_window_s: float = 1e-9,
    wide_rebin: int = 25,
    support_s: float | None = None,
) -> tuple[np.ndarray, np.ndarray, dict]:
    """Histogram the saved delta-t feather, fit the peak, and plot it.

    Pools the selected pair columns of a saved feather, converts to
    picoseconds (per the ``delta_unit`` metadata), histograms, fits the
    coincidence peak and saves the figure. Old single-column feathers
    (``delta_code`` / ``delta_ps``) are still read.

    The feather is **streamed** through 'compute_and_save_delta_histogram',
    one Arrow record batch at a time, and never held in memory - a full run's
    combined feather is far too large to load. The counts are identical to
    histogramming the whole pooled array at once; only the memory profile
    differs. The histogram is cached in ``processed/``, so re-fitting or
    changing the background model does not re-read the feather.

    Locate the feather either from ``path`` (its ``processed/``) or directly
    via ``feather_path``, which ignores ``path`` entirely - that is how a
    combined feather is plotted. Either may point at a ``*_delta_t_parts``
    folder to plot straight from the part files.

    With ``subtract_background`` the histogram that gets fitted is the *lag-0
    minus shifted-frames* residual: the accidental pedestal is measured from
    frame-shifted pairs and subtracted, rather than modelled. All the lags are
    histogrammed in the same single streaming pass. Needs a feather written
    with ``calculate_and_save_timestamp_differences(...,
    subtract_background=True)``.

    Parameters
    ----------
    path : str | None, optional
        Path to the raw-data folder (its ``processed`` holds the
        feather, ``results`` receives the plot). May be None when
        ``feather_path`` is given. Required otherwise.
    open_shutter_time_s : float | None, optional
        Open shutter time, i.e. one reference-clock period, in seconds. Used as
        the fallback for the TDC full scale, hence for the triangle half-base
        and the default streaming support. The default is None, which is fine
        only for ``background='flat'`` without ``subtract_background``; the
        triangle models need a full scale from somewhere.
    time_unit_s : float, optional
        Seconds per raw TDC code, applied to a ``delta_code`` feather. The
        default is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent.
    bin_width_s : float, optional
        Histogram bin width in seconds. The default is one TDC code.
    plot_window_s : float, optional
        Histogram/plot half-range in seconds around zero. The default is
        ``3e-9`` (3 ns).
    fit_window_s : float | None, optional
        Fit half-range in seconds around the tallest bin. Only used for
        ``background='flat'``. The default is ``1e-9`` (1 ns).
    background : str, optional
        Peak-and-background model. ``'flat'`` (default) fits
        ``Gaussian + constant`` over ``fit_window_s`` — the legacy daplis
        model. ``'triangle'`` fits ``Gaussian + triangular background`` via
        'fit_gaussian_on_triangle' over the whole histogram, which is the
        physically correct model for the Kelpie ORT free-running phase and
        does not inflate sigma. ``'two_gaussians'`` adds a second, broader
        Gaussian on the same triangle via
        'fit_two_gaussians_on_triangle' — after TDC calibration the peak shows
        a narrow core on wide tails that one Gaussian averages into a single
        inflated sigma, and this separates the two. The plot then also draws
        both components and the bare triangle. For ``'triangle'`` and
        ``'two_gaussians'`` set ``plot_window_s`` wide (ideally the whole full
        scale) so the triangle is well constrained.
    full_scale_s : float | None, optional
        TDC full scale in seconds, which **is** the triangle half-base: the
        accidental background peaks at 0 and reaches zero at
        ``+/- full_scale_s``. The default is None, which resolves it from
        ``open_shutter_time_s`` via 'core.timing.resolve_tdc_full_scale_s'.
        Only used when ``background`` is ``'triangle'`` or ``'two_gaussians'``;
        the fit refines it. Pass it explicitly when the shutter is longer than
        one clock period, because then the timestamps fold and the triangle is
        narrower than the shutter.
    feather_path : str | None, optional
        Explicit feather to read and plot directly, or a ``*_delta_t_parts``
        folder to read the parts; when given, ``path`` is ignored and the
        raw-data folder is not needed (the plot is saved beside the feather).
        The default is None (derive from ``path``).
    pairs : Sequence[str] | None, optional
        Restrict pooling to these pixel-pair column names (e.g.
        ``["16,16-20,20"]``). The default is None (pool every pair column).
        Ignored for old single-column feathers.
    label : str, optional
        Short label for the title and output file (``'spdc'`` / ``'hbt'``).
        The default is ``'spdc'``.
    cmap_title : str | None, optional
        Optional override for the plot title. The default is None.
    reuse_histogram : bool, optional
        Reuse a saved histogram whose binning and sources match, instead of
        re-streaming the feather. The default is True. Pass False to force
        the full pass.
    subtract_background : bool, optional
        Subtract the frame-shifted accidentals before fitting, and write the
        three-panel diagnostic beside the fit. The default is False. Requires a
        feather written with ``subtract_background=True``, and only makes sense
        with ``background='flat'`` - there is no triangle left to fit. Note the
        streamed histogram then spans ``support_s`` rather than
        ``plot_window_s``, since the residual out at the edge of the full scale
        is the evidence that the pedestal was accidental.
    peak_window_s : float, optional
        Half-range the model-free coincidence count is summed over, in seconds.
        ``subtract_background`` only. The default is ``1e-9`` (1 ns).
    wide_rebin : int, optional
        Histogram bins per bin of the diagnostic's two full-support panels.
        ``subtract_background`` only. The default is 25.
    support_s : float | None, optional
        Half-range the lag histograms are streamed over, in seconds.
        ``subtract_background`` only. The default is None, i.e. one TDC full
        scale - the whole physical support of a difference.

    Returns
    -------
    counts : np.ndarray
        Counts per bin - the histogram that was fitted and plotted, or the
        residual when ``subtract_background`` is set.
    centers : np.ndarray
        Bin centres in picoseconds, to re-fit without re-reading the feather.
    fit : dict
        The fit results plus ``fwhm`` (2.355 sigma), ``jitter_per_detector``
        (sigma / sqrt(2)), ``car`` (peak-to-background ratio) and ``n``. With
        ``subtract_background`` also ``model_free``, ``subtraction`` and
        ``subtraction_png``.

    Raises
    ------
    FileNotFoundError
        Raised when the delta-t '.feather' cannot be found.
    ValueError
        Raised when ``subtract_background`` is asked for on a feather with no
        frame-shifted columns, or alongside a non-flat background model.
    """
    if subtract_background:
        _check_subtraction_model(background)

    # The triangle half-base is the TDC full scale, so it has to come from the
    # calibration or the clock - never from a fixed 100 ns. Resolved only when
    # something actually needs it: a plain flat fit does not, and demanding a
    # clock for it would break every narrow-window call that never cared.
    needs_full_scale = background in ("triangle", "two_gaussians")
    if needs_full_scale or (subtract_background and support_s is None):
        full_scale_s, _ = timing.resolve_tdc_full_scale_s(
            full_scale_s=full_scale_s,
            open_shutter_time_s=open_shutter_time_s,
        )
        if support_s is None:
            support_s = full_scale_s

    # The subtraction needs the whole support in one histogram: the zoom is a
    # slice of it and the diagnostic's wide panels are a regrouping, so the
    # expensive streaming pass is still paid exactly once.
    hist_window_s = (
        _whole_bins(support_s, bin_width_s)
        if subtract_background
        else plot_window_s
    )
    counts, edges, info = compute_and_save_delta_histogram(
        path,
        feather_path=feather_path,
        time_unit_s=time_unit_s,
        bin_width_s=bin_width_s,
        plot_window_s=hist_window_s,
        pairs=pairs,
        subtract_background=subtract_background,
        reuse=reuse_histogram,
    )
    name = info["name"]
    results_dir = store.results_dir(info["root"], _KIND, create=False)
    centers = (edges[:-1] + edges[1:]) / 2
    n_in_window = int(info["n"])

    # Comb guard: warn if the data granularity is coarser than the bin.
    gran_s = info.get("granularity_s")
    if gran_s is not None and gran_s > bin_width_s * 1.5:
        print(
            f"  WARNING: delta granularity is ~{gran_s * 1e12:.0f} ps but "
            f"bin_width={bin_width_s * 1e12:.0f} ps -> the histogram will be a "
            "comb. Bin coarser, or ensure the fine timing is populated "
            "in the acquisition."
        )

    sub = None
    if subtract_background:
        sub = _subtract_lag_background(
            counts,
            centers,
            info["frames_used"],
            name=name,
            root=info["root"],
            background=background,
            bin_width_s=bin_width_s,
            zoom_window_s=plot_window_s,
            peak_window_s=peak_window_s,
            wide_rebin=wide_rebin,
            label=label,
        )
        counts = sub["residual"]
        centers = sub["centers"]
        # The model-free count, not the bin total: the residual's bins are
        # differences of two histograms, so summing them *is* the coincidence
        # number and anything outside the peak is noise around zero.
        n_in_window = int(round(sub["model_free"]["n"]))

    counts, centers, fit = _fit_and_plot_delta(
        counts,
        centers,
        name=name,
        results_dir=results_dir,
        bin_width_s=bin_width_s,
        plot_window_s=plot_window_s,
        fit_window_s=fit_window_s,
        background=background,
        full_scale_s=full_scale_s if needs_full_scale else None,
        label=label,
        cmap_title=cmap_title,
        n_in_window=n_in_window,
        stem="delta_t_sub" if subtract_background else "delta_t",
        subtracted=subtract_background,
    )
    if sub is not None:
        fit["model_free"] = sub["model_free"]
        fit["subtraction"] = sub
        fit["subtraction_png"] = sub["png_path"]
    return counts, centers, fit

collect_and_plot_delta_counts

collect_and_plot_delta_counts(path: str | None = None, *, counts_path: str | None = None, open_shutter_time_s: float | None = None, time_unit_s: float = CODE_S, bin_width_s: float = CODE_S, plot_window_s: float = 3e-09, fit_window_s: float | None = 1e-09, background: str = 'flat', full_scale_s: float | None = None, pairs: Sequence[str] | None = None, label: str = 'spdc', cmap_title: str | None = None, snap_to_grid: bool = True, subtract_background: bool = False, peak_window_s: float = 1e-09, wide_rebin: int = 25) -> tuple[ndarray, ndarray, dict]

Rebin a saved counts grid, fit the peak, and plot it.

The counterpart of 'collect_and_plot_timestamp_differences' for the feather-free path, and deliberately the same function underneath: only the histogram differs, so the two can be compared directly.

With subtract_background the histogram that gets fitted is the lag-0 minus shifted-frames residual instead of the raw one - same fit, same figure, on a distribution whose accidental pedestal has been measured away rather than modelled. Needs an artifact written with calculate_and_save_delta_counts(..., subtract_background=True).

Nothing is read from disk except one ~100 MB array, and the rebinning is in memory, so changing bin_width_s or plot_window_s costs milliseconds rather than a full pass over the data - and every pixel pair is available at once instead of whichever columns were streamed.

Parameters:

Name	Type	Description	Default
path	str | None	Data folder whose processed/ holds the counts artifact and whose results/ receives the plot. May be None when counts_path is given.	None
counts_path	str | None	Explicit counts '.npy', ignoring path. The default is None.	None
open_shutter_time_s	float | None	Open shutter time, i.e. one reference-clock period, in seconds. The fallback source for the TDC full scale, hence the triangle half-base. The default is None, which is fine only for background='flat'.	None
time_unit_s	float	Seconds per raw TDC code, for an uncalibrated grid. The default is 'core.timing.CODE_S' (~71.6 ps), which is clock-independent.	CODE_S
bin_width_s	float	Histogram bin width in seconds. The default is one TDC code. Use a whole multiple of the artifact's grid cell for an exact rebin.	CODE_S
plot_window_s	float	Histogram/plot half-range in seconds around zero. The default is 3e-9 (3 ns). Unlike the feather path, widening this is free - it is a slice of what is already loaded.	3e-09
fit_window_s	float | None	Fit half-range in seconds around the tallest bin. 'flat' background only. The default is 1e-9 (1 ns).	1e-09
background	str	'flat' (default), 'triangle' or 'two_gaussians'; see 'collect_and_plot_timestamp_differences'.	'flat'
full_scale_s	float | None	Triangle half-base seed in seconds - the TDC full scale. The default is None, which resolves it from the artifact's own open_shutter_time_s, then from this call's.	None
pairs	Sequence[str] | None	Restrict pooling to these pixel-pair labels. The default is None (pool every pair).	None
label	str	Short label for the title and output file. The default is 'spdc'.	'spdc'
cmap_title	str | None	Optional override for the plot title. The default is None.	None
snap_to_grid	bool	Round bin_width_s to a whole number of grid cells and plot_window_s to a whole number of bins. The default is True, and you almost always want it - see the note below. Pass False to take the numbers literally.	True
subtract_background	bool	Subtract the shifted-frame accidentals before fitting. The default is False. Requires a lagged artifact, and only makes sense with background='flat' - there is no triangle left to fit. Writes a second, three-panel figure whose middle panel is the evidence that the pedestal really was accidental; read it before quoting anything.	False
peak_window_s	float	Half-range the model-free coincidence count is summed over, in seconds. subtract_background only. The default is 1e-9 (1 ns).	1e-09
wide_rebin	int	Grid cells per bin in the two full-support panels of the diagnostic figure. subtract_background only. The default is 25, which smooths the triangle without hiding structure in the residual.	25

Returns:

Name	Type	Description
counts	ndarray	Counts per bin - the rebinned histogram that was fitted and plotted, or the residual when subtract_background is set.
centers	ndarray	Bin centres in picoseconds.
fit	dict	As 'collect_and_plot_timestamp_differences'. With subtract_background it also carries model_free (the 'background_subtraction.integrate_residual' summary), subtraction (the residual, background and errors) and subtraction_png.

Raises:

Type	Description
FileNotFoundError	Raised when the counts artifact cannot be found.
ValueError	Raised when subtract_background is asked for on an artifact with no lag planes, or alongside a non-flat background model.

Notes

Re-binning already-binned data onto an incommensurate width is not a rounding detail, it is a ripple. A 71 ps bin over a 7.7 ps grid holds 9 or 10 cells depending on where it lands, so a perfectly smooth distribution comes out with an ~11% bin-to-bin sawtooth - on real data that moved the peak bin by 8%.

With snap_to_grid the rebin is not merely close to the feather path's histogram, it is identical to it, because the target bins then hold a whole odd number of cells and their edges fall between cells rather than through them. Verified bin-for-bin on a 3.22 GB run: zero difference at 9, 11 and 13 cells per bin, ~1% at 8, 10 and 12. See '_snap_binning'.

The one place the two paths still disagree is the outermost bin at each end, where the feather path cuts the data at abs(delta) <= window but draws edges half a bin beyond it, so those bins can only ever be part filled. Here they are filled properly.

Source code in src/dapkel/functions/delta_t.py

def collect_and_plot_delta_counts(
    path: str | None = None,
    *,
    counts_path: str | None = None,
    open_shutter_time_s: float | None = None,
    time_unit_s: float = timing.CODE_S,
    bin_width_s: float = timing.CODE_S,
    plot_window_s: float = 3e-9,
    fit_window_s: float | None = 1e-9,
    background: str = "flat",
    full_scale_s: float | None = None,
    pairs: Sequence[str] | None = None,
    label: str = "spdc",
    cmap_title: str | None = None,
    snap_to_grid: bool = True,
    subtract_background: bool = False,
    peak_window_s: float = 1e-9,
    wide_rebin: int = 25,
) -> tuple[np.ndarray, np.ndarray, dict]:
    """Rebin a saved counts grid, fit the peak, and plot it.

    The counterpart of 'collect_and_plot_timestamp_differences' for the
    feather-free path, and deliberately the same function underneath: only the
    histogram differs, so the two can be compared directly.

    With ``subtract_background`` the histogram that gets fitted is the *lag-0
    minus shifted-frames* residual instead of the raw one - same fit, same
    figure, on a distribution whose accidental pedestal has been measured away
    rather than modelled. Needs an artifact written with
    ``calculate_and_save_delta_counts(..., subtract_background=True)``.

    Nothing is read from disk except one ~100 MB array, and the rebinning is
    in memory, so changing ``bin_width_s`` or ``plot_window_s`` costs
    milliseconds rather than a full pass over the data - and every pixel pair
    is available at once instead of whichever columns were streamed.

    Parameters
    ----------
    path : str | None, optional
        Data folder whose ``processed/`` holds the counts artifact and whose
        ``results/`` receives the plot. May be None when ``counts_path`` is
        given.
    counts_path : str | None, optional
        Explicit counts '.npy', ignoring ``path``. The default is None.
    open_shutter_time_s : float | None, optional
        Open shutter time, i.e. one reference-clock period, in seconds. The
        fallback source for the TDC full scale, hence the triangle half-base.
        The default is None, which is fine only for ``background='flat'``.
    time_unit_s : float, optional
        Seconds per raw TDC code, for an uncalibrated grid. The default is
        'core.timing.CODE_S' (~71.6 ps), which is clock-independent.
    bin_width_s : float, optional
        Histogram bin width in seconds. The default is one TDC code. Use a
        whole multiple of the artifact's grid cell for an exact rebin.
    plot_window_s : float, optional
        Histogram/plot half-range in seconds around zero. The default is
        ``3e-9`` (3 ns). Unlike the feather path, widening this is free - it is
        a slice of what is already loaded.
    fit_window_s : float | None, optional
        Fit half-range in seconds around the tallest bin. ``'flat'`` background
        only. The default is ``1e-9`` (1 ns).
    background : str, optional
        ``'flat'`` (default), ``'triangle'`` or ``'two_gaussians'``; see
        'collect_and_plot_timestamp_differences'.
    full_scale_s : float | None, optional
        Triangle half-base seed in seconds - the TDC full scale. The default is
        None, which resolves it from the artifact's own
        ``open_shutter_time_s``, then from this call's.
    pairs : Sequence[str] | None, optional
        Restrict pooling to these pixel-pair labels. The default is None
        (pool every pair).
    label : str, optional
        Short label for the title and output file. The default is ``'spdc'``.
    cmap_title : str | None, optional
        Optional override for the plot title. The default is None.
    snap_to_grid : bool, optional
        Round ``bin_width_s`` to a whole number of grid cells and
        ``plot_window_s`` to a whole number of bins. The default is True,
        and you almost always want it - see the note below. Pass False to
        take the numbers literally.
    subtract_background : bool, optional
        Subtract the shifted-frame accidentals before fitting. The default is
        False. Requires a lagged artifact, and only makes sense with
        ``background='flat'`` - there is no triangle left to fit. Writes a
        second, three-panel figure whose middle panel is the evidence that the
        pedestal really was accidental; read it before quoting anything.
    peak_window_s : float, optional
        Half-range the model-free coincidence count is summed over, in seconds.
        ``subtract_background`` only. The default is ``1e-9`` (1 ns).
    wide_rebin : int, optional
        Grid cells per bin in the two full-support panels of the diagnostic
        figure. ``subtract_background`` only. The default is 25, which smooths
        the triangle without hiding structure in the residual.

    Returns
    -------
    counts : np.ndarray
        Counts per bin - the rebinned histogram that was fitted and plotted,
        or the residual when ``subtract_background`` is set.
    centers : np.ndarray
        Bin centres in picoseconds.
    fit : dict
        As 'collect_and_plot_timestamp_differences'. With
        ``subtract_background`` it also carries ``model_free`` (the
        'background_subtraction.integrate_residual' summary), ``subtraction``
        (the residual, background and errors) and ``subtraction_png``.

    Raises
    ------
    FileNotFoundError
        Raised when the counts artifact cannot be found.
    ValueError
        Raised when ``subtract_background`` is asked for on an artifact with no
        lag planes, or alongside a non-flat background model.

    Notes
    -----
    Re-binning already-binned data onto an incommensurate width is not a
    rounding detail, it is a **ripple**. A 71 ps bin over a 7.7 ps grid holds
    9 or 10 cells depending on where it lands, so a perfectly smooth
    distribution comes out with an ~11% bin-to-bin sawtooth - on real data
    that moved the peak bin by 8%.

    With ``snap_to_grid`` the rebin is not merely close to the feather path's
    histogram, it is **identical to it**, because the target bins then hold a
    whole odd number of cells and their edges fall between cells rather than
    through them. Verified bin-for-bin on a 3.22 GB run: zero difference at 9,
    11 and 13 cells per bin, ~1% at 8, 10 and 12. See '_snap_binning'.

    The one place the two paths still disagree is the outermost bin at each
    end, where the feather path cuts the data at ``abs(delta) <= window`` but
    draws edges half a bin beyond it, so those bins can only ever be part
    filled. Here they are filled properly.
    """
    grid_counts, grid_centers, info = load_delta_counts(
        path,
        counts_path=counts_path,
        pairs=pairs,
        pool=True,
        time_unit_s=time_unit_s,
    )
    if info["missing"]:
        print(f"  no such pixel-pair row(s): {info['missing']}")

    # The artifact records the clock it was built at, so a triangle fit can get
    # its half-base without the caller restating it.
    needs_full_scale = background in ("triangle", "two_gaussians")
    if needs_full_scale:
        full_scale_s, _ = timing.resolve_tdc_full_scale_s(
            full_scale_s=full_scale_s,
            open_shutter_time_s=(
                open_shutter_time_s
                if open_shutter_time_s is not None
                else info.get("open_shutter_time_s")
            ),
        )

    cell_ps = timing.to_ps(float(info["cell_s"]))
    if snap_to_grid:
        bin_width_s, plot_window_s = _snap_binning(
            bin_width_s, plot_window_s, timing.to_s(cell_ps)
        )
    else:
        ratio = timing.to_ps(bin_width_s) / cell_ps
        if abs(ratio - round(ratio)) > 1e-9:
            print(
                f"  WARNING: bin_width={bin_width_s * 1e12:.2f} ps is "
                f"{ratio:.2f} grid cells, so bins hold {int(ratio)} or "
                f"{int(ratio) + 1} of them and the histogram carries a "
                f"~{100 / max(int(ratio), 1):.0f}% sawtooth. Leave "
                "snap_to_grid=True unless you know you want this."
            )

    sub = None
    if subtract_background:
        _check_subtraction_model(background)
        if grid_counts.ndim != 2 or not info.get("lags"):
            raise ValueError(
                "subtract_background=True needs a counts artifact written with "
                "lag planes. Recompute stage 1 with "
                "'calculate_and_save_delta_counts(..., "
                "subtract_background=True)'."
            )
        # Rebin the native grid once, over its whole support, and let the
        # helper slice the zoom out of it: one binning for the fit and the
        # diagnostic means the panels cannot disagree about a bin edge.
        lag_counts, lag_centers = _rebin_lag_planes(
            grid_counts,
            grid_centers,
            bin_width_s=bin_width_s,
            plot_window_s=_whole_bins(
                timing.to_s(float(np.abs(grid_centers).max())), bin_width_s
            ),
        )
        sub = _subtract_lag_background(
            lag_counts,
            lag_centers,
            info["frames_used"],
            name=info["name"],
            root=info["root"],
            background=background,
            bin_width_s=bin_width_s,
            zoom_window_s=plot_window_s,
            peak_window_s=peak_window_s,
            wide_rebin=wide_rebin,
            label=label,
        )
        counts = sub["residual"]
        centers = sub["centers"]
        # The model-free count, not the bin total: the residual's bins are
        # differences of two histograms, so summing them *is* the coincidence
        # number and anything outside the peak is noise around zero.
        n_in_window = int(round(sub["model_free"]["n"]))
    else:
        counts, edges = rebin_delta_counts(
            grid_counts,
            grid_centers,
            bin_width_s=bin_width_s,
            plot_window_s=plot_window_s,
        )
        centers = (edges[:-1] + edges[1:]) / 2
        n_in_window = int(counts.sum())

        print(
            f"\n> > > {int(grid_counts.sum())} difference(s) on the "
            f"{grid_centers.size}-cell grid across {len(info['labels'])} pair "
            f"row(s); {n_in_window} inside "
            f"+/-{plot_window_s * 1e9:.2f} ns, rebinned to "
            f"{counts.size} x {bin_width_s * 1e12:.0f} ps < < <"
        )

    counts, centers, fit = _fit_and_plot_delta(
        counts,
        centers,
        name=info["name"],
        results_dir=store.results_dir(info["root"], _KIND, create=False),
        bin_width_s=bin_width_s,
        plot_window_s=plot_window_s,
        fit_window_s=fit_window_s,
        background=background,
        full_scale_s=full_scale_s if needs_full_scale else None,
        label=label,
        cmap_title=cmap_title,
        n_in_window=n_in_window,
        stem="delta_counts_sub" if subtract_background else "delta_counts",
        subtracted=subtract_background,
    )
    if sub is not None:
        fit["model_free"] = sub["model_free"]
        fit["subtraction"] = sub
        fit["subtraction_png"] = sub["png_path"]
    return counts, centers, fit

calc_diff

Per-window timestamp differences for pairs of pixels.

The low-level, daplis-style counterpart of daplis.functions.calc_diff: it takes an already-structured mapping of per-pixel, per-frame timestamps and returns the differences for the requested pixel pairs. Unpacking, structuring and saving are handled one level up in 'dapkel.functions.delta_t'.

A coincidence is a within-window event: in a frame where both pixels of a pair hold a valid TDC code, the difference of their codes is one entry of the delta-t distribution. Correlated pairs (e.g. SPDC signal/idler) pile up at delta = 0.

See docs/guide/coincidences.md.

as_pixel_list

as_pixel_list(group: Pixel | Sequence[Pixel]) -> list[Pixel]

Normalise a detector spec into a list of (row, col) pixels.

Parameters:

Name	Type	Description	Default
group	tuple[int, int] | Sequence[tuple[int, int]]	A single (row, col) pixel or a sequence of them.	required

Returns:

Type	Description
list[tuple[int, int]]	List of (row, col) pixels.

Source code in src/dapkel/functions/calc_diff.py

def as_pixel_list(group: Pixel | Sequence[Pixel]) -> list[Pixel]:
    """Normalise a detector spec into a list of ``(row, col)`` pixels.

    Parameters
    ----------
    group : tuple[int, int] | Sequence[tuple[int, int]]
        A single ``(row, col)`` pixel or a sequence of them.

    Returns
    -------
    list[tuple[int, int]]
        List of ``(row, col)`` pixels.
    """
    if (
        len(group) == 2
        and isinstance(group[0], (int, np.integer))
        and isinstance(group[1], (int, np.integer))
    ):
        return [(int(group[0]), int(group[1]))]
    return [(int(r), int(c)) for r, c in group]

calculate_differences

calculate_differences(pixel_timestamps: dict[Pixel, ndarray], pixels: Sequence[Sequence[Pixel]], delta_window_native: float | None = None, mode: str = 'all_pairs') -> dict[str, ndarray]

Compute per-window timestamp differences for pairs of pixels.

For each requested pixel pair (a, b) the difference t_b - t_a is taken in every acquisition window (frame) where both pixels hold a valid timestamp. The per-pixel arrays must be aligned by frame index (same length, NaN where the pixel had no photon that frame), which is how 'delta_t' structures them.

Parameters:

Name	Type	Description	Default
pixel_timestamps	dict[tuple[int, int], ndarray]	Mapping (row, col) -> 1D array of that pixel's per-frame timestamps (in TDC code units), with NaN where the pixel did not register a photon in that frame. All arrays share the same length (the number of frames).	required
pixels	Sequence[Sequence[tuple[int, int]]]	Two groups [group_a, group_b], each a sequence of (row, col) pixels (e.g. the signal and idler blobs read off the hitmap). A single (row, col) is accepted for a one-pixel group.	required
delta_window_native	float | None	If given, keep only differences with abs(delta) <= delta_window_native. In the same units as pixel_timestamps - picoseconds when a calibration LUT has been applied upstream, raw TDC codes when it has not - which is why this one parameter is not in seconds like the rest of the package: the unit is the caller's to know. 'delta_t' converts a delta_window_s for you. The default is None (keep all).	None
mode	str	'all_pairs' (default) forms the full cross product of group_a x group_b — use this for SPDC/HBT when it is not known which pixel of a blob catches a given photon. '1v1' zips group_a[i] with group_b[i] (requires equal-length groups).	'all_pairs'

Returns:

Type	Description
dict[str, ndarray]	Mapping "ra,ca-rb,cb" -> 1D array of timestamp differences for each pixel pair.

Raises:

Type	Description
ValueError	Raised when mode is invalid, or when mode='1v1' and the two groups differ in length.

Source code in src/dapkel/functions/calc_diff.py

def calculate_differences(
    pixel_timestamps: dict[Pixel, np.ndarray],
    pixels: Sequence[Sequence[Pixel]],
    delta_window_native: float | None = None,
    mode: str = "all_pairs",
) -> dict[str, np.ndarray]:
    """Compute per-window timestamp differences for pairs of pixels.

    For each requested pixel pair ``(a, b)`` the difference ``t_b - t_a``
    is taken in every acquisition window (frame) where *both* pixels hold a
    valid timestamp. The per-pixel arrays must be aligned by frame index
    (same length, NaN where the pixel had no photon that frame), which is
    how 'delta_t' structures them.

    Parameters
    ----------
    pixel_timestamps : dict[tuple[int, int], np.ndarray]
        Mapping ``(row, col) -> 1D array`` of that pixel's per-frame
        timestamps (in TDC code units), with NaN where the pixel did not
        register a photon in that frame. All arrays share the same length
        (the number of frames).
    pixels : Sequence[Sequence[tuple[int, int]]]
        Two groups ``[group_a, group_b]``, each a sequence of
        ``(row, col)`` pixels (e.g. the signal and idler blobs read off the
        hitmap). A single ``(row, col)`` is accepted for a one-pixel group.
    delta_window_native : float | None, optional
        If given, keep only differences with
        ``abs(delta) <= delta_window_native``. **In the same units as
        ``pixel_timestamps``** - picoseconds when a calibration LUT has been
        applied upstream, raw TDC codes when it has not - which is why this one
        parameter is not in seconds like the rest of the package: the unit is
        the caller's to know. 'delta_t' converts a ``delta_window_s`` for you.
        The default is None (keep all).
    mode : str, optional
        ``'all_pairs'`` (default) forms the full cross product of
        ``group_a`` x ``group_b`` — use this for SPDC/HBT when it is not
        known which pixel of a blob catches a given photon. ``'1v1'`` zips
        ``group_a[i]`` with ``group_b[i]`` (requires equal-length groups).

    Returns
    -------
    dict[str, np.ndarray]
        Mapping ``"ra,ca-rb,cb" -> 1D array`` of timestamp differences for
        each pixel pair.

    Raises
    ------
    ValueError
        Raised when ``mode`` is invalid, or when ``mode='1v1'`` and the two
        groups differ in length.
    """
    if mode not in ("all_pairs", "1v1"):
        raise ValueError(f"mode must be 'all_pairs' or '1v1', got {mode!r}")

    group_a = as_pixel_list(pixels[0])
    group_b = as_pixel_list(pixels[1])

    if mode == "1v1":
        if len(group_a) != len(group_b):
            raise ValueError(
                "mode='1v1' needs equal-length groups "
                f"(got {len(group_a)} and {len(group_b)})."
            )
        pair_iter = list(zip(group_a, group_b, strict=True))
    else:
        pair_iter = [(a, b) for a in group_a for b in group_b]

    deltas_all: dict[str, np.ndarray] = {}
    for a, b in pair_iter:
        ta = pixel_timestamps[a]
        tb = pixel_timestamps[b]
        both = ~np.isnan(ta) & ~np.isnan(tb)
        d = tb[both] - ta[both]
        if delta_window_native is not None:
            d = d[np.abs(d) <= delta_window_native]
        deltas_all[f"{a[0]},{a[1]}-{b[0]},{b[1]}"] = d

    return deltas_all