Following scattered photons with Monte Carlo transport

A photon removed from a primary ray has not necessarily disappeared. It may turn, lose some energy and arrive at a different detector pixel. The straight-ray calculation correctly excludes it from the primary signal, but the detector may be less selective. To predict that contribution, we need paths whose next direction is decided inside the object.

We will construct those paths as samples of a stated transport model. The random history is a computational device for estimating an expected detector signal. It is not the physical exposure itself, and its sampling noise is not the detector noise described in Chapter 8. Keeping those two uncertainties separate will matter even more when the estimates become gradients.

Turn the ray integrator into a photon walker by giving each history a position, direction, energy and random identity. It must advance to a collision or material boundary, update the state after scattering, and score an accepted detector crossing. Build that walk before the image reduction: one complete contribution per source history is enough to estimate both the mean image and its sampling uncertainty. The functions in dpt.transport will supply these operations.

9.1 Beyond the uncollided photon

A photon state contains physical position $\mathbf{x}$, unit direction $\boldsymbol{\omega}$ and energy $E$. Between interactions, it travels along a straight line. At an interaction it can be absorbed or scattered into a new direction and energy. The object is stationary, photon histories do not interact with one another, and the response considered here is linear in the incident photon population.

For a reduced diagnostic-energy model, distinguish photoelectric absorption, coherent scattering and incoherent scattering. We can consider incident energies between 10 and 150 keV, but the source interval does not bound all subsequent energies: Compton scattering lowers photon energy. The coefficient tables and samplers must cover every tracked state, or the calculation must declare a low-energy termination rule and measure its effect. This range excludes pair production. It does not justify treating bound electrons as free at every scattering angle.

In the reduced model, photoelectric absorption terminates the photon. Fluorescence, atomic relaxation, electron transport and secondary bremsstrahlung are omitted. Such a model can study photon-image formation under these approximations, but it is not a dose calculation or a complete account of photons re-emitted after absorption. Add omitted processes only with matching cross sections, final-state laws and validation.

Let $\psi(\mathbf{x},\boldsymbol{\omega},E)$ be the exposure-integrated angular photon fluence, in photons per $\mathrm{mm}^2$, per steradian and per keV. Let $q$ be an internal source density in photons per $\mathrm{mm}^3$, per steradian and per keV. For each scattering process $a$, let $\mu_a(\mathbf{x},E)$ be its macroscopic coefficient and $p_a(\boldsymbol{\omega},E\mid\boldsymbol{\omega}',E')$ its normalised outgoing-state probability kernel. The stationary linear transport equation is

In equation (9.1), streaming changes fluence along a direction, total attenuation removes photons from that state, and the collision integral adds photons scattered from other states. The integration over direction is over solid angle. A kernel can contain a delta distribution in energy, as for elastic scattering, and normalisation is then understood as a probability measure. Spatial dependence of the kernel is suppressed in the notation but follows the local material.

The total coefficient is $\mu_t=\mu_{\mathrm{abs}}+\sum_a\mu_a$. Source illumination can be supplied through incoming boundary conditions instead of an internal $q$. Specify one consistent source representation. A source already launched at an entrance surface should not also be counted as a volume source.

Setting the collision source to zero leaves the uncollided solution along each characteristic, whose attenuation is the exponential developed in Chapter 2. Successive scattering events correspond to successive collision-source contributions. Monte Carlo samples those event sequences without explicitly storing a field over every position, direction and energy. The detector score is a linear functional of the resulting transport solution.

Choose the interaction law and absorption estimator through TransportSpec. The default analogue walk terminates on absorption. Continuous weighting, derived in §9.6, integrates absorption survival along scattering paths. The angular choice is idealised isotropic-elastic or unpolarised free-electron-compton, paired with compatible macroscopic coefficients. Their sum gives total extinction, and the walk interpolates those coefficients linearly between energy nodes.

@dataclass(frozen=True, slots=True)
class TransportSpec:
    """A complete physical scope and finite launch-resource policy.

    Source rays and nonnegative importance weights are supplied by the caller.
    Their sampling law must be fixed with respect to active material parameters;
    the library does not silently invent an emission distribution. Provenance
    identifies the user's coefficient source, not a fabricated material asset.
    """

    grid: MaterialGrid
    detector: PlanarDetector
    energy_kev: float
    coefficient_provenance: str
    coefficient_energies_kev: tuple[float, ...] = ()
    scattering_law: Literal["isotropic-elastic", "free-electron-compton"] = "isotropic-elastic"
    scoring: Literal["photon-count", "energy-kev"] = "photon-count"
    max_angle_trials: int = 4096
    max_events: int = 4096
    max_crossings: int = 65536
    block_dim: Literal[64, 128, 256] = 128
    estimator: Literal["analogue", "continuous-absorption"] = "analogue"

    def __post_init__(self) -> None:
        _freeze_fields(self, ("coefficient_energies_kev",))
        _physical(self.energy_kev, "source energy", positive=True)
        for node in self.coefficient_energies_kev:
            _physical(node, "coefficient energy", nonnegative=True)
        if type(self.coefficient_provenance) is not str or not self.coefficient_provenance.strip():
            raise TransportError("the supplied coefficients require a provenance identifier")
        if self.detector.z_mm <= self.grid.upper_mm[2]:
            raise TransportError("the detector plane must lie above the entire material grid")
        nodes = self.energy_nodes
        if any(right <= left for left, right in pairwise(nodes)):
            raise TransportError("coefficient energies must be strictly increasing")
        if not nodes[0] <= self.energy_kev <= nodes[-1]:
            raise TransportError("coefficient energies must contain the source energy")
        if self.scattering_law not in ("isotropic-elastic", "free-electron-compton"):
            raise TransportError("unsupported scattering law")
        if self.scattering_law == "free-electron-compton" and len(nodes) < 2:
            raise TransportError("Compton scattering needs energy-dependent coefficient tables")
        if self.scoring not in ("photon-count", "energy-kev"):
            raise TransportError("unsupported detector score")
        if self.estimator not in ("analogue", "continuous-absorption"):
            raise TransportError("unsupported transport estimator")
        _positive_integer(self.max_angle_trials, "max_angle_trials")
        _positive_integer(self.max_events, "max_events")
        _positive_integer(self.max_crossings, "max_crossings")
        if self.block_dim not in (64, 128, 256):
            raise TransportError("block_dim must be 64, 128 or 256")

    @property
    def energy_nodes(self) -> tuple[float, ...]:
        """An omitted grid declares coefficients at the one elastic source energy."""
        return self.coefficient_energies_kev or (self.energy_kev,)

MaterialGrid stores material IDs in (nz, ny, nx) order, with x contiguous, but its origin is the lower cell face and its fields are axis-aligned and piecewise constant. Chapter 4’s grid instead places its origin at the first sample centre and interpolates a continuous interior field. Converting between them changes the field model and requires an explicit decision about cell locations and coefficients.

Place PlanarDetector on an axis-aligned plane above the grid. It uses (width, height), an xy origin at the lower pixel faces, and flat index row * width + column. Convert the earlier DetectorGeometry explicitly: its shape is (height, width) and its origin is the first pixel centre. These choices determine which pixel owns a crossing.

The event, crossing and angular-trial limits protect finite resources. Reaching a limit invalidates the estimate, but it does not turn the remaining photon into an absorption event. Leaving coefficient energy support also fails explicitly. These choices keep a numerical resource guard from becoming an undocumented physical cutoff.

Figure 9.1 shows how these event sequences connect to detector contributions.

9.2 Free flights and interaction sampling

For a photon at fixed energy traversing a homogeneous region, the probability of no interaction over distance $s$ is $e^{-\mu_t s}$. Its next interaction distance $S$ has an exponential distribution when $\mu_t>0$:

Equation (9.2) comes from inverting the survival probability. The mean free path is $1/\mu_t$, in millimetres. A vacuum has $\mu_t=0$ and no finite collision distance, so traverse it to the next geometric boundary. Random-number endpoints require explicit treatment so that an accidental zero does not create an unintended infinite logarithm.

At a real collision, choose process $a$ according to its share of the total coefficient:

The probabilities in equation (9.3), including absorption, must sum to one. The same process definitions must appear in both numerator and denominator. Using a total coefficient that excludes coherent scattering while allowing a coherent event branch produces a different transport law. XCOM explicitly provides totals both with and without coherent scattering, alongside process-specific coefficients. Selecting the appropriate total is therefore part of importing the physics data. [9]

A heterogeneous region is sampled in optical depth. Draw $\tau=-\log U$ and locate the point at which the accumulated total coefficient reaches that draw:

For piecewise-constant cells, equation (9.4) can be solved by marching to boundaries and subtracting each cell’s optical thickness from the remaining $\tau$. If the next cell contributes more than the remaining draw, the collision lies inside it at remaining optical depth divided by its coefficient. Otherwise cross the boundary and continue. Preserve the finite domain and terminate or score at its declared surfaces before considering material beyond them. Geant4’s interaction-length description uses the same exponential draw in accumulated mean-free-path units and subtracts the consumed amount across steps. This supplies a useful independent specification of the flight law. [26]

Do not carry a physical distance sampled in one material unchanged into another. Carrying the remaining optical depth has the correct survival law. Resampling at a boundary can also be correct under exponential memorylessness with suitable conditioning, but mixing these two constructions without that argument invites double sampling or missed attenuation.

Implement walk with one CUDA thread carrying the photon state and residual optical depth through a complete history. Its inner loop consumes the current draw across cells. Advance the event identity only after a collision, so subdividing one material into more cells leaves the sampled flight unchanged.

@wp.func
def walk(
    initial_position: wp.vec3d,
    initial_direction: wp.vec3d,
    source_weight: wp.float64,
    source_amplitude: wp.float64,
    seed: wp.uint64,
    history: wp.uint64,
    active_material: int,
    material_ids: wp.array(dtype=wp.int32),
    energies: wp.array(dtype=wp.float64),
    absorption: wp.array(dtype=wp.float64),
    scattering: wp.array(dtype=wp.float64),
    density: wp.array(dtype=wp.float64),
    p: Parameters,
) -> Trace:
    result = Trace()
    result.pixel = MISS_PIXEL
    result.events = 0
    result.status = 0
    result.energy = p.source_energy
    result.score = wp.float64(0.0)
    result.density_score = wp.float64(0.0)
    result.absorption_depth = wp.float64(0.0)
    position = initial_position
    direction = initial_direction
    entry = enter_grid(position, direction, p)
    position = entry.position
    cell = entry.cell
    alive = entry.alive
    result.status = entry.status
    crossings = int(0)
    while alive != 0:
        if result.events >= p.max_events:
            result.status = 2
            break
        if result.energy < energies[0] or result.energy > energies[p.energy_bins - 1]:
            result.status = 5
            break
        draws = random4(seed, history, wp.uint32(result.events), wp.uint32(0))
        residual = -wp.log(draws[0])
        collision = int(0)
        material = int(0)
        interaction = wp.vec2d(wp.float64(0.0))
        while alive != 0 and collision == 0:
            flat = (cell[2] * p.shape[1] + cell[1]) * p.shape[0] + cell[0]
            material = material_ids[flat]
            interaction = coefficients(material, result.energy, energies, absorption, scattering, p)
            extinction = density[material] * (interaction[0] + interaction[1])
            absorption_rate = wp.float64(0.0)
            if p.continuous_absorption != 0:
                extinction = density[material] * interaction[1]
                absorption_rate = density[material] * interaction[0]
            face_distance = distances_to_faces(position, direction, cell, p)
            distance = wp.min(face_distance[0], wp.min(face_distance[1], face_distance[2]))
            if distance < wp.float64(0.0) or not wp.isfinite(distance):
                result.status = 4
                alive = 0
                break
            optical_distance = extinction * distance
            if (
                not wp.isfinite(extinction)
                or not wp.isfinite(optical_distance)
                or not wp.isfinite(absorption_rate)
            ):
                result.status = 4
                alive = 0
                break
            absorption_distance = wp.float64(0.0)
            if p.continuous_absorption != 0:
                segment_distance = distance
                if extinction > wp.float64(0.0) and residual < optical_distance:
                    segment_distance = residual / extinction
                # Use the realised segment, not the entire distance to the face:
                # scattering can interrupt this segment before the face is reached.
                absorption_distance = absorption_rate * segment_distance
                result.absorption_depth += absorption_distance
                if not wp.isfinite(result.absorption_depth):
                    result.status = 4
                    alive = 0
                    break
            if extinction > wp.float64(0.0) and residual < optical_distance:
                distance = residual / extinction
                position += distance * direction
                if material == active_material:
                    result.density_score += wp.float64(1.0) - residual
                    if p.continuous_absorption != 0:
                        result.density_score -= absorption_distance
                collision = 1
            else:
                if material == active_material:
                    result.density_score -= optical_distance
                    if p.continuous_absorption != 0:
                        result.density_score -= absorption_distance
                residual -= optical_distance
                position += distance * direction
                crossing = cross_faces(position, direction, cell, face_distance, distance, p)
                position = crossing.position
                cell = crossing.cell
                alive = crossing.inside
                crossings += 1
                if alive != 0 and crossings >= p.max_crossings:
                    result.status = 3
                    alive = 0
        if collision != 0:
            result.events += 1
            absorption_probability = wp.float64(0.0)
            if p.continuous_absorption == 0:
                absorption_probability = interaction[0] / (interaction[0] + interaction[1])
            if p.continuous_absorption == 0 and draws[1] < absorption_probability:
                result.status = 1
                alive = 0
            else:
                cosine = wp.float64(2.0) * draws[2] - wp.float64(1.0)
                azimuth = TWO_PI * draws[3]
                if p.compton != 0:
                    angular = compton_scatter(
                        result.energy,
                        seed,
                        history,
                        wp.uint32(result.events - 1),
                        p.max_angle_trials,
                    )
                    cosine = angular[0]
                    azimuth = angular[1]
                    result.energy = angular[2]
                    if angular[3] != wp.float64(0.0):
                        result.status = 6
                        alive = 0
                if alive != 0:
                    direction = scatter_direction(direction, cosine, azimuth)
    if result.status == 0:
        result.pixel = detector_pixel(position, direction, p)
        if result.pixel == INVALID_PIXEL:
            result.pixel = MISS_PIXEL
            result.status = 4
        if result.pixel >= 0:
            score_energy = wp.float64(1.0)
            if p.energy_score != 0:
                score_energy = result.energy
            result.score = product3(source_weight, source_amplitude, score_energy)
            if p.continuous_absorption != 0:
                result.score = attenuated_product4(
                    result.absorption_depth,
                    source_weight,
                    source_amplitude,
                    score_energy,
                    wp.float64(1.0),
                )
    if not wp.isfinite(result.score) or not wp.isfinite(result.density_score):
        result.status = 4
    return result

In the collision branch, advance by residual optical depth divided by local extinction. At a boundary, subtract the optical thickness just traversed. Use enter_grid to assign the first cell and cross_faces to cross every exactly tied face, setting coordinates from the grid faces themselves. This avoids accumulating artificial epsilon displacements. Vacuum consumes distance without consuming optical depth. Binary64 coordinates and path arithmetic help retain these distinctions across anisotropic cells and long trajectories.

The density_score accumulator will be useful in Chapter 10: for a selected material, each collision adds one and each traversed optical thickness subtracts from it. Accumulating it alongside the flight gives the derivative calculation the same segments and termination decisions as the forward walk.

A second approach uses a majorant $\bar\mu\geq\mu_t$ throughout a region. Sample candidate events at rate $\bar\mu$ and accept a real collision at position $\mathbf{x}$ with probability $\mu_t(\mathbf{x},E)/\bar\mu$. Rejected candidates are null events: position advances but energy, direction and physical weight do not change. In a homogeneous segment,

The thinning identity in equation (9.5) explains why null events preserve the target law. A loose majorant wastes work, whereas an underestimated majorant invalidates the acceptance probability. Clamping a ratio above one hides the invalid bound and biases the process. With regional or energy-dependent majorants, handle their boundaries and update the rate consistently.

Use explicit cell crossings for the piecewise-constant grid here, where each segment’s optical depth is available directly. Null-collision tracking is an alternative way to sample the same flight law when a valid majorant is available.

Figure 9.2 shows the different quantities advanced by the two flight samplers.

9.3 Scattering angle and energy change

Collision occurrence and collision outcome are separate distributions. Total coefficients determine when a process occurs. Differential laws determine the direction and energy after that process. A table of total cross sections alone cannot supply an angular sampler.

For an unpolarised photon scattering from a free electron initially at rest, let $\vartheta$ be the scattering angle and $E'$ the outgoing photon energy. Energy and momentum conservation give

Equation (9.6) is the free-electron Compton relation. Forward scattering retains the incident energy, while backscatter gives $E/(1+2E/(m_ec^2))$. The missing photon energy becomes electron kinetic energy in this model. It has not thereby been shown to deposit locally at the collision point. Geant4’s Compton final-state account gives this photon energy relation and treats the recoil electron separately. [27]

Writing $r_e$ for the classical electron radius and $k=E'/E$, the Klein–Nishina differential cross section per solid angle for that free-electron model is

To turn equation (9.7) into a sampler, normalise it over solid angle. Sampling the polar angle uniformly would omit both its physical weighting and the solid-angle measure. Use $u=\cos\vartheta$ as the polar integration variable, with $-1\leq u\leq1$, so solid-angle integration becomes integration over $\varphi$ and $u$. The azimuth is uniform for the stated unpolarised model. Rejection sampling must use a verified envelope and retain rejected draws as rejected draws, not accepted events with a convenient angle.

compton_scatter proposes a uniform cosine and accepts it with the Klein–Nishina shape divided by an envelope of two, after removing the constant cross-section prefactor. The shape is bounded above by $k^3+k\leq2$ for $0<k\leq1$, so the acceptance probability cannot exceed one under the stated model. The accepted cosine determines the outgoing energy through equation (9.6), and azimuth remains uniform.

@wp.func
def compton_scatter(
    energy: wp.float64,
    seed: wp.uint64,
    history: wp.uint64,
    event: wp.uint32,
    max_trials: int,
) -> wp.vec4d:
    """Cosine, azimuth, outgoing energy and status from the conditional free-electron law."""
    result = wp.vec4d(wp.float64(0.0), wp.float64(0.0), energy, wp.float64(6.0))
    for trial in range(max_trials):
        angular = random4(seed, history, event, wp.uint32(RANDOM_NAMESPACE_BIT) + wp.uint32(trial))
        cosine = wp.float64(2.0) * angular[0] - wp.float64(1.0)
        ratio = wp.float64(1.0) / (
            wp.float64(1.0) + energy / ELECTRON_REST_ENERGY_KEV * (wp.float64(1.0) - cosine)
        )
        # Klein-Nishina density divided by its envelope 2. This conditional
        # scattering law is independent of the active material-density scale.
        acceptance = wp.float64(0.5) * (
            ratio * ratio * ratio + ratio - ratio * ratio * (wp.float64(1.0) - cosine * cosine)
        )
        if angular[1] < acceptance:
            result = wp.vec4d(
                cosine,
                TWO_PI * angular[2],
                energy * ratio,
                wp.float64(0.0),
            )
            break
    return result

Give rejection trials their own random-counter domain so that angular retries leave the next collision’s free-flight draw untouched. Evaluate acceptance in binary64 and return a failure if the trial limit is exhausted.

Geant4 distinguishes free-electron Compton scattering from models with atomic-shell effects. Its broader models account for details such as electron binding and momentum that our Klein–Nishina sampler leaves out. Coherent scattering also needs an angular law based on atomic form factors, as described in NIST’s tabulations. An isotropic elastic law is useful for isolating the walker, while a physical Rayleigh model requires those form factors. Keep total interaction rates and conditional outcome laws as separate functions so that either can be replaced with a compatible model. [27, 8]

After sampling an angle, construct the outgoing direction relative to the incoming one. Choose a stable orthonormal pair $\mathbf{e}_1,\mathbf{e}_2$ perpendicular to $\boldsymbol{\omega}$, with $\mathbf{e}_2=\boldsymbol{\omega}\times\mathbf{e}_1$. Then

Equation (9.8) preserves unit length in exact arithmetic. Select a helper axis that is not nearly parallel to the incoming direction before forming the transverse basis. A special case that works for a beam along world $z$ is not a general rotation algorithm. Check norm and angle after the transform, including directions close to every coordinate axis.

The source and scattering model must also agree on polarisation assumptions. Adding polarisation later changes the state and the azimuthal law. It cannot be represented by relabelling an unpolarised history.

9.4 Sources, boundaries and detector scores

The straight-ray model defined $n_{0,p}$ at each detector pixel. A history calculation instead starts from a source phase-space distribution over emission position, direction and energy. Its normalisation must reproduce the same open-beam detector population when scattering and object attenuation are removed. Geometric acceptance enters through the source-to-detector trajectories, so another inverse-square factor at detector scoring would count it again.

Let $z$ denote the source state and let $F(z)$ be its non-negative physical intensity density over the exposure, with total expected emitted population $N_{\mathrm{src}}$. Sample $z$ from a normalised proposal $g(z)$ that is positive wherever the physical source can contribute. Let $S_p$ be the detector score from a history sampled according to the physical interaction laws. Then

The ratio in equation (9.9) is a source importance weight. If $g=F/N_{\mathrm{src}}$, it is simply $N_{\mathrm{src}}$. A point source or line spectrum is naturally a mixed or discrete measure, so define $F/g$ on that measure rather than dividing incompatible densities. The source integral includes its full angular acceptance, not just histories that happen to arrive at the detector.

For an absorbing detector surface, score its first accepted crossing and terminate the photon. A basic counting score is one for a hit in pixel $p$ and zero otherwise. A deterministic expected-response score uses $R_p(E)$ at the crossing. Sampling a random per-photon detector output is possible, but unnecessary if only the expected image is required and the conditional response mean is known. Averaging that response analytically can reduce computational variance.

A tally from sampled surface crossings counts photons directly. A tally formed by integrating angular fluence over a surface instead includes the projected-area factor $|\boldsymbol{\omega}\cdot\mathbf{n}|$. Mixing these two scoring descriptions adds or removes a cosine incorrectly. State whether the detector is one-sided, how its edges assign pixels, and whether a grazing or out-of-footprint crossing is accepted.

A photon leaving the object can still travel through vacuum to the detector. Leaving one object’s support is therefore not necessarily history termination. Leaving the outer simulation domain is termination only when no future contribution from outside that domain is part of the model. The detector and external domain must be placed accordingly.

For linear transport, primary and scattered scores can be separated by an interaction counter. A detector hit with zero physical scattering events belongs to the primary tally, and null events do not change that classification. Their sum must equal the total tally under the same source and detector rules.

Score the final escaped ray once for the convex grid and external absorbing detector used here. Multiply source importance weight by source_amplitude for a count score, and by arrival energy as well for an ideal incident-energy score in keV. Continuous weighting adds absorption survival along the path. Retain the same source normalisation in both estimators. The resulting image contains primary and scattered arrivals together.

9.5 Estimators, variance and uncertainty

For independent source histories indexed by $h$, collect each complete weighted contribution as $X_{hp}$. This includes all descendants of that source history if splitting or secondary production is used. The expected image and the estimated variance of its Monte Carlo estimate are

Equation (9.10) estimates the variance of the sample mean, not the variance of a physical exposure. It assumes independent, identically distributed history contributions with finite variance. Histories in a fixed stratified allocation require the corresponding stratum-wise variance calculation. Correlated descendants must be combined before estimating uncertainty at the primary-history level.

The usual standard error decreases as $N^{-1/2}$ when these assumptions hold. That scaling says nothing about bias from omitted physics, a wrong detector response, coarse geometry or a biased sampler. Increasing history count can make agreement with the wrong answer exceedingly precise.

A pure-absorption slab provides an exact sampling example. Launch identical normal rays through thickness $d$ with coefficient $\mu$, and let the ideal detector accept every survivor. Each unweighted history score is Bernoulli with success probability $T=e^{-\mu d}$:

Equation (9.11) shows why a strongly attenuated pixel needs many histories under direct survival sampling. The relative uncertainty depends on the expected number of detected histories, $NT$. Zero recorded hits do not establish a zero expected signal. A sample variance of zero in that case is also not evidence of certainty.

Figure 9.3 shows how rare surviving histories control relative precision.

A confidence interval based on a Gaussian approximation needs enough effective contributing histories and a suitable contribution distribution. Rare large importance weights can make the apparent convergence misleading. Use independent replications, inspect the distribution of history weights and report how intervals were constructed. For image-wide conclusions, pixelwise intervals do not automatically provide simultaneous coverage over every pixel.

Keep simulation count $N$ separate from physical emitted population $N_{\mathrm{src}}$. The former controls numerical precision, while the latter scales the expected exposure. Halving Monte Carlo standard error costs roughly four times as many histories under fixed variance, whereas doubling a physical exposure changes the detector’s expected counts and shot noise. They answer different questions.

history_moments implements equation (9.10) from sparse (pixel, score, status) records and requires IID source histories. Repeated launches of one deterministic source ray qualify when their conditional transport draws are independent. Deterministically stratified sources instead need independent complete-batch replications to estimate uncertainty. The caller must establish the sampling law. Realised scores alone cannot verify it.

def history_moments(
    pixel: Any,
    score: Any,
    history_status: Any,
    *,
    batch: HistoryBatch,
    workspace: EstimatorWorkspace,
    out_mean: Any,
    out_variance_of_mean: Any,
    stream: Any = None,
    validate: bool = True,
) -> None:
    """Mean and unbiased variance of that mean, for IID original-history scores.

    Source histories must be independently identically distributed; a repeated
    deterministic ray also qualifies. This sampling-law precondition cannot be
    verified from the realised score arrays.
    Deterministically stratified sources require variance across independent
    *replicated complete batches*, not this within-history formula. The two
    output images contain marginal variances; they do not claim independent
    detector pixels or supply a covariance-free loss-error bar.

    FP64 scale, mean and centred-deviation passes preserve the result range. Their
    floating summation order is not bitwise reproducible. This analogue engine
    produces at most one contribution per original history; descendant splitting
    must first combine contributions before using these moment operations.
    """
    if batch.count < 2:
        raise TransportError("variance of a mean requires at least two original histories")
    _accumulate_history_mean(
        pixel,
        score,
        history_status,
        batch=batch,
        workspace=workspace,
        out_mean=out_mean,
        out_variance_of_mean=out_variance_of_mean,
        stream=stream,
        validate=validate,
    )
    transport = workspace.transport
    pixels = transport.spec.detector.pixels
    transport._launch(
        transport._kernels.centred_moments,
        batch.count,
        [
            pixel,
            score,
            out_mean,
            workspace._scale,
            out_variance_of_mean,
            workspace._hits,
            transport._status,
        ],
    )
    transport._launch(
        transport._kernels.finish_variance,
        pixels,
        [
            out_mean,
            workspace._scale,
            workspace._hits,
            batch.count,
            out_variance_of_mean,
            transport._status,
        ],
    )
    if validate:
        transport.check_status()

Build the moment reduction in three passes: find a magnitude scale per pixel, sum normalised scores, then accumulate centred squared deviations. This avoids overflow in raw second moments and cancellation between nearly equal large sums. Include misses and absorption as zero contributions in the original-history count and in the centred deviations. Averaging only detected histories would give a conditional response instead of transmission.

history_mean shares the validation and scaled-mean passes without requesting a variance.

Allocate moment scratch per detector pixel and sparse output per history. The result gives marginal pixel-mean variances without a history-by-pixel matrix. To estimate uncertainty in an image loss, reduce each complete independent replicate to that scalar first, preserving the effects of pixel correlations. Chapter 10 uses this construction for optimisation.

9.6 Variance reduction without changing the answer

Importance sampling replaces a path distribution with a proposal that visits useful contributions more often. If a path has target density $f(H)$, proposal density $g(H)$ and score $S(H)$, its weighted estimator satisfies

Equation (9.12) requires proposal support wherever the target integrand contributes. The weight includes every altered sampling decision: source emission, flight, process choice and scattering outcome as applicable. It is not enough to correct only the final detector direction. Whether the weighted estimator has finite, smaller variance depends on the proposal and score.

Conditional expectation offers another reduction. In the pure-absorption slab, the expected survival score is known exactly as $e^{-\mu d}$. Replacing a random survival decision by that conditional expectation removes its Bernoulli noise. In scattering transport, a next-event detector estimate can similarly integrate a selected connection to the detector while retaining the random path to the scattering site. It must include the correct angular density, detector acceptance and transmittance, and avoid double counting the same contribution through another estimator.

For absorption at collisions, an implicit-capture scheme can retain only scattering continuations while reducing weight. If $\mu_s$ is the sum of scattering coefficients and $w$ the incoming statistical weight, continue with weight $w\mu_s/\mu_t$ and sample the scattering process conditionally within $\mu_s$. This preserves the expected continuation under the stated linear model. If $\mu_s=0$, no continuation remains. Weight removed as absorption is an expectation tally, not a sampled local dose.

Low-weight paths can be terminated with Russian roulette. Choose a survival probability $q$ with $0<q\leq1$, and replace weight $w$ by

The conditional identity in equation (9.13) preserves expected contributions if continuation sampling is correct. Roulette can increase variance while reducing work. Its benefit is measured at fixed computational cost. Deterministically killing every history below a weight threshold lacks the compensating survivor weight and generally biases the image.

Splitting a history into $m$ descendants of weight $w/m$ also preserves total weight. The descendants share an ancestor even when their subsequent paths are sampled independently. Combine their scores before estimating uncertainty across source histories. A buffer that silently drops excess descendants changes the estimator regardless of how sophisticated the splitting rule was.

For an estimated quantity, compare methods using uncertainty at a measured runtime and the same physical model. A method that traces fewer paths can still be slower because each path costs more, or less accurate because a few huge weights dominate. Check weight accounting before comparing runtimes.

Integrating absorption along the path

Consider what the analogue estimator asks the GPU to do in the absorption-only limit: flip enough survival coins to rediscover the exponential we derived in Chapter 2. That is a rather expensive way of expressing uncertainty about an answer we already know conditionally on the ray. TransportSpec.estimator="continuous-absorption" removes those absorption decisions from the walk. The default remains "analogue", so the two estimators can be compared under the same physical inputs.

The distinction between stored coefficients and physical rates matters here. The device arrays contain density-independent coefficients. Multiplying by the material’s density scale gives $\kappa_a=\rho_m\mu_{a,m}(E)$ and $\kappa_s=\rho_m\mu_{s,m}(E)$. The weighted walk samples scattering flights at rate $\kappa_s$ and accumulates absorption depth $\tau_a=\int_H\kappa_a\,\mathrm ds$. Its terminal score includes $\exp(-\tau_a)$. Hayakawa, Spanier and Venugopalan derive this continuous absorption construction and distinguish it from the collision-weighting scheme above. Their later erratum corrects the numerical error and efficiency values in the paper. [39, 40]

For a detector-reaching path with specified scattering events, analogue sampling supplies a scattering-rate factor at each event and the flight-survival factor $\exp[-\int_H(\kappa_s+\kappa_a)\,\mathrm ds]$. The weighted walk retains the same scattering-rate and angular factors, but its flight survival contains only $\kappa_s$. Multiplying by $\exp(-\tau_a)$ restores the absorption part of the original survival factor. The path contribution therefore has the same weight in the expected detector signal, although the two algorithms sample it differently.

Extend the same cell traversal to accumulate absorption along each travelled segment. A scattering event ends that segment at the collision. A cell crossing carries the remaining scattering optical depth onward, and the final segment to escape contributes absorption too. A zero-scattering cell therefore needs only boundary traversal and the absorption integral. The continuous_absorption branch in Listing 9.2 puts these operations in the shared walk.

We retain absorption depth as a binary64 logarithmic weight until scoring. Forming exp(-tau) first would discard a contribution when that intermediate underflows, even if a large source amplitude makes the final score representable. attenuated_product4 combines the attenuation exponent with the original factors’ binary exponents before the final rounding. Forward values and density derivatives each receive that treatment, so neither inherits the other’s lost bits.

With a fixed pencil source and zero scattering throughout the coefficient table, every history follows the same attenuation integral. Preparation can establish that this expectation is deterministic. A randomly positioned source over different material columns still produces different integrals, and scattering leaves random paths even for a fixed source. Zero empirical variance alone establishes neither case.

Continuous weighting keeps one terminal detector contribution per history, but continues walking even when absorption has made its weight tiny. Some paths therefore cost more than their analogue counterparts. Compare uncertainty at equal runtime to decide whether removing absorption noise pays for the extra tracing. This mode uses neither a weight cutoff nor roulette.

9.7 Transport histories on the GPU

A history-owned kernel assigns a source history to a thread and follows it until termination. It keeps current position, direction, energy, weight, random state and event counters nearby, but neighbouring threads can take different numbers of events. Some escape immediately while another finds yet another scattering site. The warp waits for the paths it actually has, not the average history in a planning spreadsheet.

An event-based design instead stores active states in queues and launches kernels for flights, collisions and scoring. Compaction can reduce inactive lanes and group similar work, at the cost of queue traffic, launch overhead and additional state storage. Sorting by material or energy may improve access patterns but also costs work. Neither organisation has an intrinsic speed advantage independent of the scene and hardware.

Start with a history-owned kernel so the flight, collision and detector logic stay together. Prepare the resident material IDs and coefficient tables with prepare_transport. Allocate source positions, unit directions, importance weights, positive density scales and output destinations on the device, then supply them to trace_histories. It launches the walk on the prepared stream without moving source rays through the host.

def trace_histories(
    positions: Any,
    directions: Any,
    weights: Any,
    density: Any,
    *,
    batch: HistoryBatch,
    workspace: TransportWorkspace,
    out_pixel: Any,
    out_score: Any,
    out_energy: Any,
    out_events: Any,
    out_status: Any,
    source_amplitude: float = 1.0,
    stream: Any = None,
    validate: bool = True,
) -> None:
    """Write one sparse detector score per independent original history.

    Binary64 positions are in mm; directions must be unit vectors to 1e-12 in
    squared norm. Source positions lie below the detector plane. Nonnegative
    base weights encode the caller's fixed source importance sampling, while
    ``source_amplitude`` scales the expected measurement. ``density`` contains
    strictly positive dimensionless scales multiplying both partial coefficients.

    The score is per launched history; averaging includes misses and absorption
    as zeros. It is not normalised by detected photons. `out_pixel=-1` means no
    detector hit. Score, final energy (keV), collision count and status are always
    written, including on failure; invalid outputs must be discarded as a batch.
    Continuous absorption samples scattering flights and weights each detector
    hit by absorption along its entire path; its event count counts scatterings.
    The analogue default retains sampled absorption and total-extinction flights.
    """
    wp = workspace.context.wp
    _validate_call(
        workspace,
        batch,
        positions,
        directions,
        weights,
        density,
        [
            ("out_pixel", out_pixel, wp.int32),
            ("out_score", out_score, wp.float64),
            ("out_energy", out_energy, wp.float64),
            ("out_events", out_events, wp.int32),
            ("out_status", out_status, wp.int32),
        ],
        source_amplitude,
        stream,
        validate,
    )
    workspace._launch(
        workspace._kernels.trace_histories,
        batch.count,
        [
            positions,
            directions,
            weights,
            density,
            *workspace._model_inputs(),
            wp.uint64(batch.seed),
            wp.uint64(batch.first_history),
            source_amplitude,
            out_pixel,
            out_score,
            out_energy,
            out_events,
            out_status,
            workspace._status,
        ],
    )
    if validate:
        workspace.check_status()

Use binary64 positions, directions, coefficients and scores, with integer pixel indices and terminal statuses. Store the terminal output per history and keep evolving state in the thread. Memory then grows with history count rather than total event count, although long, divergent paths can still increase register demand.

Figure 9.4 shows the state movement introduced by the two execution organisations.

Randomness must be attached to logical histories, not accidental execution order. Record a master seed and a reproducible mapping from history identifier, descendant identifier and draw counter to random values. Reusing one stream after splitting can correlate descendants, and changing scheduling should not accidentally reuse a variate. A counter-based scheme can make replay easier, but its mapping and statistical properties still require checks.

The shared random4 implements Philox4x32-10 from Salmon and colleagues’ counter-based generator family. A 64-bit seed is its key, and the counter contains a 64-bit original-history identity, event and domain. Four returned uniforms supply the ordinary flight, process, cosine and azimuth draws. Source sampling reserves a separate event range, while Compton rejection and detector observations use separate domains. HistoryBatch prevents identity wrap and records the source-draw namespace as well as the transport range. [35]

@wp.func
def random4(
    seed: wp.uint64,
    identity: wp.uint64,
    event: wp.uint32,
    domain: wp.uint32,
) -> wp.vec4d:
    """Return four open-interval uniforms; each has 32 random mantissa bits."""
    c0 = wp.uint32(identity & wp.uint64(0xFFFFFFFF))
    c1 = wp.uint32(identity >> wp.uint64(32))
    c2 = event
    c3 = domain
    k0 = wp.uint32(seed & wp.uint64(0xFFFFFFFF))
    k1 = wp.uint32(seed >> wp.uint64(32))
    for _ in range(10):
        product0 = wp.uint64(0xD2511F53) * wp.uint64(c0)
        product1 = wp.uint64(0xCD9E8D57) * wp.uint64(c2)
        low0 = wp.uint32(product0 & wp.uint64(0xFFFFFFFF))
        low1 = wp.uint32(product1 & wp.uint64(0xFFFFFFFF))
        high0 = wp.uint32(product0 >> wp.uint64(32))
        high1 = wp.uint32(product1 >> wp.uint64(32))
        c0 = high1 ^ c1 ^ k0
        c1 = low1
        c2 = high0 ^ c3 ^ k1
        c3 = low0
        k0 = k0 + wp.uint32(0x9E3779B9)
        k1 = k1 + wp.uint32(0xBB67AE85)
    scale = wp.float64(1.0 / 4294967296.0)
    return wp.vec4d(
        (wp.float64(c0) + wp.float64(0.5)) * scale,
        (wp.float64(c1) + wp.float64(0.5)) * scale,
        (wp.float64(c2) + wp.float64(0.5)) * scale,
        (wp.float64(c3) + wp.float64(0.5)) * scale,
    )

Adding half a unit before scaling places each uniform strictly between zero and one. Each still contains only 32 random bits, so extreme-tail experiments must account for that finite resolution. Chunking a launch preserves counter addresses. Reusing them can reproduce a path on the same device and build with unchanged inputs, but it does not create another independent sample, nor does it promise identical floating-point reductions across schedules.

Tally accumulation needs an explicit strategy. Direct atomic additions are simple and can contend at bright detector pixels. Per-block or tiled partial tallies reduce some contention but consume memory and require a reduction. Floating-point addition is not associative, so changing execution order can change the last bits even when the estimator is statistically unchanged. Reproducibility claims must distinguish identical random histories from bitwise-identical reductions. Table 9.1 specifies the execution evidence needed to assess these implementation choices.

The moment kernels group participating warp lanes that target the same pixel. A native CUDA helper reduces their normalised binary64 contributions through a peer-mask shuffle tree, and one leader performs the global atomic update. Scale maxima and centred deviations use the same grouping. This reduces contention at concentrated pixels without allocating a detector histogram per block.

On SM 70 and newer targets, every live lane in each tally pass reaches a ballot before participating histories combine equal pixel addresses. Misses remain outside the participant mask. Matching uses that mask, and shuffles use the resulting same-pixel peer masks. Sparse participating masks and a final partial warp must obey the same shuffle participation rules, and misses still belong in the original-history denominator. Targets below SM 70 retain individual atomics. Grouping adds work when scores are already dispersed, so its usefulness depends on where histories land, not just how many are launched.

Keep cross-section tables and working buffers on the GPU across batches. For queues or bounded history loops, propagate overflow and exhaustion to the caller or provide a continuation that preserves the estimator. Include every launched history in the accounting.

Mixed precision needs a quantity-specific justification. Positions near thin boundaries, accumulated optical depths, small statistical weights and large image tallies have different error sensitivities. Validate coordinate units and tolerances at the actual scene scale. Use a CPU reference for independently specified formulas. Check race freedom and device precision in the CUDA execution, and measure occupancy and memory traffic with the profiler.

9.8 Validate before comparing speed

Begin with laws whose answers do not depend on another transport implementation. In a homogeneous medium, check the flight survival distribution and mean free path. Across layered cells, compare survival with the exponential of the summed optical thickness. The result must not depend on subdividing a constant material into additional cells.

For a single Compton event, check the forward and backward energy limits, the angle–energy relation and outgoing direction norm. Test the angular distribution against independently integrated probabilities over bins, including the solid-angle measure. An energy histogram alone can miss an azimuth or world-frame rotation error.

At the history level, track photon energy through the reduced model. With one initial photon and no fluorescence or secondary photons, the bookkeeping identity is

In equation (9.14), detector and escaped energy are mutually exclusive terminal contributions. Absorption records the remaining photon energy at termination, recoil sums energy transferred in Compton events, and cutoff records unresolved residual photon energy. This is an energy ledger, not a spatial dose estimate. For biased sampling, the corresponding weighted tallies obey conservation in expectation rather than necessarily in each realised path.

A cutoff tally makes discarded energy visible but does not prove that its detector contribution is negligible. Lower the cutoff and compare the scored image, or establish a bound under the stated model. Likewise, conservation alone cannot validate an angular distribution: a photon can conserve energy while being sent to the wrong pixel. The tests in Table 9.2 check the transport law and estimator beyond that conservation ledger.

For an independent-code comparison, first make the physical models equal enough that agreement is meaningful. Record coefficient sources, interpolation, interaction options, geometry boundaries, source normalisation, cutoffs and scoring definitions. A reference with fluorescence enabled and a candidate with terminal photoelectric absorption need not agree in a detector region sensitive to fluorescence. Treat that as a model difference to quantify, not a tolerance to relax until the test passes.

Independent replications assess the empirical spread of estimates and the coverage of their intervals. They should also reveal whether changing a seed changes the result by an amount consistent with the reported uncertainty. Do not demand every stochastic comparison land within a fixed two-standard-error band: occasional excursions are expected, and repeated testing needs its own acceptance design.

Begin with the analytic slab-survival calculation, then exercise layered survival, the Compton angle–energy law and detector tally moments separately. To check the complete energy balance in equation (9.14), retain absorption and recoil totals in addition to the terminal energy recorded by the sparse output. Apply the null-collision and splitting checks in Table 9.2 when implementing those sampling alternatives.

Compare implementations by time to a specified uncertainty, alongside history rate, memory use and tally behaviour. These quantities expose the cost of the estimator and the GPU organisation together.

The expected image is now an integral over random histories whose probabilities depend on the object. A derivative must account for those probabilities as well as for the score of each realised path. Chapter 10 takes up that distinction: differentiating only the arithmetic left in one sampled trace can miss part of the expected derivative.