Gamma-Ray Spectroscopy

Spectrum smoothing, peak counting, energy/efficiency calibration, X-ray line algebra, SDEF decay-source cards, and the SPE readers behind Nucleide's spectroscopy toolkit.

Overview

Gamma-ray pulse-height spectra arrive as counts per channel plus a header (recording times, detector id, calibration fits). Answering “how many counts sit under this peak, and at what energy” needs a small pipeline: smooth the spectrum, delimit the peak, subtract the background, map channels to energies, correct for detector efficiency, and account for X-ray lines. Nucleide implements exactly the upstream analysis subset (E1–E8) in crates/spectroscopy, with every caller-side constant kept explicit, plus the E9 decay-source composition: caller-supplied decay lines normalized into emission probabilities and rendered as an MCNP SDEF card.

Mathematical details

E1 — Rectangular smoothing

With window width m (odd, m >= 3) and ext = (m - 1) / 2:

smooth[i] = (1/m) * sum(counts[i-ext..=i+ext]) for ext <= i < N - ext.

The first and last ext channels have no full neighbourhood and are copied unchanged.

E2 — Five-point smoothing

smooth[i] = (c[i-2] + c[i+2] + 2*c[i+1] + 2*c[i-1] + 3*c[i]) / 9 for 2 <= i < N - 2; the first and last two channels are copied. The weights are the low-statistics recommendation cited by the upstream module as “G.W. Phillips, Nucl. Instrum. Methods 153 (1978), 449” — carried here as a citation string only, with no journal-verification claim.

E3 — Background under a peak (`m == 1` only)

low = sum(counts[c1-2..c1]), high = sum(counts[c2..c2+2]), and bg = (low + high) * (c2 - c1 + 1) / 6. Only the m == 1 two-sided estimator exists; any other m is an error. There is no Compton-continuum or step model.

E4 — Gross counts

gc = sum(counts[c1..c2]), half-open: channel c2 is excluded.

E5 — Net counts

nc = gc - bg.

E3–E5 in one picture: sideband sums set the background under the c1..c2 window; net = gross − background.

E6 — Energy calibration

ebin[ch] = a0 + a1*ch + a2*ch^2 from the calib_e_fit triplet. FWHM coefficients are parsed from the file header but never evaluated.

E7 — Detector-efficiency evaluation (energy in MeV)

With caller-supplied coefficients a_i:

fit == 1: eff = exp(sum a_i * (ln E)^i);
fit == 2: eff = exp(sum a_i * (1/E)^i);
anything else is an error.

E7-fit — Efficiency-coefficient fit (extension of E7, not a new equation)

The upstream module ships no efficiency-calibration fitting routine, so the fit has no upstream parity — it is pinned by synthetic closed-form gates only. Over caller points (E_i, eff_i, w_i) with E_i > 0, eff_i > 0, and caller-supplied weights w_i >= 0, the fit solves the log-space linear least squares

y_i = ln eff_i, X[i][j] = (ln E_i)^j (fit == 1) or (1/E_i)^j (fit == 2), min_β sum w_i (X_i β − y_i)^2,

returning order + 1 coefficients that feed the E7 evaluation above. The solve runs through the workspace least-squares kernel (nucleide-linalg lstsq, faer column-pivoted QR); a zero weight skips its row, all-zero weights and fewer points than coefficients are errors. This is the one coefficient-fitting path in the toolkit — every other coefficient vector (energy fits, X-ray constants, decay lines) stays a caller input.

E8 — X-ray line algebra

From caller-supplied atomic constants (K/L fluorescence yields, the kb_to_ka / ka2_to_ka1 ratios, the K-fill probability, and the four line energies) plus K/L conversion coefficients:

with k_conv: xk = k_fluor * k_conv, xka = xk / (1 + kb_to_ka), xka1 = xka / (1 + ka2_to_ka1), xka2 = xka - xka1, xkb = xk - xka, and — when l_conv is also present — xl = (l_conv + k_conv * prob) * l_fluor;
with only l_conv: xl = l_conv * l_fluor;
with neither, all intensities are zero.

Output order is always Ka1, Ka2, Kb, L. No atomic table is vendored, so every constant is an explicit input; an absent conversion is None (the upstream NaN sentinel).

E9 — Decay-source normalization and SDEF cards

The upstream point-source emitter renders one monoenergetic particle at a point: SDEF POS=<x> <y> <z>, VEC=<u> <v> <w> DIR=1 on a continuation line only when a direction component is nonzero, then ERG=<E>, WGT=<w>, and PAR=<designator>. E9 generalizes it to a caller-supplied line list (E_i, I_i) — no line data is vendored, every pair is an input. With duplicates merged by summing and bins sorted ascending:

p_i = I_i / sum_j I_j, so the emitted probabilities sum to 1.0.

Behavior is specified, not inherited: empty input, negative energies or intensities, and non-finite values are errors; zero-intensity lines are dropped (a zero-probability SP1 entry has unverified transport meaning); input order is irrelevant.

Line lists also arrive from files: the runtime decay-lines TSV interchange (energy_MeV intensity pairs, # comments, blank lines skipped) parses to the same (E_i, I_i) pairs E9 consumes — same validation, no vendored data. No legacy line-list reader ships: upstream line families query a compiled data store with no text format to mirror, ENSDF text is out of proportion, ICRP-107 exports are license-blocked, and ENDF mean energies are not lines.

The card keeps the upstream field order. One surviving line keeps the upstream inline ERG=<E>; several lines switch to the discrete-distribution form ERG=D1, with the energies on an SI1 L card and the probabilities on a paired SP1 D card, wrapped at 80 columns onto 5-space continuations. Numbers render in C++ default six-significant-digit form, matching the upstream C++ emitter. The PAR= designator comes from the nucleide-nuclei particle dialect (MCNP5 n/p/e; MCNP6 adds h for protons).

The distribution syntax is verified surface — the card shapes lex, parse, and round-trip through an MCNP input parser. The sampling semantics behind them (what the L/D flags select, SPn-omitted defaults, zero/negative probability handling, MCNP renormalization) are not verified from local files and are the caller’s responsibility against the MCNP manual. Re-checked against the parser suite: it exercises only this surface, so the semantics stay unverified.

Algorithm

Smoothing and counting are single passes over plain vectors. Channel bounds c1/c2 index the counts vector positionally: a nonzero start_chan_num never shifts the slice (the upstream quirk), while the c2 <= max(channels) range check runs first and out-of-range slice stops clamp like Python slices.

The .spe readers reproduce the upstream grammars quirk-for-quirk. Dollar format requires the $SPEC_ID: magic (else an unsupported-format error), resolves duplicate tags first-wins, raises on missing tags, reads $MEAS_TIM: as live-then-real, takes $DATA: start last with last + 1 one-float lines after it, and reads the $MCA_CAL:/$SHAPE_CAL: triplets two lines below each tag (a trailing keV token is ignored). $ROI:, $PRESETS:, and $ENER_FIT: are never read. Plain format rejects the dollar magic, reads key: value headers (unknown keys ignored, later keys overwrite), parses channel: counts lines after the SPECTRUM marker, splits double-spaced Energy Fit:/FWHM Fit: values at single-space indices 0/2/4, and re-joins the three colon parts of the acquisition time.

Worked example

Counts [2, 5, 1, 6, 3, 8, 4], peak window c1 = 2, c2 = 5:

E1 (m = 5): [2, 5, 17/5, 23/5, 22/5, 8, 4].
E2: [2, 5, 30/9, 39/9, 42/9, 8, 4].
E3: low = 7, high = 12, bg = 19 * 4 / 6 = 76/6.
E4: gc = 1 + 6 + 3 = 10. E5: nc = 10 - 76/6 = -16/6.

Assumptions and limitations

Counts, channels, fits, atomic constants, decay lines, and efficiency-fit points/weights are all caller inputs; nothing is fetched. Coefficient fitting exists on exactly one path (E7-fit above); every other coefficient vector is an input, never fitted.
FWHM fits are parsed, never evaluated: no evaluation form for the triplet is pinned by the upstream analysis subset, so no evaluator exists here.
X-ray energies are in keV; efficiency energies are in MeV; SDEF decay-line energies are in MeV.
Smoothing edges copy rather than extrapolate.
E9 renders the SDEF/distribution card shapes only; MCNP sampling semantics (distribution flags, defaults, renormalization) are the caller’s responsibility against the MCNP manual.

Peak search, peak fitting, activity calculations, and plotting all sit outside this toolkit: none exist in the upstream analysis subset, and plotting stays user-side (the upstream plot helpers are an isolated GUI module). The E7-fit efficiency-coefficient fit is not peak fitting — it fits the detector-response curve to caller calibration points. Decay-spectrum line data also stays out: no evaluated line library is vendored, so E9 sources are built purely from caller-supplied lines. The validation/spectroscopy_vs_pyne.py harness checks the toolkit two ways: synthetic closed-form checks in-repo and a live PyNE cross-check in the validation container.