From 0c42eeee76cd76f32c35c1c94af935f5fdea98ef Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Tue, 26 May 2026 13:43:41 -0700 Subject: [PATCH 1/2] Add extrapolation predicate processing for non-bfield (ie line) trajectory tracks --- Fit/Track.hh | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Fit/Track.hh b/Fit/Track.hh index 237f2982..521a78b3 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -814,6 +814,15 @@ namespace KinKal { } retval = true; } else { + // geometric extrapolation. Extend to the maximum + auto& endpiece = tdir == TimeDir::forwards ? fittraj_->backPtr() : fittraj_->frontPtr(); + TimeRange newrange = tdir == TimeDir::forwards ? + TimeRange(endpiece->range().begin(),endpiece->range().end()+xtest.maxDt()) + : + TimeRange(endpiece->range().begin()-xtest.maxDt(),endpiece->range().end()); + endpiece->setRange(newrange); + // invoke the test + xtest.needsExtrapolation(*fittraj_,tdir); retval = true; } } From 5ca92df188d919ed551eddc4d9b6bd762be05893 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Tue, 26 May 2026 18:16:21 -0700 Subject: [PATCH 2/2] Fixes for extrapolation --- Fit/Track.hh | 25 +++++++++++++++---------- General/TimeRange.hh | 15 +++++++++++++-- 2 files changed, 28 insertions(+), 12 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index 521a78b3..aeb7d83e 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -800,7 +800,7 @@ 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 = bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()); // always positive + double dt = std::min(bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()),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()))); addDomain(domain,tdir,true); // use exact transport @@ -814,16 +814,21 @@ namespace KinKal { } retval = true; } else { - // geometric extrapolation. Extend to the maximum + // geometric extrapolation of the end piece; no need to protect auto& endpiece = tdir == TimeDir::forwards ? fittraj_->backPtr() : fittraj_->frontPtr(); - TimeRange newrange = tdir == TimeDir::forwards ? - TimeRange(endpiece->range().begin(),endpiece->range().end()+xtest.maxDt()) - : - TimeRange(endpiece->range().begin()-xtest.maxDt(),endpiece->range().end()); - endpiece->setRange(newrange); - // invoke the test - xtest.needsExtrapolation(*fittraj_,tdir); - retval = true; + double time = tdir == TimeDir::forwards ? endpiece->range().end() : endpiece->range().begin(); + 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()) + : + 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) < xtest.maxDt()); } } return retval; diff --git a/General/TimeRange.hh b/General/TimeRange.hh index 46875a82..c733efdc 100644 --- a/General/TimeRange.hh +++ b/General/TimeRange.hh @@ -11,8 +11,15 @@ namespace KinKal { class TimeRange { public: TimeRange() : range_{0.0,0.0} {} // default to a null range (matches no times) - TimeRange(double begin, double end) : range_{begin,end} { - if(begin > end)throw std::invalid_argument("Invalid Time Range"); } + TimeRange(double begin, double end,bool ordered=true) : range_{begin,end} { + if(begin > end){ + if(ordered){ + throw std::invalid_argument("Invalid Time Range"); + } else { + range_[0] = end; range_[1] = begin; + } + } + } double begin() const { return range_[0]; } double end() const { return range_[1]; } double rbegin() const { return range_[1]; } @@ -27,6 +34,10 @@ namespace KinKal { return (begin() <= other.begin() && end() >= other.end()); } // force time to be in range void forceRange(double& time) const { time = std::min(std::max(time,begin()),end()); } + // test if a given time is beyond the range in the given direction + bool beyond(double time, TimeDir tdir) const { + return tdir == TimeDir::forwards ? time > range_[1] : time < range_[0]; + } // augment another range void combine(TimeRange const& other ) { range_[0] = std::min(begin(),other.begin());