Support KinematicLine extrapolation

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;

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());