Fixes

Fit/Status.hh

@@ -16,6 +16,7 @@ namespace KinKal {
     Chisq chisq_; // current chisquared
     std::string comment_; // further information about the status
     bool usable() const { return status_ > unfit && status_ < lowNDOF; }
+    bool unusable() const { return !usable(); }
     bool needsFit() const { return status_ == unfit || status_ == unconverged; }
     Status(unsigned miter,unsigned iter=0,status stat=unfit,const char* comment="") : miter_(miter), iter_(iter), status_(stat), chisq_(NParams()),comment_(comment){}
     static std::string statusName(status stat);

Fit/Track.hh

@@ -159,7 +159,7 @@ namespace KinKal {
       KKEFFCOL effects_; // effects used in this fit, sorted by time
       HITCOL hits_; // hits used in this fit
       EXINGCOL exings_; // material xings used in this fit
-      DOMAINCOL domains_; // DomainWall domains used in this fit, contiguous and sorted by time
+      DOMAINCOL domains_; // domains used in this fit, contiguous and sorted by time
   };
   // minimal constructor for subclasses. The resulting object has no fit.
   template <class KTRAJ> Track<KTRAJ>::Track(Config const& cfg, BFieldMap const& bfield) :
@@ -248,12 +248,13 @@ namespace KinKal {
   template <class KTRAJ> void Track<KTRAJ>::extend(Config const& cfg, HITCOL& hits, EXINGCOL& exings) {
     // require the existing fit to be usable
     if(!fitStatus().usable())return;
-    // remember the previous config
-    auto const& oldconfig = config_.back();
-    // update the configuration
-    config_.push_back(cfg);
     // configuation check
-    if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
+    if(cfg.schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
+    // save the new configuration
+    config_.push_back(cfg);
+    // remember the previous config
+    auto oldciter = config_.rbegin(); oldciter++;
+    auto const& oldconfig = *oldciter;
     // find the range of the added information, and extend as needed
     TimeRange exrange = fittraj_->range();
     if(hits.size() >0 || exings.size() > 0){
@@ -302,7 +303,7 @@ 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 pieces correspond to the new domains. The new traj is geometrically equivalent, but not parametrically equal.
   template <class KTRAJ> void Track<KTRAJ>::replaceDomains(DOMAINCOL const& domains) {
     // if domains exist, clear them and remove all DomainWall effects
     if(domains_.size() > 0){
@@ -319,17 +320,21 @@ namespace KinKal {
       }
     }
     auto newtraj = std::make_unique<PKTRAJ>();
-    // loop over domains
+    // 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<KTRAJ>;
       using DKTRAJ = std::deque<KTRAJPTR>;
       using DKTRAJCITER = typename DKTRAJ::const_iterator;
       DKTRAJCITER first,last;
       fittraj_->pieceRange(domain->range(),first,last);
-      // loop over these pieces
+      // loop over these pieces; first and last can be the same!
       auto olditer = first;
-      while(olditer != last){
+      do {
         auto const& oldpiece = **olditer;
         // copy this piece, translating bnom to this domain's field
         KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid());
@@ -340,14 +345,14 @@ namespace KinKal {
           newpiece.range() = TimeRange(tstart,tend);
           newtraj->append(newpiece);
         }
-        ++olditer;
-      }
+        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
+    // swap out the fit trajectory; this will be used as reference for the next iterations
     fittraj_.swap(newtraj);
   }
 
@@ -694,7 +699,10 @@ namespace KinKal {
     }
     KKEFFFWDBND fwdbnds; // bounding sites used for fitting
     KKEFFREVBND revbnds;
-    setBounds(fwdbnds,revbnds);
+    if(!setBounds(fwdbnds,revbnds)){
+      status().status_ = Status::lowNDOF;
+      return;
+    }
     // initialize the fit state where we left off processing
     FitStateArray states;
     TimeRange fitrange(fwdbnds[0]->get()->time(),revbnds[0]->get()->time());

General/TimeRange.cc

@@ -1,5 +1,24 @@
 #include "KinKal/General/TimeRange.hh"
 namespace KinKal {
+
+  bool TimeRange::restrict(TimeRange const& other ) {
+    bool retval = overlaps(other);
+    if(retval){
+      range_[0] = std::max(begin(),other.begin());
+      range_[1] = std::min(end(),other.end());
+    }
+    return retval;
+  }
+
+  void TimeRange::extend(double time, TimeDir tdir) {
+    // require the resulting range to be >0
+    if(tdir == TimeDir::forwards){
+      if(time > range_[0])range_[1] = time;
+    } else {
+      if(time < range_[1])range_[0] = time;
+    }
+  }
+
   std::ostream& operator <<(std::ostream& ost, TimeRange const& trange) {
     ost << " Range [" << trange.begin() << "," << trange.end() << "]";
     return ost;

General/TimeRange.hh

@@ -5,6 +5,8 @@
 #include <ostream>
 #include <array>
 #include <stdexcept>
+#include "KinKal/General/TimeDir.hh"
+
 namespace KinKal {
   class TimeRange {
     public:
@@ -32,14 +34,9 @@ namespace KinKal {
       }
       // restrict the range to the overlap with another range. If there is no overlap
       // leave the object unchanged and return 'false'
-      bool restrict(TimeRange const& other ) {
-        bool retval = overlaps(other);
-        if(retval){
-          range_[0] = std::max(begin(),other.begin());
-          range_[1] = std::min(end(),other.end());
-        }
-        return retval;
-      }
+      bool restrict(TimeRange const& other );
+      // extend in a given direction
+      void extend(double time, TimeDir tdir);
     private:
       std::array<double,2> range_; // range of times
   };