Skip to content

Fix dedx_get_inverse_stp() branch selection for monotone STP curves; add Bragg-peak STP tool #121

Open

Description

@grzanka

Part of #118 (dedx_extra.c migration) · Phase C. This is the highest-value item — a real bug that affects every consumer, currently worked around only in dedx_web.

Bug

dedx_get_inverse_stp() in src/dedx_tools.c uses find_min() to locate the Bragg-peak energy. For a monotonically decreasing stopping-power curve (e.g. proton on water, ICRU49), find_min() returns -1, which becomes x2 = -1 in the bisection — producing negative intermediate energies and a spurious dedx_get_stp() error instead of a valid inverse.

Separately, the side semantics are inconsistent: the documented API is side = 0 (low-energy branch) / side = 1 (high-energy branch), but the implementation maps side < 0 → low branch, so neither 0 nor 1 ever selects the ascending branch — both currently return the descending-branch result.

The web project already worked around both in wasm/dedx_extra.c (dedx_get_inverse_stp_flat) with a robust strategy:

  1. Sample STP on a log-spaced energy grid to locate the peak.
  2. If the peak is at the leftmost sample → monotone curve → bisect the single descending branch over the full range.
  3. If the peak is interior → side == 0 selects [emin, e_peak] (ascending), side == 1 selects [e_peak, emax] (descending), with a guard for requested STP below STP(emin).

Proposed change

  1. Port the sampled-peak + branch-aware bisection into the core dedx_get_inverse_stp() so all consumers (CLI, Python, WASM) get correct results, and fix the side ∈ {0,1} mapping to match the documented contract.
  2. Add a Bragg-peak STP tool to dedx_tools.h (the web dedx_get_bragg_peak_stp samples ~300 log-spaced energies and returns the max), e.g.:
    double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err);

Acceptance criteria

  • Regression test: proton/water ICRU49 (monotone) returns a valid inverse for a STP between STP(emax) and STP(emin) — no error, no negative energy.
  • Test that side = 0 and side = 1 return the ascending vs descending branch respectively for a curve with an interior Bragg peak.
  • dedx_get_bragg_peak_stp() declared in dedx_tools.h + test.
  • Docstrings state units (MeV cm²/g, MeV/nucl).

Reference

dedx_web wasm/dedx_extra.c dedx_get_inverse_stp_flat (lines ~264–396) — the proven implementation to port.


Filed via Claude Code as part of the dedx_extra.c → libdedx migration plan.

Activity

  1. grzanka commented on Jul 24, 2026

    @grzanka
    ContributorAuthor

    Fixed in #148.

    While reproducing this, the bug turned out to be more severe than described above: find_min()'s hardcoded x ∈ [0.01, 10] search domain plus the broken side mapping meant dedx_get_inverse_stp() failed (err=203) for essentially every ordinary proton/water query I tried (PSTAR and ICRU49, both sides) — not just the monotone-curve case. #148 ports the sampled-peak + branch-aware bisection from dedx_web's dedx_get_inverse_stp_flat as proposed, adds dedx_get_bragg_peak_stp(), and fixes an adjacent bug where the function reloaded an already-loaded config unconditionally (corrupting the workspace), unlike dedx_get_csda() which already guards on config->loaded.

  2. grzanka commented on Jul 24, 2026

    @grzanka
    ContributorAuthor

    One naming adjustment to the spec above: the new function in #148 is dedx_get_max_stp(ws, config, err), not dedx_get_bragg_peak_stp() as proposed here (and as named in dedx_web's reference implementation).

    Reasoning: "Bragg peak" properly refers to the maximum of a depth-dose curve — a quantity that also depends on range straggling and isn't something this library computes. What this function actually returns is the maximum of the stopping-power-vs-energy curve, which is a related but distinct quantity. dedx_get_max_stp() names it precisely instead of borrowing dosimetry terminology; the header docstring calls out the distinction explicitly so nobody assumes it accounts for straggling.

    wasm/dedx_extra.c is slated for removal once this migration completes (per #118), so it'll adopt the new name rather than the other way around.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions