Skip to content

Add correlated K/Kpri scoring and scatter-dominant mode to egs_kerma - #1442

Open
mainegra wants to merge 3 commits into
feature-egs_kerma-remove-rho_cvfrom
feature-egs_kerma-scatter-mode
Open

Add correlated K/Kpri scoring and scatter-dominant mode to egs_kerma#1442
mainegra wants to merge 3 commits into
feature-egs_kerma-remove-rho_cvfrom
feature-egs_kerma-scatter-mode

Conversation

@mainegra

Copy link
Copy Markdown
Contributor

Summary

  • Adds correlated $K/K_\text{pri}$ output with statistically correct uncertainty,
    via a per-region cross-term array scg_kp[j][ir] $= \sum_i K_i \cdot q_i$
    accumulated once per history in endHistory.
  • Adds a second scoring mode (scatter correction = no) that scores
    $K_\text{scat}$ (kerma from latch $\neq$ 0 photons) instead of $K_\text{pri}$,
    and derives $K/K_\text{pri} = K/(K - K_\text{scat})$.
    Default (scatter correction = yes) scores $K_\text{pri}$ (latch $=$ 0
    photons) directly.

Physics

Both modes report $R = K / K_\text{pri}$. They differ in what secondary
quantity $q$ is scored and how $K_\text{pri}$ is obtained from it.

Scatter-correction mode (scatter correction = yes, default)

$q = K_\text{pri}$ is scored directly (latch $=$ 0 photons).
$K_\text{pri}$ is read straight from the scoring array, so

$$R = \frac{K}{q}$$

Applying the delta method to $R = K/q$:

$$\left(\frac{\sigma_R}{R}\right)^2 = \delta_K^2 + \delta_q^2 - 2 \cdot r_c$$

where $\delta_X = \sigma_X / \bar X$ is the relative standard deviation
and $r_c = cov(K,q)/(\bar K \bar q)$ is the
relative covariance (not the Pearson coefficient; unbounded).

Scatter-dominant mode (scatter correction = no)

$q = K_\text{scat}$ is scored (latch $\neq$ 0 photons).
$K_\text{pri} = K - K_\text{scat}$, so

$$R = \frac{K}{K - q}$$

Applying the delta method to $R = K/(K-q)$:

$$\left(\frac{\sigma_R}{R}\right)^2 = (R-1)^2 \left(\delta_K^2 + \delta_q^2 - 2 \cdot r_c\right)$$

where $r_c = cov(K, K_\text{scat})/(\bar K \bar K_\text{scat})$.

The $(R-1)^2$ prefactor makes this mode dramatically less efficient at
depth: at 30 mfp in iron $R \approx 86$, so even near-perfect
correlation between $K$ and $K_\text{scat}$ is insufficient to overcome
the $\approx 85\times$ amplification of any residual decorrelation.
Benchmark runs show scatter-dominant mode is ${\sim}5.7\times$ less
efficient per history than scatter-correction mode at 30 mfp.

Note on $r_c$

$r_c$ is the relative covariance $cov(K,q)/(\bar K \bar q)$,
estimated from the cross-term array as

$$r_c = \frac{\overline{Kq} - \bar K \bar q}{\bar K \bar q} = \frac{\sum_i K_i q_i / N - \bar K \bar q}{\bar K \bar q}$$

It is not the Pearson correlation coefficient and is not bounded to
$[-1,1]$. Do not clamp it.

Input

:start scoring options:

    score primaries    = yes
    scatter correction = no    # default yes; set no for Kscat mode
    ...
:stop scoring options:

Implementation notes

  • scg_kp is allocated only when score primaries = yes; fully
    serialized in outputData / readData / addState / resetCounter
    for correct parallel job combination.
  • kerma_p stores either $K_\text{pri}$ or $K_\text{scat}$ depending
    on is_scatter_correction; the latch condition at scoring time is
    !latch (correction mode) or latch != 0 (scatter mode).

mainegra and others added 2 commits July 30, 2026 09:29
  Introduce is_scatter_correction flag (default true = score Kpri) and
  input key "scatter correction = no" to switch to scatter-dominant mode
  (score Kscat instead of Kpri), computing K/Kpri = K/(K-Kscat).

  Both modes accumulate the cross-term Sum(K*q) per region in scg_kp[][],
  enabling a correlated ratio estimate that eliminates the large positive
  covariance between K and Kpri/Kscat. The correlated uncertainty formula
  is σ_BUF/BUF = sqrt(δ_K² + δ_q² - 2·rc) for scatter-correction mode,
  and |BUF-1|·sqrt(δ_K² + δ_S² - 2·rc) for scatter-dominant mode (delta
  method applied to BUF = 1/(1-ρ) where ρ = Kscat/K).

  scg_kp is fully serialized in outputData/readData/addState/resetCounter
  for correct parallel job combination. Column header "K/Kpri" replaces
  "BUF" in all output strings; code variable names BUF/dBUF are unchanged.

  Co-Authored-By: Claude Sonnet 4.6
rc = cov(K,q)/(K*q) — it is unbounded and should never be clamped
to [-1,1].  Add a comment at both rc sites so a future reader does
not "fix" it by treating it as a Pearson coefficient.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
@mainegra
mainegra requested a review from a team as a code owner July 30, 2026 14:41
@mainegra
mainegra requested review from ftessier and rtownson and removed request for a team July 30, 2026 14:41
Completes propagation of 409e460 to the top of the stack.  This is the
branch the cluster builds from, so results produced before this merge
carry the splitting-rate-dependent bias and are not comparable with
results produced after it.

Verified after the fix: the map dependence at 15 mfp drops from 12.8
sigma to 1.7 sigma, and both maps agree with two independent
calculations that use no importance sampling (30.671 and 30.217 against
30.304 and 30.328).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant