Skip to content

Many thicknesses and extrapolation

Running Monte Carlo at every thickness is expensive. A more practical approach is to run Monte Carlo at a few thin shields, then interpolate and extrapolate the build-up factor to all thicknesses.

Why Gaussian Processes?

Build-up factors don't follow a simple functional form; they depend on material, energy, geometry, and thickness in complex ways. Fitting a polynomial or exponential would impose assumptions about the shape that may not hold.

A Gaussian Process (GP) makes no assumption about the functional form. Instead it:

  • Learns the shape from the data: whatever the Monte Carlo points show, the GP follows
  • Gives uncertainty estimates: tight near data, wider far away, honestly reflecting what we don't know
  • Accounts for Monte Carlo noise: each data point has a statistical uncertainty that the GP uses to avoid overfitting

This makes it robust across different materials and geometries without needing to choose a fitting function.

BuildupTable wraps the inference-tools GP with adaptive kernel length scales to prevent overfitting.

Basic example

Run Monte Carlo at 4 thicknesses, then extrapolate to 10:

import rad_point_kernel as rpk

concrete = rpk.Material(
    composition={
        "H": 0.01, "O": 0.53, "Si": 0.34,
        "Ca": 0.04, "Al": 0.03, "Fe": 0.01,
    },
    density=2.3,
    fraction="mass",
)

SOURCE = 1e12
source = rpk.Source(particle="photon", energy=1e6)

# Step 1: Monte Carlo at 4 thicknesses
mc_thicknesses = [5, 10, 15, 20]
mc_geometries = [
    [rpk.Layer(thickness=t, material=concrete)]
    for t in mc_thicknesses
]

mc_results = rpk.compute_buildup(
    geometries=mc_geometries,
    source=source,
    quantities=["dose-AP"],
)

for t, r in zip(mc_thicknesses, mc_results):
    print(f"  {t:>2d} cm: B = {r.buildup['dose-AP']}")

# Step 2: Build interpolation table
table = rpk.BuildupTable(
    points=[{"thickness": t} for t in mc_thicknesses],
    results=mc_results,
)

# Step 3: Predict build-up at any thickness
all_thicknesses = mc_thicknesses + [30, 50, 75, 100, 150, 200]
for t in all_thicknesses:
    layers = [rpk.Layer(thickness=t, material=concrete)]
    bi = table.interpolate(thickness=t)

    result = rpk.calculate_dose(
        source_strength=SOURCE,
        layers=layers,
        source=source,
        geometry="AP",
        buildup=bi,
    )

    status = "EXTRAPOLATED" if bi.is_extrapolated else "interpolated"
    print(f"  {t:>3d} cm: dose = {result.dose_rate} Sv/hr, "
          f"B = {bi.value} +/- {bi.sigma} ({status})")

   5 cm: B = 1.503618772779222
  10 cm: B = 2.128954332208891
  15 cm: B = 2.83890281280358
  20 cm: B = 3.6120520881534257
    5 cm: dose = 37.02471193455802 Sv/hr, B = 1.503629832611645 +/- 0.0028376455152916086 (interpolated)
   10 cm: dose = 6.27082541694349 Sv/hr, B = 2.128541692543146 +/- 0.006592418865115662 (interpolated)
   15 cm: dose = 1.780317494068822 Sv/hr, B = 2.841094264402011 +/- 0.008540907206003022 (interpolated)
   20 cm: dose = 0.6085888766755234 Sv/hr, B = 3.607762018475563 +/- 0.01972956034191378 (interpolated)
   30 cm: dose = 0.08843746674627623 Sv/hr, B = 5.150259511245304 +/- 0.13181979189344095 (EXTRAPOLATED)
   50 cm: dose = 0.0023670647790297265 Sv/hr, B = 7.299503474915477 +/- 0.8977759872977322 (EXTRAPOLATED)
   75 cm: dose = 2.5906605180319583e-05 Sv/hr, B = 7.160070088664549 +/- 2.353390760619516 (EXTRAPOLATED)
  100 cm: dose = 2.816291604352237e-07 Sv/hr, B = 5.511920354207172 +/- 3.221983208615853 (EXTRAPOLATED)
  150 cm: dose = 5.7471100471521575e-11 Sv/hr, B = 4.015498700434222 +/- 3.4556898298339447 (EXTRAPOLATED)
  200 cm: dose = 1.996624316552037e-14 Sv/hr, B = 3.9350110032853722 +/- 3.456496825511351 (EXTRAPOLATED)

InterpolationResult

table.interpolate() returns an InterpolationResult with:

  • value - predicted build-up factor
  • sigma - GP uncertainty (1 standard deviation)
  • is_extrapolated - True if the query point is outside the range of Monte Carlo data
  • extrapolated_axes - dict showing which axes are extrapolated and by how much

Uncertainty grows with distance from data points:

r_near = table.interpolate(thickness=15)
r_far = table.interpolate(thickness=200)
print(f"Near MC data (t=15):  sigma = {r_near.sigma}")
print(f"Far from data (t=200): sigma = {r_far.sigma}")

Near MC data (t=15):  sigma = 0.008540907206003022
Far from data (t=200): sigma = 3.456496825511351

Using InterpolationResult as build-up

You can pass an InterpolationResult directly to any calculation function:

bi = table.interpolate(thickness=50)
result = rpk.calculate_dose(
    source_strength=1e12,
    layers=[rpk.Layer(thickness=50, material=concrete)],
    source=source,
    geometry="AP",
    buildup=bi,
)
print(f"Dose at 50 cm concrete: {result.dose_rate} Sv/hr (B = {bi.value})")

Dose at 50 cm concrete: 0.0023670647790297265 Sv/hr (B = 7.299503474915477)

Using your own build-up values

You don't have to use BuildupTable - if you have build-up factors from another source (literature, your own fitting, a different interpolation method), you can apply them directly:

# Your own B value from any source
my_B = 3.2
layers = [rpk.Layer(thickness=50, material=concrete)]
result = rpk.calculate_dose(
    source_strength=1e12,
    layers=layers,
    source=source,
    geometry="AP",
    buildup=rpk.BuildupModel.constant(my_B),
)
print(f"Dose with B={my_B}: {result.dose_rate} Sv/hr")

Dose with B=3.2: 0.0010376880179487598 Sv/hr

This means you can use any interpolation or fitting method you prefer (scipy, scikit-learn, a lookup table, or even a hand-drawn curve) and feed the result into the point-kernel calculation.

Multi-dimensional extrapolation

BuildupTable supports N-dimensional parameter spaces. For example, with water and concrete thicknesses as two axes:

import rad_point_kernel as rpk

water = rpk.Material(composition={"H2O": 1.0}, density=1.0)
concrete = rpk.Material(
    composition={
        "H": 0.01, "O": 0.53, "Si": 0.34,
        "Ca": 0.04, "Al": 0.03, "Fe": 0.01,
    },
    density=2.3,
    fraction="mass",
)

# Monte Carlo at a grid of (water, concrete) thicknesses
mc_water = [0, 10, 20]
mc_conc = [10, 20, 30]
points = []
geometries = []
for w in mc_water:
    for c in mc_conc:
        points.append({"water": w, "conc": c})
        layers = [rpk.Layer(thickness=1000)]
        if w > 0:
            layers.append(rpk.Layer(thickness=w, material=water))
        layers.append(rpk.Layer(thickness=c, material=concrete))
        geometries.append(layers)

source = rpk.Source(particle="neutron", energy=14.06e6)
mc_results = rpk.compute_buildup(
    geometries=geometries,
    source=source,
    quantities=["dose-AP"],
)

table = rpk.BuildupTable(points=points, results=mc_results)

# Query at any (water, concrete) combination
bi = table.interpolate(water=15, conc=25)
print(f"B = {bi.value} +/- {bi.sigma}")

B = 1.1165191179805656 +/- 0.008030358886349747

Plotting build-up factors with uncertainty

from io import StringIO
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

# After running Monte Carlo and building a table (as above)
thicknesses = np.linspace(5, 200, 100)
b_values = []
b_lo = []
b_hi = []

for t in thicknesses:
    bi = table.interpolate(thickness=float(t), warn=False)
    b_values.append(bi.value)
    b_lo.append(bi.value - bi.sigma)
    b_hi.append(bi.value + bi.sigma)

fig, ax = plt.subplots()
ax.plot(thicknesses, b_values, "b-", label="GP prediction")
ax.fill_between(thicknesses, b_lo, b_hi, color="blue", alpha=0.15, label="1-sigma")

# Monte Carlo points
mc_bs = [r.buildup["dose-AP"] for r in mc_results]
ax.plot(mc_thicknesses, mc_bs, "ko", markersize=7, label="Monte Carlo")

ax.set_xlabel("Thickness (cm)")
ax.set_ylabel("Build-up factor B")
ax.legend()

buf = StringIO()
fig.savefig(buf, format="svg")
plt.close(fig)
print(buf.getvalue())
2026-04-22T09:09:53.293694 image/svg+xml Matplotlib v3.10.8, https://matplotlib.org/

Plotting dose with uncertainty

from io import StringIO
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

doses = []
doses_lo = []
doses_hi = []

for t in all_thicknesses:
    layers = [rpk.Layer(thickness=t, material=concrete)]
    bi = table.interpolate(thickness=t, warn=False)
    pk = rpk.calculate_dose(
        source_strength=SOURCE,
        layers=layers,
        source=source,
        geometry="AP",
    )
    doses.append(pk.dose_rate * bi.value)
    doses_lo.append(pk.dose_rate * (bi.value - bi.sigma))
    doses_hi.append(pk.dose_rate * (bi.value + bi.sigma))

fig, ax = plt.subplots()
ax.plot(all_thicknesses, doses, "b-", label="PK with build-up")
ax.fill_between(all_thicknesses, doses_lo, doses_hi, color="blue", alpha=0.15)

# Monte Carlo reference points
mc_doses = [r.mc["dose-AP"] * SOURCE for r in mc_results]
ax.plot(mc_thicknesses, mc_doses, "ko", markersize=7, label="Monte Carlo")

ax.set_xlabel("Thickness (cm)")
ax.set_ylabel("Dose rate (Sv/hr)")
ax.set_yscale("log")
ax.legend()

buf = StringIO()
fig.savefig(buf, format="svg")
plt.close(fig)
print(buf.getvalue())
2026-04-22T09:09:53.369072 image/svg+xml Matplotlib v3.10.8, https://matplotlib.org/