Merge pull request #209 from brownd1978/regrow

Support for regrowing fits from seeds.

Author: brownd1978

Detector/ElementXing.hh

@@ -33,8 +33,7 @@ namespace KinKal {
       virtual KTRAJ const& referenceTrajectory() const =0; // trajectory WRT which the xing is defined
       virtual std::vector<MaterialXing>const&  matXings() const =0; // Effect of each physical material component of this detector element on this trajectory
       virtual void print(std::ostream& ost=std::cout,int detail=0) const =0;
-      // crossings  without material are inactive
-      bool active() const { return matXings().size() > 0; }
+      virtual bool active() const =0;
       // momentum change and variance increase associated with crossing this element forwards in time, in spatial basis
       void momentumChange(SVEC3& dmom, SMAT& dmomvar) const;
       // parameter change associated with crossing this element forwards in time

Examples/StrawXing.hh

@@ -57,6 +57,7 @@ namespace KinKal {
       KTRAJ const& referenceTrajectory() const override { return tpca_.particleTraj(); }
       std::vector<MaterialXing>const&  matXings() const override { return mxings_; }
       void print(std::ostream& ost=std::cout,int detail=0) const override;
+      bool active() const override { return mxings_.size() > 0; }
       // accessors
       auto const& closestApproach() const { return tpca_; }
       auto const& strawMaterial() const { return smat_; }

Fit/Domain.hh

@@ -11,8 +11,8 @@
 namespace KinKal {
   class Domain : public TimeRange {
     public:
-      Domain(double lowtime, double range, VEC3 const& bnom, double tol) : range_(lowtime,lowtime+range), bnom_(bnom), tol_(tol) {}
-      Domain(TimeRange const& range, VEC3 const& bnom, double tol) : range_(range), bnom_(bnom), tol_(tol) {}
+      Domain(double lowtime, double range, VEC3 const& bnom) : range_(lowtime,lowtime+range), bnom_(bnom) {}
+      Domain(TimeRange const& range, VEC3 const& bnom) : range_(range), bnom_(bnom) {}
       // clone op for reinstantiation
       Domain(Domain const&);
       std::shared_ptr< Domain > clone(CloneContext&) const;
@@ -22,21 +22,17 @@ namespace KinKal {
       double begin() const { return range_.begin();}
       double end() const { return range_.end();}
       double mid() const { return range_.mid();}
-      double tolerance() const { return tol_; }
       auto const& bnom() const { return bnom_; }
       void updateBNom( VEC3 const& bnom) { bnom_ = bnom; }; // used in DomainWall updating
     private:
       TimeRange range_; // range of this domain
       VEC3 bnom_; // nominal BField for this domain
-      double tol_; // tolerance used to create this domain
   };
 
   // clone op for reinstantiation
   Domain::Domain(Domain const& rhs):
       range_(rhs.range()),
-      bnom_(rhs.bnom()),
-      tol_(rhs.tolerance()){
-    /**/
+      bnom_(rhs.bnom()){
   }
 
   std::shared_ptr<Domain> Domain::clone(CloneContext& context) const{

Fit/Status.hh

@@ -17,7 +17,7 @@ namespace KinKal {
     std::string comment_; // further information about the status
     bool usable() const { return status_ > unfit && status_ < lowNDOF; }
     bool needsFit() const { return status_ == unfit || status_ == unconverged; }
-    Status(unsigned miter,unsigned iter=0) : miter_(miter), iter_(iter), status_(unfit), chisq_(NParams()){}
+    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);
   };
   std::ostream& operator <<(std::ostream& os, Status const& fitstatus );

Fit/Track.hh

@@ -20,7 +20,7 @@ namespace KinKal {
       bool inRange(double t) const {return t >= begin() && t < end(); }
       bool null() const { return end() == begin(); }
       bool overlaps(TimeRange const& other ) const {
-        return inRange(other.begin()) || inRange(other.end()) || contains(other) || other.contains(*this); }
+        return inRange(other.begin()) || other.inRange(begin()) || contains(other) || other.contains(*this); }
       bool contains(TimeRange const& other) const {
         return (begin() <= other.begin() && end() >= other.end()); }
       // force time to be in range

General/TimeRange.hh

@@ -464,7 +464,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) {
     }
   }
   // create and fit the track
-  KKTRK kktrk(config,*BF,seedtraj,thits,dxings);
+  KKTRK kktrk(config,*BF,straj,thits,dxings);
   if(extend && kktrk.fitStatus().usable())kktrk.extend(exconfig,exthits, exdxings);
   if(!printbad)kktrk.print(cout,detail);
   TFile fitfile((KTRAJ::trajName() + string("FitTest") + tfname + string(".root")).c_str(),"RECREATE");

Tests/FitTest.hh

@@ -10,6 +10,7 @@
 #include "KinKal/General/MomBasis.hh"
 #include "KinKal/General/TimeRange.hh"
 #include <deque>
+#include <array>
 #include <memory>
 #include <ostream>
 #include <stdexcept>
@@ -20,6 +21,8 @@ namespace KinKal {
     public:
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       using DKTRAJ = std::deque<KTRAJPTR>;
+      using DKTRAJITER = DKTRAJ::iterator;
+      using DKTRAJCITER = DKTRAJ::const_iterator;
       // forward calls to the pieces
       VEC3 position3(double time) const { return nearestPiece(time).position3(time); }
       VEC3 velocity(double time) const { return nearestPiece(time).velocity(time); }
@@ -55,6 +58,8 @@ namespace KinKal {
       KTRAJ& back() { return *pieces_.back(); }
       KTRAJPTR const& frontPtr() const { return pieces_.front(); }
       KTRAJPTR const& backPtr() const { return pieces_.back(); }
+      void pieceRange(TimeRange const& range, DKTRAJCITER& first, DKTRAJCITER& last  ) const;
+      void pieceRange(TimeRange const& range,DKTRAJITER& first, DKTRAJITER& last );
       size_t nearestIndex(double time) const;
       DKTRAJ const& pieces() const { return pieces_; }
       // test for spatial gaps
@@ -68,13 +73,17 @@ namespace KinKal {
   template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::setRange(TimeRange const& trange, bool trim) {
     // trim pieces as necessary
     if(trim){
-      while(pieces_.size() > 1 && trange.begin() > pieces_.front()->range().end() ) pieces_.pop_front();
-      while(pieces_.size() > 1 && trange.end() < pieces_.back()->range().begin() ) pieces_.pop_back();
+      auto ipiece = pieces_.begin();
+      while (ipiece != pieces_.end() && !(trange.overlaps((*ipiece)->range())))++ipiece;
+      if(ipiece != pieces_.begin())pieces_.erase(pieces_.begin(),--ipiece);
+      auto jpiece=pieces_.rbegin();
+      while(jpiece != pieces_.rend() && !(trange.overlaps((*jpiece)->range())))++jpiece;
+      pieces_.erase(jpiece.base(),pieces_.end());
     } else if(trange.begin() > pieces_.front()->range().end() || trange.end() < pieces_.back()->range().begin())
-      throw std::invalid_argument("Invalid Range");
+      throw std::invalid_argument("PiecewiseTrajectory::setRange; Invalid Range");
     // update piece range
-    pieces_.front()->setRange(TimeRange(trange.begin(),pieces_.front()->range().end()));
-    pieces_.back()->setRange(TimeRange(pieces_.back()->range().begin(),trange.end()));
+    pieces_.front()->range().restrict(trange);
+    pieces_.back()->range().restrict(trange);
   }
 
   template <class KTRAJ> PiecewiseTrajectory<KTRAJ>::PiecewiseTrajectory(KTRAJ const& piece) : pieces_(1,std::make_shared<KTRAJ>(piece))
@@ -89,13 +98,13 @@ namespace KinKal {
         prepend(newpiece,allowremove);
         break;
       default:
-        throw std::invalid_argument("Invalid direction");
+        throw std::invalid_argument("PiecewiseTrajectory::add; Invalid direction");
     }
   }
 
   template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::prepend(KTRAJ const& newpiece, bool allowremove) {
     // new piece can't have null range
-    if(newpiece.range().null())throw std::invalid_argument("Can't prepend null range traj");
+    if(newpiece.range().null())throw std::invalid_argument("PiecewiseTrajectory::prepend; Can't prepend null range traj");
     if(pieces_.empty()){
       pieces_.push_back(std::make_shared<KTRAJ>(newpiece));
     } else {
@@ -104,7 +113,7 @@ namespace KinKal {
         if(allowremove)
           *this = PiecewiseTrajectory(newpiece);
         else
-          throw std::invalid_argument("range overlap");
+          throw std::invalid_argument("PiecewiseTrajector::prepend; range overlap");
       } else {
         // find the piece that needs to be modified
         size_t ipiece = nearestIndex(newpiece.range().end());
@@ -122,15 +131,15 @@ namespace KinKal {
           pieces_.push_front(std::make_shared<KTRAJ>(newpiece));
           pieces_.front()->range() = TimeRange(tmin,pieces_.front()->range().end());
         } else {
-          throw std::invalid_argument("range error");
+          throw std::invalid_argument("PiecewiseTrajectory::prepend; range error");
         }
       }
     }
   }
 
   template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::append(KTRAJ const& newpiece, bool allowremove) {
     // new piece can't have null range
-    if(newpiece.range().null())throw std::invalid_argument("Can't append null range traj");
+    if(newpiece.range().null())throw std::invalid_argument("PiecewiseTrajectory::append; Can't append null range traj");
     if(pieces_.empty()){
       pieces_.push_back(std::make_shared<KTRAJ>(newpiece));
     } else {
@@ -139,7 +148,7 @@ namespace KinKal {
         if(allowremove)
           *this = PiecewiseTrajectory(newpiece);
         else
-          throw std::invalid_argument("range overlap");
+          throw std::invalid_argument("PiecewiseTrajectory::append; range overlap");
       } else {
         // find the piece that needs to be modified
         size_t ipiece = nearestIndex(newpiece.range().begin());
@@ -159,7 +168,7 @@ namespace KinKal {
           pieces_.push_back(std::make_shared<KTRAJ>(newpiece));
           pieces_.back()->range() = TimeRange(pieces_.back()->range().begin(),tmax);
         } else {
-          throw std::invalid_argument("range error");
+          throw std::invalid_argument("PiecewiseTrajectory::append; range error");
         }
       }
     }
@@ -253,6 +262,35 @@ namespace KinKal {
     return ost;
   }
 
+  template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::pieceRange(TimeRange const& range,
+      std::deque<std::shared_ptr<KTRAJ>>::const_iterator& first,
+      std::deque<std::shared_ptr<KTRAJ>>::const_iterator& last ) const {
+    first = last = pieces_.end();
+    // check for no overlap
+    if(this->range().overlaps(range)){
+      // find the first and last pieces which overlap with the range. They can be the same piece.
+      first = pieces_.cbegin();
+      while(first != pieces_.cend() && !((*first)->range().overlaps(range))) ++first;
+      auto rlast = pieces_.crbegin();
+      while(rlast != pieces_.crend() && !((*rlast)->range().overlaps(range))) ++rlast;
+      last = (rlast+1).base(); // convert back to forwards-iterator
+    }
+  }
+
+  template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::pieceRange(TimeRange const& range,
+      std::deque<std::shared_ptr<KTRAJ>>::iterator& first,
+      std::deque<std::shared_ptr<KTRAJ>>::iterator& last) {
+    first = last = pieces_.end();
+    // check for no overlap
+    if(this->range().overlaps(range)){
+      first = pieces_.begin();
+      while(first != pieces_.end() && !((*first)->range().overlaps(range))) ++first;
+      auto rlast = pieces_.rbegin();
+      while(rlast != pieces_.rend() && !((*rlast)->range().overlaps(range))) ++rlast;
+      last= (rlast+1).base();
+    }
+  }
+
   // clone op for reinstantiation
   template<class KTRAJ>
   PiecewiseTrajectory<KTRAJ>::PiecewiseTrajectory(PiecewiseTrajectory<KTRAJ> const& rhs){