From cd84ae325b60c8dded6764cef39ba39af7109697 Mon Sep 17 00:00:00 2001 From: robmina Date: Tue, 30 Jun 2026 11:42:01 -0500 Subject: [PATCH 1/6] Fit: domain step-floor, field-break, and empty-seed guard for bfcorr cosmics Robustness fixes to the bfcorr=true CentralHelix domain machinery for near-vertical cosmic tracks that exit the DS bore (|B| -> 0): - Config: add mindtstep_ (min domain step), domainmargin_ (clamp domains to the active range +/- margin), and minfield_ (stop extrapolation before the dPardB pole). Defaults preserve legacy behaviour (margin = max(), minfield = 0). - createDomains/extendDomains: floor the domain step at mindtstep_ so a vanishing rangeInTolerance in a high-gradient/low-momentum region can no longer spawn an unbounded number of micro-domains (CPU hang / OOM). Clamp domain bounds to the active range +/- domainmargin_ so a uniform-field track just past the tracker doesn't sample Bz->0 at a domain midpoint and blow the CentralHelix reference off by metres via the singular bfrac/(1+bfrac) correction. - createEffects: connect only domains that actually abut, instead of throwing "Invalid domains" across a gap between separate contiguous domain blocks. - extrapolate: break before |B| < minfield_, tested at the new domain midpoint (not the frontier), handing the field-free region to a straight line tail. - convertSeed: guard against an empty fittraj_ when a degenerate active range yields zero domains -- createEffects would otherwise build a Measurement against an empty PiecewiseTrajectory and throw std::length_error, aborting the whole art event. Treat it as an unfittable track (outsidemap) so the module drops it cleanly. Includes temporary domain-count / empty-seed diagnostics (marked "remove before PR"). Co-Authored-By: Claude Opus 4.8 --- Fit/Config.cc | 3 ++ Fit/Config.hh | 12 ++++++++ Fit/Track.hh | 83 +++++++++++++++++++++++++++++++++++++++++++-------- 3 files changed, 85 insertions(+), 13 deletions(-) diff --git a/Fit/Config.cc b/Fit/Config.cc index 4b3370ca..d564f24f 100644 --- a/Fit/Config.cc +++ b/Fit/Config.cc @@ -8,6 +8,9 @@ namespace KinKal { << " diverge dpar chisq " << kkconfig.pdchisq_ << " diverge traj gap (mm) " << kkconfig.divgap_ << " fractional momentum tolerance " << kkconfig.tol_ + << " min domain step (ns) " << kkconfig.mindtstep_ + << " min field (T) " << kkconfig.minfield_ + << " domain margin (ns) " << kkconfig.domainmargin_ << " min NDOF " << kkconfig.minndof_ << " BField correction " << kkconfig.bfcorr_ << " with " << kkconfig.schedule().size() diff --git a/Fit/Config.hh b/Fit/Config.hh index a73ff62a..ebcf81d2 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -29,6 +30,17 @@ namespace KinKal { double pdchisq_ = 1.0e6; // maximum allowed parameter change (units of chisqred) WRT previous reference double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm) double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps + double mindtstep_ = 1.0e-3; // ns: hard floor on the BField domain step, bounding the domain count in + // high-gradient/low-momentum regions where rangeInTolerance -> ~0 (else the + // bfcorr domain walk takes ~MaxDt/dt micro-steps -> OOM / time-budget truncation) + double minfield_ = 0.0; // T: if >0, stop the bfcorr extrapolation once |B| drops below this, so the + // ill-conditioned CentralHelix is never driven into B->0 (the physical runaway origin) + double domainmargin_ = std::numeric_limits::max(); // ns: max time a BField fit domain may extend + // beyond the active (hit) range. Default = unclamped (legacy). Set small (e.g. 0) + // to confine the fit's domains to the measurement region, so it never samples the + // field outside the solenoid bore -- where, for tracks that leave the magnet just + // past the tracker (cosmics), Bz->0 makes the CentralHelix dPardB ~ 1/(1+dB/B) + // correction singular and destroys the fit. Affects only createDomains/extendDomains. unsigned minndof_ = 5; // minimum number of DOFs to continue fit bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit bool ends_ = true; // process the passive effects at each end of the track after schedule completion diff --git a/Fit/Track.hh b/Fit/Track.hh index aeb7d83e..c43b20db 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -56,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -391,6 +392,18 @@ namespace KinKal { newpiece.range() = domain->range(); fittraj_->append(newpiece); } + // A degenerate active range (an inverted/zero-width detectorRange for a pathological seed) can leave + // createDomains with zero domains, so the loop above appends nothing and fittraj_ is empty. createEffects + // would then build a Measurement against this empty PiecewiseTrajectory and throw + // std::length_error("Empty PiecewiseTrajectory!"), which -- being neither std::invalid_argument nor + // caught -- aborts the whole art event (and corrupts the output -> TTree::SetEntries). Treat it as an + // unfittable track: outsidemap stops the fit cleanly (needsFit()==false skips createEffects) and the + // module drops it via goodFit()==false. fittraj_ is left as a valid (empty) object, never null. + if(fittraj_->pieces().empty()){ + std::cout << "CONVERTSEED EMPTY range=[" << range.begin() << "," << range.end() << "] ndomains=" << domains.size() << std::endl; // DIAG - remove before PR + history_.emplace_back(0,0,Status::outsidemap, "Empty seed trajectory (no domains)"); + return; + } } else { // use the middle of the range as the nominal BField for this fit: double tref = range.mid(); @@ -418,8 +431,12 @@ namespace KinKal { auto prevdom = nextdom; ++nextdom; while( nextdom != domains.cend() ){ - if(fabs(prevdom->get()->end()-nextdom->get()->begin())>1e-10)throw std::invalid_argument("Invalid domains"); - effects_.emplace_back(std::make_unique(*prevdom,*nextdom ,*fittraj_)); + // only connect domains that actually abut. A gap means prevdom and nextdom belong to different + // contiguous blocks (e.g. separate low/high fit extensions either side of the existing core domains); + // bridging across that gap would create a spurious DomainWall spanning the whole core, so skip it + // (this previously threw "Invalid domains" and aborted the fit). + if(fabs(prevdom->get()->end()-nextdom->get()->begin())<=1e-10) + effects_.emplace_back(std::make_unique(*prevdom,*nextdom ,*fittraj_)); prevdom = nextdom; ++nextdom; } @@ -461,6 +478,8 @@ namespace KinKal { status().comment_ = status().comment_ + error.what(); } } + // CALIBRATION PROBE – remove before PR + std::cout << "KinKal::Track::fit ndomains=" << domains_.size() << " status=" << fitStatus().status_ << std::endl; if(config().plevel_ > Config::none)print(std::cout, config().plevel_); } @@ -666,9 +685,12 @@ namespace KinKal { double time = drange.begin(); while(time > fitrange.begin()){ auto const& ktraj = fittraj_->nearestPiece(time); - double dt = bfield_.rangeInTolerance(ktraj,time,config().tol_); - TimeRange range(time-dt,time); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); + double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); + // clamp the domain low bound to the active range minus domainmargin_ (see createDomains; default + // margin = max() leaves this unclamped) + double dlo = std::max(time-dt, fitrange.begin() - config().domainmargin_); + TimeRange range(dlo,time); + Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); addDomain(domain,TimeDir::backwards); time = domain.begin(); } @@ -678,14 +700,19 @@ namespace KinKal { double time = drange.end(); while(time < fitrange.end()){ auto const& ktraj = fittraj_->nearestPiece(time); - double dt = bfield_.rangeInTolerance(ktraj,time,config().tol_); - TimeRange range(time,time+dt); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); + double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); + // clamp the domain high bound to the active range plus domainmargin_ (see createDomains; default + // margin = max() leaves this unclamped) + double dhi = std::min(time+dt, fitrange.end() + config().domainmargin_); + TimeRange range(time,dhi); + Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); addDomain(domain,TimeDir::forwards); time = domain.end(); } } } + // CALIBRATION PROBE – remove before PR + if(retval) std::cout << "KinKal::Track::extendDomains ndomains=" << domains_.size() << " fitrange=[" << fitrange.begin() << "," << fitrange.end() << "]" << std::endl; return retval; } @@ -750,15 +777,33 @@ namespace KinKal { auto const& ktraj = ptraj.nearestPiece(range.begin()); // catch exceptions if the fit extends beyond the range of the field map try { - double trange = bfield_.rangeInTolerance(ktraj,range.begin(),config().tol_); + // floor the step at config().mindtstep_: in high-gradient/low-momentum regions rangeInTolerance + // can return a vanishing step, which otherwise produces an unbounded number of domains (CPU hang/OOM) + double trange = std::max(bfield_.rangeInTolerance(ktraj,range.begin(),config().tol_),config().mindtstep_); // define 1st domain to have the 1st effect in the middle. This avoids effects having exactly the same time double tstart = range.begin() - 0.5*trange; do { // see how far we can go on the current traj before the DomainWall change causes the momentum estimate to go out of tolerance // note this assumes the trajectory is accurate (geometric extrapolation only) auto const& ktraj = ptraj.nearestPiece(tstart); - trange = bfield_.rangeInTolerance(ktraj,tstart,config().tol_); - domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)))); + trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); + // Clamp the domain BOUNDS to [range - domainmargin_, range + domainmargin_]. The walk steps + // ~0.5*trange beyond each end to bracket the edge effects, and in a uniform field (e.g. the + // tracker) trange is large, so an unclamped domain spans metres past the active region -- for a + // track that exits the solenoid bore there (a cosmic, just past the tracker), Bz->0 and the + // CentralHelix dPardB ~ bfrac/(1+bfrac) correction (bfrac=ΔBz/|Bnom|) is singular at bfrac=-1, + // blowing the reference off by metres and killing the fit. The fit has no measurements outside + // the active range, so a finite domainmargin_ confines the domains (and the field sampling at + // their midpoints, here and in convertSeed) to it; in a uniform field that collapses to a single + // in-bore domain (bfcorr=true reduces to the bfcorr=false behaviour the tracker needs). The + // default domainmargin_ = max() leaves the legacy unclamped behaviour unchanged. + double dlo = std::max(tstart, range.begin() - config().domainmargin_); + double dhi = std::min(tstart+trange, range.end() + config().domainmargin_); + if(dhi > dlo){ + TimeRange drange(dlo,dhi); + auto const& straj = ptraj.nearestPiece(drange.mid()); + domains.emplace(std::make_shared(drange,bfield_.fieldVect(straj.position3(drange.mid())))); + } // start the next domain at the end of this one tstart += trange; } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect @@ -800,9 +845,21 @@ namespace KinKal { while(fabs(time-tstart) < xtest.maxDt() && xtest.needsExtrapolation(*fittraj_,tdir) ){ // create a domain for this extrapolation auto const& ktraj = fittraj_->nearestPiece(time); - double dt = std::min(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),xtest.maxDtStep()); // always positive + // stop cleanly if the helix has gone singular + if( !std::isfinite(ktraj.momentum(time)) ) break; + // clamp the step: the floor (mindtstep_) bounds the domain count so a vanishing rangeInTolerance in a + // high-gradient/low-momentum region can no longer exhaust the MaxDt budget in micro-steps (OOM / truncation) + double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); // always positive TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); + // stop cleanly once we reach the (near) field-free region: the CentralHelix is ill-conditioned + // as |B|->0 (dPardB pole), the physical origin of the domain-walk runaway. Test the field where + // the NEW domain's BNom is sampled (its midpoint), not the current frontier -- a single coarse + // domain wall can otherwise jump the bnom straight from the tracker field to the field-free DS- + // shell region, hitting the pole before a frontier-only check sees it. Leaving from a valid state + // lets the caller (the line tail) continue as a straight line. + auto domainfield = bfield_.fieldVect(ktraj.position3(range.mid())); + if( config().minfield_ > 0.0 && domainfield.R() < config().minfield_ ) break; + Domain domain(range,domainfield); addDomain(domain,tdir,true); // use exact transport time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); } From a9861f3bb0fc393df71f4feabaeb12de0712d27e Mon Sep 17 00:00:00 2001 From: robmina Date: Tue, 30 Jun 2026 13:32:58 -0500 Subject: [PATCH 2/6] Fit: add MaxDomains cap on the BField domain walk; drop temporary diagnostics - Config: new maxdomains_ (default unlimited, preserving legacy behaviour) printed in operator<<. A hard cap on the number of BField domains a single fit may accumulate. - Track.hh createDomains/extendDomains: when domains exceed maxdomains_, throw -- caught by the existing fit()/processEnds() handlers, recording the fit as failed so the unusable track is dropped. A diverging low-momentum track can otherwise build ~1e5 domains (each a KKDW effect + traj piece) -> wasted CPU + ~GB memory before being dropped anyway. Validated over 112 cosmic files: usable tracks peak at 368 domains, runaways at 1e4-1e5; a cap of 1000 (set in the reco config) cleanly separates them. - Remove the temporary domain-count / empty-seed cout diagnostics added during the investigation (the empty-seed outsidemap guard itself is retained). Co-Authored-By: Claude Opus 4.8 --- Fit/Config.cc | 1 + Fit/Config.hh | 5 +++++ Fit/Track.hh | 15 ++++++++++----- 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/Fit/Config.cc b/Fit/Config.cc index d564f24f..8d3ee93f 100644 --- a/Fit/Config.cc +++ b/Fit/Config.cc @@ -11,6 +11,7 @@ namespace KinKal { << " min domain step (ns) " << kkconfig.mindtstep_ << " min field (T) " << kkconfig.minfield_ << " domain margin (ns) " << kkconfig.domainmargin_ + << " max domains " << kkconfig.maxdomains_ << " min NDOF " << kkconfig.minndof_ << " BField correction " << kkconfig.bfcorr_ << " with " << kkconfig.schedule().size() diff --git a/Fit/Config.hh b/Fit/Config.hh index ebcf81d2..21f53d6b 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -41,6 +41,11 @@ namespace KinKal { // field outside the solenoid bore -- where, for tracks that leave the magnet just // past the tracker (cosmics), Bz->0 makes the CentralHelix dPardB ~ 1/(1+dB/B) // correction singular and destroys the fit. Affects only createDomains/extendDomains. + unsigned maxdomains_ = std::numeric_limits::max(); // hard cap on the number of BField domains a + // single fit may accumulate. Default = unlimited (legacy). Set finite (e.g. 1000) to + // abort a fit whose iterative domain walk runs away: a diverging low-momentum track can + // build ~1e5 domains (wasted CPU + ~GB memory) before it is dropped anyway. The over-cap + // fit is failed cleanly. Affects createDomains/extendDomains. unsigned minndof_ = 5; // minimum number of DOFs to continue fit bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit bool ends_ = true; // process the passive effects at each end of the track after schedule completion diff --git a/Fit/Track.hh b/Fit/Track.hh index c43b20db..11a2fac9 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -400,7 +400,6 @@ namespace KinKal { // unfittable track: outsidemap stops the fit cleanly (needsFit()==false skips createEffects) and the // module drops it via goodFit()==false. fittraj_ is left as a valid (empty) object, never null. if(fittraj_->pieces().empty()){ - std::cout << "CONVERTSEED EMPTY range=[" << range.begin() << "," << range.end() << "] ndomains=" << domains.size() << std::endl; // DIAG - remove before PR history_.emplace_back(0,0,Status::outsidemap, "Empty seed trajectory (no domains)"); return; } @@ -478,8 +477,6 @@ namespace KinKal { status().comment_ = status().comment_ + error.what(); } } - // CALIBRATION PROBE – remove before PR - std::cout << "KinKal::Track::fit ndomains=" << domains_.size() << " status=" << fitStatus().status_ << std::endl; if(config().plevel_ > Config::none)print(std::cout, config().plevel_); } @@ -693,6 +690,11 @@ namespace KinKal { Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); addDomain(domain,TimeDir::backwards); time = domain.begin(); + // hard cap: a diverging low-momentum track can otherwise accumulate ~1e5 domains here (each adds a + // KKDW effect + traj piece) -> wasted CPU + ~GB memory before it is dropped anyway. Abort the runaway; + // iterate()'s caller catches this and records the fit as failed (the track is unusable -> dropped). + if(domains_.size() > config().maxdomains_) + throw std::runtime_error("Fit exceeded MaxDomains (BField domain walk runaway)"); } } // then forwards @@ -708,11 +710,11 @@ namespace KinKal { Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); addDomain(domain,TimeDir::forwards); time = domain.end(); + if(domains_.size() > config().maxdomains_) + throw std::runtime_error("Fit exceeded MaxDomains (BField domain walk runaway)"); } } } - // CALIBRATION PROBE – remove before PR - if(retval) std::cout << "KinKal::Track::extendDomains ndomains=" << domains_.size() << " fitrange=[" << fitrange.begin() << "," << fitrange.end() << "]" << std::endl; return retval; } @@ -806,6 +808,9 @@ namespace KinKal { } // start the next domain at the end of this one tstart += trange; + // hard cap (backup to extendDomains): bail if the initial domain build itself runs away + if(domains.size() > config().maxdomains_) + throw std::runtime_error("createDomains exceeded MaxDomains"); } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect } catch (std::exception const& error) { retval = false; From 6f32e49fc651dce3b6a006c4d8044cf8dbe67080 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Fri, 24 Jul 2026 12:57:49 -0500 Subject: [PATCH 3/6] Remove robust domains guard and shorten long comments. --- Fit/Config.hh | 20 +++--------- Fit/Track.hh | 84 ++++++++++++++++++--------------------------------- 2 files changed, 34 insertions(+), 70 deletions(-) diff --git a/Fit/Config.hh b/Fit/Config.hh index 21f53d6b..a162823d 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -30,22 +30,10 @@ namespace KinKal { double pdchisq_ = 1.0e6; // maximum allowed parameter change (units of chisqred) WRT previous reference double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm) double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps - double mindtstep_ = 1.0e-3; // ns: hard floor on the BField domain step, bounding the domain count in - // high-gradient/low-momentum regions where rangeInTolerance -> ~0 (else the - // bfcorr domain walk takes ~MaxDt/dt micro-steps -> OOM / time-budget truncation) - double minfield_ = 0.0; // T: if >0, stop the bfcorr extrapolation once |B| drops below this, so the - // ill-conditioned CentralHelix is never driven into B->0 (the physical runaway origin) - double domainmargin_ = std::numeric_limits::max(); // ns: max time a BField fit domain may extend - // beyond the active (hit) range. Default = unclamped (legacy). Set small (e.g. 0) - // to confine the fit's domains to the measurement region, so it never samples the - // field outside the solenoid bore -- where, for tracks that leave the magnet just - // past the tracker (cosmics), Bz->0 makes the CentralHelix dPardB ~ 1/(1+dB/B) - // correction singular and destroys the fit. Affects only createDomains/extendDomains. - unsigned maxdomains_ = std::numeric_limits::max(); // hard cap on the number of BField domains a - // single fit may accumulate. Default = unlimited (legacy). Set finite (e.g. 1000) to - // abort a fit whose iterative domain walk runs away: a diverging low-momentum track can - // build ~1e5 domains (wasted CPU + ~GB memory) before it is dropped anyway. The over-cap - // fit is failed cleanly. Affects createDomains/extendDomains. + double mindtstep_ = 0.0; // ns: hard floor on the BField domain step (0 = legacy; >0 bounds the domain count where rangeInTolerance->0) + double minfield_ = 0.0; // T: if >0, stop the bfcorr extrapolation once |B| drops below this (keeps CentralHelix out of the B->0 runaway) + double domainmargin_ = std::numeric_limits::max(); // ns: max time a fit domain may extend beyond the active range (max = unclamped/legacy) + unsigned maxdomains_ = std::numeric_limits::max(); // hard cap on BField domains per fit; over-cap fits fail cleanly (max = unlimited/legacy) unsigned minndof_ = 5; // minimum number of DOFs to continue fit bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit bool ends_ = true; // process the passive effects at each end of the track after schedule completion diff --git a/Fit/Track.hh b/Fit/Track.hh index 11a2fac9..3a564e36 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -392,17 +392,6 @@ namespace KinKal { newpiece.range() = domain->range(); fittraj_->append(newpiece); } - // A degenerate active range (an inverted/zero-width detectorRange for a pathological seed) can leave - // createDomains with zero domains, so the loop above appends nothing and fittraj_ is empty. createEffects - // would then build a Measurement against this empty PiecewiseTrajectory and throw - // std::length_error("Empty PiecewiseTrajectory!"), which -- being neither std::invalid_argument nor - // caught -- aborts the whole art event (and corrupts the output -> TTree::SetEntries). Treat it as an - // unfittable track: outsidemap stops the fit cleanly (needsFit()==false skips createEffects) and the - // module drops it via goodFit()==false. fittraj_ is left as a valid (empty) object, never null. - if(fittraj_->pieces().empty()){ - history_.emplace_back(0,0,Status::outsidemap, "Empty seed trajectory (no domains)"); - return; - } } else { // use the middle of the range as the nominal BField for this fit: double tref = range.mid(); @@ -430,12 +419,9 @@ namespace KinKal { auto prevdom = nextdom; ++nextdom; while( nextdom != domains.cend() ){ - // only connect domains that actually abut. A gap means prevdom and nextdom belong to different - // contiguous blocks (e.g. separate low/high fit extensions either side of the existing core domains); - // bridging across that gap would create a spurious DomainWall spanning the whole core, so skip it - // (this previously threw "Invalid domains" and aborted the fit). - if(fabs(prevdom->get()->end()-nextdom->get()->begin())<=1e-10) - effects_.emplace_back(std::make_unique(*prevdom,*nextdom ,*fittraj_)); + // must be contiguous + if(fabs(prevdom->get()->end()-nextdom->get()->begin())>1e-10)throw std::invalid_argument("Invalid domains"); + effects_.emplace_back(std::make_unique(*prevdom,*nextdom ,*fittraj_)); prevdom = nextdom; ++nextdom; } @@ -683,16 +669,15 @@ namespace KinKal { while(time > fitrange.begin()){ auto const& ktraj = fittraj_->nearestPiece(time); double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); - // clamp the domain low bound to the active range minus domainmargin_ (see createDomains; default - // margin = max() leaves this unclamped) + // clamp the domain low bound to the active range minus domainmargin_ (max = unclamped/legacy) double dlo = std::max(time-dt, fitrange.begin() - config().domainmargin_); TimeRange range(dlo,time); - Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); + // sample BNom at the domain-midpoint piece only when confined (domainmargin_ set); else legacy piece (nearest time) + auto const& straj = (config().domainmargin_ < std::numeric_limits::max()) ? fittraj_->nearestPiece(range.mid()) : ktraj; + Domain domain(range,bfield_.fieldVect(straj.position3(range.mid()))); addDomain(domain,TimeDir::backwards); time = domain.begin(); - // hard cap: a diverging low-momentum track can otherwise accumulate ~1e5 domains here (each adds a - // KKDW effect + traj piece) -> wasted CPU + ~GB memory before it is dropped anyway. Abort the runaway; - // iterate()'s caller catches this and records the fit as failed (the track is unusable -> dropped). + // abort a runaway domain walk (caught by iterate()'s caller -> fit failed -> track dropped) if(domains_.size() > config().maxdomains_) throw std::runtime_error("Fit exceeded MaxDomains (BField domain walk runaway)"); } @@ -703,11 +688,12 @@ namespace KinKal { while(time < fitrange.end()){ auto const& ktraj = fittraj_->nearestPiece(time); double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); - // clamp the domain high bound to the active range plus domainmargin_ (see createDomains; default - // margin = max() leaves this unclamped) + // clamp the domain high bound to the active range plus domainmargin_ (max = unclamped/legacy) double dhi = std::min(time+dt, fitrange.end() + config().domainmargin_); TimeRange range(time,dhi); - Domain domain(range,bfield_.fieldVect(fittraj_->nearestPiece(range.mid()).position3(range.mid()))); + // sample BNom at the domain-midpoint piece only when confined (domainmargin_ set); else legacy piece (nearest time) + auto const& straj = (config().domainmargin_ < std::numeric_limits::max()) ? fittraj_->nearestPiece(range.mid()) : ktraj; + Domain domain(range,bfield_.fieldVect(straj.position3(range.mid()))); addDomain(domain,TimeDir::forwards); time = domain.end(); if(domains_.size() > config().maxdomains_) @@ -779,8 +765,7 @@ namespace KinKal { auto const& ktraj = ptraj.nearestPiece(range.begin()); // catch exceptions if the fit extends beyond the range of the field map try { - // floor the step at config().mindtstep_: in high-gradient/low-momentum regions rangeInTolerance - // can return a vanishing step, which otherwise produces an unbounded number of domains (CPU hang/OOM) + // floor the step at mindtstep_ so a vanishing rangeInTolerance can't spawn unbounded domains (CPU hang/OOM) double trange = std::max(bfield_.rangeInTolerance(ktraj,range.begin(),config().tol_),config().mindtstep_); // define 1st domain to have the 1st effect in the middle. This avoids effects having exactly the same time double tstart = range.begin() - 0.5*trange; @@ -789,22 +774,19 @@ namespace KinKal { // note this assumes the trajectory is accurate (geometric extrapolation only) auto const& ktraj = ptraj.nearestPiece(tstart); trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); - // Clamp the domain BOUNDS to [range - domainmargin_, range + domainmargin_]. The walk steps - // ~0.5*trange beyond each end to bracket the edge effects, and in a uniform field (e.g. the - // tracker) trange is large, so an unclamped domain spans metres past the active region -- for a - // track that exits the solenoid bore there (a cosmic, just past the tracker), Bz->0 and the - // CentralHelix dPardB ~ bfrac/(1+bfrac) correction (bfrac=ΔBz/|Bnom|) is singular at bfrac=-1, - // blowing the reference off by metres and killing the fit. The fit has no measurements outside - // the active range, so a finite domainmargin_ confines the domains (and the field sampling at - // their midpoints, here and in convertSeed) to it; in a uniform field that collapses to a single - // in-bore domain (bfcorr=true reduces to the bfcorr=false behaviour the tracker needs). The - // default domainmargin_ = max() leaves the legacy unclamped behaviour unchanged. - double dlo = std::max(tstart, range.begin() - config().domainmargin_); - double dhi = std::min(tstart+trange, range.end() + config().domainmargin_); - if(dhi > dlo){ - TimeRange drange(dlo,dhi); - auto const& straj = ptraj.nearestPiece(drange.mid()); - domains.emplace(std::make_shared(drange,bfield_.fieldVect(straj.position3(drange.mid())))); + // domainmargin_ confines the domains to the active range +/- margin, keeping the field sampling out of the B->0 region past the tracker where the CentralHelix dPardB correction is singular + if(config().domainmargin_ < std::numeric_limits::max()){ + // confined: clamp the bounds, drop a fully clamped-out domain, sample BNom at the clamped midpoint + double dlo = std::max(tstart, range.begin() - config().domainmargin_); + double dhi = std::min(tstart+trange, range.end() + config().domainmargin_); + if(dhi > dlo){ + TimeRange drange(dlo,dhi); + auto const& straj = ptraj.nearestPiece(drange.mid()); + domains.emplace(std::make_shared(drange,bfield_.fieldVect(straj.position3(drange.mid())))); + } + } else { + // legacy (default): always emplace (even a zero-width boundary domain), sampling BNom at ktraj, so this is bit-identical to upstream + domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)))); } // start the next domain at the end of this one tstart += trange; @@ -835,6 +817,8 @@ namespace KinKal { tmax = std::max(tmax,exing->time()); } } + // no (active) effects leaves tmin>tmax; return a null range instead of an invalid one that would throw + if(tmax < tmin) return TimeRange(); return TimeRange(tmin,tmax); } @@ -850,18 +834,10 @@ namespace KinKal { while(fabs(time-tstart) < xtest.maxDt() && xtest.needsExtrapolation(*fittraj_,tdir) ){ // create a domain for this extrapolation auto const& ktraj = fittraj_->nearestPiece(time); - // stop cleanly if the helix has gone singular - if( !std::isfinite(ktraj.momentum(time)) ) break; - // clamp the step: the floor (mindtstep_) bounds the domain count so a vanishing rangeInTolerance in a - // high-gradient/low-momentum region can no longer exhaust the MaxDt budget in micro-steps (OOM / truncation) + // floor the step at mindtstep_ so a vanishing rangeInTolerance can't exhaust MaxDt in micro-steps (OOM) double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); // always positive TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); - // stop cleanly once we reach the (near) field-free region: the CentralHelix is ill-conditioned - // as |B|->0 (dPardB pole), the physical origin of the domain-walk runaway. Test the field where - // the NEW domain's BNom is sampled (its midpoint), not the current frontier -- a single coarse - // domain wall can otherwise jump the bnom straight from the tracker field to the field-free DS- - // shell region, hitting the pole before a frontier-only check sees it. Leaving from a valid state - // lets the caller (the line tail) continue as a straight line. + // stop before the near-field-free region (dPardB pole), testing |B| at the new domain's midpoint not the frontier, so the line tail can take over auto domainfield = bfield_.fieldVect(ktraj.position3(range.mid())); if( config().minfield_ > 0.0 && domainfield.R() < config().minfield_ ) break; Domain domain(range,domainfield); From 7a5b6bd8941658800cd094e40cc6cce552e6c965 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Tue, 28 Jul 2026 13:37:07 -0500 Subject: [PATCH 4/6] Option to handoff helical->geometric linear extrapolation when B field ~= 0. --- Fit/Config.cc | 1 + Fit/Config.hh | 5 +- Fit/Track.hh | 287 ++++++++++++++++++++++++++++++------------- General/BFieldMap.hh | 5 + 4 files changed, 211 insertions(+), 87 deletions(-) diff --git a/Fit/Config.cc b/Fit/Config.cc index 8d3ee93f..eefcadb8 100644 --- a/Fit/Config.cc +++ b/Fit/Config.cc @@ -14,6 +14,7 @@ namespace KinKal { << " max domains " << kkconfig.maxdomains_ << " min NDOF " << kkconfig.minndof_ << " BField correction " << kkconfig.bfcorr_ + << " zero-field extrap handoff " << kkconfig.zerofield_extrap_ << " with " << kkconfig.schedule().size() << " Meta-iterations:" << std::endl; for(auto const& miconfig : kkconfig.schedule() ) { diff --git a/Fit/Config.hh b/Fit/Config.hh index a162823d..701bb533 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -31,11 +31,12 @@ namespace KinKal { double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm) double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps double mindtstep_ = 0.0; // ns: hard floor on the BField domain step (0 = legacy; >0 bounds the domain count where rangeInTolerance->0) - double minfield_ = 0.0; // T: if >0, stop the bfcorr extrapolation once |B| drops below this (keeps CentralHelix out of the B->0 runaway) - double domainmargin_ = std::numeric_limits::max(); // ns: max time a fit domain may extend beyond the active range (max = unclamped/legacy) + double minfield_ = 0.0; // T: if >0, and zerofield_extrap_ is enabled, hand bfcorr extrapolation off to free-particle continuation once |B| drops below this + double domainmargin_ = std::numeric_limits::max(); // ns: max time a fit domain may extend beyond the active range (max = unclamped/legacy overhang; finite = confine walk+sampling to range±margin) unsigned maxdomains_ = std::numeric_limits::max(); // hard cap on BField domains per fit; over-cap fits fail cleanly (max = unlimited/legacy) unsigned minndof_ = 5; // minimum number of DOFs to continue fit bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit + bool zerofield_extrap_ = false; // if true: (1) bfcorr extrapolate() hands off to geometric free-particle continuation outside the map / below minfield_; (2) createDomains stops DomainWalls at that edge instead of failing Extension; (3) replaceDomains charge/mass mismatch soft-keeps the prior usable fit (CH cosmic CRV). Default false preserves legacy LH/CH behaviour. bool ends_ = true; // process the passive effects at each end of the track after schedule completion printLevel plevel_ = none; // print level // schedule of meta-iterations. These will be executed sequentially until completion or failure diff --git a/Fit/Track.hh b/Fit/Track.hh index 3a564e36..fea0356d 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -209,6 +209,12 @@ namespace KinKal { auto jdom= domains.rbegin(); while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom; domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator + // Hit/ParameterHit times can fall outside every saved domain (e.g. DomainMargin=0 CHTruthSeed + // with a short domainBounds span vs a longer traj piece). Soft-fail instead of deref empty. + if(domains.empty()){ + history_.emplace_back(0,0,Status::outsidemap, "Empty domains after detector-range trim"); + return; + } // trim the trajectory to this range detrange.combine((*domains.begin())->range()); detrange.combine((*domains.rbegin())->range()); @@ -270,7 +276,19 @@ namespace KinKal { // create domains for the whole range dok &= createDomains(*fittraj_,exrange, domains); // replace previous domains with these. This replaces the trajectory and bfield-related effects - if(dok)replaceDomains(domains); + if(config().zerofield_extrap_ && dok && domains.empty()){ + // Map-edge stop before any domain: do not call replaceDomains on an empty set. + dok = false; + } else if(dok){ + // CH rebuild under Extension's tighter BCorrTolerance can flip omega/charge near the + // map edge → ParticleTrajectory::append throws. Keep the usable construction fit. + try { + replaceDomains(domains); + } catch (std::invalid_argument const&) { + if(!config().zerofield_extrap_) throw; + dok = false; + } + } } else { // create domains just for the extensions TimeRange exlow(exrange.begin(),fittraj_->range().begin()); @@ -288,7 +306,9 @@ namespace KinKal { } } if(!dok){ - // domain calculation failed: abort the fit + // domain calculation failed. With ZeroFieldExtrap, keep a previously usable fit + // (map-edge truncation is preferred inside createDomains; this is a safety net). + if(config().zerofield_extrap_ && fitStatus().usable()) return; history_.push_back(Status(0)); status().status_ = Status::outsidemap; status().comment_ = std::string("Extension error"); @@ -305,11 +325,48 @@ namespace KinKal { // replace domains when DomainWall correction is added or changed. the traj must also be replaced, so that // the pieces correspond to the new domains. The new traj is geometrically equivalent, but not parametrically equal. + // Build the replacement traj first so a failed append (e.g. CH charge flip) leaves domains_/effects_/fittraj_ intact. template void Track::replaceDomains(DOMAINCOL const& domains) { - // if domains exist, clear them and remove all DomainWall effects + auto newtraj = std::make_unique(); + // loop over domains, splitting the overlapping traj pieces at the domain walls, and transforming them to reference the domain's field + // This increases the number of traj pieces. + // extend the existing traj to the domain range (restored on failure) + TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end()); + TimeRange front_range = fittraj_->front().range(); + TimeRange back_range = fittraj_->back().range(); + fittraj_->setRange(drange); + try { + for(auto const& domain : domains) { + // find the range of existing ptraj pieces that overlaps with this domain's range + using KTRAJPTR = std::shared_ptr; + using DKTRAJ = std::deque; + using DKTRAJCITER = typename DKTRAJ::const_iterator; + DKTRAJCITER first,last; + fittraj_->pieceRange(domain->range(),first,last); + // loop over these pieces; first and last can be the same! + auto olditer = first; + do { + auto const& oldpiece = **olditer; + // copy this piece, translating bnom to this domain's field + KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); + // set the range for this piece, making sure it is non-zero + double tstart = std::max(domain->begin(), oldpiece.range().begin()); + double tend = std::min(domain->end(),oldpiece.range().end()); + if(tstart < tend){ + newpiece.range() = TimeRange(tstart,tend); + newtraj->append(newpiece); + } + if(olditer != last)++olditer; + } while(olditer != last); + } + } catch (...) { + fittraj_->front().setRange(front_range); + fittraj_->back().setRange(back_range); + throw; + } + // commit: clear old domains / DomainWall effects, retarget remaining effects, swap traj if(domains_.size() > 0){ domains_.clear(); - // remove all existing DomainWall effects auto ieff = effects_.begin(); while(ieff != effects_.end()){ const KKDW* kkbf = dynamic_cast(ieff->get()); @@ -320,40 +377,9 @@ namespace KinKal { } } } - auto newtraj = std::make_unique(); - // loop over domains, splitting the overlapping traj pieces at the domain walls, and transforming them to reference the domain's field - // This increases the number of traj pieces. - // extend the existing traj to the domain range - TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end()); - fittraj_->setRange(drange); - for(auto const& domain : domains) { - // find the range of existing ptraj pieces that overlaps with this domain's range - using KTRAJPTR = std::shared_ptr; - using DKTRAJ = std::deque; - using DKTRAJCITER = typename DKTRAJ::const_iterator; - DKTRAJCITER first,last; - fittraj_->pieceRange(domain->range(),first,last); - // loop over these pieces; first and last can be the same! - auto olditer = first; - do { - auto const& oldpiece = **olditer; - // copy this piece, translating bnom to this domain's field - KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); - // set the range for this piece, making sure it is non-zero - double tstart = std::max(domain->begin(), oldpiece.range().begin()); - double tend = std::min(domain->end(),oldpiece.range().end()); - if(tstart < tend){ - newpiece.range() = TimeRange(tstart,tend); - newtraj->append(newpiece); - } - if(olditer != last)++olditer; - } while(olditer != last); - } - // switch over any existing effects to reference this traj (could be none) for (auto& eff : effects_) { eff->updateReference(*newtraj); } - // swap out the fit trajectory; this will be used as reference for the next iterations fittraj_.swap(newtraj); } @@ -392,6 +418,13 @@ namespace KinKal { newpiece.range() = domain->range(); fittraj_->append(newpiece); } + // Degenerate active range (or DomainMargin confinement dropping every domain) can leave + // createDomains with zero domains → empty fittraj_. createEffects would then throw + // std::length_error("Empty PiecewiseTrajectory!") and abort the art event. Soft-fail instead. + if(fittraj_->pieces().empty()){ + history_.emplace_back(0,0,Status::outsidemap, "Empty seed trajectory (no domains)"); + return; + } } else { // use the middle of the range as the nominal BField for this fit: double tref = range.mid(); @@ -762,40 +795,63 @@ namespace KinKal { template bool Track::createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const { bool retval(true); if(config().bfcorr_ ) { - auto const& ktraj = ptraj.nearestPiece(range.begin()); + // With ZeroFieldExtrap: stop DomainWalls at the map / B≈0 edge instead of throwing. + // Keeps domains built so far and lets Extension proceed (same gate as extrapolate handoff). + auto atZeroField = [&](VEC3 const& pos) -> bool { + if(!config().zerofield_extrap_) return false; + bool outside = !bfield_.inRange(pos); + VEC3 b = outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(pos); + return outside || BFieldMap::isZeroField(b) + || (config().minfield_ > 0.0 && b.R() < config().minfield_); + }; // catch exceptions if the fit extends beyond the range of the field map try { - // floor the step at mindtstep_ so a vanishing rangeInTolerance can't spawn unbounded domains (CPU hang/OOM) - double trange = std::max(bfield_.rangeInTolerance(ktraj,range.begin(),config().tol_),config().mindtstep_); - // define 1st domain to have the 1st effect in the middle. This avoids effects having exactly the same time - double tstart = range.begin() - 0.5*trange; - do { - // see how far we can go on the current traj before the DomainWall change causes the momentum estimate to go out of tolerance - // note this assumes the trajectory is accurate (geometric extrapolation only) - auto const& ktraj = ptraj.nearestPiece(tstart); - trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); - // domainmargin_ confines the domains to the active range +/- margin, keeping the field sampling out of the B->0 region past the tracker where the CentralHelix dPardB correction is singular - if(config().domainmargin_ < std::numeric_limits::max()){ - // confined: clamp the bounds, drop a fully clamped-out domain, sample BNom at the clamped midpoint - double dlo = std::max(tstart, range.begin() - config().domainmargin_); - double dhi = std::min(tstart+trange, range.end() + config().domainmargin_); - if(dhi > dlo){ - TimeRange drange(dlo,dhi); + if(config().domainmargin_ < std::numeric_limits::max()){ + // Confined (DomainMargin set): walk ONLY within active range ± margin. Do not use the + // legacy half-domain overhang past that window — for near-uniform B (cosmic CentralHelix) + // rangeInTolerance is huge, so begin-0.5*trange geometrically extrapolates far outside + // the hits / BField maps even though every hit is inside. + double const tlo = range.begin() - config().domainmargin_; + double const thi = range.end() + config().domainmargin_; + double tstart = tlo; + while(tstart < thi){ + auto const& ktraj = ptraj.nearestPiece(tstart); + if(atZeroField(ktraj.position3(tstart))) break; + double trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); + double dhi = std::min(tstart + trange, thi); + if(dhi > tstart){ + TimeRange drange(tstart,dhi); auto const& straj = ptraj.nearestPiece(drange.mid()); - domains.emplace(std::make_shared(drange,bfield_.fieldVect(straj.position3(drange.mid())))); + VEC3 midpos = straj.position3(drange.mid()); + if(atZeroField(midpos)) break; + domains.emplace(std::make_shared(drange,bfield_.fieldVect(midpos))); } - } else { - // legacy (default): always emplace (even a zero-width boundary domain), sampling BNom at ktraj, so this is bit-identical to upstream - domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)))); + tstart = dhi; + if(domains.size() > config().maxdomains_) + throw std::runtime_error("createDomains exceeded MaxDomains"); } - // start the next domain at the end of this one - tstart += trange; - // hard cap (backup to extendDomains): bail if the initial domain build itself runs away - if(domains.size() > config().maxdomains_) - throw std::runtime_error("createDomains exceeded MaxDomains"); - } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect + } else { + // Legacy (default DomainMargin = max): half-domain overhang so the first/last effect sits + // mid-domain. Bit-identical to upstream when zerofield_extrap_ is false. + auto const& ktraj0 = ptraj.nearestPiece(range.begin()); + if(atZeroField(ktraj0.position3(range.begin()))) return true; + double trange = std::max(bfield_.rangeInTolerance(ktraj0,range.begin(),config().tol_),config().mindtstep_); + double tstart = range.begin() - 0.5*trange; + do { + auto const& ktraj = ptraj.nearestPiece(tstart); + if(atZeroField(ktraj.position3(tstart))) break; + trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); + VEC3 sample = ktraj.position3(tstart+0.5*trange); + if(atZeroField(sample)) break; + domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(sample))); + tstart += trange; + if(domains.size() > config().maxdomains_) + throw std::runtime_error("createDomains exceeded MaxDomains"); + } while(tstart < range.end() + 0.5*trange); + } } catch (std::exception const& error) { - retval = false; + // ZeroFieldExtrap: treat unexpected map samples as a soft stop (keep domains so far). + if(!config().zerofield_extrap_) retval = false; } } return retval; @@ -826,31 +882,93 @@ namespace KinKal { bool retval = fitStatus().usable(); if(retval){ if(config().bfcorr_){ - // test for extrapolation outside the bfield map range - try { - // iterate until the extrapolation condition is met + // Opt-in zero-field handoff (zerofield_extrap_): for CH cosmic CRV extrapolation that must + // leave the map. Default false → legacy bfcorr domain walk unchanged (LH-safe). + if(config().zerofield_extrap_){ + auto geometricExtend = [&](double tmax_remaining) { + if(tmax_remaining <= 0.0) return; + auto& endpiece = tdir == TimeDir::forwards ? fittraj_->backPtr() : fittraj_->frontPtr(); + double time = tdir == TimeDir::forwards ? endpiece->range().end() : endpiece->range().begin(); + double tstart = time; + bool needsext(true); + do { + TimeRange newrange = tdir == TimeDir::forwards ? + TimeRange(endpiece->range().begin(),endpiece->range().end()+xtest.maxDtStep()) + : + TimeRange(endpiece->range().begin()-xtest.maxDtStep(),endpiece->range().end()); + endpiece->setRange(newrange); + time = tdir == TimeDir::forwards ? endpiece->range().end() : endpiece->range().begin(); + needsext = xtest.needsExtrapolation(*fittraj_,tdir); + } while(needsext && fabs(time-tstart) < tmax_remaining); + }; + + bool handed_off = false; double time = tdir == TimeDir::forwards ? domains_.crbegin()->get()->end() : domains_.cbegin()->get()->begin(); double tstart = time; - while(fabs(time-tstart) < xtest.maxDt() && xtest.needsExtrapolation(*fittraj_,tdir) ){ - // create a domain for this extrapolation - auto const& ktraj = fittraj_->nearestPiece(time); - // floor the step at mindtstep_ so a vanishing rangeInTolerance can't exhaust MaxDt in micro-steps (OOM) - double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); // always positive - TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); - // stop before the near-field-free region (dPardB pole), testing |B| at the new domain's midpoint not the frontier, so the line tail can take over - auto domainfield = bfield_.fieldVect(ktraj.position3(range.mid())); - if( config().minfield_ > 0.0 && domainfield.R() < config().minfield_ ) break; - Domain domain(range,domainfield); - addDomain(domain,tdir,true); // use exact transport - time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); + try { + while(fabs(time-tstart) < xtest.maxDt() && xtest.needsExtrapolation(*fittraj_,tdir) ){ + auto const& ktraj = fittraj_->nearestPiece(time); + if( !std::isfinite(ktraj.momentum(time)) ) break; + + VEC3 frontier = ktraj.position3(time); + bool outside = !bfield_.inRange(frontier); + // Only sample the map when inside it — never call fieldDeriv out of range + VEC3 bfront = outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(frontier); + bool zerofield = outside || BFieldMap::isZeroField(bfront) + || (config().minfield_ > 0.0 && bfront.R() < config().minfield_); + if(zerofield){ + // Leave the bfcorr / DomainWall path. Free-particle continuation is geometric + // range-extend of the current end piece (no CH rebuild at B≈0). + handed_off = true; + break; + } + + double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); + TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); + VEC3 midpos = ktraj.position3(range.mid()); + bool mid_outside = !bfield_.inRange(midpos); + VEC3 domainfield = mid_outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(midpos); + if(mid_outside || BFieldMap::isZeroField(domainfield) + || (config().minfield_ > 0.0 && domainfield.R() < config().minfield_)){ + handed_off = true; + break; + } + Domain domain(range,domainfield); + addDomain(domain,tdir,true); + time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); + } + } catch (std::exception const& error) { + history_.push_back(Status(0)); + status().status_ = Status::outsidemap; + status().comment_ = std::string("Extrapolation error"); + retval = false; } - } catch (std::exception const& error) { - history_.push_back(Status(0)); - status().status_ = Status::outsidemap; - status().comment_ = std::string("Extrapolation error"); - retval = false; + if(retval && handed_off && xtest.needsExtrapolation(*fittraj_,tdir)){ + geometricExtend(xtest.maxDt() - fabs(time-tstart)); + } + } else { + // Legacy bfcorr extrapolation (default): unchanged domain walk + try { + double time = tdir == TimeDir::forwards ? domains_.crbegin()->get()->end() : domains_.cbegin()->get()->begin(); + double tstart = time; + while(fabs(time-tstart) < xtest.maxDt() && xtest.needsExtrapolation(*fittraj_,tdir) ){ + auto const& ktraj = fittraj_->nearestPiece(time); + double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); + TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); + auto domainfield = bfield_.fieldVect(ktraj.position3(range.mid())); + if( config().minfield_ > 0.0 && domainfield.R() < config().minfield_ ) break; + Domain domain(range,domainfield); + addDomain(domain,tdir,true); + time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); + } + } catch (std::exception const& error) { + history_.push_back(Status(0)); + status().status_ = Status::outsidemap; + status().comment_ = std::string("Extrapolation error"); + retval = false; + } + retval = true; } - retval = true; } else { // geometric extrapolation of the end piece; no need to protect auto& endpiece = tdir == TimeDir::forwards ? fittraj_->backPtr() : fittraj_->frontPtr(); @@ -858,7 +976,6 @@ namespace KinKal { double tstart = time; bool needsext(true); do { - // extend the range by the step dt TimeRange newrange = tdir == TimeDir::forwards ? TimeRange(endpiece->range().begin(),endpiece->range().end()+xtest.maxDtStep()) : diff --git a/General/BFieldMap.hh b/General/BFieldMap.hh index 8ed493db..0e9b1603 100644 --- a/General/BFieldMap.hh +++ b/General/BFieldMap.hh @@ -33,6 +33,11 @@ namespace KinKal { BFieldMap& operator =(BFieldMap const& ) = delete; // speed of light in units to convert Tesla to mm (bending radius) static double constexpr cbar() { return CLHEP::c_light/1000.0; } + // |B| below this (T) is treated as physically zero for ZeroFieldExtrap handoff decisions. + // Fit paths must not sample here; with zerofield_extrap_, extrapolation leaves the bfcorr + // domain walk and continues by geometric range-extend (no CH rebuild at B≈0). + static double constexpr zeroField() { return 1.0e-6; } + static bool isZeroField(VEC3 const& bvec) { return bvec.R() < zeroField(); } // templated interface for interacting with kinematic trajectory classes // how far can you go along the given kinematic trajectory till BField inhomogeneity makes the momentum accuracy out of (fractional) tolerance template double rangeInTolerance(KTRAJ const& ktraj, double tstart, double tol) const; From ed008025959b987de0c225872f212e1c8cfec040 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Thu, 6 Aug 2026 20:21:25 -0500 Subject: [PATCH 5/6] Remove MaxDomains config, add new status code for invalid trajectory piece, and move B field valid region testing into constructor so that it is precomputed at run start rather than testing for each seed/track. --- Fit/Config.cc | 3 - Fit/Config.hh | 7 +- Fit/Status.cc | 2 + Fit/Status.hh | 2 +- Fit/Track.hh | 275 +++++++++++++++++-------------- General/BFieldMap.hh | 129 ++++++++++++++- Trajectory/ParticleTrajectory.hh | 12 +- 7 files changed, 292 insertions(+), 138 deletions(-) diff --git a/Fit/Config.cc b/Fit/Config.cc index eefcadb8..ab92005a 100644 --- a/Fit/Config.cc +++ b/Fit/Config.cc @@ -9,12 +9,9 @@ namespace KinKal { << " diverge traj gap (mm) " << kkconfig.divgap_ << " fractional momentum tolerance " << kkconfig.tol_ << " min domain step (ns) " << kkconfig.mindtstep_ - << " min field (T) " << kkconfig.minfield_ << " domain margin (ns) " << kkconfig.domainmargin_ - << " max domains " << kkconfig.maxdomains_ << " min NDOF " << kkconfig.minndof_ << " BField correction " << kkconfig.bfcorr_ - << " zero-field extrap handoff " << kkconfig.zerofield_extrap_ << " with " << kkconfig.schedule().size() << " Meta-iterations:" << std::endl; for(auto const& miconfig : kkconfig.schedule() ) { diff --git a/Fit/Config.hh b/Fit/Config.hh index 701bb533..0756d10b 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -30,13 +30,10 @@ namespace KinKal { double pdchisq_ = 1.0e6; // maximum allowed parameter change (units of chisqred) WRT previous reference double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm) double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps - double mindtstep_ = 0.0; // ns: hard floor on the BField domain step (0 = legacy; >0 bounds the domain count where rangeInTolerance->0) - double minfield_ = 0.0; // T: if >0, and zerofield_extrap_ is enabled, hand bfcorr extrapolation off to free-particle continuation once |B| drops below this - double domainmargin_ = std::numeric_limits::max(); // ns: max time a fit domain may extend beyond the active range (max = unclamped/legacy overhang; finite = confine walk+sampling to range±margin) - unsigned maxdomains_ = std::numeric_limits::max(); // hard cap on BField domains per fit; over-cap fits fail cleanly (max = unlimited/legacy) + double mindtstep_ = 0.0; // ns: minimum domain range. >0 bounds the domain count at (walk window)/mindtstep_; 0 leaves it unbounded + double domainmargin_ = std::numeric_limits::max(); // ns: max time a domain may extend beyond the active range; finite confines the walk and its field sampling to range +/- margin unsigned minndof_ = 5; // minimum number of DOFs to continue fit bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit - bool zerofield_extrap_ = false; // if true: (1) bfcorr extrapolate() hands off to geometric free-particle continuation outside the map / below minfield_; (2) createDomains stops DomainWalls at that edge instead of failing Extension; (3) replaceDomains charge/mass mismatch soft-keeps the prior usable fit (CH cosmic CRV). Default false preserves legacy LH/CH behaviour. bool ends_ = true; // process the passive effects at each end of the track after schedule completion printLevel plevel_ = none; // print level // schedule of meta-iterations. These will be executed sequentially until completion or failure diff --git a/Fit/Status.cc b/Fit/Status.cc index dda58614..7edef0f3 100644 --- a/Fit/Status.cc +++ b/Fit/Status.cc @@ -21,6 +21,8 @@ namespace KinKal { return "OutsideBFieldMap "; case Status::failed: return "Failed "; + case Status::incompatiblepiece: + return "IncompatiblePiece "; } } diff --git a/Fit/Status.hh b/Fit/Status.hh index f94fa954..e769afb1 100644 --- a/Fit/Status.hh +++ b/Fit/Status.hh @@ -9,7 +9,7 @@ namespace KinKal { // struct to define fit status struct Status { - enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status + enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed,incompatiblepiece}; // fit status unsigned miter_; // meta-iteration number; unsigned iter_; // iteration number; status status_; // current status diff --git a/Fit/Track.hh b/Fit/Track.hh index fea0356d..f070b9bc 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -56,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -117,6 +118,7 @@ namespace KinKal { std::vector const& history() const { return history_; } Status const& fitStatus() const { return history_.back(); } // most recent status PKTRAJ const& fitTraj() const { return *fittraj_; } + bool hasTraj() const { return static_cast(fittraj_); } // false for a fit that failed before the traj was built KKEFFCOL const& effects() const { return effects_; } Config const& config() const { return config_.back(); } CONFIGCOL const& configs() const { return config_; } @@ -137,6 +139,12 @@ namespace KinKal { void convertSeed(KTRAJ const& seedtraj,TimeRange const& refrange, DOMAINCOL& domains); void fit(); // process the effects and create the trajectory. This executes the current schedule bool createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const; + // build one usable domain starting at tstart, or nothing if the field there can't support one + std::optional createDomain(PKTRAJ const& ptraj, double tstart, double tend) const; + // input preconditions that don't depend on the BField. Checking them here keeps unusable input out + // of the domain walk and out of the trajectory; records the reason and returns false on failure. + bool validInput(TimeRange const& detrange); + bool validInput(TimeRange const& detrange, KTRAJ const& seedtraj); bool setBounds(KKEFFFWDBND& fwdbnds,KKEFFREVBND& revbnds); // set the bounds. Returns false if the bounds are empty bool extendDomains(TimeRange const& fitrange); // extend domains if the fit range changes. Return value says if domains were added void updateDomains(PKTRAJ const& ptraj); // Update domains between iterations @@ -145,7 +153,7 @@ namespace KinKal { void initFitState(FitStateArray& states, TimeRange const& fitrange, double dwt=1.0); PKTRAJPTR initTraj(FitState& state, TimeRange const& fitrange); bool canIterate() const; - void replaceDomains(DOMAINCOL const& domains); + bool replaceDomains(DOMAINCOL const& domains); void extendTraj(DOMAINCOL const& domains); void processEnds(); // add a single domain within the tolerance and extend the fit in the specified direction. @@ -186,6 +194,7 @@ namespace KinKal { template void Track::fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj) { auto detrange = detectorRange(hits,exings,true); + if(!validInput(detrange,seedtraj)) return; // convert the seed traj to a piecewaise traj. This creates the domains DOMAINCOL domains; convertSeed(seedtraj,detrange,domains); @@ -201,6 +210,7 @@ namespace KinKal { fittraj_ = std::move(fittraj); // steal the underlying object // truncate the domains and fit trajectory to be within the detector range auto detrange = detectorRange(hits,exings,true); + if(!validInput(detrange)) return; if(domains.size() > 0){ auto idom = domains.begin(); // stop at the 1st domain overlaping the detector range, and erase all elements up to that point @@ -209,8 +219,8 @@ namespace KinKal { auto jdom= domains.rbegin(); while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom; domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator - // Hit/ParameterHit times can fall outside every saved domain (e.g. DomainMargin=0 CHTruthSeed - // with a short domainBounds span vs a longer traj piece). Soft-fail instead of deref empty. + // hit times can fall outside every saved domain when the domain span is shorter than the + // trajectory piece; soft-fail rather than dereference an empty set if(domains.empty()){ history_.emplace_back(0,0,Status::outsidemap, "Empty domains after detector-range trim"); return; @@ -276,17 +286,17 @@ namespace KinKal { // create domains for the whole range dok &= createDomains(*fittraj_,exrange, domains); // replace previous domains with these. This replaces the trajectory and bfield-related effects - if(config().zerofield_extrap_ && dok && domains.empty()){ + if(bfield_.protecting() && dok && domains.empty()){ // Map-edge stop before any domain: do not call replaceDomains on an empty set. dok = false; } else if(dok){ - // CH rebuild under Extension's tighter BCorrTolerance can flip omega/charge near the - // map edge → ParticleTrajectory::append throws. Keep the usable construction fit. - try { - replaceDomains(domains); - } catch (std::invalid_argument const&) { - if(!config().zerofield_extrap_) throw; - dok = false; + // a tighter extension tolerance can flip omega near a collapsing field, which the + // parameterization reads as a charge change; leave the track untouched and record why + if(!replaceDomains(domains)){ + history_.push_back(Status(0)); + status().status_ = Status::incompatiblepiece; + status().comment_ = std::string("Domain replacement: incompatible piece"); + return; } } } else { @@ -306,9 +316,9 @@ namespace KinKal { } } if(!dok){ - // domain calculation failed. With ZeroFieldExtrap, keep a previously usable fit + // domain calculation failed. Under low-field protection, keep a previously usable fit // (map-edge truncation is preferred inside createDomains; this is a safety net). - if(config().zerofield_extrap_ && fitStatus().usable()) return; + if(bfield_.protecting() && fitStatus().usable()) return; history_.push_back(Status(0)); status().status_ = Status::outsidemap; status().comment_ = std::string("Extension error"); @@ -325,24 +335,27 @@ namespace KinKal { // replace domains when DomainWall correction is added or changed. the traj must also be replaced, so that // the pieces correspond to the new domains. The new traj is geometrically equivalent, but not parametrically equal. - // Build the replacement traj first so a failed append (e.g. CH charge flip) leaves domains_/effects_/fittraj_ intact. - template void Track::replaceDomains(DOMAINCOL const& domains) { + // Two-phase: build the replacement without touching any member, then commit. Returns false, leaving + // the track untouched, if a transformed piece cannot describe the same particle as those before it. + template bool Track::replaceDomains(DOMAINCOL const& domains) { + if(domains.size() == 0) return false; + TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end()); + TimeRange oldrange = fittraj_->range(); auto newtraj = std::make_unique(); // loop over domains, splitting the overlapping traj pieces at the domain walls, and transforming them to reference the domain's field // This increases the number of traj pieces. - // extend the existing traj to the domain range (restored on failure) - TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end()); - TimeRange front_range = fittraj_->front().range(); - TimeRange back_range = fittraj_->back().range(); - fittraj_->setRange(drange); - try { - for(auto const& domain : domains) { + for(auto const& domain : domains) { + using KTRAJPTR = std::shared_ptr; + using DKTRAJ = std::deque; + using DKTRAJCITER = typename DKTRAJ::const_iterator; + // clamp the lookup to the existing traj: domains may extend past it, and extending fittraj_ to + // reach them is the pre-mutation that used to need rollback. drange is restored on newtraj below. + double tlo = std::max(domain->begin(),oldrange.begin()); + double thi = std::min(domain->end(),oldrange.end()); + if(tlo < thi){ // find the range of existing ptraj pieces that overlaps with this domain's range - using KTRAJPTR = std::shared_ptr; - using DKTRAJ = std::deque; - using DKTRAJCITER = typename DKTRAJ::const_iterator; DKTRAJCITER first,last; - fittraj_->pieceRange(domain->range(),first,last); + fittraj_->pieceRange(TimeRange(tlo,thi),first,last); // loop over these pieces; first and last can be the same! auto olditer = first; do { @@ -353,17 +366,26 @@ namespace KinKal { double tstart = std::max(domain->begin(), oldpiece.range().begin()); double tend = std::min(domain->end(),oldpiece.range().end()); if(tstart < tend){ + // test only where the old code would actually have appended (and thrown) + if(!newtraj->compatible(newpiece)) return false; newpiece.range() = TimeRange(tstart,tend); newtraj->append(newpiece); } if(olditer != last)++olditer; } while(olditer != last); + } else { + // domain lies entirely outside the current traj: transform the nearest piece to cover it + auto const& oldpiece = fittraj_->nearestPiece(domain->range().mid()); + KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); + if(!newtraj->compatible(newpiece)) return false; + newpiece.range() = domain->range(); + newtraj->append(newpiece); } - } catch (...) { - fittraj_->front().setRange(front_range); - fittraj_->back().setRange(back_range); - throw; } + if(newtraj->pieces().size() == 0) return false; + // restore the full domain span on the replacement, matching the range the old code reached by + // extending fittraj_ up front + newtraj->setRange(drange); // commit: clear old domains / DomainWall effects, retarget remaining effects, swap traj if(domains_.size() > 0){ domains_.clear(); @@ -381,6 +403,7 @@ namespace KinKal { eff->updateReference(*newtraj); } fittraj_.swap(newtraj); + return true; } template void Track::extendTraj(DOMAINCOL const& domains ) { @@ -397,6 +420,30 @@ namespace KinKal { fittraj_->setRange(temprange); } + // no active hits or material crossings: there is nothing to fit, and a null range would otherwise be + // walked for domains and then set on the trajectory + template bool Track::validInput(TimeRange const& detrange) { + if(detrange.null()){ + history_.emplace_back(0,0,Status::lowNDOF, "No active hits or material crossings"); + return false; + } + return true; + } + + // as above, plus the seed itself must be finite: a NaN seed otherwise propagates into the + // parameterization and is only caught much later, if at all + template bool Track::validInput(TimeRange const& detrange, KTRAJ const& seedtraj) { + if(!validInput(detrange)) return false; + double tref = detrange.mid(); + auto spos = seedtraj.position3(tref); + if(!std::isfinite(seedtraj.momentum(tref)) || !std::isfinite(spos.X()) || + !std::isfinite(spos.Y()) || !std::isfinite(spos.Z())){ + history_.emplace_back(0,0,Status::failed, "Non-finite seed trajectory"); + return false; + } + return true; + } + template void Track::convertSeed(KTRAJ const& seedtraj,TimeRange const& range, DOMAINCOL& domains) { // if we're making local DomainWall corrections, divide the trajectory into domain pieces. Each will have equivalent parameters, but relative // to the local field @@ -416,11 +463,16 @@ namespace KinKal { auto bf = bfield_.fieldVect(seedtraj.position3(domain->mid())); KTRAJ newpiece(seedtraj,bf,domain->mid()); newpiece.range() = domain->range(); + // same routine incompatibility as replaceDomains; this append was previously unguarded, so a + // CentralHelix omega sign flip here threw std::invalid_argument out of the Track constructor + if(!fittraj_->compatible(newpiece)){ + history_.emplace_back(0,0,Status::incompatiblepiece, "Seed conversion: incompatible piece"); + return; + } fittraj_->append(newpiece); } - // Degenerate active range (or DomainMargin confinement dropping every domain) can leave - // createDomains with zero domains → empty fittraj_. createEffects would then throw - // std::length_error("Empty PiecewiseTrajectory!") and abort the art event. Soft-fail instead. + // zero domains leaves fittraj_ empty, which createEffects would report by throwing + // std::length_error out of the constructor; soft-fail instead if(fittraj_->pieces().empty()){ history_.emplace_back(0,0,Status::outsidemap, "Empty seed trajectory (no domains)"); return; @@ -702,17 +754,14 @@ namespace KinKal { while(time > fitrange.begin()){ auto const& ktraj = fittraj_->nearestPiece(time); double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); - // clamp the domain low bound to the active range minus domainmargin_ (max = unclamped/legacy) + // clamp the domain low bound to the active range minus domainmargin_ (max = unclamped) double dlo = std::max(time-dt, fitrange.begin() - config().domainmargin_); TimeRange range(dlo,time); - // sample BNom at the domain-midpoint piece only when confined (domainmargin_ set); else legacy piece (nearest time) + // sample BNom at the domain-midpoint piece when confined (domainmargin_ set), else at the nearest piece auto const& straj = (config().domainmargin_ < std::numeric_limits::max()) ? fittraj_->nearestPiece(range.mid()) : ktraj; Domain domain(range,bfield_.fieldVect(straj.position3(range.mid()))); addDomain(domain,TimeDir::backwards); time = domain.begin(); - // abort a runaway domain walk (caught by iterate()'s caller -> fit failed -> track dropped) - if(domains_.size() > config().maxdomains_) - throw std::runtime_error("Fit exceeded MaxDomains (BField domain walk runaway)"); } } // then forwards @@ -721,16 +770,14 @@ namespace KinKal { while(time < fitrange.end()){ auto const& ktraj = fittraj_->nearestPiece(time); double dt = std::max(bfield_.rangeInTolerance(ktraj,time,config().tol_),config().mindtstep_); - // clamp the domain high bound to the active range plus domainmargin_ (max = unclamped/legacy) + // clamp the domain high bound to the active range plus domainmargin_ (max = unclamped) double dhi = std::min(time+dt, fitrange.end() + config().domainmargin_); TimeRange range(time,dhi); - // sample BNom at the domain-midpoint piece only when confined (domainmargin_ set); else legacy piece (nearest time) + // sample BNom at the domain-midpoint piece when confined (domainmargin_ set), else at the nearest piece auto const& straj = (config().domainmargin_ < std::numeric_limits::max()) ? fittraj_->nearestPiece(range.mid()) : ktraj; Domain domain(range,bfield_.fieldVect(straj.position3(range.mid()))); addDomain(domain,TimeDir::forwards); time = domain.end(); - if(domains_.size() > config().maxdomains_) - throw std::runtime_error("Fit exceeded MaxDomains (BField domain walk runaway)"); } } } @@ -781,6 +828,12 @@ namespace KinKal { for(auto const& stat : history_) ost << stat << endl; } ost << " Fit Result "; + // convertSeed returns before fittraj_ is built when domain initialization fails, so a soft-failed + // fit has no trajectory to print; dereferencing it here segfaulted + if(!hasTraj()){ + ost << "(no trajectory: " << fitStatus().comment_ << ")" << endl; + return; + } fitTraj().print(ost,detail); if(detail > Config::basic) { ost << " Reference "; @@ -791,70 +844,58 @@ namespace KinKal { for(auto const& eff : effects()) eff.get()->print(ost,detail-3); } } + // build one domain starting at tstart, clipped to end no later than tend. Returns nothing when the + // field can't support a domain here; the map decides that, this just reports it. + template std::optional Track::createDomain(PKTRAJ const& ptraj, double tstart, double tend) const { + auto const& ktraj = ptraj.nearestPiece(tstart); + if(!bfield_.usable(ktraj.position3(tstart))) return std::nullopt; + double trange = bfield_.domainStep(ktraj,tstart,config().tol_,config().mindtstep_); + double dhi = std::min(tstart+trange,tend); + if(dhi <= tstart) return std::nullopt; + TimeRange drange(tstart,dhi); + // the domain carries the field sampled at its midpoint. Requiring that sample to be usable is what + // protection buys; unprotected, fieldVect just reports null outside the map, as it always did. + // Sample the midpoint on the midpoint's own piece only when confined (domainmargin_ set); the + // unconfined walk samples it on the start piece, as extendDomains does. + bool confined = config().domainmargin_ < std::numeric_limits::max(); + auto const& straj = confined ? ptraj.nearestPiece(drange.mid()) : ktraj; + VEC3 midpos = straj.position3(drange.mid()); + if(bfield_.protecting() && !bfield_.usable(midpos)) return std::nullopt; + return Domain(drange,bfield_.fieldVect(midpos)); + } + // divide a trajectory into magnetic 'domains' used to apply the DomainWall corrections template bool Track::createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const { - bool retval(true); - if(config().bfcorr_ ) { - // With ZeroFieldExtrap: stop DomainWalls at the map / B≈0 edge instead of throwing. - // Keeps domains built so far and lets Extension proceed (same gate as extrapolate handoff). - auto atZeroField = [&](VEC3 const& pos) -> bool { - if(!config().zerofield_extrap_) return false; - bool outside = !bfield_.inRange(pos); - VEC3 b = outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(pos); - return outside || BFieldMap::isZeroField(b) - || (config().minfield_ > 0.0 && b.R() < config().minfield_); - }; - // catch exceptions if the fit extends beyond the range of the field map - try { - if(config().domainmargin_ < std::numeric_limits::max()){ - // Confined (DomainMargin set): walk ONLY within active range ± margin. Do not use the - // legacy half-domain overhang past that window — for near-uniform B (cosmic CentralHelix) - // rangeInTolerance is huge, so begin-0.5*trange geometrically extrapolates far outside - // the hits / BField maps even though every hit is inside. - double const tlo = range.begin() - config().domainmargin_; - double const thi = range.end() + config().domainmargin_; - double tstart = tlo; - while(tstart < thi){ - auto const& ktraj = ptraj.nearestPiece(tstart); - if(atZeroField(ktraj.position3(tstart))) break; - double trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); - double dhi = std::min(tstart + trange, thi); - if(dhi > tstart){ - TimeRange drange(tstart,dhi); - auto const& straj = ptraj.nearestPiece(drange.mid()); - VEC3 midpos = straj.position3(drange.mid()); - if(atZeroField(midpos)) break; - domains.emplace(std::make_shared(drange,bfield_.fieldVect(midpos))); - } - tstart = dhi; - if(domains.size() > config().maxdomains_) - throw std::runtime_error("createDomains exceeded MaxDomains"); - } - } else { - // Legacy (default DomainMargin = max): half-domain overhang so the first/last effect sits - // mid-domain. Bit-identical to upstream when zerofield_extrap_ is false. - auto const& ktraj0 = ptraj.nearestPiece(range.begin()); - if(atZeroField(ktraj0.position3(range.begin()))) return true; - double trange = std::max(bfield_.rangeInTolerance(ktraj0,range.begin(),config().tol_),config().mindtstep_); - double tstart = range.begin() - 0.5*trange; - do { - auto const& ktraj = ptraj.nearestPiece(tstart); - if(atZeroField(ktraj.position3(tstart))) break; - trange = std::max(bfield_.rangeInTolerance(ktraj,tstart,config().tol_),config().mindtstep_); - VEC3 sample = ktraj.position3(tstart+0.5*trange); - if(atZeroField(sample)) break; - domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(sample))); - tstart += trange; - if(domains.size() > config().maxdomains_) - throw std::runtime_error("createDomains exceeded MaxDomains"); - } while(tstart < range.end() + 0.5*trange); - } - } catch (std::exception const& error) { - // ZeroFieldExtrap: treat unexpected map samples as a soft stop (keep domains so far). - if(!config().zerofield_extrap_) retval = false; + if(!config().bfcorr_) return true; + // No usable domain means: soft stop when protection is on (keep what we have, let the caller + // proceed), otherwise the failure the map used to signal by throwing out of fieldDeriv. + if(config().domainmargin_ < std::numeric_limits::max()){ + // confined (domainmargin_ set): walk only within the active range +/- margin, clipping the last domain + double const tlo = range.begin() - config().domainmargin_; + double const thi = range.end() + config().domainmargin_; + double tstart = tlo; + while(tstart < thi){ + auto domain = createDomain(ptraj,tstart,thi); + if(!domain) return bfield_.protecting(); + domains.emplace(std::make_shared(*domain)); + tstart = domain->end(); } + } else { + // Unconfined (default): half-domain overhang so the first/last effect sits mid-domain. The loop + // bound tracks the current step, as upstream, so domains are unclipped and full length. + auto const& ktraj0 = ptraj.nearestPiece(range.begin()); + if(!bfield_.usable(ktraj0.position3(range.begin()))) return bfield_.protecting(); + double trange = bfield_.domainStep(ktraj0,range.begin(),config().tol_,config().mindtstep_); + double tstart = range.begin() - 0.5*trange; + do { + auto domain = createDomain(ptraj,tstart,std::numeric_limits::max()); + if(!domain) return bfield_.protecting(); + trange = domain->range().range(); + domains.emplace(std::make_shared(*domain)); + tstart = domain->end(); + } while(tstart < range.end() + 0.5*trange); } - return retval; + return true; } template TimeRange Track::detectorRange(HITCOL& hits, EXINGCOL& exings,bool active) { @@ -882,9 +923,9 @@ namespace KinKal { bool retval = fitStatus().usable(); if(retval){ if(config().bfcorr_){ - // Opt-in zero-field handoff (zerofield_extrap_): for CH cosmic CRV extrapolation that must - // leave the map. Default false → legacy bfcorr domain walk unchanged (LH-safe). - if(config().zerofield_extrap_){ + // opt-in low-field handoff, for extrapolation that must leave the field map; with minfield_ 0 + // the unprotected domain walk below is used unchanged + if(bfield_.protecting()){ auto geometricExtend = [&](double tmax_remaining) { if(tmax_remaining <= 0.0) return; auto& endpiece = tdir == TimeDir::forwards ? fittraj_->backPtr() : fittraj_->frontPtr(); @@ -910,15 +951,11 @@ namespace KinKal { auto const& ktraj = fittraj_->nearestPiece(time); if( !std::isfinite(ktraj.momentum(time)) ) break; - VEC3 frontier = ktraj.position3(time); - bool outside = !bfield_.inRange(frontier); - // Only sample the map when inside it — never call fieldDeriv out of range - VEC3 bfront = outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(frontier); - bool zerofield = outside || BFieldMap::isZeroField(bfront) - || (config().minfield_ > 0.0 && bfront.R() < config().minfield_); - if(zerofield){ - // Leave the bfcorr / DomainWall path. Free-particle continuation is geometric - // range-extend of the current end piece (no CH rebuild at B≈0). + // the map decides whether this point can carry field-corrected transport; asking it first + // also keeps rangeInTolerance from sampling fieldDeriv out of range + if(!bfield_.usable(ktraj.position3(time))){ + // leave the bfcorr / DomainWall path; free-particle continuation is a geometric + // range-extend of the current end piece, with no parameter rebuild at B ~ 0 handed_off = true; break; } @@ -926,14 +963,11 @@ namespace KinKal { double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); VEC3 midpos = ktraj.position3(range.mid()); - bool mid_outside = !bfield_.inRange(midpos); - VEC3 domainfield = mid_outside ? VEC3(0.0,0.0,0.0) : bfield_.fieldVect(midpos); - if(mid_outside || BFieldMap::isZeroField(domainfield) - || (config().minfield_ > 0.0 && domainfield.R() < config().minfield_)){ + if(!bfield_.usable(midpos)){ handed_off = true; break; } - Domain domain(range,domainfield); + Domain domain(range,bfield_.fieldVect(midpos)); addDomain(domain,tdir,true); time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); } @@ -947,7 +981,7 @@ namespace KinKal { geometricExtend(xtest.maxDt() - fabs(time-tstart)); } } else { - // Legacy bfcorr extrapolation (default): unchanged domain walk + // no low-field protection: the unprotected bfcorr domain walk, unchanged try { double time = tdir == TimeDir::forwards ? domains_.crbegin()->get()->end() : domains_.cbegin()->get()->begin(); double tstart = time; @@ -956,7 +990,6 @@ namespace KinKal { double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); auto domainfield = bfield_.fieldVect(ktraj.position3(range.mid())); - if( config().minfield_ > 0.0 && domainfield.R() < config().minfield_ ) break; Domain domain(range,domainfield); addDomain(domain,tdir,true); time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); diff --git a/General/BFieldMap.hh b/General/BFieldMap.hh index 0e9b1603..4db89430 100644 --- a/General/BFieldMap.hh +++ b/General/BFieldMap.hh @@ -11,6 +11,7 @@ #include #include #include +#include #include namespace KinKal { @@ -27,24 +28,144 @@ namespace KinKal { virtual bool inRange(VEC3 const& position) const = 0; virtual ~BFieldMap(){} virtual void print(std::ostream& os ) const = 0; - BFieldMap(){} + // smallest |B| (T) usable for field-corrected transport; 0 (default) disables low-field + // protection, leaving an unusable sample a hard failure for the caller + BFieldMap(double minfield=0.0) : minfield_(minfield) {} + double minField() const { return minfield_; } + bool protecting() const { return minfield_ > 0.0; } // is low-field protection enabled? + // Can this position support field-corrected transport? Answered from the region computed once at + // construction, so this costs no field evaluation; maps registering no region use the direct test. + bool usable(VEC3 const& position) const { + if(!region_.built_) return directlyUsable(position); + return region_.usable(position); + } + // pre-registered value returned for positions outside the usable region; no lookup needed + VEC3 const& invalidField() const { return invalidfield_; } + void printUsableRegion(std::ostream& os) const { + if(!region_.built_){ + os << "BField usable region: none precomputed (direct test, minField " << minfield_ << " T)" << std::endl; + } else { + size_t ncell = region_.usable_.size(); + os << "BField usable region: minField " << minfield_ << " T, cylindrical r<" << region_.rmax_ + << " z[" << region_.zlow_ << "," << region_.zlow_+region_.zdim_ << "] mm, cell " << region_.cell_ + << " mm x " << region_.nphi_ << " phi bins, " << region_.nr_ << "x" << region_.nphi_ << "x" + << region_.nz_ << " = " << ncell << " cells, " << region_.nsampled_ << " field samples, " + << region_.nusable_ << " usable (" << (ncell>0 ? 100.0*region_.nusable_/ncell : 0.0) << "%)" + << std::endl; + } + } // disallow copy and equivalence BFieldMap(BFieldMap const& ) = delete; BFieldMap& operator =(BFieldMap const& ) = delete; // speed of light in units to convert Tesla to mm (bending radius) static double constexpr cbar() { return CLHEP::c_light/1000.0; } - // |B| below this (T) is treated as physically zero for ZeroFieldExtrap handoff decisions. - // Fit paths must not sample here; with zerofield_extrap_, extrapolation leaves the bfcorr - // domain walk and continues by geometric range-extend (no CH rebuild at B≈0). + // |B| below this (T) is treated as physically zero when deciding to hand off extrapolation; + // fit paths must not sample here static double constexpr zeroField() { return 1.0e-6; } static bool isZeroField(VEC3 const& bvec) { return bvec.R() < zeroField(); } // templated interface for interacting with kinematic trajectory classes // how far can you go along the given kinematic trajectory till BField inhomogeneity makes the momentum accuracy out of (fractional) tolerance template double rangeInTolerance(KTRAJ const& ktraj, double tstart, double tol) const; + // the domain step at tstart: rangeInTolerance with a floor applied. Callers must have checked + // usable(position) first -- rangeInTolerance samples fieldDeriv, which is only defined in range. + template double domainStep(KTRAJ const& ktraj, double tstart, double tol, double mindtstep) const { + return std::max(rangeInTolerance(ktraj,tstart,tol),mindtstep); + } // integrate the residual magentic force over the given kinematic trajectory and range due to the difference between the true field and the nominal field in the template VEC3 integrate(KTRAJ const& ktraj, TimeRange const& trange) const; + protected: + // Map the usable region once, at construction, on a cylindrical lattice about the z axis; called by + // subclasses once their extents are known. nphi need only follow the departure from axial symmetry. + void buildUsableRegion(double rmax, double zlow, double zhigh, double cellsize, unsigned nphi); + // True only when the precomputed region positively excludes this point, so a map can return its + // pre-registered invalid value without a lookup. Must not call usable(): that would recurse. + bool preRegisteredInvalid(VEC3 const& position) const { + return region_.built_ && !region_.usable(position); + } + // the direct, per-sample test. Used to build the region, and as the fallback for maps that + // register none. + bool directlyUsable(VEC3 const& position) const { + if(!inRange(position)) return false; // outside: fieldDeriv is undefined here + if(!protecting()) return true; // unprotected: in-range was the only requirement + return fieldVect(position).R() >= std::max(minfield_,zeroField()); + } + private: + double minfield_; // smallest usable |B| (T); 0 disables low-field protection + VEC3 invalidfield_ = VEC3(0.0,0.0,0.0); // pre-registered return outside the usable region + // Precomputed occupancy of the usable region, in cylindrical coordinates about the z axis: for a + // near-axially-symmetric field the boundary is close to r = rmax(phi,z), so azimuth needs few bins. + struct UsableRegion { + bool built_ = false; + double rmax_ = 0.0, zlow_ = 0.0, zdim_ = 0.0; + double cell_ = 0.0; // radial and axial bin size + double dphi_ = 0.0; // azimuthal bin size + unsigned nr_ = 0, nphi_ = 0, nz_ = 0; + std::vector usable_; + size_t nsampled_ = 0; // field samples taken while building (the run-start cost) + size_t nusable_ = 0; + bool usable(VEC3 const& p) const { + double z = p.Z()-zlow_; + if(z < 0.0 || z >= zdim_) return false; + double r = std::sqrt(p.X()*p.X()+p.Y()*p.Y()); + if(r >= rmax_) return false; + double phi = std::atan2(p.Y(),p.X()); + if(phi < 0.0) phi += 2.0*M_PI; + unsigned ir = static_cast(r/cell_); + unsigned iz = static_cast(z/cell_); + unsigned ip = static_cast(phi/dphi_); + if(ip >= nphi_) ip = nphi_-1; // guard the wrap at exactly 2pi + return usable_[(static_cast(ir)*nphi_ + ip)*nz_ + iz]; + } + }; + UsableRegion region_; }; + // Corners are sampled once each and a cell is marked usable only if all eight are, so the recorded + // region is conservative: it never claims space the field cannot support. + inline void BFieldMap::buildUsableRegion(double rmax, double zlow, double zhigh, double cellsize, unsigned nphi) { + auto& r = region_; + // a non-positive cell (or radius) cannot describe a region: leave the direct test in place rather + // than building a degenerate one that would report everything unusable + if(cellsize <= 0.0 || rmax <= 0.0 || zhigh <= zlow) return; + r.rmax_ = rmax; r.zlow_ = zlow; r.zdim_ = zhigh-zlow; r.cell_ = cellsize; + r.nr_ = std::max(1u,static_cast(std::ceil(rmax/cellsize))); + r.nz_ = std::max(1u,static_cast(std::ceil(r.zdim_/cellsize))); + r.nphi_ = std::max(1u,nphi); + r.dphi_ = 2.0*M_PI/r.nphi_; + r.rmax_ = r.nr_*cellsize; r.zdim_ = r.nz_*cellsize; + // corner lattice: (nr+1) x nphi x (nz+1), azimuth periodic so no extra plane + unsigned mr = r.nr_+1, mz = r.nz_+1; + std::vector node(static_cast(mr)*r.nphi_*mz,false); + for(unsigned ir=0; ir(ir)*r.nphi_ + ip)*mz + iz] = directlyUsable(pos); + } + } + } + r.nsampled_ = node.size(); + r.usable_.assign(static_cast(r.nr_)*r.nphi_*r.nz_,false); + for(unsigned ir=0; ir(ir+dr)*r.nphi_ + (dp ? ipn : ip))*mz + (iz+dz)]; + r.usable_[(static_cast(ir)*r.nphi_ + ip)*r.nz_ + iz] = ok; + } + } + } + r.nusable_ = std::count(r.usable_.begin(),r.usable_.end(),true); + r.built_ = true; + } + template VEC3 BFieldMap::integrate(KTRAJ const& ktraj, TimeRange const& trange) const { // take a fixed number of steps. This may fail for long ranges FIXME! unsigned nsteps(10); diff --git a/Trajectory/ParticleTrajectory.hh b/Trajectory/ParticleTrajectory.hh index d4743f86..7b89d83c 100644 --- a/Trajectory/ParticleTrajectory.hh +++ b/Trajectory/ParticleTrajectory.hh @@ -18,15 +18,19 @@ namespace KinKal { // construct from an initial piece, which also provides kinematic information ParticleTrajectory(KTRAJ const& piece) : PTTRAJ(piece) {} ParticleTrajectory() : PTTRAJ() {} + // does this piece describe the same particle (mass and charge) as the existing ones? Callers that + // expect incompatibility as a routine outcome test this instead of catching append's throw + bool compatible(KTRAJ const& newpiece) const { + return PTTRAJ::pieces().size() == 0 || + (fabs(newpiece.mass()-mass()) <= 1e-6 && newpiece.charge() == charge()); + } // append and prepend to check mass and charge consistency void append(KTRAJ const& newpiece, bool allowremove=false) { - if(PTTRAJ::pieces().size() > 0){ - if(fabs(newpiece.mass()-mass())>1e-6 || newpiece.charge() != charge()) throw std::invalid_argument("Invalid particle parameters"); - } + if(!compatible(newpiece)) throw std::invalid_argument("Invalid particle parameters"); PTTRAJ::append(newpiece,allowremove); } void prepend(KTRAJ const& newpiece, bool allowremove=false) { - if(fabs(newpiece.mass()-mass())>1e-6 || newpiece.charge() != charge()) throw std::invalid_argument("Invalid particle parameters"); + if(!compatible(newpiece)) throw std::invalid_argument("Invalid particle parameters"); PTTRAJ::prepend(newpiece,allowremove); } // kinematic interface From f68cb765b463c090eca9bd7f6c62d2ea36076a90 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Thu, 20 Aug 2026 13:39:46 -0500 Subject: [PATCH 6/6] Replace voxel field pre-registering with simpler check. Add guards for invalid field in CentralHelix ctor. --- Fit/Track.hh | 108 ++++++++++++++++++-------------- General/BFieldMap.hh | 123 ++++--------------------------------- Trajectory/CentralHelix.cc | 9 +++ 3 files changed, 81 insertions(+), 159 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index f070b9bc..f3652817 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -293,6 +293,9 @@ namespace KinKal { // a tighter extension tolerance can flip omega near a collapsing field, which the // parameterization reads as a charge change; leave the track untouched and record why if(!replaceDomains(domains)){ + // A fit that already converged is not invalidated by the EXTENSION failing to re-domain it: + // keep it, as the low-field handoff did before this returned a status instead of throwing. + if(bfield_.protecting() && fitStatus().usable()) return; history_.push_back(Status(0)); status().status_ = Status::incompatiblepiece; status().comment_ = std::string("Domain replacement: incompatible piece"); @@ -340,52 +343,45 @@ namespace KinKal { template bool Track::replaceDomains(DOMAINCOL const& domains) { if(domains.size() == 0) return false; TimeRange drange(domains.begin()->get()->begin(),domains.rbegin()->get()->end()); - TimeRange oldrange = fittraj_->range(); auto newtraj = std::make_unique(); + // Split from a COPY extended to the domain span. The original code extended fittraj_ itself and + // rolled its end-piece ranges back on failure; copying keeps the piece splitting identical to that + // (the copy ctor deep-copies every piece) while leaving fittraj_ untouched until the commit below. + PKTRAJ srctraj(*fittraj_); + // setRange throws when the domain span is disjoint from the trajectory. Here that is a routine + // outcome, not an error -- report it like any other failure to re-domain, so the caller can keep a + // fit that already converged, instead of the exception unwinding out to the module and losing it. + if(drange.begin() > srctraj.front().range().end() || + drange.end() < srctraj.back().range().begin()) return false; + srctraj.setRange(drange); // loop over domains, splitting the overlapping traj pieces at the domain walls, and transforming them to reference the domain's field // This increases the number of traj pieces. for(auto const& domain : domains) { using KTRAJPTR = std::shared_ptr; using DKTRAJ = std::deque; using DKTRAJCITER = typename DKTRAJ::const_iterator; - // clamp the lookup to the existing traj: domains may extend past it, and extending fittraj_ to - // reach them is the pre-mutation that used to need rollback. drange is restored on newtraj below. - double tlo = std::max(domain->begin(),oldrange.begin()); - double thi = std::min(domain->end(),oldrange.end()); - if(tlo < thi){ - // find the range of existing ptraj pieces that overlaps with this domain's range - DKTRAJCITER first,last; - fittraj_->pieceRange(TimeRange(tlo,thi),first,last); - // loop over these pieces; first and last can be the same! - auto olditer = first; - do { - auto const& oldpiece = **olditer; - // copy this piece, translating bnom to this domain's field - KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); - // set the range for this piece, making sure it is non-zero - double tstart = std::max(domain->begin(), oldpiece.range().begin()); - double tend = std::min(domain->end(),oldpiece.range().end()); - if(tstart < tend){ - // test only where the old code would actually have appended (and thrown) - if(!newtraj->compatible(newpiece)) return false; - newpiece.range() = TimeRange(tstart,tend); - newtraj->append(newpiece); - } - if(olditer != last)++olditer; - } while(olditer != last); - } else { - // domain lies entirely outside the current traj: transform the nearest piece to cover it - auto const& oldpiece = fittraj_->nearestPiece(domain->range().mid()); + // find the range of existing ptraj pieces that overlaps with this domain's range + DKTRAJCITER first,last; + srctraj.pieceRange(domain->range(),first,last); + // loop over these pieces; first and last can be the same! + auto olditer = first; + do { + auto const& oldpiece = **olditer; + // copy this piece, translating bnom to this domain's field KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); - if(!newtraj->compatible(newpiece)) return false; - newpiece.range() = domain->range(); - newtraj->append(newpiece); - } + // set the range for this piece, making sure it is non-zero + double tstart = std::max(domain->begin(), oldpiece.range().begin()); + double tend = std::min(domain->end(),oldpiece.range().end()); + if(tstart < tend){ + // test only where the old code would actually have appended (and thrown) + if(!newtraj->compatible(newpiece)) return false; + newpiece.range() = TimeRange(tstart,tend); + newtraj->append(newpiece); + } + if(olditer != last)++olditer; + } while(olditer != last); } if(newtraj->pieces().size() == 0) return false; - // restore the full domain span on the replacement, matching the range the old code reached by - // extending fittraj_ up front - newtraj->setRange(drange); // commit: clear old domains / DomainWall effects, retarget remaining effects, swap traj if(domains_.size() > 0){ domains_.clear(); @@ -460,8 +456,15 @@ namespace KinKal { fittraj_ = std::make_unique(); for(auto const& domain : domains) { // Set the DomainWall to the start of this domain - auto bf = bfield_.fieldVect(seedtraj.position3(domain->mid())); - KTRAJ newpiece(seedtraj,bf,domain->mid()); + // the domain walk should never hand back a domain whose midpoint has no usable field, but a + // null bnom makes the parameterization degenerate (omega -> signed zero, momentum 0/0), so + // never construct a piece from an unchecked sample + auto bf = bfield_.usableField(seedtraj.position3(domain->mid())); + if(!bf){ + history_.emplace_back(0,0,Status::outsidemap, "Seed conversion: unusable field at domain"); + return; + } + KTRAJ newpiece(seedtraj,*bf,domain->mid()); newpiece.range() = domain->range(); // same routine incompatibility as replaceDomains; this append was previously unguarded, so a // CentralHelix omega sign flip here threw std::invalid_argument out of the Track constructor @@ -480,9 +483,13 @@ namespace KinKal { } else { // use the middle of the range as the nominal BField for this fit: double tref = range.mid(); - VEC3 bf = bfield_.fieldVect(seedtraj.position3(tref)); + auto bf = bfield_.usableField(seedtraj.position3(tref)); + if(!bf){ + history_.emplace_back(0,0,Status::outsidemap, "Seed conversion: unusable field at reference"); + return; + } // create the first piece. Note this constructor adjusts the parameters according to the local field - KTRAJ firstpiece(seedtraj,bf,tref); + KTRAJ firstpiece(seedtraj,*bf,tref); firstpiece.range() = range; // create the piecewise trajectory from this fittraj_ = std::make_unique(firstpiece); @@ -738,7 +745,10 @@ namespace KinKal { template void Track::updateDomains(PKTRAJ const& ptraj) { for(auto& domain : domains_) { - domain->updateBNom(bfield_.fieldVect(ptraj.position3(domain->mid()))); + // an unusable sample would zero the domain's BNom and make every piece built from it degenerate; + // keep the field the domain was created with instead + auto bf = bfield_.usableField(ptraj.position3(domain->mid())); + if(bf) domain->updateBNom(*bf); } } @@ -858,10 +868,13 @@ namespace KinKal { // Sample the midpoint on the midpoint's own piece only when confined (domainmargin_ set); the // unconfined walk samples it on the start piece, as extendDomains does. bool confined = config().domainmargin_ < std::numeric_limits::max(); - auto const& straj = confined ? ptraj.nearestPiece(drange.mid()) : ktraj; - VEC3 midpos = straj.position3(drange.mid()); - if(bfield_.protecting() && !bfield_.usable(midpos)) return std::nullopt; - return Domain(drange,bfield_.fieldVect(midpos)); + double tmid = confined ? drange.mid() : tstart + 0.5*trange; + auto const& straj = confined ? ptraj.nearestPiece(tmid) : ktraj; + VEC3 midpos = straj.position3(tmid); + // one interpolation gives both the test and the domain's BNom + auto midfield = bfield_.usableField(midpos); + if(bfield_.protecting() && !midfield) return std::nullopt; + return Domain(drange,midfield ? *midfield : bfield_.fieldVect(midpos)); } // divide a trajectory into magnetic 'domains' used to apply the DomainWall corrections @@ -963,11 +976,12 @@ namespace KinKal { double dt = std::clamp(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),config().mindtstep_,xtest.maxDtStep()); TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); VEC3 midpos = ktraj.position3(range.mid()); - if(!bfield_.usable(midpos)){ + auto midfield = bfield_.usableField(midpos); // one interpolation: test AND the BNom below + if(!midfield){ handed_off = true; break; } - Domain domain(range,bfield_.fieldVect(midpos)); + Domain domain(range,*midfield); addDomain(domain,tdir,true); time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); } diff --git a/General/BFieldMap.hh b/General/BFieldMap.hh index 4db89430..f152ef40 100644 --- a/General/BFieldMap.hh +++ b/General/BFieldMap.hh @@ -11,7 +11,7 @@ #include #include #include -#include +#include #include namespace KinKal { @@ -33,27 +33,15 @@ namespace KinKal { BFieldMap(double minfield=0.0) : minfield_(minfield) {} double minField() const { return minfield_; } bool protecting() const { return minfield_ > 0.0; } // is low-field protection enabled? - // Can this position support field-corrected transport? Answered from the region computed once at - // construction, so this costs no field evaluation; maps registering no region use the direct test. - bool usable(VEC3 const& position) const { - if(!region_.built_) return directlyUsable(position); - return region_.usable(position); - } - // pre-registered value returned for positions outside the usable region; no lookup needed - VEC3 const& invalidField() const { return invalidfield_; } - void printUsableRegion(std::ostream& os) const { - if(!region_.built_){ - os << "BField usable region: none precomputed (direct test, minField " << minfield_ << " T)" << std::endl; - } else { - size_t ncell = region_.usable_.size(); - os << "BField usable region: minField " << minfield_ << " T, cylindrical r<" << region_.rmax_ - << " z[" << region_.zlow_ << "," << region_.zlow_+region_.zdim_ << "] mm, cell " << region_.cell_ - << " mm x " << region_.nphi_ << " phi bins, " << region_.nr_ << "x" << region_.nphi_ << "x" - << region_.nz_ << " = " << ncell << " cells, " << region_.nsampled_ << " field samples, " - << region_.nusable_ << " usable (" << (ncell>0 ? 100.0*region_.nusable_/ncell : 0.0) << "%)" - << std::endl; - } + // The field here, if this position can support field-corrected transport; nothing if it cannot. + std::optional usableField(VEC3 const& position) const { + if(!inRange(position)) return std::nullopt; // outside: fieldDeriv is undefined here + VEC3 bf = fieldVect(position); + if(protecting() && bf.R() < std::max(minfield_,zeroField())) return std::nullopt; + return bf; } + // predicate form, for the callers that only decide and never use the field + bool usable(VEC3 const& position) const { return usableField(position).has_value(); } // disallow copy and equivalence BFieldMap(BFieldMap const& ) = delete; BFieldMap& operator =(BFieldMap const& ) = delete; @@ -66,106 +54,17 @@ namespace KinKal { // templated interface for interacting with kinematic trajectory classes // how far can you go along the given kinematic trajectory till BField inhomogeneity makes the momentum accuracy out of (fractional) tolerance template double rangeInTolerance(KTRAJ const& ktraj, double tstart, double tol) const; - // the domain step at tstart: rangeInTolerance with a floor applied. Callers must have checked - // usable(position) first -- rangeInTolerance samples fieldDeriv, which is only defined in range. + // the domain step at tstart: rangeInTolerance with a floor applied. Callers must have confirmed + // the position is in range first -- rangeInTolerance samples fieldDeriv, undefined outside it. template double domainStep(KTRAJ const& ktraj, double tstart, double tol, double mindtstep) const { return std::max(rangeInTolerance(ktraj,tstart,tol),mindtstep); } // integrate the residual magentic force over the given kinematic trajectory and range due to the difference between the true field and the nominal field in the template VEC3 integrate(KTRAJ const& ktraj, TimeRange const& trange) const; - protected: - // Map the usable region once, at construction, on a cylindrical lattice about the z axis; called by - // subclasses once their extents are known. nphi need only follow the departure from axial symmetry. - void buildUsableRegion(double rmax, double zlow, double zhigh, double cellsize, unsigned nphi); - // True only when the precomputed region positively excludes this point, so a map can return its - // pre-registered invalid value without a lookup. Must not call usable(): that would recurse. - bool preRegisteredInvalid(VEC3 const& position) const { - return region_.built_ && !region_.usable(position); - } - // the direct, per-sample test. Used to build the region, and as the fallback for maps that - // register none. - bool directlyUsable(VEC3 const& position) const { - if(!inRange(position)) return false; // outside: fieldDeriv is undefined here - if(!protecting()) return true; // unprotected: in-range was the only requirement - return fieldVect(position).R() >= std::max(minfield_,zeroField()); - } private: double minfield_; // smallest usable |B| (T); 0 disables low-field protection - VEC3 invalidfield_ = VEC3(0.0,0.0,0.0); // pre-registered return outside the usable region - // Precomputed occupancy of the usable region, in cylindrical coordinates about the z axis: for a - // near-axially-symmetric field the boundary is close to r = rmax(phi,z), so azimuth needs few bins. - struct UsableRegion { - bool built_ = false; - double rmax_ = 0.0, zlow_ = 0.0, zdim_ = 0.0; - double cell_ = 0.0; // radial and axial bin size - double dphi_ = 0.0; // azimuthal bin size - unsigned nr_ = 0, nphi_ = 0, nz_ = 0; - std::vector usable_; - size_t nsampled_ = 0; // field samples taken while building (the run-start cost) - size_t nusable_ = 0; - bool usable(VEC3 const& p) const { - double z = p.Z()-zlow_; - if(z < 0.0 || z >= zdim_) return false; - double r = std::sqrt(p.X()*p.X()+p.Y()*p.Y()); - if(r >= rmax_) return false; - double phi = std::atan2(p.Y(),p.X()); - if(phi < 0.0) phi += 2.0*M_PI; - unsigned ir = static_cast(r/cell_); - unsigned iz = static_cast(z/cell_); - unsigned ip = static_cast(phi/dphi_); - if(ip >= nphi_) ip = nphi_-1; // guard the wrap at exactly 2pi - return usable_[(static_cast(ir)*nphi_ + ip)*nz_ + iz]; - } - }; - UsableRegion region_; }; - // Corners are sampled once each and a cell is marked usable only if all eight are, so the recorded - // region is conservative: it never claims space the field cannot support. - inline void BFieldMap::buildUsableRegion(double rmax, double zlow, double zhigh, double cellsize, unsigned nphi) { - auto& r = region_; - // a non-positive cell (or radius) cannot describe a region: leave the direct test in place rather - // than building a degenerate one that would report everything unusable - if(cellsize <= 0.0 || rmax <= 0.0 || zhigh <= zlow) return; - r.rmax_ = rmax; r.zlow_ = zlow; r.zdim_ = zhigh-zlow; r.cell_ = cellsize; - r.nr_ = std::max(1u,static_cast(std::ceil(rmax/cellsize))); - r.nz_ = std::max(1u,static_cast(std::ceil(r.zdim_/cellsize))); - r.nphi_ = std::max(1u,nphi); - r.dphi_ = 2.0*M_PI/r.nphi_; - r.rmax_ = r.nr_*cellsize; r.zdim_ = r.nz_*cellsize; - // corner lattice: (nr+1) x nphi x (nz+1), azimuth periodic so no extra plane - unsigned mr = r.nr_+1, mz = r.nz_+1; - std::vector node(static_cast(mr)*r.nphi_*mz,false); - for(unsigned ir=0; ir(ir)*r.nphi_ + ip)*mz + iz] = directlyUsable(pos); - } - } - } - r.nsampled_ = node.size(); - r.usable_.assign(static_cast(r.nr_)*r.nphi_*r.nz_,false); - for(unsigned ir=0; ir(ir+dr)*r.nphi_ + (dp ? ipn : ip))*mz + (iz+dz)]; - r.usable_[(static_cast(ir)*r.nphi_ + ip)*r.nz_ + iz] = ok; - } - } - } - r.nusable_ = std::count(r.usable_.begin(),r.usable_.end(),true); - r.built_ = true; - } - template VEC3 BFieldMap::integrate(KTRAJ const& ktraj, TimeRange const& trange) const { // take a fixed number of steps. This may fail for long ranges FIXME! unsigned nsteps(10); diff --git a/Trajectory/CentralHelix.cc b/Trajectory/CentralHelix.cc index 3bf09854..d61238bc 100644 --- a/Trajectory/CentralHelix.cc +++ b/Trajectory/CentralHelix.cc @@ -32,6 +32,12 @@ namespace KinKal { CentralHelix::CentralHelix(VEC4 const &pos0, MOM4 const &mom0, int charge, VEC3 const &bnom, TimeRange const &trange) : trange_(trange), mass_(mom0.M()), bnom_(bnom) { + // A null nominal field has no valid parameterization here: radToMom below is 0, so omega comes out + // as signed zero and momentum() as 0/0 = NaN, while charge() -- which reads the sign of omega -- + // becomes meaningless. Refuse rather than return a silently degenerate object. Callers that can meet + // a field-free region test BFieldMap::usableField() first; this is the backstop behind them. + if(bnom_.R() < BFieldMap::zeroField()) + throw std::invalid_argument("CentralHelix::CentralHelix; null BNom"); // Transform into the system where Z is along the Bfield. This is a pure rotation about the origin VEC4 pos(pos0); MOM4 mom(mom0); @@ -85,6 +91,9 @@ namespace KinKal { } void CentralHelix::resetBNom(VEC3 const& bnom) { + // same degeneracy as construction: refuse to move an existing piece onto a null field + if(bnom.R() < BFieldMap::zeroField()) + throw std::invalid_argument("CentralHelix::resetBNom; null BNom"); bnom_ = bnom; setTransforms(); }