Stability controls for the bfcorr BField domain walk

[RobMina]

Summary

Bringing up bfcorr=true CentralHelix→CRV extrapolation for Mu2e cosmics exposed four ways the BField
"domain walk" can fail: an OOM, a multi-minute CPU spin, a lost fit, and an aborted event. They share
one root cause — the domain step has no lower bound and the domain set no upper bound — compounded by
the fact that the CentralHelix parameterization is singular as |B|→0 (p_T ∝ B/ω; the dPardB ~ bfrac/(1+bfrac) correction has a pole at bfrac = −1). Near-vertical cosmics leaving the DS bore just
past the tracker drive the walk straight into that region.

The failure modes

#	Symptom	Mechanism
1	OOM — RSS 0 → 6.7 GB in ~20 s, then killed	Backward (CRV-bound) extrapolation of a near-vertical cosmic degenerates climbing into the DS bore: traj pieces shrink to ~4e-5 ns, reported momentum collapses unphysically (11.46 GeV → 11 MeV), rangeInTolerance returns a vanishing dt, and addDomain appends an effect + traj piece per step without bound.
2	Multi-minute CPU hang (flat memory)	A tight BCorrTolerance drives createDomains' trange→0, so tstart stops advancing and the loop spins emitting deduped domains into a std::set.
3	Fit lost — throw std::invalid_argument("Invalid domains")	extend builds extension domains in independent createDomains calls, each with its own ±0.5·trange margin; the blocks don't tile contiguously and createEffects requires abutting domains.
4	Whole event aborted — throw std::length_error("Empty PiecewiseTrajectory!")	A degenerate seed yields zero domains, so the seed trajectory comes out empty; the throw is uncaught and art drops the event, corrupting the output entry count.

New knobs

knob	default	effect at default	effect when set
Config::mindtstep_	0.0 ns	max(rangeInTolerance, 0) == rangeInTolerance	floor on the domain step; bounds the domain count where rangeInTolerance→0. Must not exceed the extrapolation's maxDtStep() (checked, throws)
Config::domainmargin_	DBL_MAX ns	unconfined walk, unchanged	confine domains — and the field sampling at their midpoints — to the active range ± margin
BFieldMap(double minfield)	0.0 T	no low-field protection	domains are only built where |B| ≥ minfield; extrapolation below it hands off to geometric free-particle continuation, before the dPardB pole

Always-on changes

These are not behind a knob. Each only fires where the previous code had undefined behavior or an
uncaught throw, so a well-behaved track is unaffected:

• extend grows existing domains with extendDomains instead of independent createDomains calls, so extension domains always abut (Update2 #3).
• detectorRange returns a null TimeRange when there are no active hits/xings, instead of an inverted range that throws.
• Seed conversion and replaceDomains soft-fail with a status (no domains, domains not spanning the trajectory, empty trajectory — Mat #4) instead of throwing or dereferencing an empty container. replaceDomains builds the replacement before touching domains_/effects_, so a failure leaves the track intact.
• A piece that can't join the trajectory (charge/mass mismatch, e.g. an omega sign flip) is reported as the new Status::unphysical, via ParticleTrajectory::compatible(), instead of throwing out of append.
• CentralHelix::constructible() (trivially true for LoopHelix/KinematicLine) refuses a null field or a circle centered on the axis before construction; resetBNom refuses a null BNom. Domains are never built with a null BNom.
• Extrapolation: a non-finite momentum or a caught error fails the extrapolation with a recorded status.

Verification

No effect when unconfigured. Mu2e Offline built twice on the same source base, differing only by this
branch versus its base 2011f34, identical input, no knob set, compared exactly (rtol = 0):

fit type	sample	result
LoopHelix (production OnSpill reco)	20,000 events, 3,349 tracks	IDENTICAL — 243 leaves
CentralHelix (cosmic reco)	10,000 events	IDENTICAL — 527 leaves
KinematicLine (extracted cosmics reco)	20,000 events	IDENTICAL — 478 leaves

Stability. Event 5160 under a 4 GB address-space cap, knobs off versus on, with nothing else changed:

	knobs off (= current head)	knobs on
peak RSS	3.2 GB (at the cap; ~6.7 GB uncapped)	1.0 GB
wall time	43.5 s (27% CPU — 672k page faults)	4.0 s
outcome	fails — FatalRootError finalizing a 367,596-domain runaway	completes, exit 0; degenerate track bounded and dropped cleanly

With the knobs set (mindtstep_ 0.03 ns, domainmargin_ 0, minfield 1e-6 T), bfcorr=true CentralHelix→CRV
reco runs over 112 Run-2 cosmic files (~98k events) with no crashes, hangs or domain errors.

The Mu2e Offline plumbing for these knobs will follow in a separate Offline PR.

[brownd1978]

Thanks for this PR. You've clearly understood a lot about how the fit works and how it can fail, and made important improvements.
These fixes make CreateDomains quite long, which suggests it should be refactored. I see big blocks that look very similar: maybe those need to be encapsulated in new functions or lambdas. I wonder if some protection could be moved to BFieldMap. Perhaps BFieldMap could define 'boundary' and '0-field' regions on construction, so that those wouldn't have to be test for by each track. It might help to add a createDomain function that returns a guaranteed usable domain, which could deal with things like minimum field.
Protections connected with sign flipping could maybe be moved to CentralHelix; that should allow omega to change sign while keeping the same particle charge, at least within the same piecetraj.
I suspect the config can be simplified; perhaps a min and max domain range are sufficient. Limits on physical ranges (time) have clearer interpretation than limits on container size.
I believe the minimum field value could go into the BFieldMap constructor: LH and KL won't ever care, so it's fine to set it to a single, global value.
I'm also curious if unusable inputs could be identified and failed before any field calculations, to avoid burdening downstream code with testing.
It's always good to catch failure conditions, but some of the throw/catch blocks seem designed to handle common/routine conditions, which might be clearer to implement as normal functions.
The term 'legacy' in comments won't mean much in a few years: referring to conditions that make the block relevant will be more useful.
I'm happy to discuss if any points above aren't clear.

[RobMina]

Thanks Dave! Sorry for the delay.

I think I've addressed most of your comments and concerns.

The one that I deferred is allowing omega to change sign, because that would require significant rewrites in CentralHelix (the derivatives dMdPar would flip sign incorrectly without adding extra logic). The sign flips I saw in testing were associated with a null bnom value passed into the constructor, so I added guards for that instead.

[brownd1978]

I believe this PR addresses 2 underlying failure cases when extrapolating CH tracks: extrapolating through the rapidly varying field inside the DS coil, and extrapolating poor fits with 'unphysical' interpretation. I believe the intent is to prevent these from causing an (uncaught) exception, which stops event processing, or from causing an 'infinite' loop.

I believe the min and max domain sizes address most of these problems.

I'd like to better understand the origin of the remaining problems, to see if there's a way to pre-identify them and abort the processing before getting into detailed algorithms.

In any event, the new low-level protections are still a good addition.

[brownd1978]

extend builds extension domains in independent createDomains calls, each with its own ±0.5·trange margin

New domains should always be created adjacent to existing domains. the low-level (private) createDomains function could be given a range which introduces a domain range gap; that would be a bug in the calling function.

[RobMina]

Hi Dave, thanks for the followup!

The createDomains calls that had caused noncontiguous domains were in extend. I replaced that block with single call to extendDomains. However, that did re-introduce another throw-catch to handle the out-of-map case, let me know if you want that part reworked.

In A/B testing I found that the production KL output was dropping ~25% of tracks relative to the pre-PR baseline because of out-of-map B field checks, which are not even needed for the KL fit. So I reworked that. The current state of the PR reproduces baseline outputs for all production configs.

I also reworded the comments to refer only to current behavior, rather than pre-PR.

Lastly, I reworked replaceDomains and extend to remove the assumption about the failure reason from the latter.

[brownd1978]

I just realized that many of my comments are still pending, so I'm releasing this review now.
Overall this is great work, thanks for digging out these error cases, I just have a few quibbles.

[brownd1978]

Thanks Rob, I realized that many of my GitHub comments were pending releasing a review, a lot of what I said yesterday must have been cryptic! I just released the review. As long as the logic in extend for replacing domains when the bfield config changes remains, updating it to use extendDomains when extending the range makes sense. As long as the catch-throws all end up with the fit status correct adding a layer is fine. Cheers, Dave

…

On Wed, Sep 16, 2026 at 10:03 AM Rob Mina ***@***.***> wrote: *RobMina* left a comment (KFTrack/KinKal#224) <#224 (comment)> Hi Dave, thanks for the followup! The createDomains calls that had caused noncontiguous domains were in extend. I replaced that block with single call to extendDomains. However, that did re-introduce another throw-catch to handle the out-of-map case, let me know if you want that part reworked. In A/B testing I found that the production KL output was dropping ~25% of tracks relative to the pre-PR baseline because of out-of-map B field checks, which are not even needed for the KL fit. So I reworked that. The current state of the PR reproduces baseline outputs for all production configs. I also reworded the comments to refer only to current behavior, rather than pre-PR. Lastly, I reworked replaceDomains and extend to remove the assumption about the failure reason from the latter. — Reply to this email directly, view it on GitHub <#224?email_source=notifications&email_token=ABAH5754XLB6Y6APDRPDFHT5PLBV5A5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTKNZQGEZTONZVGQ42M4TFMFZW63VHMNXW23LFNZ2KKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-5701377549>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABAH572XFMQF3LEGEPEOTCD5PLBV5AVCNFSNUABFKJSXA33TNF2G64TZHMZDIMJZHE2TENJYHNEXG43VMU5TKMBQGI3TOMBZGE22C5QC> . You are receiving this because you commented.Message ID: ***@***.***>

-- David Brown ***@***.*** Office Phone (510) 486-7261 Lawrence Berkeley National Lab M/S 50R5008 (50-6026C) Berkeley, CA 94720

[RobMina]

Hi Dave, thanks for the comments! I pushed a commit that I think addresses all of the concerns.

All production configs produce identical output. CentralHelix extrapolates to the CRV with bfcorr=true reliably whereas before it failed on most jobs.

[brownd1978]

Apart a minor quibble this now looks very good. Thanks for taking the time to go through my comments.

[oksuzian]

PR Review Summary — #224

Reviewed at head 134126b540c92a208c7622b8759ab33476ceec24. First pass, second-opinion review after @brownd1978's APPROVED review at the same head; his outstanding request (name the 1.0 in CentralHelix::constructible) is not repeated here.

Decision

• 🟡 comment only. Finding 1 is worth resolving before merge; the rest are minor.

Scope understood

• Bounds the bfcorr BField domain walk so a degenerate CentralHelix track cannot run the domain count away (OOM), spin on a vanishing step, or escape as an uncaught std::invalid_argument / std::length_error.
• Two new Config knobs (mindtstep_, domainmargin_), both no-ops at their defaults; low-field protection moved into the BFieldMap constructor (minfield_), so protecting() is off unless the experiment passes it.
• New soft-fail paths (Status::unphysical, validInput, hasTraj) plus a two-phase replaceDomains.

Findings

1. 🟠 [S1] extendDomains was left out of the hardening, and one of the new throws now escapes through it

   • Evidence: Fit/Track.hh:745-781. The two domains it builds sample the field with a bare bfield_.fieldVect(straj.position3(range.mid())) (lines 762 and 778) and take the step as std::max(rangeInTolerance(...), mindtstep_) with no positivity guard. createDomain (line 845) guards both cases — bfield_.usable() at 851, usableField() at 865, and if(dhi <= tstart) return std::nullopt at 853 — and updateDomains (line 740) keeps the previous field when the sample is unusable, so extendDomains is the one domain-producing path with neither guard. Its two call sites, iterate() at line 573 and processEnds() at line 793, are not inside a try, unlike the new try in extend().
     • Null BNom: Domain with bnom = (0,0,0) reaches addDomain → DomainWall::append → newpiece.resetBNom(next_->bnom()) (Fit/DomainWall.hh:109), which now throws std::invalid_argument("CentralHelix::resetBNom; null BNom") (Trajectory/CentralHelix.cc:93-95). Nothing between there and the caller catches it, so it leaves Track::fit() — the uncaught-throw failure mode #4 that this PR removes elsewhere, recreated by the new guard on a path that does not pre-check. initTraj has the same shape at line 628 (front.resetBNom(domain->bnom())), fed by domains that createDomain builds from the raw fieldVect fallback at line 867 whenever protecting() is false.
     • Reachability in Mu2e: KKBField::fieldVect returns exactly (0,0,0) for a point outside all maps (Offline/Mu2eKinKal/src/KKBField.cc:14-19), and KKBField::inRange tests only the inner maps' bounding boxes, so inRange can be true where the lookup fails. rangeInTolerance returns the whole piece length when d2pdt2 <= 1e-10 (General/BFieldMap.hh), so range.mid() can sit far outside the map even when time is inside it. Only CentralHelix has the new guard — LoopHelix::resetBNom (Trajectory/LoopHelix.cc:84) and KinematicLine do not — which is consistent with the identical LoopHelix/KinematicLine A/B output, and leaves the exposure on the CentralHelix cosmic path this PR targets.
     • Non-advancing step: rangeInTolerance returns sqrt(dp/d2pdt2) with dp = momentum(tstart)*tol, and the collapsing momentum described in failure mode #1 drives that to zero. With mindtstep_ at its default 0.0, dt == 0 gives TimeRange(time,time), which is a legal null range (General/TimeRange.hh:14-31), Domain::operator< dedupes it in the std::set, and time = domain.begin() never advances — failure mode #2, in the function extend() now routes through.
   • Impact: an event aborted by an uncaught exception, or a CPU spin, in the configuration Offline runs today (minfield_ unset, so protecting() is false and mindtstep_ is 0.0).
   • Suggested fix: build these two domains with createDomain as well, so extendDomains inherits the field check and the guaranteed-positive step; failing that, replicate both guards and decide what extendDomains should return when it cannot advance.

2. 🟡 [S2] std::clamp with an unordered pair of bounds

   • Evidence: Fit/Track.hh:966 and 1005, std::clamp(rangeInTolerance(...), config().mindtstep_, xtest.maxDtStep()). std::clamp requires !(hi < lo); nothing relates mindtstep_ (a KinKal Config value) to maxDtStep() (from the experiment's extrapolation predicate — Mu2e sets 1.0 ns for three of four and 10.0 ns for one, Offline/Mu2eKinKal/fcl/prolog.fcl:377,390,403,416). Any mindtstep_ above 1.0 ns — a plausible value for the failure mode #1 fix — is undefined behavior; libstdc++ returns lo, silently exceeding the extrapolation's own maximum step. The pre-PR line was std::min(rangeInTolerance(...), xtest.maxDtStep()), which could not violate a precondition.
   • Suggested fix: apply the two bounds in a defined order (std::max(std::min(dt, maxDtStep), std::min(mindtstep_, maxDtStep))), or validate mindtstep_ <= maxDtStep() once and say which bound wins.

3. 🟡 [S2] The two extrapolate branches report a caught exception differently

   • Evidence: in the protecting branch the catch at Fit/Track.hh:982 leaves retval = false; in the unprotected branch the identical catch at line 1016 is immediately overwritten by retval = true at line 1018. That line predates the PR, but the restructuring is what makes the two configurations disagree, and the return value is live in Offline: KKExtrap.hh:156 does if(!ktrk.extrapolate(tdir,toPerim)) return;, so once Offline passes a minfield the same out-of-map exception will abandon the rest of the extrapolation chain where today it continues.
   • Suggested fix: drop line 1018, or comment why a caught out-of-map exception is a success only for the unprotected walk.

4. 🟡 [S2] Non-finite momentum stops the extrapolation with no trace

   • Evidence: Fit/Track.hh:956, if( !std::isfinite(ktraj.momentum(time)) ) break; — no handed_off, no history_ entry, and no geometric continuation, so the fit stays usable() with a silently truncated extrapolation. Every other early exit in that branch records something. This is the same objection you raised on earlier revisions at Track.hh:324 and :377.
   • Suggested fix: record a Status before breaking, or set handed_off so the existing handoff bookkeeping covers it.

5. ⚪ [S3] Two small items

   • The PR description still documents minfield_, maxdomains_ and zerofield_extrap_ as new Config knobs; the final diff adds only mindtstep_ and domainmargin_, with minfield_ moved to the BFieldMap constructor and the other two dropped. Worth refreshing before merge, since the body is the merge record.
   • ParticleTrajectory::compatible returns true on an empty piece list, so prepend to an empty trajectory no longer throws from mass() (Trajectory/ParticleTrajectory.hh:44-45, "this will throw for empty"). That looks like an improvement, but it is a contract change for prepend that the PR does not mention.

Verified

• 🟢 CentralHelix::constructible faithfully mirrors the constructor it stands in for. The rotation matches setTransforms() (CentralHelix.cc:134) term for term, and momToRad, radius, amsign, phirm and lcent match CentralHelix.cc:48-63; the two refusals correspond to the constructor's own throws — null BNom at line 39 and rcent < 1.0 at line 65. The ParticleState constructor delegates to the pos/mom constructor (line 114), so the predicate covers the constructor actually reached. It is therefore not stricter than the constructor and cannot refuse a seed that previously fitted. The duplication and the unnamed 1.0 remain as you noted.
• 🟢 Status::unphysical is appended last, so no existing enumerator changes value; usable() (status_ > unfit && status_ < lowNDOF) classifies it unusable, and the statusName switch covers it. Mu2e does not persist the KinKal enum — KKFit.hh:687,699 maps it into TrkFitFlag — so there is no backward-compatibility question for existing files.
• 🟢 replaceDomains is behavior-preserving on the success path: the copy is built and committed before domains_/effects_ are touched, and the ordering of updateReference and the traj swap is unchanged. The added guard matches the setRange failure it describes. It does cost one deep copy of the trajectory per call; that is on the domain-replacement path only, not per iteration, so it is unlikely to matter.

Validation check

• Build/tests run: no. Nothing was compiled or run for this review; the PR's own A/B evidence (identical output for LoopHelix, CentralHelix cosmics and KalSeed domain bounds; 3.2 GB/43.5 s → 1.0 GB/4.0 s on event 5160) was read, not reproduced.
• Config contract check: partial. The new knobs are no-ops at their defaults, and protecting() is off until an experiment passes minfield to the BFieldMap constructor.
• Cross-repo consistency: needs follow-up. No open Mu2e PR sets minfield, mindtstep_ or domainmargin_, so in Offline as it stands the low-field handoff branch is unreachable and the unprotected paths — including the one in finding 1 — are the ones that run.

Residual risk

• Finding 1 is a code-path argument, not a reproduced crash; it needs a domain whose midpoint samples a null field, or a zero step, on the extendDomains path.
• The bfcorr=true CentralHelix→CRV configuration this PR is for is not exercised by the A/B samples, which ran production configs with the knobs off.

Author follow-ups

1. Give extendDomains the same field and positive-step guards as createDomain, ideally by calling it (finding 1).
2. Make the std::clamp bounds well-ordered, or validate mindtstep_ <= maxDtStep() (finding 2).
3. Decide whether a caught extrapolation exception should fail the extrapolation in both branches (finding 3), and record the non-finite-momentum stop in history_ (finding 4).
4. Refresh the PR description to the knobs that survived, and say which Offline PR will pass minfield (finding 5).

[RobMina]

Thanks for the review! I agree with all the suggestions. All are now implemented.

[brownd1978]

Thanks for all your work on this Rob!