diff --git a/Fit/Track.hh b/Fit/Track.hh index 237f2982..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,7 +814,21 @@ namespace KinKal { } retval = true; } else { - retval = true; + // geometric extrapolation of the end piece; no need to protect + 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 { + // 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());