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.

The new package boundary is dpt.transport. It consumes a declared cell grid, partial interaction coefficients and source states, then returns one complete detector contribution per original history. A separate estimator turns those contributions into image means and uncertainty. This division keeps the transport law inspectable before a tally or inverse objective compresses its outputs. The deterministic primary operators remain useful as independent limiting cases and as the implementation used when uncollided transmission is the intended model.

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.

TransportSpec declares the implemented subset of this model. Absorption terminates a history under the default analogue estimator, while continuous weighting integrates its survival factor instead (§9.6). Scattering uses either the explicitly idealised isotropic-elastic law or unpolarised free-electron-compton. Neither option supplies coherent form factors, bound-electron corrections or fluorescence. The caller provides compatible macroscopic absorption and scattering coefficients. Their sum defines total extinction, and linear interpolation between the declared energy nodes is part of this transport model.

@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.

There is a second layout boundary to observe: this package’s PlanarDetector.shape is (width, height), with an xy origin at the lower pixel faces and flat index row * width + column. The deterministic DetectorGeometry uses (height, width) and a first-pixel-centre origin. The transport detector is an axis-aligned plane above the grid. Passing the earlier detector metadata through unchanged would therefore alter both pixel placement and indexing, even if all the arrays still had the expected length.

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.

The library uses residual optical depth. In walk, one CUDA thread carries position, direction, energy and the current draw through the whole history. The inner loop consumes that draw across cells, and only a completed collision advances the event identity used for the next draw. An arbitrary subdivision of a material region therefore does not introduce another random flight sample.

@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

The collision branch advances by residual optical depth divided by local extinction. The boundary branch subtracts the optical thickness actually traversed. enter_grid handles initial cell ownership. cross_faces crosses every exactly tied face and assigns its coordinate from the grid face itself. Both are shared device functions in the same canonical file. There is no accumulated epsilon displacement to walk the photon away from a boundary. Vacuum consumes distance without consuming optical depth. Binary64 coordinates and path arithmetic retain these distinctions across anisotropic cells and long trajectories.

The same walk also accumulates density_score for a selected active material. For now it is diagnostic state with a specific derivative purpose: Chapter 10 derives why a collision adds one and each traversed optical thickness subtracts from it. Forward tracing and derivative replay will therefore use the same flight and termination decisions, with no separate implementation of the physics to drift out of agreement.

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.

The walk engine uses explicit cell crossings. Null-collision tracking provides the alternative derived here.

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

Rejection trials occupy their own counter domains. A rejected candidate cannot consume the next collision’s free-flight draw, and exhausting the trial limit returns a failure status. Acceptance arithmetic uses binary64. This routine samples a conditional outcome after scattering has been selected, and it does not supply a material’s macroscopic scattering coefficient or turn an arbitrary coefficient table into a free-electron model. That compatibility remains the caller’s physical assumption.

Bound-electron scattering modifies this account through scattering functions, shell effects and Doppler broadening. Coherent scattering changes direction with negligible energy change in the elastic approximation, but its angular distribution depends on atomic form factors. An isotropic coherent sampler is an illustrative mathematical model, not a validated Rayleigh model for diagnostic imaging. Combining tabulated interaction probabilities with a free-electron final-state approximation is a deliberate approximation whose consequences require checking. The Geant4 manual distinguishes its free-electron treatment from a model with atomic-shell effects. NIST likewise describes coherent form factors and incoherent scattering functions in its tabulations. A total attenuation value cannot replace those differential models. [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.

Both supported transport estimators write at most one detector hit per history. A convex cell grid followed by its external detector permits scoring the final escaped straight ray once. The analogue count score is the source importance weight times source_amplitude, while the ideal energy score also multiplies by the arrival energy in keV. Continuous absorption weighting, introduced in §9.6, includes the absorption survival weight along the realised path. Both use the same original source normalisation. An energy-kev score remains incident photon energy, with no electronics calibration or local dose interpretation. The total tally already includes uncollided survivors, so it must not be added to a separately calculated primary image.

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. Its caller must supply IID source histories. A repeated deterministic source ray also qualifies because only the conditional transport draws are random. Deterministically stratified sources require independent complete-batch replications for their uncertainty calculation. This is a precondition on the sampling law, not a property the function can establish from the realised scores. There is no boolean argument that certifies 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()

The implementation first finds a magnitude scale for each pixel, sums normalised scores, then accumulates centred squared deviations. This avoids forming a potentially overflowing raw second moment and subtracting two nearly equal large quantities. Misses and absorption contribute zero but remain in batch.count, and the finishing pass includes their centred deviations even though they have no detector address. Averaging only the hits would estimate a conditional response to detection and give the wrong transmission.

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

Persistent moment scratch scales with detector pixels, while history outputs scale with launched histories, and neither contains a history-by-pixel matrix. The two output images give marginal variances of the pixel means. They do not make pixels independent, and summing those variances would not produce the uncertainty of a general image loss. Chapter 10 will estimate that scalar uncertainty across independent complete replicates.

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.

The existing cell traversal performs both integrals. A scattering event interrupts the current segment, so only its travelled portion contributes absorption. Crossing a cell face carries the unused scattering optical depth into the next material, and the final segment to escape contributes absorption too. A cell with zero scattering therefore needs neither a fictitious collision nor division by a zero rate. Its absorption integral still counts. The continuous_absorption branch in Listing 9.2 implements these choices inside 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 the existing one-history, one-terminal-pixel accounting and 28-byte forward output record. It also keeps walking after absorption has made a contribution tiny, so some scattering histories cost more than their analogue counterparts. There is no weight cutoff or roulette concealed in this mode. Event or crossing exhaustion invalidates the estimate. Whether the changed sampling measure buys useful accuracy per second is an execution question, not a consequence of having written “variance reduction” above the code.

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.

Our implementation uses the history-owned organisation. prepare_transport binds resident material IDs and coefficient arrays to one CUDA stream. The caller supplies resident source positions, unit directions, importance weights and positive density scales, and owns every output destination. trace_histories validates the call and launches the history kernel. It never allocates a queue or transfers a source ray inside the launch path.

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()

Positions and directions use binary64 vectors, coefficients and scores use binary64, and pixel indices and terminal statuses use integers. The five forward outputs require 28 bytes per history, independent of its event count. Thread-local state grows with the chosen history algorithm rather than with an allocated record of every event. This avoids an event tape, but variable history lengths still cause divergence and can raise 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 active state resident on the GPU across a batch. Allocate working buffers before the hot path, and make overflow a reported failure or a correct continuation mechanism. A hard event limit can protect a run from hanging, but reaching it invalidates an otherwise unqualified estimate unless a justified residual treatment exists. Reporting only the histories that finished is selection bias.

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.

experiments/transport-validation/run.py checks the engine against an independently specified slab expectation and records the source and configuration. More detailed CUDA tests check layered survival, Compton angle–energy laws, history identity and tally moments. The sparse output records final energy and terminal status. It does not include the separate absorption and recoil totals needed for the energy ledger in equation (9.14). The null-collision and splitting checks in Table 9.2 apply when those samplers are used.

Only after these checks should a performance report compare history rate, time to a specified uncertainty, memory consumption and tally behaviour on a named GPU. A fast run with a biased termination rule has measured the cost of a different problem.

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.