GridBased/Notebooks/XRF_AgMo_1_factor_10.ipynb

XRF AgMo_1_factor_10 -- polychromatic GridBased model on measured lines¶

2603 MaxIV, MOENCH014 at 150 V, full 400x400 chip flood-illuminated by the fluorescence of an Ag/Mo sample; slit factor 10, 0.67 % occupancy, 1.0e8 photons in 100 000 frames. Mo Ka 17.44 / Kb 19.64 and Ag Ka 22.10 keV sit in the middle of the 3-30 keV training range; Compton + elastic scatter forms a bump at 28-30 keV. The clean case -- XRF_CuCrCo_1_factor_10.ipynb has the low-energy lines.

Model Results/260919_poly3to30keV_Moench014_c96_00, the default since 2026-09-19: [32, 64, 96, 96], flat 3-30 keV, chosen from four sizes in Results/overnight_poly3to30keV/report. Reconstruction from scripts/05_infer_measurement.py --config Configs/infer_measurement_poly_Moench014.yaml --run AgMo_1_factor_10 --save-photons (FP16, masked column common mode). The previous 5-30 keV model's reconstructions are kept as <run>_poly5to30. Compared throughout:

• grid direct -- the model output;
• grid corrected -- $E = g^{-1}(E_\mathrm{pred})$, where $g(E)$ is the MC position of a line of true energy $E$ (scripts/07_energy_response.py);
• classical 3x3 -- 3x3 sum around a local maximum on the same calibrated frames (first 10 000 frames), no model.

In [1]:

import json
import sys
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

REPO = Path('/home/xie_x1/MLXID/GridBased')
sys.path.insert(0, str(REPO / 'src'))
from xrf_analysis import (EnergyCorrection, classical_and_model, fit_lines,
                          line_centroid, load_reconstruction)

CONFIG = REPO / 'Configs/infer_measurement_poly_Moench014.yaml'
RUN = 'AgMo_1_factor_10'
EXPERIMENT = REPO / 'Results/260919_poly3to30keV_Moench014_c96_00'
RESPONSE = EXPERIMENT / 'energy_response.npz'
GRID = Path('/mnt/sls_det_storage/moench_data/MLXID/InferenceResults/GridBased/XRF_2603MaxIV')
EDGES = np.arange(0.0, 40.025, 0.05)             # 50 eV bins
CENTRES = 0.5 * (EDGES[1:] + EDGES[:-1])

FIT_GROUPS = [       # each line fitted on its own window
    dict(label='Mo Ka', alpha=['Mo_Ka'], window=(16.0, 18.9)),
    dict(label='Mo Kb', alpha=['Mo_Kb'], window=(18.75, 20.7)),
    dict(label='Ag Ka', alpha=['Ag_Ka'], window=(20.6, 23.6)),
]   # Ag Kb sits on the rising Compton bump and is left out
SPECTRUM_ZOOM = (14.5, 27.5)
LINES = [n for g in FIT_GROUPS for n in g['alpha']]
PEAKS = LINES + [e[0] for g in FIT_GROUPS for e in g.get('extra', ())]
COLOURS = {'grid direct': '#1f6fb4', 'grid corrected': '#10695f',
           'classical 3x3': '#6b6b6b'}

In [2]:

reconstruction = load_reconstruction(GRID, RUN)
summary = reconstruction['summary']
n_frames = summary['n_frames']
energy = np.array(reconstruction['photons'][:, 2])
correction = EnergyCorrection(RESPONSE)                    # core-based g
correction_mean = EnergyCorrection(RESPONSE, estimator='mean')

reference = classical_and_model(CONFIG, RUN, n_frames=10_000)   # cached
classical, grid_subset = reference['classical'], reference['grid']
n_subset = int(reference['n_frames'])

energies = {'grid direct': energy, 'grid corrected': correction(energy),
            'classical 3x3': classical[:, 3]}
frames_of = {'grid direct': n_frames, 'grid corrected': n_frames,
             'classical 3x3': n_subset}
print(f"{RUN}: {n_frames} frames, {summary['photons_per_frame']:.0f} photons/frame, "
      f"occupancy {summary['occupancy_valid_region']*100:.2f} %, "
      f"{(energy < 3).mean()*100:.1f} % of the photons below 3 keV")

AgMo_1_factor_10: 100000 frames, 1054 photons/frame, occupancy 0.69 %, 1.1 % of the photons below 3 keV

1. The MC energy response¶

The energy head is a regression trained on a flat 3-30 keV spectrum. On MC (occupancy <= 2.5 %) the core of the predicted energy -- what a line fit measures -- stays within ~25 meV of the truth from 4 to 28 keV and moves only within ~0.75 keV of the edges (+98 meV at 3.2 keV, -93 meV at 29.6 keV), where the energy scale is also compressed. The mean is pulled towards the centre by ~0.1 keV at 5-12 keV and more near the edges, because the few per cent of photons whose energy is off by more than 1 keV (the tails) are cut by the bounded range. That mean is the energy_bias_keV of eval_summary.json. It is not a line shift, and correcting with it moves lines at 5-12 keV by ~0.1 keV the wrong way. The correction used here therefore inverts the core position $g(E)$; with the 3-30 keV model it matters only below ~4 keV.

In [3]:

response = np.load(RESPONSE)
centre = response['table_centre']
figure, axes = plt.subplots(1, 2, figsize=(13, 4.2))
for curve, colour, label in ((correction, COLOURS['grid corrected'], 'core of E_pred (used)'),
                             (correction_mean, '#c0392b', 'mean of E_pred')):
    axes[0].plot(curve.grid, curve.g_grid - curve.grid, color=colour, lw=1.5, label=label)
for name in LINES:
    axes[0].axvline(line_centroid(name), color='k', lw=.5, ls=':')
axes[0].axhline(0, color='k', lw=.6); axes[0].set_ylim(-0.6, 0.3)
axes[0].set_xlabel('E_true [keV]'); axes[0].set_ylabel('g(E) - E [keV]')
axes[0].set_title('MC: where a line of true energy E appears'); axes[0].legend(fontsize=8)
axes[0].grid(alpha=.3)
axes[1].plot(centre, response['table_core_sigma'], 'k.', ms=3, label='core sigma, raw')
axes[1].plot(correction.grid, correction.mc_sigma(correction.grid) / correction.slope(correction.grid),
             color=COLOURS['grid corrected'], lw=1.2, label="core sigma / g'  (true keV)")
axes[1].set_ylim(0, 0.4); axes[1].set_xlabel('E_true [keV]'); axes[1].set_ylabel('keV')
axes[1].set_title('MC: energy resolution'); axes[1].legend(fontsize=8); axes[1].grid(alpha=.3)
figure.tight_layout()
for name in LINES:
    e = line_centroid(name)
    print(f"{name}: g - E = {(correction.forward(np.array([e]))[0] - e)*1000:+.0f} meV, "
          f"g' = {correction.slope(e):.2f}, MC core sigma {correction.mc_sigma(e):.3f} keV")

Mo_Ka: g - E = -12 meV, g' = 1.00, MC core sigma 0.275 keV
Mo_Kb: g - E = -8 meV, g' = 1.00, MC core sigma 0.290 keV
Ag_Ka: g - E = -6 meV, g' = 1.00, MC core sigma 0.301 keV

2. Spectra¶

Per frame, 50 eV bins; grey = outside the 3-30 keV training range.

In [4]:

histograms = {name: np.histogram(values, EDGES)[0] for name, values in energies.items()}
figure, axes = plt.subplots(2, 1, figsize=(14, 8.5))
for axis, log in zip(axes, (True, False)):
    for name, counts in histograms.items():
        axis.plot(CENTRES, counts / frames_of[name], lw=1, color=COLOURS[name], label=name)
    for name in LINES:
        axis.axvline(line_centroid(name), color='k', lw=.4, ls='--')
    axis.axvspan(0, 3, color='0.9'); axis.axvspan(30, 40, color='0.9')
    axis.set_xlabel('energy [keV]'); axis.set_ylabel('photons / frame / 50 eV')
    axis.grid(alpha=.3)
    if log:
        axis.set_yscale('log'); axis.set_xlim(2, 33); axis.set_ylim(1e-4, None)
        axis.legend(fontsize=9); axis.set_title(RUN)
    else:
        axis.set_xlim(*SPECTRUM_ZOOM)
figure.tight_layout()

3. Line positions and energy resolution¶

Grid lines: Gaussian core + a wide Gaussian at the same position (the tails of section 1) + an erfc low-energy step, on an exp(polynomial) background. Classical lines: Gaussian + step. Ka is fitted as the tabulated doublet, so the position is the intensity-weighted centroid. sigma is the core width; FWHM is that of the whole line shape.

In [5]:

fits = {}
for variant, counts in histograms.items():
    for group in FIT_GROUPS:
        fits[variant, group['label']] = fit_lines(
            CENTRES, counts, group['alpha'], group.get('beta'), window=group['window'],
            beta_ratio_fixed=group.get('ratio_fixed'),
            fit_stretch=group.get('fit_stretch', True), extra_peaks=group.get('extra', ()),
            tails=variant != 'classical 3x3')


def peak(variant, name):
    group = next(g for g in FIT_GROUPS
                 if name in g['alpha'] or name in [e[0] for e in g.get('extra', ())])
    return fits[variant, group['label']]['lines'][name]


def fwhm_of(variant, name):
    """FWHM of the line shape; for the extra peaks (no doublet) from sigma."""
    line = peak(variant, name)
    return line.get('fwhm', 2.3548 * line['sigma'])


print(f"{'':8s} {'table':>6s} | {'position [keV]':^26s} | {'core sigma [keV]':^26s} | "
      f"{'FWHM [keV]':^20s} | MC sigma")
print(f"{'':8s} {'':>6s} | {'direct':>8s} {'corr.':>8s} {'3x3':>8s} | "
      f"{'direct':>8s} {'corr.':>8s} {'3x3':>8s} | {'grid':>6s} {'3x3':>6s} {'ratio':>6s} |")
for name in PEAKS:
    table = f"{line_centroid(name):6.2f}" if name in LINES else f"{'':6s}"
    positions = [peak(v, name)['mu'] for v in energies]
    sigmas = [peak(v, name)['sigma'] for v in energies]
    grid_fwhm = fwhm_of('grid direct', name)
    classical_fwhm = fwhm_of('classical 3x3', name)
    mc = correction.mc_sigma(line_centroid(name)) if name in LINES else \
        correction.mc_sigma(peak('grid corrected', name)['mu'])
    print(f"{name:8s} {table} | " + ' '.join(f'{p:8.2f}' for p in positions) + ' | ' +
          ' '.join(f'{s:8.3f}' for s in sigmas) +
          f" | {grid_fwhm:6.2f} {classical_fwhm:6.2f} {grid_fwhm/classical_fwhm:6.2f} | {mc:.3f}")
print('grid FWHM: direct output')

          table |       position [keV]       |      core sigma [keV]      |      FWHM [keV]      | MC sigma
                |   direct    corr.      3x3 |   direct    corr.      3x3 |   grid    3x3  ratio |
Mo_Ka     17.44 |    17.47    17.48    17.52 |    0.327    0.327    0.464 |   0.82   1.09   0.75 | 0.275
Mo_Kb     19.64 |    19.66    19.67    19.74 |    0.353    0.352    0.492 |   0.83   1.16   0.72 | 0.290
Ag_Ka     22.10 |    22.10    22.11    22.17 |    0.378    0.378    0.489 |   0.95   1.15   0.82 | 0.301
grid FWHM: direct output

In [6]:

shown = ('grid direct', 'classical 3x3')
figure, axes = plt.subplots(2, len(FIT_GROUPS) * 2, figsize=(4.4 * len(FIT_GROUPS) * 2, 6),
                            gridspec_kw=dict(height_ratios=[3, 1]), squeeze=False)
for column, (group, variant) in enumerate((g, v) for g in FIT_GROUPS for v in shown):
    fit = fits[variant, group['label']]
    x = fit['centres']
    axes[0, column].plot(x, fit['counts'], 'k.', ms=2.5)
    axes[0, column].plot(x, fit['model'](x), color=COLOURS[variant], lw=1.3)
    for part, curve in fit['components'](x).items():
        axes[0, column].plot(x, curve, lw=.7, ls='--', label=part)
    axes[0, column].set_yscale('log')
    axes[0, column].set_ylim(max(fit['counts'].min() * 0.5, 1), fit['counts'].max() * 2)
    axes[0, column].set_title(f"{variant}: {group['label']}", fontsize=9)
    axes[0, column].legend(fontsize=6)
    axes[1, column].plot(x, (fit['counts'] - fit['model'](x)) / np.sqrt(np.maximum(fit['counts'], 1)),
                         '.', ms=2.5, color=COLOURS[variant])
    axes[1, column].axhline(0, color='k', lw=.6)
    axes[1, column].set_xlabel('energy [keV]'); axes[1, column].set_ylabel('pull')
figure.tight_layout()

In [7]:

figure, axis = plt.subplots(figsize=(7.5, 4.4))
for offset, variant in zip((-0.08, 0.0, 0.08), energies):
    axis.errorbar([line_centroid(n) + offset for n in LINES],
                  [peak(variant, n)['sigma'] for n in LINES],
                  [peak(variant, n)['sigma_err'] for n in LINES],
                  fmt='o', ms=5, capsize=3, color=COLOURS[variant], label=variant)
e = np.linspace(min(line_centroid(n) for n in LINES) - 1, max(line_centroid(n) for n in LINES) + 1, 200)
e = e[(e >= 5) & (e <= 30)]
axis.plot(e, correction.mc_sigma(e), 'k--', lw=1, label='MC core sigma (raw)')
axis.plot(e, correction.mc_sigma(e) / correction.slope(e), 'k:', lw=1, label="MC / g' (true keV)")
axis.set_xlabel('line energy [keV]'); axis.set_ylabel('core sigma [keV]'); axis.set_ylim(0, None)
axis.set_title('energy resolution'); axis.legend(fontsize=8); axis.grid(alpha=.3)
figure.tight_layout()

4. Reconstruction efficiency¶

Without truth, three handles on the 10 000 frames of the classical reference:

• isolated photons: classical clusters with no other cluster within 4 px that have a grid photon within 1.5 px. Above 5 keV the 2.5 keV classical threshold finds every isolated photon, so this is the absolute efficiency. The MC equivalent is photons with no neighbour within 4 px.
• pile-up: classical clusters holding two or more grid photons within 2 px.
• energy balance: summed reconstructed energy / deposited energy per frame (calibrated frame minus its baseline).

In [8]:

match = reference['classical_match']
found = match >= 0
isolated = reference['isolation_classical'] > 4.0
above5 = classical[:, 3] >= 5
edges_e = np.arange(2.5, 30.01, 0.5)
which = np.digitize(classical[:, 3], edges_e) - 1
fraction = np.array([found[isolated & (which == i)].mean() if (isolated & (which == i)).sum() > 200
                     else np.nan for i in range(len(edges_e) - 1)])

mc = json.load(open(EXPERIMENT / 'eval_summary.json'))['1%']['pair']
far = [i for i, edge in enumerate(mc['edges']) if edge >= 4.0]
mc_isolated = sum(mc['n_match'][i] for i in far) / sum(mc['n_true'][i] for i in far)

figure, axis = plt.subplots(figsize=(7.5, 4.2))
axis.plot(0.5 * (edges_e[1:] + edges_e[:-1]), fraction * 100, 'o-', ms=3,
          color=COLOURS['grid direct'], label='isolated classical clusters found by the grid')
axis.axhline(mc_isolated * 100, color='k', ls='--', lw=1, label='MC, isolated photons, 1 %')
axis.axvspan(0, 3, color='0.9'); axis.set_ylim(90, 100.5); axis.set_xlim(2.5, 30)
axis.set_xlabel('classical 3x3 energy [keV]'); axis.set_ylabel('%')
axis.set_title('detection efficiency'); axis.legend(fontsize=8); axis.grid(alpha=.3)
figure.tight_layout()

multiplicity = reference['grid_per_classical']
occupancy = summary['occupancy_valid_region']
balance = reference['balance'][:, 1].sum()
low = (classical[:, 3] >= 3) & ~above5
print(f"isolated photons >= 5 keV found: {found[isolated & above5].mean()*100:.2f} %   "
      f"(3-5 keV: {found[isolated & low].mean()*100:.1f} %;  MC: {mc_isolated*100:.2f} %)")
print(f"grid / classical photon count: {len(grid_subset) / len(classical):.3f}")
print(f"classical clusters with >= 2 grid photons: {(multiplicity >= 2).mean()*100:.1f} %   "
      f"(Poisson, second photon within 2 px: {(1 - np.exp(-occupancy * np.pi * 4))*100:.1f} %)")
print(f"energy balance: grid {grid_subset[:, 3].sum() / balance:.4f}, "
      f"classical {classical[:, 3].sum() / balance:.4f}")

isolated photons >= 5 keV found: 99.97 %   (3-5 keV: 92.9 %;  MC: 99.96 %)
grid / classical photon count: 1.091
classical clusters with >= 2 grid photons: 8.5 %   (Poisson, second photon within 2 px: 8.3 %)
energy balance: grid 0.9968, classical 0.9980

Conclusions¶

• Lines. All resolved; the grid puts them within 0.03 keV of the tables, the classical 3x3 sum 0.06-0.10 keV high (the MOENCH014 LUT reads 0.3-0.5 % high here). The MC-based correction changes nothing at these energies (<= 12 meV).
• Energy resolution. FWHM 0.82 / 0.83 / 0.95 keV for Mo Ka / Mo Kb / Ag Ka against 1.09 / 1.16 / 1.15 keV for the classical sum, 18-28 % narrower. The core widths (0.33-0.38 keV) are 19-26 % above the MC (0.28-0.30 keV): at these energies the data carry an extra term of ~1 % of E (pixel-to-pixel gain residuals of the calibration, not in the MC). The same holds for the previous 5-30 keV model.
• Efficiency. 99.97 % of the isolated photons above 5 keV (MC 99.96 %). The grid finds 9.1 % more photons than the classical finder: pile-up the 3x3 finder merges (8.5 % of its clusters hold two grid photons, 8.3 % expected from Poisson), plus low-energy photons that the 3-30 keV model now reports. Energy balance 0.997.

In [9]:

output = dict(
    run=RUN, epoch=summary['epoch'], n_frames=n_frames,
    photons_per_frame=summary['photons_per_frame'], occupancy=occupancy,
    lines={name: {variant: dict(mu=float(peak(variant, name)['mu']),
                                sigma=float(peak(variant, name)['sigma']),
                                fwhm=float(fwhm_of(variant, name)))
                  for variant in energies}
           for name in PEAKS},
    efficiency_isolated_above5keV=float(found[isolated & above5].mean()),
    efficiency_isolated_mc=mc_isolated,
    grid_over_classical_count=len(grid_subset) / len(classical),
    energy_balance_grid=float(grid_subset[:, 3].sum() / balance),
)
with open(GRID / RUN / 'xrf_energy_analysis.json', 'w') as handle:
    json.dump(output, handle, indent=1)