Skip to content

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

@grzanka

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.

Metadata

Metadata

Assignees

No one assigned

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions