Calculate build-up with MC

The point-kernel method computes uncollided flux: the fraction of particles that travel through the shield without scattering. In reality, scattered particles also contribute to dose behind the shield. The build-up factor B corrects for this:

Corrected = Point-kernel * B

compute_buildup runs OpenMC Monte Carlo simulations on a list of geometries and returns the ratio of Monte Carlo result to point-kernel result for each.

Note

This requires OpenMC. See Installation for instructions.

Basic example

Compute the photon dose build-up factor for 10 cm of iron:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
layers = [rpk.Layer(thickness=10, material=iron)]

source = rpk.Source(particle="photon", energy=1e6)
results = rpk.compute_buildup(
    geometries=[layers],
    source=source,
    quantities=["dose-AP-photon"],
)

r = results[0]
print(f"MC dose:     {r.mc['dose-AP-photon']} (per source particle)")
print(f"PK dose:     {r.pk['dose-AP-photon']}")
print(f"Build-up B:  {r.buildup['dose-AP-photon']}")

MC dose:     1.642166554122132e-16 (per source particle)
PK dose:     3.2606076549874897e-17
Build-up B:  5.036381950493926

Multiple thicknesses example

Pass a list of layer lists to run Monte Carlo on several shield configurations in one call:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
thicknesses = [5, 10, 15, 20]
geometries = [[rpk.Layer(thickness=t, material=iron)] for t in thicknesses]

source = rpk.Source(particle="photon", energy=1e6)
results = rpk.compute_buildup(
    geometries=geometries,
    source=source,
    quantities=["dose-AP-photon"],
)

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

 5 cm: B = 2.733900428287289
10 cm: B = 5.03638195049393
15 cm: B = 7.657670772820534
20 cm: B = 10.883381763549531

Quantity strings

The quantities argument specifies what to tally. The particle is part of the quantity name so a result's keys never depend on which Source produced them:

Quantity string	Description
`"flux-neutron"`	Neutron flux (neutron source)
`"flux-photon"`	Photon flux (photon source)
`"dose-AP-neutron"`	Neutron dose, AP geometry (neutron source)
`"dose-AP-photon"`	Photon dose, AP geometry (photon source)
`"dose-AP-coupled-photon"`	Secondary photon dose, AP geometry (neutron source only)
`"dose-AP-total"`	Neutron + secondary-photon dose at AP (neutron source only)

The same -{particle} / -coupled-photon suffix rules apply to all six dose geometries (AP, PA, RLAT, LLAT, ROT, ISO).

You can request multiple quantities in a single call:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
layers = [rpk.Layer(thickness=10, material=iron)]

source = rpk.Source(particle="neutron", energy=14.1e6)
results = rpk.compute_buildup(
    geometries=[layers],
    source=source,
    quantities=["dose-AP-neutron", "flux-neutron"],
)

r = results[0]
print(f"Dose B: {r.buildup['dose-AP-neutron']}")
print(f"Flux B: {r.buildup['flux-neutron']}")

Dose B: 2.596163329059408
Flux B: 3.9010053493146293

BuildupResult

Each element in the returned list is a BuildupResult with these dictionaries, keyed by quantity string:

r.mc - Monte Carlo tally value (per source particle)
r.mc_std_dev - Monte Carlo standard deviation
r.pk - point-kernel reference value
r.buildup - build-up factor (mc / pk)

It also has r.optical_thickness (the total optical thickness of the geometry).

Applying the build-up factor

Pass BuildupResult directly

The simplest approach. The calculate_* functions auto-detect the correct quantity key:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
layers = [rpk.Layer(thickness=10, material=iron)]
PARTICLES_PER_HOUR = 1e12 * 3600  # 1e12 photons/sec activity

source = rpk.Source(particle="photon", energy=1e6)
results = rpk.compute_buildup(
    geometries=[layers],
    source=source,
    quantities=["dose-AP-photon"],
)
r = results[0]

corrected = rpk.calculate_dose(
    layers=layers,
    source=source,
    geometry="AP",
    buildup=r,
).scale(strength=PARTICLES_PER_HOUR)
print(f"Dose with build-up: {corrected.dose} Sv/hr")

Dose with build-up: 0.5911799594839677 Sv/hr

Use BuildupModel.constant() manually

If you have a known build-up factor from another source:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
layers = [rpk.Layer(thickness=10, material=iron)]

# Apply a fixed buildup factor of 2.5 - this multiplies the PK dose
# by 2.5 regardless of thickness. Useful when you have B from an
# external source (literature, previous calculation, etc.)
buildup = rpk.BuildupModel.constant(2.5)

source = rpk.Source(particle="photon", energy=1e6)
result = rpk.calculate_dose(
    layers=layers,
    source=source,
    geometry="AP",
    buildup=buildup,
).scale(strength=1e12 * 3600)  # photons/hour, dose in Sv/hr
print(f"Dose with B=2.5: {result.dose} Sv/hr")

Dose with B=2.5: 0.2934546889488741 Sv/hr

Cross sections path

If the OPENMC_CROSS_SECTIONS environment variable is not set, pass the path directly:

source = rpk.Source(particle="photon", energy=1e6)
results = rpk.compute_buildup(
    geometries=[layers],
    source=source,
    quantities=["dose-AP-photon"],
    cross_sections="/path/to/cross_sections.xml",
)

Secondary photons (coupled transport)

When using a neutron source, nuclear reactions in the shield produce secondary gamma rays. To include these, add "dose-AP-coupled-photon" to your quantities. The Monte Carlo runs with coupled neutron-photon transport automatically:

import rad_point_kernel as rpk

iron = rpk.Material(composition={"Fe": 1.0}, density=7.874)
layers = [rpk.Layer(thickness=10, material=iron)]

source = rpk.Source(particle="neutron", energy=14.1e6)
results = rpk.compute_buildup(
    geometries=[layers],
    source=source,
    quantities=["dose-AP-total"],
)

r = results[0].scale(strength=1e16)  # neutrons per shot, dose in Sv/shot
print(f"Neutron dose:    {r.mc['dose-AP-neutron']} Sv/shot")
print(f"Secondary gamma: {r.mc['dose-AP-coupled-photon']} Sv/shot")
print(f"Total dose:      {r.mc['dose-AP-total']} Sv/shot")

Neutron dose:    2915.568278678489 Sv/shot
Secondary gamma: 26.382481529374573 Sv/shot
Total dose:      2941.950760207864 Sv/shot

When both "dose-AP-neutron" and "dose-AP-coupled-photon" are requested for the same geometry, a synthetic "dose-AP-total" quantity is added automatically (sum of the two; standard deviation combined in quadrature; buildup factor referenced to the neutron PK). The same rule applies to PA, RLAT, LLAT, ROT, and ISO.

Since you're already running Monte Carlo for neutron build-up, coupled transport adds some extra compute cost and gives the secondary photon dose.

Simulation parameters

compute_buildup accepts optional parameters to control the Monte Carlo simulation:

Parameter	Default	Description
`particles_per_batch`	10,000	Particles per batch
`max_batches`	100	Maximum number of batches
`trigger_rel_err`	0.05	Target relative error on tallies