diff --git a/.gitignore b/.gitignore index a5e39f1..652c889 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,5 @@ simspad *.vscode __pycache__ spice/*_tran.txt +sipm_fast.json +fast_kernel.npy diff --git a/README.md b/README.md index 6bbc2c7..57e76d1 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,32 @@ The input files can be created using examples in the examples directory (see the Python helpers in `examples/python/simspad.py`). The simulation streams in bounded memory, so arbitrarily long traces can be run. +#### Output pulse shaping + +By default SimSPAD returns the bare avalanche charge per step. `-S/--shape` +applies a length-preserving output filter: + +| mode | output | +|---|---| +| `none` | bare avalanche charge train (default) | +| `gaussian` | Gaussian FIR, FWHM = `tauFwhm` (an idealised amplifier) | +| `fast` | bipolar AC-coupled fast-output terminal, two-pole (`tauLoad`, `tauRecovery`) | +| `bench` | `fast` then `gaussian` — what a scope capture of the real device looks like | +| `kernel` | convolve with a **tabulated** impulse response (the actual device pulse) named by `kernelFile` in the params | + +`--shape kernel` is the faithful path for an arbitrary SiPM: supply the device's +real fast-output pulse (measured, or from a SPICE deck) as a 1-D `.npy` of +current per unit avalanche charge, point `kernelFile`/`kernelDt` at it, and the +simulator convolves it directly rather than assuming a functional form. + +The easiest way to build that kernel is **`examples/python/pulse_to_kernel.py`**, +which turns a single-photoelectron fast-output pulse — from a bench scope capture +**or** a SPICE run — into the kernel `.npy` plus an augmented params file (see +`examples/python/README.md`). For the shipped MicroFJ-30020 there is also a +one-shot C++ tool, `spice_kernel`, that runs the bundled ngspice equivalent +circuit and additionally fits the cheaper two-pole `fast` mode as a cross-check +(see `spice/README.md`). + ### Web Application Once SimSPAD server is running, you are able to send a POST request to `http://localhost:33232/simspad`. diff --git a/examples/python/README.md b/examples/python/README.md index c7f2b3d..0b708be 100644 --- a/examples/python/README.md +++ b/examples/python/README.md @@ -1,3 +1,118 @@ # Python Examples -An importable .py exists with simspad.py, and example usage is shown in example.py \ No newline at end of file +Helpers for SimSPAD's self-describing IO and for building output-shaping +kernels. + +| File | What it is | +|---|---| +| `simspad.py` | Importable helpers: write/read device parameters (JSON) and waveforms (`.npy`), and a small web client. | +| `example.py` | End-to-end usage: configure a device, simulate, save the params + waveforms. | +| `pulse_to_kernel.py` | Turn a **measured or simulated single-PE fast-output pulse** into a `--shape kernel` kernel (see below). | + +## Bring your own SiPM: the `kernel` shape mode + +SimSPAD's core Monte-Carlo models the **avalanche statistics** of *your* device +from its parameters (`numMicrocell`, `vBr`, `tauRecovery`, `pdeMax`, `vChr`, +`cCell`). On top of that, `--shape` filters the output to mimic a real readout. +The `kernel` mode is the faithful, device-specific one: it convolves the +avalanche charge train with a **tabulated fast-output impulse response** — the +actual shape of your terminal — instead of a built-in approximation. + +`pulse_to_kernel.py` builds that kernel from one single-photoelectron (1 p.e.) +fast-output pulse. The pulse can come from **a bench capture** or **a SPICE +run** — the tool only needs a `time, signal` table. + +### What the kernel is + +The kernel `g(t)` is the fast-output **current per unit avalanche charge** +(units 1/s): + +``` +g(t) = i_f(t) / Q_avalanche +``` + +`i_f` is the fast-output current (`= v / R_load` if you captured a voltage), and +`Q_avalanche` is the charge one fired microcell deposits. Convolving `g` with +SimSPAD's charge-per-step train then gives physical fast-output current. One +captured 1-PE pulse *is* one avalanche, so dividing by `Q_avalanche` is the +entire normalisation — the measured amplitude already carries the AC-coupling +fraction and the rail RC of the real terminal. + +You supply `Q_avalanche` one of three ways (most → least self-consistent): + +* `--device params.json` → `Q = (vBias − vBr) · cCell`. **Recommended:** this is + the *same* charge SimSPAD assigns a fired cell, so the convolution units match + the simulation exactly. +* `--gain G` → `Q = G · e` (datasheet gain × electron charge). +* `--charge C` → `Q = C` directly (Coulombs). + +> Note: normalising by the device's own `(vBias−vBr)·cCell` rather than a +> circuit's avalanche charge matters when the two disagree (e.g. a SPICE deck +> whose junction capacitance differs from the model `cCell`). The `--device` +> path keeps the kernel consistent with the simulation it feeds. + +### Example A — bench capture + +Export an oscilloscope trace of a single dark-count fast-output pulse as CSV +(`time, voltage`, here in nanoseconds across a 50 Ω load): + +```bash +python pulse_to_kernel.py \ + --pulse single_pe.csv --time-unit ns \ + --signal-units voltage --load 50 \ + --device j30020.json --params-out j30020_kernel.json +simspad -p j30020_kernel.json -i light.npy -o resp.npy --shape kernel +``` + +This writes `fast_kernel.npy` and a copy of `j30020.json` with `kernelFile` / +`kernelDt` added. `--dt` defaults to the device's `dt`; the onset is +auto-detected (override with `--t0`). + +### Example B — SPICE deck + +Run a fast-output equivalent circuit (see `spice/README.md` for the shipped +J-Series deck and how to adapt it to your device), then feed its `wrdata` dump +straight in — `pulse_to_kernel.py` reads whitespace tables and skips headers, so +no conversion step is needed: + +```bash +ngspice -b spice/jseries_fastout.cir # -> spice/fastout_tran.txt (t v(f) t v(a)) +python pulse_to_kernel.py \ + --pulse spice/fastout_tran.txt --time-col 0 --signal-col 1 \ + --signal-units voltage --load 50 \ + --device j30020.json --params-out j30020_kernel.json +simspad -p j30020_kernel.json -i light.npy -o resp.npy --shape kernel +``` + +### Options + +`--help` lists everything. The common ones: + +| flag | meaning | +|---|---| +| `--pulse PATH` | trace file: CSV / whitespace / ngspice `wrdata` | +| `--time-col N`, `--signal-col N` | 0-based column indices (default 0, 1) | +| `--time-unit {s,ms,us,ns,ps}` | unit of the time column (default `s`) | +| `--signal-units {current,voltage}` | `voltage` needs `--load OHMS` | +| `--device / --gain / --charge` | the `Q_avalanche` normalisation (pick one) | +| `--dt SECONDS` | kernel sample step (default: device `dt`) | +| `--t0 SECONDS` | avalanche onset (default: auto-detect) | +| `--kernel-out`, `--params-out` | output paths | + +### Verify the maths + +```bash +python pulse_to_kernel.py --selftest +``` + +Round-trips an analytic bipolar pulse through the tool and checks the kernel +reproduces it (and integrates to ~0). This mirrors the C++ `make test` shaping +checks. + +## The C++ `spice_kernel` tool + +`src/spice_kernel.cpp` (`make spice_kernel`) is a one-shot C++ pipeline +specialised to the shipped J-Series deck: it runs ngspice, emits the kernel, and +*also* fits the cheaper two-pole `--shape fast` approximation as a cross-check. +`pulse_to_kernel.py` is the general, source-agnostic path — prefer it for an +arbitrary device or a bench measurement. diff --git a/examples/python/pulse_to_kernel.py b/examples/python/pulse_to_kernel.py new file mode 100755 index 0000000..83a82f1 --- /dev/null +++ b/examples/python/pulse_to_kernel.py @@ -0,0 +1,377 @@ +#!/usr/bin/env python3 +"""Turn a measured (or simulated) single-photon fast-output pulse into a +SimSPAD ``--shape kernel`` kernel. + +SimSPAD's ``kernel`` shape mode convolves the avalanche charge train with a +*tabulated* fast-output impulse response, so the simulated output carries the +real shape of *your* device's terminal instead of a two-pole approximation. +This tool builds that kernel from a single-photoelectron (1 p.e.) pulse, taken +from any source: + + * a **bench capture** -- an oscilloscope trace of one dark-count / single-PE + fast-output pulse, exported as CSV; + * a **SPICE run** -- ``ngspice -b your_device.cir`` writing a ``wrdata`` table + of the fast-output node (see ``spice/README.md``); + * any 2-column ``time, signal`` table. + +What the kernel is +------------------ +The kernel ``g(t)`` is the fast-output **current per unit avalanche charge** +(units 1/s):: + + g(t) = i_f(t) / Q_avalanche + +where ``i_f`` is the fast-output current (``= v / R_load`` if you captured a +voltage across a load) and ``Q_avalanche`` is the charge one fired microcell +deposits. Convolving ``g`` with SimSPAD's avalanche charge-per-step train +(Coulombs/step) then yields physical fast-output current (A). One captured 1-PE +pulse *is* one avalanche, so dividing by ``Q_avalanche`` is the whole +normalisation -- the measured amplitude already carries the AC-coupling fraction +and the rail RC of the real terminal. + +``Q_avalanche`` can be given three ways (most to least self-consistent): + + * ``--device params.json`` -> ``Q = (vBias - vBr) * cCell`` (the *same* charge + SimSPAD assigns a fired cell, so the convolution units match exactly); + * ``--gain G`` -> ``Q = G * e`` (datasheet gain x electron charge); + * ``--charge C`` -> ``Q = C`` (Coulombs, if you know it directly). + +Output +------ +Writes ``fast_kernel.npy`` (1-D float64, the kernel sampled at the simulation +``dt``) and, when ``--device`` is given, an augmented parameter file with +``kernelFile``/``kernelDt`` added, ready to run:: + + python pulse_to_kernel.py --pulse single_pe.csv --signal-units voltage \ + --load 50 --device j30020.json --params-out j30020_kernel.json + simspad -p j30020_kernel.json -i light.npy -o resp.npy --shape kernel + +Run ``python pulse_to_kernel.py --selftest`` to verify the maths against an +analytic bipolar pulse, or ``--help`` for the full option list. +""" +import argparse +import sys + +import numpy as np + +try: + from simspad import read_params # augment an existing device parameter file +except ImportError: # allow running from outside examples/python/ + read_params = None + +ELECTRON_CHARGE = 1.602176634e-19 # C + +# Time-unit suffixes accepted for --time-unit and recognised in column headers. +TIME_UNITS = {"s": 1.0, "ms": 1e-3, "us": 1e-6, "ns": 1e-9, "ps": 1e-12} + + +def load_trace(path, time_col=0, signal_col=1): + """Read two numeric columns from a CSV / whitespace / ngspice-wrdata table. + + Lines that are blank or whose first column is non-numeric (headers, SPICE + comment lines beginning ``*`` or ``#``) are skipped, so the same reader + handles a scope CSV with a header row and a raw ngspice ``wrdata`` dump. + Comma-separated files are detected automatically; everything else is split + on whitespace. + """ + times, signals = [], [] + with open(path) as f: + for line in f: + line = line.strip() + if not line or line[0] in "#*;": + continue + fields = line.split(",") if "," in line else line.split() + if max(time_col, signal_col) >= len(fields): + continue + try: + t = float(fields[time_col]) + s = float(fields[signal_col]) + except ValueError: + continue # header / non-numeric row + times.append(t) + signals.append(s) + if len(times) < 4: + raise ValueError( + f"{path}: found only {len(times)} numeric rows in columns " + f"({time_col}, {signal_col}) -- check --time-col/--signal-col" + ) + return np.asarray(times, dtype=float), np.asarray(signals, dtype=float) + + +def detect_onset(t, i_f, frac=0.01): + """Estimate the avalanche onset time: the foot of the rising edge. + + Finds the positive peak and walks back to the last sample below + ``frac * peak``. Override with ``--t0`` if your trace has a pre-trigger + offset or a non-standard baseline. + """ + peak_idx = int(np.argmax(i_f)) + peak = i_f[peak_idx] + if peak <= 0: + return t[0] + thresh = frac * peak + onset_idx = 0 + for k in range(peak_idx, -1, -1): + if i_f[k] < thresh: + onset_idx = k + break + return t[onset_idx] + + +def pulse_to_kernel(t, signal, q_avalanche, dt, t0=None, load=None, onset_frac=0.01): + """Build a SimSPAD kernel from a single-PE fast-output pulse. + + Parameters + ---------- + t, signal : 1-D arrays + Time [s] and the captured fast-output signal (current, or voltage if + ``load`` is given). + q_avalanche : float + Charge one fired microcell deposits [C]; the normalisation. + dt : float + Simulation time step [s] to sample the kernel at (== ``kernelDt``). + t0 : float or None + Avalanche onset time [s]; auto-detected when None. + load : float or None + Load resistance [ohm]; if given, ``signal`` is a voltage and the + current is ``signal / load``. + onset_frac : float + Fraction-of-peak threshold for onset auto-detection. + + Returns + ------- + kernel : 1-D float64 array + ``g(t) = i_f(t) / q_avalanche`` resampled onto the ``dt`` grid from the + onset (units 1/s). + info : dict + Diagnostics: ``t0``, ``n_taps``, ``net_charge`` (per unit avalanche + charge; ~0 for an AC-coupled terminal), ``fwhm`` of the positive lobe. + """ + t = np.asarray(t, dtype=float) + i_f = np.asarray(signal, dtype=float) + if load is not None: + if load <= 0: + raise ValueError("--load must be > 0 ohm") + i_f = i_f / load + order = np.argsort(t) # tolerate unsorted / non-monotone exports + t, i_f = t[order], i_f[order] + + if t0 is None: + t0 = detect_onset(t, i_f, onset_frac) + if t0 >= t[-1]: + raise ValueError(f"onset t0={t0:g}s is at/after the trace end {t[-1]:g}s") + if q_avalanche <= 0: + raise ValueError("avalanche charge must be > 0 C") + if dt <= 0: + raise ValueError("dt must be > 0 s") + + # Resample i_f onto a uniform grid at the simulation dt, measured from the + # onset; samples past the trace end hold the last value (the tail). + span = t[-1] - t0 + n = int(np.floor(span / dt)) + 1 + grid = np.arange(n) * dt + i_grid = np.interp(grid + t0, t, i_f, left=0.0, right=i_f[-1]) + kernel = i_grid / q_avalanche + + net = float(np.sum(kernel) * dt) # net charge per unit avalanche charge + info = { + "t0": float(t0), + "n_taps": int(n), + "net_charge": net, + "fwhm": _positive_fwhm(grid, kernel), + "peak": float(np.max(kernel)), + } + return kernel.astype("= pk / 2)[0] + if above.size < 2: + return 0.0 + return float(t[above[-1]] - t[above[0]]) + + +def resolve_avalanche_charge(args, device_params): + """Pick Q_avalanche [C] from --charge / --gain / --device (in that order).""" + if args.charge is not None: + return args.charge, "explicit --charge" + if args.gain is not None: + return args.gain * ELECTRON_CHARGE, f"--gain {args.gain:g} x e" + if device_params is not None: + v_over = device_params["vBias"] - device_params["vBr"] + q = v_over * device_params["cCell"] + return q, "(vBias - vBr) * cCell from --device" + raise SystemExit( + "error: need an avalanche charge -- pass one of --device, --gain or --charge" + ) + + +def _read_device(path): + """Load a device parameter JSON as a plain dict (no SimSPAD import needed).""" + import json + + with open(path) as f: + return json.load(f) + + +def main(argv=None): + p = argparse.ArgumentParser( + description="Build a SimSPAD --shape kernel from a single-PE fast-output pulse.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + epilog="See examples/python/README.md and spice/README.md for worked examples.", + ) + p.add_argument("--pulse", help="trace file: CSV / whitespace / ngspice wrdata") + p.add_argument("--time-col", type=int, default=0, help="0-based time column") + p.add_argument("--signal-col", type=int, default=1, help="0-based signal column") + p.add_argument("--time-unit", choices=sorted(TIME_UNITS), default="s", + help="unit of the time column") + p.add_argument("--signal-units", choices=["current", "voltage"], default="current", + help="captured signal: current [A], or voltage [V] (needs --load)") + p.add_argument("--load", type=float, default=None, + help="load resistance [ohm], required for --signal-units voltage") + + norm = p.add_argument_group("normalisation (pick one)") + norm.add_argument("--device", help="device params JSON; derives Q = (vBias-vBr)*cCell and dt") + norm.add_argument("--gain", type=float, default=None, help="device gain G (Q = G*e)") + norm.add_argument("--charge", type=float, default=None, help="avalanche charge [C] directly") + + p.add_argument("--dt", type=float, default=None, + help="simulation dt [s] to sample the kernel at (default: device dt)") + p.add_argument("--t0", type=float, default=None, + help="avalanche onset time [s] (default: auto-detect)") + p.add_argument("--onset-frac", type=float, default=0.01, + help="fraction-of-peak threshold for onset auto-detection") + p.add_argument("--kernel-out", default="fast_kernel.npy", help="kernel .npy output path") + p.add_argument("--params-out", default=None, + help="write an augmented device params JSON here (requires --device)") + p.add_argument("--selftest", action="store_true", + help="verify against an analytic bipolar pulse and exit") + args = p.parse_args(argv) + + if args.selftest: + return _selftest() + + if not args.pulse: + p.error("--pulse is required (or use --selftest)") + if args.signal_units == "voltage" and args.load is None: + p.error("--signal-units voltage requires --load OHMS") + + device_params = _read_device(args.device) if args.device else None + q_av, q_src = resolve_avalanche_charge(args, device_params) + + dt = args.dt + if dt is None: + if device_params is None: + p.error("--dt is required when --device is not given") + dt = device_params["dt"] + + t_raw, sig = load_trace(args.pulse, args.time_col, args.signal_col) + t = t_raw * TIME_UNITS[args.time_unit] + load = args.load if args.signal_units == "voltage" else None + + kernel, info = pulse_to_kernel( + t, sig, q_av, dt, t0=args.t0, load=load, onset_frac=args.onset_frac + ) + np.save(args.kernel_out, kernel) + if not args.kernel_out.endswith(".npy"): + # np.save appends .npy; keep our reported name honest. + args.kernel_out += ".npy" + + print("=== pulse -> --shape kernel ===") + print(f" source pulse : {args.pulse} ({len(t)} samples)") + print(f" Q_avalanche : {q_av:.4e} C [{q_src}]") + print(f" onset t0 : {info['t0'] * 1e9:.3f} ns" + + (" (auto)" if args.t0 is None else " (given)")) + print(f" kernel : {info['n_taps']} taps @ dt = {dt * 1e12:.1f} ps") + print(f" positive FWHM: {info['fwhm'] * 1e9:.3f} ns") + # Net charge tells you which terminal you built: a fast output is AC-coupled + # (bipolar, integrates to ~0); a standard/slow output is DC-coupled (unipolar, + # carries the whole avalanche charge, so ~1). A standard output landing well + # off 1 usually means the device cCell disagrees with the circuit's cell cap. + net = info["net_charge"] + kind = "bipolar / AC-coupled fast output, expect ~0" if abs(net) < 0.1 \ + else "unipolar / DC-coupled standard (slow) output, expect ~1" + print(f" net charge/avalanche = {net:+.3e} [{kind}]") + print(f" wrote {args.kernel_out}") + + if args.params_out: + if args.device is None: + p.error("--params-out requires --device (the params to augment)") + _write_augmented_params(args.device, args.params_out, args.kernel_out, dt) + print(f" wrote {args.params_out}") + print(f" run: simspad -p {args.params_out} -i light.npy -o resp.npy --shape kernel") + elif args.device: + print(f' add to {args.device}: "kernelFile": "{args.kernel_out}", "kernelDt": {dt:g}') + return 0 + + +def _write_augmented_params(device_in, params_out, kernel_file, kernel_dt): + """Copy a device params JSON, adding kernelFile/kernelDt.""" + if read_params is not None: + sipm = read_params(device_in) + sipm.kernelFile = kernel_file + sipm.kernelDt = kernel_dt + sipm.write_params(params_out) + else: # standalone fallback: edit the JSON dict directly + import json + + with open(device_in) as f: + d = json.load(f) + d["kernelFile"] = kernel_file + d["kernelDt"] = kernel_dt + with open(params_out, "w") as f: + json.dump(d, f, indent=2) + + +def _selftest(): + """Round-trip an analytic bipolar pulse through pulse_to_kernel(). + + Samples h(t) = A*(e^{-t/tauL}/tauL - e^{-t/tauR}/tauR) at fine resolution, + treats it as a measured current trace, and checks the produced kernel + reproduces h(t) (after the dt resampling) and integrates to ~0. + """ + dt = 1e-11 + tauL, tauR = 2.0e-9, 24.5e-9 + A = tauR / (tauR - tauL) + q_av = 1.0 # normalise to 1 so the kernel equals h directly + + # Fine "measured" trace with a 5 ns pre-trigger baseline. + t = np.arange(0, 200e-9, dt / 3) + t0 = 5e-9 + u = t - t0 + h = np.where(u >= 0, A * (np.exp(-u / tauL) / tauL - np.exp(-u / tauR) / tauR), 0.0) + + kernel, info = pulse_to_kernel(t, h, q_av, dt) # onset auto-detected + + # Expected kernel: h at the SAME absolute sample times the tool used + # (grid measured from the *detected* onset), so the comparison is not thrown + # off by the sub-sample onset-detection offset at the sharp peak. + grid = np.arange(info["n_taps"]) * dt + want_u = (grid + info["t0"]) - t0 # time relative to the true onset + want = np.where( + want_u >= 0, A * (np.exp(-want_u / tauL) / tauL - np.exp(-want_u / tauR) / tauR), 0.0 + ) + rel = np.max(np.abs(kernel - want)) / np.max(np.abs(want)) + + onset_err = abs(info["t0"] - t0) + ok_shape = rel < 1e-2 + ok_onset = onset_err < 0.2e-9 + ok_net = abs(info["net_charge"]) < 5e-3 + print("=== pulse_to_kernel self-test ===") + print(f" onset detected : {info['t0'] * 1e9:.3f} ns (want 5.000) " + f"{'PASS' if ok_onset else 'FAIL'}") + print(f" shape rel err : {rel:.2e} (want < 1e-2) " + f"{'PASS' if ok_shape else 'FAIL'}") + print(f" net charge : {info['net_charge']:+.2e} (want ~0) " + f"{'PASS' if ok_net else 'FAIL'}") + ok = ok_shape and ok_onset and ok_net + print(" RESULT:", "PASS" if ok else "FAIL") + return 0 if ok else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/examples/python/simspad.py b/examples/python/simspad.py index 427069d..6538439 100644 --- a/examples/python/simspad.py +++ b/examples/python/simspad.py @@ -19,8 +19,14 @@ # Optional parameters with their defaults (absent keys keep the default on # both sides, so files without them remain valid). tauLoad is the fast-output # rail RC time constant in seconds, used by the `fast`/`bench` shape modes. +# kernelFile/kernelDt point the `kernel` shape mode at a tabulated fast-output +# impulse response (a 1-D .npy of current per unit avalanche charge, units 1/s, +# sampled every kernelDt seconds); both default to None and are emitted only +# when set, so plain device files are unchanged. OPTIONAL_PARAM_DEFAULTS = { "tauLoad": 2.0e-9, + "kernelFile": None, + "kernelDt": None, } @@ -40,7 +46,12 @@ def __init__(self, *args, **optional): # -- parameters (JSON) -------------------------------------------------- def params_dict(self): d = {k: getattr(self, k) for k in PARAM_KEYS} - d.update({k: getattr(self, k) for k in OPTIONAL_PARAM_DEFAULTS}) + # Optional keys are emitted only when set (None -> omitted), so device + # files without a kernel stay byte-for-byte as before. + for k in OPTIONAL_PARAM_DEFAULTS: + v = getattr(self, k) + if v is not None: + d[k] = v d["numMicrocell"] = int(d["numMicrocell"]) return d diff --git a/makefile b/makefile index c755604..5ad2888 100644 --- a/makefile +++ b/makefile @@ -10,7 +10,8 @@ TARGET_SERVER := server TARGET_TEST := test INCLUDE := -Ilib/ SRC_ALL := $(wildcard src/*.cpp) -SRC := $(filter-out src/server.cpp, $(SRC_ALL)) +SRC := $(filter-out src/server.cpp src/spice_kernel.cpp, $(SRC_ALL)) +TARGET_SPICE := spice_kernel VERSION := $(shell git describe --tags --always --dirty) CXXFLAGS += -DVERSION=\"$(VERSION)\" @@ -30,7 +31,7 @@ $(APP_DIR)/$(TARGET): $(OBJECTS) -include $(DEPENDENCIES) -.PHONY: all test build clean debug release info +.PHONY: all test build clean debug release info spice_kernel build: @mkdir -p $(APP_DIR) @@ -68,3 +69,7 @@ server: ./src/server.cpp ./src/sipm.cpp ./src/utilities.cpp ./src/pages.cpp ./sr simspad: ./src/main.cpp ./src/sipm.cpp ./src/utilities.cpp $(CXX) $(CXXFLAGS) -o $(APP_DIR)/$(TARGET) ./src/main.cpp ./src/sipm.cpp ./src/utilities.cpp + +spice_kernel: ./src/spice_kernel.cpp ./src/sipm.cpp ./src/utilities.cpp + @mkdir -p $(APP_DIR) + $(CXX) $(CXXFLAGS) $(INCLUDE) -o $(APP_DIR)/$(TARGET_SPICE) ./src/spice_kernel.cpp ./src/sipm.cpp ./src/utilities.cpp $(LDFLAGS) diff --git a/spice/README.md b/spice/README.md index 73c739b..a37597a 100644 --- a/spice/README.md +++ b/spice/README.md @@ -2,10 +2,12 @@ `jseries_fastout.cir` is an [ngspice](http://ngspice.sourceforge.net/) deck for the Onsemi **J-Series** SiPM fast-output terminal (values for **MicroFJ-30020**, -overvoltage 2.5 V). It exists to give the `--shape fast` / `--shape bench` pulse -kernels a physical provenance: rather than hand-picking the fast-rail time -constant, you can derive it from a circuit whose component values trace to the -datasheet. +overvoltage 2.5 V). It exists to give the output-shaping modes a physical +provenance: rather than hand-picking a pulse shape, you derive it from a circuit +whose component values trace to the datasheet. The `spice_kernel` tool turns +this deck into a tabulated impulse response the simulator convolves directly +(`--shape kernel`), and also fits the cheaper two-pole `--shape fast` +approximation as a cross-check. ## What it models @@ -24,8 +26,23 @@ The avalanche is triggered by a switch at t = 5 ns; the deck runs a 120 ns transient and writes node voltages `v(f)` (fast) and `v(a)` (anode) to `spice/fastout_tran.txt` (currents are `v/RL`). -## Relation to `--shape fast` +## Two ways to use this pulse +**1. `--shape kernel` — the actual curve (recommended, faithful).** +SimSPAD convolves the avalanche charge train with the *tabulated* fast-output +impulse response taken straight from this deck: + +``` +g(t) = i_f(t) / Q_av, Q_av = (Cd+Cq)·OV (units 1/s) +``` + +i.e. the fast-terminal current per unit avalanche charge, resampled onto the +simulation `dt` and stored in `fast_kernel.npy`. No functional form is assumed, +so whatever shape your circuit (or a measured pulse) produces is what the +simulator uses. This is the right path for *any* SiPM: change the deck, rerun, +get that device's real fast output. + +**2. `--shape fast` — a two-pole reduction (cheaper, approximate).** The bipolar shaper added in #28 implements, per unit avalanche charge, ``` @@ -33,17 +50,52 @@ h(t) = A·( e^{-t/tauLoad}/tauLoad − e^{-t/tauRecovery}/tauRecovery ), A = tauRecovery/(tauRecovery − tauLoad), ∫ h dt = 0 ``` -The two time constants map onto the circuit as: +with the two time constants mapped onto the circuit as: * `tauLoad = RLFAST · N · Cf` — the fast-rail RC (≈ 2.5 ns for these values); * `tauRecovery` — the cell's *loaded* recovery constant (array loading slows the bare `Rq·(Cd+Cf)`; the fit gives ≈ 24.5 ns here). -So `--shape fast` is a two-pole reduction of this circuit. Issue #29 tracks the -pipeline that runs this deck and emits the fitted `tauLoad` as a ready-to-use -`params.json` value. +So `--shape fast` is a two-pole reduction of `--shape kernel`. It's a recursive +one-pole pair (O(1)/sample, streaming), whereas `kernel` is an FIR convolution +with the full pulse (O(taps)/sample, buffers the trace) — pick `fast` when you +need speed and the two-pole shape is good enough, `kernel` when you want the +real curve. -## Running it +## Deriving the kernel (`spice_kernel`) + +`src/spice_kernel.cpp` is the spice → kernel pipeline. It runs this deck, +extracts the single-photon fast pulse, normalises it by the avalanche charge, +resamples it onto the simulation `dt`, and writes `fast_kernel.npy` plus a +ready-to-run `params.json` that points at it (`kernelFile`/`kernelDt`). As a +cross-check it also fits the two-pole `tauLoad`/`tauRecovery` and reports their +RMS error against the SPICE pulse (and writes `tauLoad`, so `--shape fast` works +off the same file). Build and run from the repository root (needs `ngspice` on +PATH): + +``` +make spice_kernel +./build/apps/spice_kernel +# -> sipm_fast.json + fast_kernel.npy +simspad -p sipm_fast.json -i light.npy -o resp.npy --shape kernel # the actual curve +simspad -p sipm_fast.json -i light.npy -o resp.npy --shape fast # two-pole fit +``` + +`tauLoad` is read from the deck's `.param` values (`R_Lfast·N·C_f`), so it can +never drift from the circuit. Sample run: + +``` +tauLoad = R_Lfast*N*C_f = 2.500 ns +tauRecovery = loaded fit = 24.539 ns +peak-normalised shape RMS error (kernel vs spice, <40 ns): 0.088 [PASS] +kernel: 11500 taps @ dt = 10.0 ps, net charge/avalanche = +6.400e-04 (want ~0) +``` + +The kernel's net charge is near zero (the fast terminal is AC-coupled); the +small residual is the recovery tail still decaying when the 120 ns transient +ends (≈4.7·tauRecovery), and shrinks if you lengthen the `tran` window. + +## Running the deck directly ``` ngspice -b spice/jseries_fastout.cir # writes spice/fastout_tran.txt @@ -53,3 +105,95 @@ Tested with ngspice-45. The transient dump is a build artifact (gitignored). Single-photon checks against the datasheet (run from the deck): FWHM ≈ 1.9 ns into 50 Ω, positive-lobe charge ≈ `Cf·OV`, near-zero net charge (AC-coupled), anode charge ≈ `G·e`. + +## Included decks + +All are Onsemi J-Series (MICROJ-SERIES/D rev.7), with inputs from the datasheet +in `docs/datasheets/microj-series-datasheet.pdf` (Tables 1–3, gain quoted at +OV = +2.5 V). Each deck writes both the fast (`v(f)`) and standard/anode (`v(a)`) +outputs, so one deck serves both terminals. + +| Deck | Device | Cells | Recharge τ | Fast cap | Provenance | +|---|---|---|---|---|---| +| `jseries_fastout.cir` | MicroFJ-30020 (3 mm, 20 µm) | 14410 | — | 50 pF | Bench-calibrated; the validation anchor | +| `microfj_30035_fastout.cir` | MicroFJ-30035 (3 mm, 35 µm) | 5676 | 45 ns | 40 pF | **Datasheet** | +| `microfj_40035_fastout.cir` | MicroFJ-40035 (4 mm, 35 µm) | 9260 | 48 ns | 70 pF | **Datasheet** | +| `microfj_60035_fastout.cir` | MicroFJ-60035 (6 mm, 35 µm) | 22292 | 50 ns | 160 pF | **Datasheet** | +| `template_fastout.cir` | **fast-output template** | — | — | — | Fill in datasheet inputs; internals derived | +| `template_stdout.cir` | **standard-output template** | — | — | — | For devices with NO fast pin (e.g. Hamamatsu MPPCs) | + +The datasheet's own consistency is a nice check on the model: with Cμ = G·e/OV = +185.9 fF, the array capacitance N·Cμ reproduces the datasheet anode capacitances +to 96–100 % (60035: 4143 vs 4140 pF). + +### Fidelity: fast output vs standard output + +The **standard (anode) output** is governed by the datasheet recharge τ and is +well-matched. The **fast output** is harder: the simple `N·Cf` fast-rail lumping +*over-broadens* it, and the error grows with array size — the decks give FWHM +~1.8 / 2.7 / 5.4 ns for the 30035 / 40035 / 60035 vs the datasheet's 1.5 / 1.7 / +3.0 ns. So for an **exact fast output**, don't rely on the lumped circuit: +**digitise the datasheet's measured pulse** (Figs. 5 and 6 give the fast and +standard shapes for all three) and feed that curve straight to +`pulse_to_kernel.py`. That is the whole advantage of the tabulated-kernel +approach — a measured curve always beats a model. + +## Your own device: run SPICE → kernel (the general path) + +For any device, the recommended route is the source-agnostic +`examples/python/pulse_to_kernel.py`, which turns *any* 1-PE fast-output trace +(a SPICE `wrdata` dump **or** a bench oscilloscope capture) into a kernel. The +SPICE half is just "run a deck that writes the fast node": + +1. **Make the deck.** Copy `template_fastout.cir` and fill in its *datasheet + inputs* block — gain, overvoltage, breakdown, cell count `NCELL`, recharge τ, + total fast capacitance, and the load `RLFAST`. Everything internal (`Cd`, + `Cq`, `Rq`, `Cf`) is **derived** from those via `.param` expressions, so you + never hand-pick a component value. Point the `wrdata v(f)` line at your + own output path. Fire one cell with the `VFIRE` switch; the exact trigger + time doesn't matter (the tool auto-detects the onset). Set the `tran` window + to ≳ 8 × recharge τ so the AC-coupled pulse fully returns to zero. + +2. **Run it** and feed the dump in: + + ```bash + ngspice -b spice/mydevice.cir # -> writes your wrdata file + python examples/python/pulse_to_kernel.py \ + --pulse spice/mydevice_tran.txt --time-col 0 --signal-col 1 \ + --signal-units voltage --load 50 \ + --device mydevice.json --params-out mydevice_kernel.json + simspad -p mydevice_kernel.json -i light.npy -o resp.npy --shape kernel + ``` + +`mydevice.json` is your device's parameter file (the avalanche-model parameters +SimSPAD simulates — SPICE only supplies the *terminal shape*). See +`examples/python/README.md` for the full option list, the bench-capture path, +and the normalisation choices. + +### Standard ("slow") output instead of the fast output + +The standard/anode output is what most SiPMs are read on (the fast terminal is a +specialised extra pin). Every fast-output deck already writes the anode node +`v(a)` as its second `wrdata` column, so you build a **slow kernel** from the +*same* run — just point the tool at the anode column and load: + +```bash +python examples/python/pulse_to_kernel.py \ + --pulse spice/mydevice_tran.txt --time-col 0 --signal-col 3 \ + --signal-units voltage --load 10 \ + --device mydevice.json --params-out mydevice_slow.json +``` + +(`wrdata v(f) v(a)` lays the columns out as `t v(f) t v(a)`, so the anode is +column 3; `--load 10` is `RLANODE`.) For a device with **no fast pin** at all, +start from `template_stdout.cir`, which models the standard output directly +(single `wrdata v(a)` column → `--signal-col 1`). + +Unlike the AC-coupled fast output (bipolar, net charge ≈ 0), the standard output +is **DC-coupled** — unipolar, carrying the full avalanche charge, so its kernel's +net charge ≈ 1. `pulse_to_kernel.py` reports which regime it detected. + +> Security note: an ngspice deck can run arbitrary shell commands (via +> `.control` `shell`/`system`), so only run decks you trust — and never wire +> "upload a deck" into the web server. The bench-capture path takes only numbers +> and has no such surface. diff --git a/spice/jseries_fastout.cir b/spice/jseries_fastout.cir index 89dcd81..dc214be 100644 --- a/spice/jseries_fastout.cir +++ b/spice/jseries_fastout.cir @@ -59,7 +59,11 @@ VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) .options reltol=1e-4 abstol=1e-15 chgtol=1e-18 .control -tran 2p 120n 0 10p +* 350 ns: the loaded fast tail (tau ~ 24.5 ns) must decay to ~1e-6 of the +* lobe, or a kernel tabulated from this run carries a net-charge residual +* (truncation at 120 ns left +1.3% of the lobe unbalanced, which a pure- +* integrator observer downstream turns into volts of drift). +tran 2p 350n 0 10p * load voltages; the python driver converts to currents via RLFAST/RLANODE wrdata spice/fastout_tran.txt v(f) v(a) quit diff --git a/spice/microfj_30035_fastout.cir b/spice/microfj_30035_fastout.cir new file mode 100644 index 0000000..556fa80 --- /dev/null +++ b/spice/microfj_30035_fastout.cir @@ -0,0 +1,74 @@ +* MicroFJ-30035 (Onsemi J-Series, 3 mm, 35 um) fast-output equivalent circuit. +* Same topology as template_fastout.cir / jseries_fastout.cir. +* +* ALL INPUTS BELOW ARE FROM THE J-SERIES DATASHEET (MICROJ-SERIES/D rev.7, +* docs/datasheets/microj-series-datasheet.pdf), Tables 1-3, at OV = +2.5 V: +* NCELL = 5676 (Table 2, 3 mm 35 um) +* GAIN = 2.9e6 (Table 3, anode-cathode, +2.5 V) +* VBR = 24.5 (Table 1, 24.2-24.7 V) +* TAU_RECHARGE = 45 ns (Table 3, microcell recharge time constant) +* CF_TOTAL = 40 pF (Table 3, fast-output capacitance) +* anode capacitance 1070 pF (Table 3) -- reproduced as N*(Cd+Cq): the deck's +* Cmicro = G*e/OV = 185.9 fF gives N*Cmicro = 1055 pF (99% of datasheet), +* an independent check that gain/cells/capacitance are mutually consistent. +* Datasheet sanity targets: fast-output FWHM 1.5 ns, anode rise time 90 ps. +* +* MODEL FIDELITY: the standard (anode) output is set by the datasheet recharge +* tau and is well-matched. The lumped N*Cf fast-rail over-broadens the FAST +* output -- this deck gives ~1.8 ns FWHM vs the datasheet 1.5 ns (the error +* grows with array size). For an exact fast output, digitise the datasheet's +* MEASURED pulse (Fig. 5/6) and feed it to pulse_to_kernel.py instead. +* +* ===================== DATASHEET INPUTS ===================== +.param GAIN = 2.9e6 ; Table 3, anode-cathode gain at +2.5 V +.param OV = 2.5 ; overvoltage the gain is quoted at +.param VBR = 24.5 ; Table 1 (24.2-24.7 V) +.param NCELL = 5676 ; Table 2 (3 mm, 35 um) +.param TAU_RECHARGE = 45n ; Table 3, microcell recharge time constant +.param CF_TOTAL = 40p ; Table 3, fast-output capacitance (N*Cf) +.param RLFAST = 50 +.param RLANODE = 10 +.param FQ = 0.086 ; Cq/(Cd+Cq), from the 30020 deck's 5.5/64 split +.param RD = 900 ; avalanche series R; sets the sub-ns rise edge +* ============================================================= +* +.param E = 1.602176634e-19 +.param CMICRO = {GAIN*E/OV} +.param CQ = {FQ*CMICRO} +.param CD = {CMICRO-CQ} +.param RQ = {TAU_RECHARGE/CMICRO} +.param CF = {CF_TOTAL/NCELL} +.param NPAS = {NCELL-1} + +VBIAS k 0 DC {VBR+OV} + +CD1 k nf {CD} +S1 k nav fire 0 SWMOD +RD1 nav nbr {RD} +VBR1 nbr nf DC {VBR} +RQ1 nf a {RQ} +CQ1 nf a {CQ} +CF1 nf f {CF} + +CDP k np {NPAS*CD} +RQP np a {RQ/NPAS} +CQP np a {NPAS*CQ} +CFP np f {NPAS*CF} + +RLF f 0 {RLFAST} +RLA a 0 {RLANODE} + +VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) +.model SWMOD SW(RON=1u ROFF=1G VT=0.5 VH=0.1) + +.options reltol=1e-4 abstol=1e-15 chgtol=1e-18 + +.control +* >= ~8 * TAU_RECHARGE so the AC-coupled fast pulse fully returns to zero. +tran 2p 600n 0 10p +* fast (v(f)) and standard/anode (v(a)) outputs; pulse_to_kernel converts to +* current via the load R (RLFAST for v(f), RLANODE for v(a)). +wrdata spice/microfj_30035_tran.txt v(f) v(a) +quit +.endc +.end diff --git a/spice/microfj_40035_fastout.cir b/spice/microfj_40035_fastout.cir new file mode 100644 index 0000000..1582f81 --- /dev/null +++ b/spice/microfj_40035_fastout.cir @@ -0,0 +1,72 @@ +* MicroFJ-40035 (Onsemi J-Series, 4 mm, 35 um) fast-output equivalent circuit. +* Same topology as template_fastout.cir / jseries_fastout.cir. +* +* ALL INPUTS BELOW ARE FROM THE J-SERIES DATASHEET (MICROJ-SERIES/D rev.7, +* docs/datasheets/microj-series-datasheet.pdf), Tables 1-3, at OV = +2.5 V: +* NCELL = 9260 (Table 2, 4 mm 35 um) +* GAIN = 2.9e6 (Table 3, anode-cathode, +2.5 V) +* VBR = 24.5 (Table 1, 24.2-24.7 V) +* TAU_RECHARGE = 48 ns (Table 3, microcell recharge time constant) +* CF_TOTAL = 70 pF (Table 3, fast-output capacitance) +* anode capacitance 1800 pF (Table 3) -- N*Cmicro = 1721 pF (96%), consistent. +* Datasheet sanity targets: fast-output FWHM 1.7 ns, anode rise time 90 ps. +* +* MODEL FIDELITY: the standard (anode) output is set by the datasheet recharge +* tau and is well-matched. The lumped N*Cf fast-rail over-broadens the FAST +* output -- this deck gives ~2.7 ns FWHM vs the datasheet 1.7 ns (the error +* grows with array size). For an exact fast output, digitise the datasheet's +* MEASURED pulse (Fig. 5/6) and feed it to pulse_to_kernel.py instead. +* +* ===================== DATASHEET INPUTS ===================== +.param GAIN = 2.9e6 ; Table 3, anode-cathode gain at +2.5 V +.param OV = 2.5 ; overvoltage the gain is quoted at +.param VBR = 24.5 ; Table 1 (24.2-24.7 V) +.param NCELL = 9260 ; Table 2 (4 mm, 35 um) +.param TAU_RECHARGE = 48n ; Table 3, microcell recharge time constant +.param CF_TOTAL = 70p ; Table 3, fast-output capacitance (N*Cf) +.param RLFAST = 50 +.param RLANODE = 10 +.param FQ = 0.086 ; Cq/(Cd+Cq), from the 30020 deck's 5.5/64 split +.param RD = 900 ; avalanche series R; sets the sub-ns rise edge +* ============================================================= +* +.param E = 1.602176634e-19 +.param CMICRO = {GAIN*E/OV} +.param CQ = {FQ*CMICRO} +.param CD = {CMICRO-CQ} +.param RQ = {TAU_RECHARGE/CMICRO} +.param CF = {CF_TOTAL/NCELL} +.param NPAS = {NCELL-1} + +VBIAS k 0 DC {VBR+OV} + +CD1 k nf {CD} +S1 k nav fire 0 SWMOD +RD1 nav nbr {RD} +VBR1 nbr nf DC {VBR} +RQ1 nf a {RQ} +CQ1 nf a {CQ} +CF1 nf f {CF} + +CDP k np {NPAS*CD} +RQP np a {RQ/NPAS} +CQP np a {NPAS*CQ} +CFP np f {NPAS*CF} + +RLF f 0 {RLFAST} +RLA a 0 {RLANODE} + +VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) +.model SWMOD SW(RON=1u ROFF=1G VT=0.5 VH=0.1) + +.options reltol=1e-4 abstol=1e-15 chgtol=1e-18 + +.control +* >= ~8 * TAU_RECHARGE so the AC-coupled fast pulse fully returns to zero. +tran 2p 600n 0 10p +* fast (v(f)) and standard/anode (v(a)) outputs; pulse_to_kernel converts to +* current via the load R (RLFAST for v(f), RLANODE for v(a)). +wrdata spice/microfj_40035_tran.txt v(f) v(a) +quit +.endc +.end diff --git a/spice/microfj_60035_fastout.cir b/spice/microfj_60035_fastout.cir new file mode 100644 index 0000000..2152e71 --- /dev/null +++ b/spice/microfj_60035_fastout.cir @@ -0,0 +1,76 @@ +* MicroFJ-60035 (Onsemi J-Series, 6 mm, 35 um) fast-output equivalent circuit. +* Same topology as template_fastout.cir / jseries_fastout.cir. +* +* ALL INPUTS BELOW ARE FROM THE J-SERIES DATASHEET (MICROJ-SERIES/D rev.7, +* docs/datasheets/microj-series-datasheet.pdf), Tables 1-3, at OV = +2.5 V: +* NCELL = 22292 (Table 2, 6 mm 35 um) +* GAIN = 2.9e6 (Table 3, anode-cathode, +2.5 V) +* VBR = 24.5 (Table 1, 24.2-24.7 V) +* TAU_RECHARGE = 50 ns (Table 3, microcell recharge time constant) +* CF_TOTAL = 160 pF (Table 3, fast-output capacitance) +* anode capacitance 4140 pF (Table 3) -- N*Cmicro = 4143 pF (100%), consistent. +* Datasheet sanity targets: fast-output FWHM 3.0 ns, anode rise time 180 ps. +* +* The 6 mm part has ~4x the cells of the 3 mm, hence ~4x the fast-rail RC, so its +* fast output is correspondingly broader (datasheet FWHM 3.0 ns vs 1.5 ns) -- a +* real, modelled size effect. +* +* MODEL FIDELITY: the standard (anode) output is set by the datasheet recharge +* tau and is well-matched. The lumped N*Cf fast-rail over-broadens the FAST +* output most at this size -- this deck gives ~5.4 ns FWHM vs the datasheet +* 3.0 ns. For an exact fast output, digitise the datasheet's MEASURED pulse +* (Fig. 5/6) and feed it to pulse_to_kernel.py instead. +* +* ===================== DATASHEET INPUTS ===================== +.param GAIN = 2.9e6 ; Table 3, anode-cathode gain at +2.5 V +.param OV = 2.5 ; overvoltage the gain is quoted at +.param VBR = 24.5 ; Table 1 (24.2-24.7 V) +.param NCELL = 22292 ; Table 2 (6 mm, 35 um) +.param TAU_RECHARGE = 50n ; Table 3, microcell recharge time constant +.param CF_TOTAL = 160p ; Table 3, fast-output capacitance (N*Cf) +.param RLFAST = 50 +.param RLANODE = 10 +.param FQ = 0.086 ; Cq/(Cd+Cq), from the 30020 deck's 5.5/64 split +.param RD = 900 ; avalanche series R; sets the sub-ns rise edge +* ============================================================= +* +.param E = 1.602176634e-19 +.param CMICRO = {GAIN*E/OV} +.param CQ = {FQ*CMICRO} +.param CD = {CMICRO-CQ} +.param RQ = {TAU_RECHARGE/CMICRO} +.param CF = {CF_TOTAL/NCELL} +.param NPAS = {NCELL-1} + +VBIAS k 0 DC {VBR+OV} + +CD1 k nf {CD} +S1 k nav fire 0 SWMOD +RD1 nav nbr {RD} +VBR1 nbr nf DC {VBR} +RQ1 nf a {RQ} +CQ1 nf a {CQ} +CF1 nf f {CF} + +CDP k np {NPAS*CD} +RQP np a {RQ/NPAS} +CQP np a {NPAS*CQ} +CFP np f {NPAS*CF} + +RLF f 0 {RLFAST} +RLA a 0 {RLANODE} + +VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) +.model SWMOD SW(RON=1u ROFF=1G VT=0.5 VH=0.1) + +.options reltol=1e-4 abstol=1e-15 chgtol=1e-18 + +.control +* >= ~8 * TAU_RECHARGE; the larger fast-rail RC needs a longer window. +tran 2p 800n 0 10p +* fast (v(f)) and standard/anode (v(a)) outputs; pulse_to_kernel converts to +* current via the load R (RLFAST for v(f), RLANODE for v(a)). +wrdata spice/microfj_60035_tran.txt v(f) v(a) +quit +.endc +.end diff --git a/spice/template_fastout.cir b/spice/template_fastout.cir new file mode 100644 index 0000000..5c5ce01 --- /dev/null +++ b/spice/template_fastout.cir @@ -0,0 +1,83 @@ +* SiPM fast-output equivalent circuit -- PARAMETERISED TEMPLATE. +* +* Copy this file, edit the "datasheet inputs" block for your device, and run it +* to get a 1-PE fast-output pulse, then feed that into +* examples/python/pulse_to_kernel.py to build a SimSPAD `--shape kernel` kernel. +* Everything below the inputs block is DERIVED, so you only ever supply numbers +* you can read off a datasheet -- nothing is hand-picked. +* +* Topology (one fired microcell + lumped passive array), per the Onsemi +* C-/J-Series fast-output structure, datasheet Fig. 9 (MICROJ-SERIES/D rev.5): +* each microcell is an APD junction capacitance Cd from the cathode rail to an +* internal node, a quench resistor Rq (with parasitic Cq) from that node to the +* anode rail, and a fast-output coupling cap Cf from the node to the fast rail. +* The fired cell self-terminates at exactly the overvoltage OV, so its avalanche +* charge is (Cd+Cq)*OV = G*e, then it recharges through Rq. +* +* NOTE: this models devices that HAVE a dedicated fast output (Onsemi C-/J-Series +* style). Devices without one (e.g. most Hamamatsu MPPCs) have no `f` node -- +* read the standard/anode output v(a) instead. +* +* ===================== DATASHEET INPUTS (edit these) ===================== +.param GAIN = 1.0e6 ; microcell gain at OV (datasheet) +.param OV = 2.5 ; overvoltage [V] the gain is quoted at (datasheet) +.param VBR = 24.5 ; breakdown voltage [V] (datasheet) +.param NCELL = 14410 ; number of microcells (datasheet) +.param TAU_RECHARGE = 15n ; microcell recharge time constant [s] (datasheet) +.param CF_TOTAL = 50p ; total fast-output capacitance N*Cf [F] (datasheet, or +* ~ a few % of the terminal capacitance Ct) +.param RLFAST = 50 ; fast-output load [ohm] (your readout, often 50) +.param RLANODE = 10 ; anode/standard-output load [ohm] (your readout) +* --- second-order knobs (defaults are fine for most devices) --- +.param FQ = 0.086 ; quench parasitic Cq as a fraction of (Cd+Cq) +.param RD = 900 ; avalanche series R; sets the sub-ns rise edge +* ========================================================================== +* +* ---------------------------- DERIVED ----------------------------------- +.param E = 1.602176634e-19 ; electron charge [C] +.param CMICRO = {GAIN*E/OV} ; total microcell cap Cd+Cq = G*e/OV +.param CQ = {FQ*CMICRO} ; quench parasitic capacitance +.param CD = {CMICRO-CQ} ; junction capacitance +.param RQ = {TAU_RECHARGE/CMICRO} ; quench resistor (recharge tau / cap) +.param CF = {CF_TOTAL/NCELL} ; per-cell fast coupling cap +.param NPAS = {NCELL-1} ; the remaining passive cells, lumped + +* cathode bias (AC ground) +VBIAS k 0 DC {VBR+OV} + +* ---- fired microcell ---------------------------------------------------- +CD1 k nf {CD} +S1 k nav fire 0 SWMOD +RD1 nav nbr {RD} +VBR1 nbr nf DC {VBR} +RQ1 nf a {RQ} +CQ1 nf a {CQ} +CF1 nf f {CF} + +* ---- remaining NPAS passive microcells, lumped --------------------------- +CDP k np {NPAS*CD} +RQP np a {RQ/NPAS} +CQP np a {NPAS*CQ} +CFP np f {NPAS*CF} + +* ---- output loads --------------------------------------------------------- +RLF f 0 {RLFAST} +RLA a 0 {RLANODE} + +* avalanche trigger: fire one cell at t = 5 ns (the exact time is uncritical; +* the discharge completes in ~Rd*(Cd+Cq), and pulse_to_kernel auto-detects the +* onset). +VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) +.model SWMOD SW(RON=1u ROFF=1G VT=0.5 VH=0.1) + +.options reltol=1e-4 abstol=1e-15 chgtol=1e-18 + +.control +* Run for >= ~8 * TAU_RECHARGE so the AC-coupled fast pulse fully returns to +* zero (a short window leaves a small net-charge residual in the kernel). +tran 2p 120n 0 10p +* load-node voltages; pulse_to_kernel converts to currents via the load R. +wrdata spice/template_tran.txt v(f) v(a) +quit +.endc +.end diff --git a/spice/template_stdout.cir b/spice/template_stdout.cir new file mode 100644 index 0000000..f43ad0f --- /dev/null +++ b/spice/template_stdout.cir @@ -0,0 +1,75 @@ +* SiPM STANDARD ("slow") output equivalent circuit -- PARAMETERISED TEMPLATE. +* +* The standard/anode output is the readout most SiPMs use (the fast output is a +* specialised third terminal that only some devices -- Onsemi C-/J-Series -- +* expose). Use THIS template for the standard output, and in particular for +* devices with NO fast pin at all (e.g. most Hamamatsu MPPCs). +* +* Copy this file, edit the "datasheet inputs" block, run it, and feed the dump +* into examples/python/pulse_to_kernel.py to build a SimSPAD `--shape kernel` +* kernel for the standard output. Everything below the inputs is DERIVED. +* +* Topology: one fired microcell + lumped passive array (the standard SiPM +* equivalent circuit). Each microcell is an APD junction capacitance Cd from the +* cathode rail to an internal node, and a quench resistor Rq (with parasitic Cq) +* from that node to the anode rail, which is the readout. There is NO fast-output +* coupling cap. The fired cell self-terminates at the overvoltage OV (avalanche +* charge (Cd+Cq)*OV = G*e) and recharges through Rq, so the standard output is a +* fast leading edge followed by a slow recharge tail -- UNIPOLAR (DC-coupled, it +* carries the full avalanche charge), unlike the bipolar fast output. +* +* ===================== DATASHEET INPUTS (edit these) ===================== +.param GAIN = 1.0e6 ; microcell gain at OV (datasheet) +.param OV = 2.5 ; overvoltage [V] the gain is quoted at (datasheet) +.param VBR = 24.5 ; breakdown voltage [V] (datasheet) +.param NCELL = 14410 ; number of microcells (datasheet) +.param TAU_RECHARGE = 15n ; microcell recharge time constant [s] (datasheet) +.param RLOAD = 50 ; standard-output readout load [ohm] (your readout) +* --- second-order knobs (defaults are fine for most devices) --- +.param FQ = 0.086 ; quench parasitic Cq as a fraction of (Cd+Cq) +.param RD = 900 ; avalanche series R; sets the sub-ns rise edge +* ========================================================================== +* +* ---------------------------- DERIVED ----------------------------------- +.param E = 1.602176634e-19 ; electron charge [C] +.param CMICRO = {GAIN*E/OV} ; total microcell cap Cd+Cq = G*e/OV +.param CQ = {FQ*CMICRO} ; quench parasitic capacitance +.param CD = {CMICRO-CQ} ; junction capacitance +.param RQ = {TAU_RECHARGE/CMICRO} ; quench resistor (recharge tau / cap) +.param NPAS = {NCELL-1} ; the remaining passive cells, lumped + +* cathode bias (AC ground) +VBIAS k 0 DC {VBR+OV} + +* ---- fired microcell (no fast-output cap) -------------------------------- +CD1 k nf {CD} +S1 k nav fire 0 SWMOD +RD1 nav nbr {RD} +VBR1 nbr nf DC {VBR} +RQ1 nf a {RQ} +CQ1 nf a {CQ} + +* ---- remaining NPAS passive microcells, lumped --------------------------- +CDP k np {NPAS*CD} +RQP np a {RQ/NPAS} +CQP np a {NPAS*CQ} + +* ---- standard-output readout load (on the anode rail) -------------------- +RLA a 0 {RLOAD} + +* avalanche trigger: fire one cell at t = 5 ns (exact time uncritical; the tool +* auto-detects the onset). +VFIRE fire 0 PULSE(0 1 5n 20p 20p 2n 1u) +.model SWMOD SW(RON=1u ROFF=1G VT=0.5 VH=0.1) + +.options reltol=1e-4 abstol=1e-15 chgtol=1e-18 + +.control +* Run for >= ~8 * TAU_RECHARGE so the slow recharge tail fully decays. +tran 2p 200n 0 10p +* anode (standard-output) node voltage; pulse_to_kernel converts to current via +* the load R. The standard output is column 2 here (single wrdata column). +wrdata spice/template_stdout_tran.txt v(a) +quit +.endc +.end diff --git a/src/main.cpp b/src/main.cpp index b6eecf0..26f4376 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -82,19 +82,26 @@ void print_info(chrono::duration elapsed, SiPM &sipm, size_t inputSize, // bench - fast then gaussian: the physical terminal followed by a // ~Gaussian amplifier, i.e. what an oscilloscope capture of the // real device looks like +// kernel - tabulated fast-output kernel: convolve with the actual measured +// or SPICE-derived impulse response named by "kernelFile" in the +// params JSON (the real device pulse, not a two-pole fit) // The fast shaper is a pair of recursive one-pole filters, so it streams -// chunk-by-chunk with its state carried across chunks. The Gaussian is a FIR -// convolution; rather than handle inter-chunk overlap, modes that need it -// buffer the full response vector and shape it in one pass (8 bytes per -// sample of extra memory -- the streaming bound is given up for those modes). +// chunk-by-chunk with its state carried across chunks. The Gaussian and the +// tabulated kernel are FIR convolutions; rather than handle inter-chunk +// overlap, modes that need them buffer the full response vector and shape it in +// one pass (8 bytes per sample of extra memory -- the streaming bound is given +// up for those modes). void simulate(string params_file, string fname_in, string fname_out, bool silence, string shapeMode) { const size_t chunk = 1u << 16; // 65536 samples per block const bool wantFast = (shapeMode == "fast") || (shapeMode == "bench"); const bool wantGaussian = (shapeMode == "gaussian") || (shapeMode == "bench"); + const bool wantKernel = (shapeMode == "kernel"); SiPM sipm = load_params_json(params_file); + if (wantKernel && sipm.fastKernel.empty()) + throw runtime_error("--shape kernel requires a \"kernelFile\" in the params JSON"); NpyReader reader(fname_in); size_t N = reader.count(); NpyWriter writer(fname_out, N); @@ -121,8 +128,8 @@ void simulate(string params_file, string fname_in, string fname_out, bool silenc reader.rewind(); auto start = chrono::steady_clock::now(); double outSum = 0.0; - vector full; // buffered response, only used when wantGaussian - if (wantGaussian) + vector full; // buffered response, used by the FIR modes (gaussian/kernel) + if (wantGaussian || wantKernel) { full.reserve(N); } @@ -149,7 +156,7 @@ void simulate(string params_file, string fname_in, string fname_out, bool silenc sipm.shape_fast_chunk(outbuf.data(), shapebuf.data(), got); res = shapebuf.data(); } - if (wantGaussian) + if (wantGaussian || wantKernel) { full.insert(full.end(), res, res + got); } @@ -174,6 +181,12 @@ void simulate(string params_file, string fname_in, string fname_out, bool silenc vector shaped = conv1d(full, get_gaussian(sipm.dt, sipm.tauFwhm)); writer.write(shaped.data(), shaped.size()); } + else if (wantKernel) + { + // Causal FIR with the tabulated fast-output impulse response. + vector shaped = sipm.shape_kernel(full); + writer.write(shaped.data(), shaped.size()); + } writer.close(); auto end = chrono::steady_clock::now(); @@ -207,8 +220,10 @@ static void show_usage(string name) << "\t-o,--output OUTPUT\tResponse output path (.npy) [required]\n" << "\t-S,--shape MODE\t\tOutput pulse shaping: none (default),\n" << "\t\t\t\tgaussian (FWHM = tauFwhm), fast (bipolar AC-coupled\n" - << "\t\t\t\tfast-output terminal, tau = tauLoad), or bench\n" - << "\t\t\t\t(fast then gaussian, like a real scope capture)" + << "\t\t\t\tfast-output terminal, tau = tauLoad), bench\n" + << "\t\t\t\t(fast then gaussian, like a real scope capture), or\n" + << "\t\t\t\tkernel (convolve with the tabulated impulse response\n" + << "\t\t\t\tin the params' kernelFile, e.g. a SPICE-derived pulse)" << endl; } @@ -298,10 +313,11 @@ int main(int argc, char *argv[]) if (!a) return EXIT_FAILURE; shapeMode = a; - if (shapeMode != "none" && shapeMode != "gaussian" && shapeMode != "fast" && shapeMode != "bench") + if (shapeMode != "none" && shapeMode != "gaussian" && shapeMode != "fast" && + shapeMode != "bench" && shapeMode != "kernel") { cerr << "error: unknown shape mode '" << shapeMode - << "' (expected none, gaussian, fast or bench)." << endl; + << "' (expected none, gaussian, fast, bench or kernel)." << endl; return EXIT_FAILURE; } } diff --git a/src/sipm.cpp b/src/sipm.cpp index 779a87d..9564de8 100644 --- a/src/sipm.cpp +++ b/src/sipm.cpp @@ -276,6 +276,73 @@ vector SiPM::shape_fast(vector inputVec) return out; } +// Linear-resample the tabulated kernel from its native kernelDt onto a grid of +// spacing targetDt. Returns fastKernel unchanged when the spacings already +// match (the spice_kernel pipeline writes the kernel at the simulation dt, so +// this is normally a no-op); the interpolation path exists so a hand-supplied +// kernel sampled at any spacing still works. +vector SiPM::resample_kernel(double targetDt) const +{ + if (fastKernel.empty() || kernelDt <= 0.0 || targetDt <= 0.0) + { + return {}; + } + if (fabs(kernelDt - targetDt) <= 1e-18) + { + return fastKernel; + } + const size_t K = fastKernel.size(); + const double span = (double)(K - 1) * kernelDt; // total supported time [s] + const size_t M = (size_t)floor(span / targetDt) + 1; + vector out(M, 0.0); + for (size_t mi = 0; mi < M; mi++) + { + const double p = ((double)mi * targetDt) / kernelDt; // position in source samples + const size_t k = (size_t)p; + if (k >= K - 1) + { + out[mi] = fastKernel[K - 1]; + } + else + { + const double frac = p - (double)k; + out[mi] = fastKernel[k] * (1.0 - frac) + fastKernel[k + 1] * frac; + } + } + return out; +} + +// Tabulated-kernel (FIR) shaping: out[i] = sum_k g[k] * charge[i-k], a causal +// convolution of the avalanche charge-per-step train with the resampled +// fast-output impulse response g (units 1/s). g already carries the AC-coupling +// charge fraction C_f/(C_d+C_q) and the rail RC of the real terminal, so the +// output is the physical fast-output current per step and, like the analytic +// `fast` shaper, integrates to ~zero over the (zero-net-charge) pulse. The +// convolution is causal (no centring): a detection at step i can only affect +// outputs at i, i+1, ... Returns all-zero if no kernel is loaded. +vector SiPM::shape_kernel(const vector &charge) +{ + const vector g = resample_kernel(dt); + const size_t N = charge.size(); + vector out(N, 0.0); + if (g.empty()) + { + return out; + } + const size_t K = g.size(); + for (size_t i = 0; i < N; i++) + { + double acc = 0.0; + const size_t kmax = (i + 1 < K) ? (i + 1) : K; // clamp: charge[i-k] needs i-k >= 0 + for (size_t k = 0; k < kmax; k++) + { + acc += g[k] * charge[i - k]; + } + out[i] = acc; + } + return out; +} + // Seed Random Engines // TODO improve this code - appears to give the same result for all runs within the same second void SiPM::seed_engines() diff --git a/src/sipm.hpp b/src/sipm.hpp index a0580de..a827fb8 100644 --- a/src/sipm.hpp +++ b/src/sipm.hpp @@ -63,6 +63,18 @@ class SiPM // files keep working. double tauLoad = 2.0e-9; + // Tabulated fast-output impulse response (the `kernel` shape mode). This is + // the device's actual single-photon fast-output current per unit avalanche + // charge (units 1/s), sampled uniformly at `kernelDt` seconds. Convolving + // it causally with the avalanche charge-per-step train reproduces the real + // terminal shape directly, with no two-pole (`fast`) approximation -- so + // any device whose fast pulse can be measured or SPICE-simulated can be + // dropped in. Populated from the JSON "kernelFile"/"kernelDt"; empty + // otherwise, in which case the `fast`/`bench`/`gaussian` modes are unaffected. + std::vector fastKernel; + double kernelDt = 0.0; // sample spacing of fastKernel [s]; 0 => unset + std::string kernelFile; // source .npy path, round-tripped through the JSON + SiPM(unsigned long numMicrocell_in, double vBias_in, double vBr_in, double tauRecovery_in, double tauFwhm_in, double digitalThreshold_in, double ccell_in, double Vchr_in, double PDE_max_in); SiPM(unsigned long numMicrocell_in, double vBias_in, double vBr_in, double tauRecovery_in, double digitalThreshold_in, double ccell_in, double Vchr_in, double PDE_max_in); @@ -108,6 +120,13 @@ class SiPM std::vector shape_fast(std::vector inputVec); + // Tabulated-kernel shaping (the `kernel` mode): causal FIR convolution of + // the avalanche charge-per-step train with fastKernel, resampled from its + // native kernelDt onto the simulation dt. Returns the physical fast-output + // current per step; length-preserving. Returns all-zero if no kernel is + // loaded. O(N * kernelTaps), so unlike `fast` it buffers the whole trace. + std::vector shape_kernel(const std::vector &charge); + private: std::vector microcellTimes; double simClock = 0.0; // running simulation time, carried across chunks @@ -124,6 +143,10 @@ class SiPM void seed_engines(void); + // Linear-resample fastKernel from kernelDt onto a grid of spacing targetDt + // (a copy when the two already match, as the spice_kernel pipeline writes). + std::vector resample_kernel(double targetDt) const; + double unif_rand_double(double a, double b); unsigned long unif_rand_int(unsigned long a, unsigned long b); diff --git a/src/spice_kernel.cpp b/src/spice_kernel.cpp new file mode 100644 index 0000000..f40450c --- /dev/null +++ b/src/spice_kernel.cpp @@ -0,0 +1,284 @@ +// spice_kernel: derive the --shape kernel impulse response from the ngspice +// equivalent circuit (spice/jseries_fastout.cir). +// +// This is the spice -> kernel pipeline. It runs the J-Series fast-output +// equivalent circuit, extracts the single-photon fast-output pulse, and writes +// it out as a *tabulated impulse response* the simulator convolves directly: +// +// g(t) = i_f(t) / Q_av, Q_av = (C_d + C_q) * OV (units 1/s) +// +// i.e. the fast-terminal current per unit avalanche charge, resampled onto the +// simulation dt and saved to fast_kernel.npy. The emitted parameter file points +// at it ("kernelFile"/"kernelDt"), so +// simspad -p sipm_fast.json -i light.npy -o resp.npy --shape kernel +// shapes the trace with the *actual* circuit pulse -- no two-pole fit. Drop in +// any device's deck (or a measured pulse) and you get its real fast output. +// +// As a cross-check it also fits the two-pole approximation the analytic +// shape_fast() uses, +// h(t) = A * ( e^{-t/tauLoad}/tauLoad - e^{-t/tauRecovery}/tauRecovery ), +// A = tauRecovery / (tauRecovery - tauLoad), +// and reports its RMS error against the SPICE pulse (tauLoad is also written to +// the params file, so the cheaper --shape fast still works off the same file). +// +// tauLoad = R_Lfast * N * C_f (fast-rail RC, read from the deck so it +// cannot drift from the circuit) +// tauRecovery = loaded recovery constant, fitted on the fast tail +// +// Prereq: ngspice on PATH (tested with ngspice-45). Run from the repository +// root so the deck's relative wrdata path resolves: +// make spice_kernel && ./build/apps/spice_kernel + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "sipm.hpp" +#include "utilities.hpp" + +using namespace std; + +namespace +{ +const double T0 = 5e-9; // avalanche trigger time set in the deck (VFIRE ... 5n) + +// Parse a SPICE numeric literal with an optional engineering suffix. +double spice_value(string tok) +{ + for (auto &c : tok) + c = (char)tolower((unsigned char)c); + // "meg" must be tested before "m"; suffixes are case-folded already. + static const vector> suffix = { + {"meg", 1e6}, {"f", 1e-15}, {"p", 1e-12}, {"n", 1e-9}, {"u", 1e-6}, + {"m", 1e-3}, {"k", 1e3}, {"g", 1e9}, {"t", 1e12}}; + for (const auto &s : suffix) + { + if (tok.size() > s.first.size() && + tok.compare(tok.size() - s.first.size(), s.first.size(), s.first) == 0) + { + return stod(tok.substr(0, tok.size() - s.first.size())) * s.second; + } + } + return stod(tok); +} + +// Pull the `.param NAME=VALUE` lines we need straight from the deck. +map read_deck_params(const string &path) +{ + ifstream f(path); + if (!f) + throw runtime_error("cannot open deck: " + path); + map p; + regex re(R"(^\.param\s+(\w+)\s*=\s*\{?([^}\s]+)\}?)", regex::icase); + string line; + while (getline(f, line)) + { + smatch m; + if (regex_search(line, m, re)) + { + string name = m[1].str(); + for (auto &c : name) + c = (char)toupper((unsigned char)c); + try + { + p[name] = spice_value(m[2].str()); // skip braced expressions (e.g. NPAS={NCELL-1}) + } + catch (const invalid_argument &) + { + } + } + } + return p; +} + +double fwhm(const vector &t, const vector &y) +{ + double pk = *max_element(y.begin(), y.end()); + size_t lo = 0, hi = 0; + for (size_t i = 0; i < y.size(); i++) + if (y[i] >= pk / 2) + { + if (!lo) + lo = i; + hi = i; + } + return t[hi] - t[lo]; +} +} // namespace + +int main(int argc, char **argv) +{ + const string deck = argc > 1 ? argv[1] : "spice/jseries_fastout.cir"; + const string tran = "spice/fastout_tran.txt"; + const string out_params = argc > 2 ? argv[2] : "sipm_fast.json"; + + // 1. Run the circuit. + const string cmd = "ngspice -b " + deck + " > /dev/null 2>&1"; + if (system(cmd.c_str()) != 0) + { + cerr << "ngspice failed (is it on PATH? run from the repo root). cmd: " << cmd << "\n"; + return 2; + } + + // 2. Read circuit constants from the deck (so this can't drift from it). + map cir = read_deck_params(deck); + const double CD = cir.at("CD"), CQ = cir.at("CQ"), CF = cir.at("CF"); + const double RLFAST = cir.at("RLFAST"), OV = cir.at("OV"); + const unsigned long NCELL = (unsigned long)llround(cir.at("NCELL")); + + // 3. Read the transient (wrdata columns: t v(f) t v(a)); fast current = v(f)/R. + ifstream tf(tran); + if (!tf) + { + cerr << "cannot open transient output: " << tran << "\n"; + return 2; + } + vector t, i_f; + { + double a, b, c, d; + while (tf >> a >> b >> c >> d) + { + t.push_back(a); + i_f.push_back(b / RLFAST); + } + } + if (t.size() < 16) + { + cerr << "transient too short (" << t.size() << " rows)\n"; + return 2; + } + + // 4. Derive the kernel constants. + const double tau_load = RLFAST * (double)NCELL * CF; // fast-rail RC, from the circuit + const double kappa = CF / (CD + CQ); // fast/avalanche charge fraction + // Loaded recovery constant: log-linear least squares on the negative tail, + // well clear of the fast edge. + double sx = 0, sy = 0, sxx = 0, sxy = 0; + long n = 0; + for (size_t k = 0; k < t.size(); k++) + { + if (t[k] > T0 + 15e-9 && t[k] < T0 + 75e-9 && i_f[k] < 0) + { + double x = t[k], y = log(-i_f[k]); + sx += x; + sy += y; + sxx += x * x; + sxy += x * y; + n++; + } + } + const double slope = (n * sxy - sx * sy) / (n * sxx - sx * sx); + const double tau_rec = -1.0 / slope; + + // 5. Validate the analytic kernel against the SPICE pulse (peak-normalised + // shape, so we test the shape the device makes, not the gain). + const double A = tau_rec / (tau_rec - tau_load); + vector tt, hk, sp; + for (size_t k = 0; k < t.size(); k++) + { + double u = t[k] - T0; + if (u < 0) + continue; + tt.push_back(u); + hk.push_back(A * (exp(-u / tau_load) / tau_load - exp(-u / tau_rec) / tau_rec)); + sp.push_back(i_f[k]); + } + const double hk_pk = *max_element(hk.begin(), hk.end()); + const double sp_pk = *max_element(sp.begin(), sp.end()); + double se = 0; + long m = 0; + for (size_t k = 0; k < tt.size(); k++) + if (tt[k] < 40e-9) + { + double e = hk[k] / hk_pk - sp[k] / sp_pk; + se += e * e; + m++; + } + const double rms = sqrt(se / m); + const double fwhm_spice = fwhm(t, i_f), fwhm_kernel = fwhm(tt, hk); + + printf("=== spice -> --shape fast kernel ===\n"); + printf(" circuit: N=%lu, C_f=%.2f fF, R_Lfast=%.0f ohm, OV=%.1f V\n", + NCELL, CF * 1e15, RLFAST, OV); + printf(" tauLoad = R_Lfast*N*C_f = %7.3f ns (-> params 'tauLoad')\n", tau_load * 1e9); + printf(" tauRecovery = loaded fit = %7.3f ns (fast-tail recovery)\n", tau_rec * 1e9); + printf(" kappa = C_f/(C_d+C_q) = %.4e\n", kappa); + printf(" FWHM: spice %.2f ns vs analytic kernel %.2f ns\n", fwhm_spice * 1e9, fwhm_kernel * 1e9); + printf(" peak-normalised shape RMS error (kernel vs spice, <40 ns): %.3f\n", rms); + const double tol = 0.10; + const bool ok = rms < tol; + printf(" [%s] shape RMS %.3f %s %.2f\n", ok ? "PASS" : "FAIL", rms, ok ? "<" : ">=", tol); + + // 6. Build the tabulated kernel straight from the SPICE curve. This is the + // actual device pulse -- the --shape kernel mode convolves it with the + // avalanche charge train, with no two-pole approximation. The kernel is + // g(t) = i_f(t) / Q_av, the fast-terminal current per unit avalanche + // charge (units 1/s), resampled from the (non-uniform) SPICE grid onto a + // uniform grid at the simulation dt and zeroed before the avalanche (T0). + const double dt_sim = 1e-11; // simulation dt (== svars[0] below) + const double Q_av = (CD + CQ) * OV; // avalanche charge per fired cell [C] + vector kernel; + { + const double tEnd = t.back() - T0; // supported time past the avalanche + const size_t M = (size_t)floor(tEnd / dt_sim) + 1; + kernel.assign(M, 0.0); + size_t j = 0; // monotone cursor into the SPICE samples (tt increases) + for (size_t mi = 0; mi < M; mi++) + { + const double tt = (double)mi * dt_sim + T0; // absolute SPICE time + while (j + 1 < t.size() && t[j + 1] < tt) + j++; + if (j + 1 >= t.size()) + { + kernel[mi] = i_f.back() / Q_av; + continue; + } + const double frac = (tt - t[j]) / (t[j + 1] - t[j]); + const double v = i_f[j] * (1.0 - frac) + i_f[j + 1] * frac; + kernel[mi] = v / Q_av; + } + } + double knet = 0.0; + for (double g : kernel) + knet += g; + knet *= dt_sim; // net charge per unit avalanche charge (want ~0: AC-coupled) + printf(" kernel: %zu taps @ dt = %.1f ps, net charge/avalanche = %+.3e (want ~0)\n", + kernel.size(), dt_sim * 1e12, knet); + + const string kernel_file = "fast_kernel.npy"; + { + NpyWriter kw(kernel_file, kernel.size()); + kw.write(kernel.data(), kernel.size()); + kw.close(); + } + + // 7. Emit a ready-to-run J30020 parameter file pointing at the kernel (the + // remaining device parameters match examples/python/example.py). tauLoad + // is still written so the analytic `fast` mode also works off this file. + vector svars = { + dt_sim, // dt + (double)NCELL, // numMicrocell + OV + 24.5, // vBias = OV + vBr + 24.5, // vBr + 2.2 * 14e-9, // tauRecovery (device) + 0.46, // pdeMax + 2.04, // vChr + 4.6e-14, // cCell + 1.5e-9, // tauFwhm + 0.0}; // digitalThreshold + SiPM sipm(svars); + sipm.tauLoad = tau_load; + sipm.kernelFile = kernel_file; + sipm.kernelDt = dt_sim; + save_params_json(out_params, sipm); + printf(" wrote %s + %s\n", out_params.c_str(), kernel_file.c_str()); + printf(" run: simspad -p %s -i light.npy -o resp.npy --shape kernel\n", out_params.c_str()); + + return ok ? 0 : 1; +} diff --git a/src/utilities.cpp b/src/utilities.cpp index e666daf..bba1af2 100644 --- a/src/utilities.cpp +++ b/src/utilities.cpp @@ -392,6 +392,26 @@ void NpyWriter::close() fout.close(); } +// Slurp a whole 1-D float64 .npy into a vector. Kernels are small (the J-Series +// fast pulse is ~tens of ns / a few thousand taps), so a non-streaming read is +// fine here. +vector read_npy_vector(const string &filename) +{ + NpyReader reader(filename); + size_t n = reader.count(); + vector v(n, 0.0); + size_t got = 0; + while (got < n) + { + size_t k = reader.read(v.data() + got, n - got); + if (k == 0) + break; + got += k; + } + v.resize(got); + return v; +} + // =========================================================================== // Flat-JSON device parameters // =========================================================================== @@ -427,14 +447,47 @@ map parse_flat_json(const string &s) m[key] = val; i = p + (size_t)(endp - start); } + else if (p < s.size() && s[p] == '"') + { + // Quoted string value (e.g. "kernelFile"): skip past its closing + // quote, otherwise the next key search lands inside the value and + // every key after it is misparsed -- the value string becomes a + // key and steals the *next* key's number ("kernelDt" was lost + // this way whenever it followed "kernelFile"). + size_t close = s.find('"', p + 1); + i = (close == string::npos) ? s.size() : close + 1; + } else { - i = colon + 1; // non-numeric value: skip (not expected in our schema) + i = colon + 1; // other non-numeric value: skip } } return m; } +// Extract a string value ("key": "value") from the flat JSON. Returns "" when +// the key is absent or its value is not a quoted string (parse_flat_json only +// handles numeric values, so string-valued keys need this companion). +static string parse_flat_json_string(const string &s, const string &key) +{ + const string needle = "\"" + key + "\""; + size_t k = s.find(needle); + if (k == string::npos) + return ""; + size_t colon = s.find(':', k + needle.size()); + if (colon == string::npos) + return ""; + size_t p = colon + 1; + while (p < s.size() && isspace((unsigned char)s[p])) + ++p; + if (p >= s.size() || s[p] != '"') // value is not a string + return ""; + size_t q2 = s.find('"', p + 1); + if (q2 == string::npos) + return ""; + return s.substr(p + 1, q2 - p - 1); +} + // Parameter key order matches SiPM::dump_configuration(). static const char *kParamKeys[10] = { "dt", "numMicrocell", "vBias", "vBr", "tauRecovery", @@ -447,7 +500,8 @@ SiPM load_params_json(const string &filename) throw runtime_error("cannot open params file: " + filename); stringstream ss; ss << f.rdbuf(); - map m = parse_flat_json(ss.str()); + const string text = ss.str(); + map m = parse_flat_json(text); vector svars(10); for (int i = 0; i < 10; i++) @@ -470,6 +524,23 @@ SiPM load_params_json(const string &filename) throw runtime_error("params JSON: tauLoad must differ from tauRecovery"); sipm.tauLoad = tl->second; } + + // Optional tabulated fast-output kernel (the `kernel` shape mode). The + // "kernelFile" is a path to a 1-D float64 .npy holding the fast-output + // impulse response per unit avalanche charge (units 1/s); "kernelDt" is its + // sample spacing in seconds (defaults to the simulation dt, which is what + // the spice_kernel pipeline writes). The path is resolved relative to the + // current working directory. + string kf = parse_flat_json_string(text, "kernelFile"); + if (!kf.empty()) + { + sipm.kernelFile = kf; + auto kd = m.find("kernelDt"); + sipm.kernelDt = (kd != m.end() && kd->second > 0) ? kd->second : sipm.dt; + sipm.fastKernel = read_npy_vector(kf); + if (sipm.fastKernel.empty()) + throw runtime_error("params JSON: kernelFile loaded no samples: " + kf); + } return sipm; } @@ -490,8 +561,15 @@ string sipm_to_json(SiPM &sipm) } // tauLoad is appended explicitly (not via dump_configuration, whose // 10-double layout the legacy .bin format depends on). - o << " \"tauLoad\": " << sipm.tauLoad << "\n"; - o << "}\n"; + o << " \"tauLoad\": " << sipm.tauLoad; + // The tabulated-kernel keys are emitted only when a kernel is set, so plain + // device files stay unchanged. + if (!sipm.kernelFile.empty()) + { + o << ",\n \"kernelFile\": \"" << sipm.kernelFile << "\""; + o << ",\n \"kernelDt\": " << sipm.kernelDt; + } + o << "\n}\n"; return o.str(); } diff --git a/src/utilities.hpp b/src/utilities.hpp index 2eacbfe..3076838 100644 --- a/src/utilities.hpp +++ b/src/utilities.hpp @@ -74,6 +74,11 @@ SiPM load_params_json(const std::string &filename); std::string sipm_to_json(SiPM &sipm); void save_params_json(const std::string &filename, SiPM &sipm); +// Read a whole 1-D little-endian float64 .npy file into a vector (used for the +// small tabulated kernels of the `kernel` shape mode; for bulk waveforms use +// the streaming NpyReader above). +std::vector read_npy_vector(const std::string &filename); + std::vector conv1d(std::vector inputVec, std::vector kernel); std::vector get_gaussian(double dt, double tauFwhm); diff --git a/test/shaping.hpp b/test/shaping.hpp index 052c5fa..790b094 100644 --- a/test/shaping.hpp +++ b/test/shaping.hpp @@ -214,6 +214,79 @@ bool TEST_shaping() passed_all = passed_all & passed; } + // --- 9. Tabulated kernel: causal FIR, identity and two-tap ------------ + { + // Identity kernel {1} reproduces the input with no time shift. + sipm.fastKernel = {1.0}; + sipm.kernelDt = dt; + vector probe = {0.0, 1.0, 2.0, 3.0, 0.5, 0.0}; + vector out = sipm.shape_kernel(probe); + passed = (out.size() == probe.size()); + for (size_t i = 0; passed && i < probe.size(); i++) + { + passed = passed && (out[i] == probe[i]); + } + + // Causal two-tap kernel {2, 3}: out[i] = 2*x[i] + 3*x[i-1]. + sipm.fastKernel = {2.0, 3.0}; + sipm.kernelDt = dt; + vector out2 = sipm.shape_kernel(probe); + bool ok2 = (out2.size() == probe.size()); + for (size_t i = 0; ok2 && i < probe.size(); i++) + { + double want = 2.0 * probe[i] + (i ? 3.0 * probe[i - 1] : 0.0); + ok2 = ok2 && (fabs(out2[i] - want) < 1e-12); + } + passed = passed && ok2; + cout << "Kernel FIR identity + 2-tap: causal, no shift\t\t\t\t" + << (passed ? "\033[32;49;1mPASS\033[0m" : "\033[31;49;1mFAIL\033[0m") << endl; + passed_all = passed_all & passed; + } + + // --- 10. Tabulated kernel: linear resampling (kernelDt != dt) ---------- + { + // Source kernel {0, 10} at twice the sim dt resamples (linearly) onto + // the dt grid as {0, 5, 10}; feeding a unit impulse reads that back. + sipm.fastKernel = {0.0, 10.0}; + sipm.kernelDt = 2.0 * dt; + vector impulse = {1.0, 0.0, 0.0, 0.0}; + vector out = sipm.shape_kernel(impulse); + passed = (out.size() == impulse.size()) && (fabs(out[0] - 0.0) < 1e-12) && + (fabs(out[1] - 5.0) < 1e-12) && (fabs(out[2] - 10.0) < 1e-12); + cout << "Kernel resample 2*dt -> dt: {0,10} -> {0,5,10}\t\t\t" + << (passed ? "\033[32;49;1mPASS\033[0m" : "\033[31;49;1mFAIL\033[0m") << endl; + passed_all = passed_all & passed; + } + + // --- 11. Tabulated kernel reproduces the analytic fast shaper ---------- + // Sampling the analytic bipolar response h(t) at dt and using it as the + // kernel should reproduce shape_fast() (which is the same h binned by a + // one-pole pair) to within the discretisation difference. Ties the new + // faithful path back to the validated analytic one. + { + const double A = tauRecovery / (tauRecovery - tauLoad); + vector g(n, 0.0); + for (int i = 0; i < n; i++) + { + double t = (i + 0.5) * dt; // bin midpoint, matching the one-pole binning + g[i] = A * (exp(-t / tauLoad) / tauLoad - exp(-t / tauRecovery) / tauRecovery) * dt; + } + sipm.fastKernel = g; + sipm.kernelDt = dt; + vector kshaped = sipm.shape_kernel(q); // q is the unit impulse from above + double hMax = 0.0, errMax = 0.0; + for (int i = 0; i < n; i++) + { + hMax = max(hMax, fabs(shaped[i])); + errMax = max(errMax, fabs(kshaped[i] - shaped[i])); + } + double relErr = errMax / hMax; + passed = relErr < 1e-3; + cout << "Kernel vs analytic fast: max rel err " << relErr << " (want < 1e-3)\t" + << (passed ? "\033[32;49;1mPASS\033[0m" : "\033[31;49;1mFAIL\033[0m") << endl; + passed_all = passed_all & passed; + } + string prefix = passed_all ? "\033[32;49;1m" : "\033[31;49;1m"; string outStatus = passed_all ? "PASS\n" : "FAIL\a\n"; cout << prefix << BAR_STRING << endl;