Regrow fixes

Detector/ElementXing.hh

@@ -54,29 +54,30 @@ namespace KinKal {
   }
 
   template <class KTRAJ> void ElementXing<KTRAJ>::momentumChange(SVEC3& dmom, SMAT& dmomvar) const {
-    // compute the parameter effect for forwards time
+    // compute the momentum change moving forwards time in global cartesian basis
+    // first, compute the change and variance in the momentum basis (along and perp to the momentum direction)
     double dm, paramomvar, perpmomvar;
     materialEffects(dm, paramomvar,perpmomvar);
-    // momentum change in forward time direction /mdue to energy loss; this is along the momentum
+    // momentum change in forward time direction is due to energy loss; this is along the momentum
     auto momdir = referenceTrajectory().direction(time());
     dmom = dm*SVEC3(momdir.X(),momdir.Y(),momdir.Z());
-    // now update the covariance; this includes smearing from energy straggling and multiple scattering
-    // move variances into a matrix
+    // now compute the covariance; this includes smearing from energy straggling and multiple scattering.
+    // This is a diagonal matrix in momentum basis; no physical correlation between stragling and scattering, or orthogonal scattering.
     ROOT::Math::SVector<double, 6>  varvec(paramomvar, 0, perpmomvar, 0, 0, perpmomvar);
-    SMAT mmvar(varvec);
+    SMAT mmvar(varvec); // construct matrix from upper-diagonal elements
     // loop over the momentum change basis directions and create the transform matrix between that and global Cartesian basis
     SSMAT dmdxyz; // momentum basis -> Cartesian conversion matrix
     for(int idir=0;idir<MomBasis::ndir; idir++) {
       auto mdir = referenceTrajectory().direction(time(),static_cast<MomBasis::Direction>(idir));
       SVEC3 vmdir(mdir.X(), mdir.Y(), mdir.Z());
       dmdxyz.Place_in_col(vmdir,0,idir);
     }
-    // return variance in global Cartesian coordinates
+    // return covariance in global Cartesian coordinates
     dmomvar = ROOT::Math::Similarity(dmdxyz,mmvar);
   }
 
   template <class KTRAJ> Parameters ElementXing<KTRAJ>::parameterChange(double varscale) const {
-    // compute this xing's effect on momentum in global Cartesian
+    // compute this xing's effect on parameters and covariance. First compute the momentum change and variance
     SVEC3 dmom;
     SMAT dmomvar;
     momentumChange(dmom,dmomvar);
@@ -85,6 +86,7 @@ namespace KinKal {
     auto dmomp = dPdM*dmom;
     // scale covariance as needed
     dmomvar*= varscale;
+    // jacobean transform of covariance into parameter space
     auto dmompvar = ROOT::Math::Similarity(dPdM,dmomvar);
     return Parameters(dmomp,dmompvar);
   }

Detector/Residual.hh

@@ -18,21 +18,22 @@ namespace KinKal {
       double parameterVariance() const  { return pvar_; }
       double variance() const { return mvar_ + pvar_; }
       DVEC const& dRdP() const { return dRdP_; } // derivative of this residual WRT parameters
-      double chisq() const { return active_ ? (value_*value_)/variance() : 0.0; }
-      double chi() const { return active_ ? value_/sqrt(variance()): 0.0; }
-      double pull() const { return chi(); }
-      unsigned nDOF() const { return active_ ? 1 : 0; }
       bool active() const { return active_; }
+      double chisq() const { return active() ? (value_*value_)/variance() : 0.0; }
+      double chi() const { return active() ? value_/sqrt(variance()): 0.0; }
+      double pull() const { return chi(); }
+      unsigned nDOF() const { return active() ? 1 : 0; }
       // calculate the weight WRT some parameters implied by this residual.  Optionally scale the variance
       Weights weight(DVEC const& params, double varscale=1.0) const;
-      Residual(double value, double mvar, double pvar, bool active, DVEC const& dRdP) : value_(value), mvar_(mvar), pvar_(pvar), active_(active), dRdP_(dRdP){}
+      Residual(double value, double mvar, double pvar, DVEC const& dRdP,bool active=true) :
+        value_(value), mvar_(mvar), pvar_(pvar), dRdP_(dRdP), active_(active){}
       Residual() : value_(0.0), mvar_(-1.0), pvar_(-1.0), active_(false) {}
     private:
       double value_;  // value for this residual
       double mvar_; // estimated variance due to measurement uncertainty
       double pvar_; // estimated variance due to parameter uncertainty
-      bool active_; // whether this residual is active or not
       DVEC dRdP_; // derivative of this residual WRT the trajectory parameters, evaluated at the reference parameters
+      bool active_;
   };
   std::ostream& operator <<(std::ostream& ost, Residual const& res);
 }

Detector/ResidualHit.hh

@@ -47,8 +47,12 @@ namespace KinKal {
     // project the parameter differnce to residual space and 'correct' the reference residual to be WRT these parameters
     double uresid = resid.value() - ROOT::Math::Dot(dpvec,resid.dRdP());
     double pvar = ROOT::Math::Similarity(resid.dRdP(),params.covariance());
-    if(pvar<0) throw std::runtime_error("Covariance projection inconsistency");
-    return Residual(uresid,resid.variance(),pvar,resid.active(),resid.dRdP());
+    bool active = resid.active();
+    if(pvar<0.0){
+      active = false;
+      pvar = resid.variance();
+    }
+    return Residual(uresid,resid.variance(),pvar,resid.dRdP(),active);
 
   }
 

Examples/ScintHit.hh

@@ -112,7 +112,7 @@ namespace KinKal {
     // Might want to do more updating (set activity) based on DOCA in future: TODO
     double dd2 = tpca_.dirDot()*tpca_.dirDot();
     double totvar = tvar_ + wvar_*dd2/(saxis_.speed(sdist)*saxis_.speed(sdist)*(1.0-dd2));
-    rresid_ = Residual(tpca_.deltaT(),totvar,0.0,true,tpca_.dTdP());
+    rresid_ = Residual(tpca_.deltaT(),totvar,0.0,tpca_.dTdP());
     this->updateWeight(config);
   }
 

Examples/SimpleWireHit.hh

@@ -148,16 +148,16 @@ namespace KinKal {
     }
     rresid_[tresid] = rresid_[dresid] = Residual();
     if(whstate_.active()){
-      rresid_[tresid] = Residual(ca_.deltaT() - tot_, totvar_,0.0,true,ca_.dTdP()); // always constrain to TOT; this stabilizes the fit
+      rresid_[tresid] = Residual(ca_.deltaT() - tot_, totvar_,0.0,ca_.dTdP()); // always constrain to TOT; this stabilizes the fit
       if(whstate_.useDrift()){
         // translate PCA to residual. Use ambiguity assignment to convert drift time to a drift radius
         double dr = dvel_*whstate_.lrSign()*ca_.deltaT() -ca_.doca();
         DVEC dRdP = dvel_*whstate_.lrSign()*ca_.dTdP() -ca_.dDdP();
-        rresid_[dresid] = Residual(dr,tvar_*dvel_*dvel_,0.0,true,dRdP);
+        rresid_[dresid] = Residual(dr,tvar_*dvel_*dvel_,0.0,dRdP);
       } else {
         // interpret DOCA against the wire directly as a residuals
         double nulldvar = dvel_*dvel_*(ca_.deltaT()*ca_.deltaT()+0.8);
-        rresid_[dresid] = Residual(ca_.doca(),nulldvar,0.0,true,ca_.dDdP());
+        rresid_[dresid] = Residual(ca_.doca(),nulldvar,0.0,ca_.dDdP());
       }
     }
     // now update the weight

Fit/Material.hh

@@ -99,7 +99,7 @@ namespace KinKal {
   }
 
   template<class KTRAJ> void Material<KTRAJ>::extrapolate(PTRAJ& ptraj,TimeDir tdir) {
-  // make sure the traj can be extrapolated
+  // make sure the traj can be appended
     if((tdir == TimeDir::forwards && ptraj.back().range().begin() > time()) ||
         (tdir == TimeDir::backwards && ptraj.front().range().end() < time()) )
       throw std::invalid_argument("Material: Can't append piece");

Fit/Track.hh

@@ -164,10 +164,10 @@ namespace KinKal {
   // minimal constructor for subclasses. The resulting object has no fit.
   template <class KTRAJ> Track<KTRAJ>::Track(Config const& cfg, BFieldMap const& bfield) :
     config_{cfg}, bfield_(bfield)
-  {
-    if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
-    history_.push_back(Status(0,0,Status::unfit, "Construction"));
-  }
+    {
+      if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule");
+      history_.push_back(Status(0,0,Status::unfit, "Construction"));
+    }
 
   // construct from configuration, reference (seed) fit, hits,and materials specific to this fit. This will compute the domains according to the configuration before fitting.
   //
@@ -198,26 +198,28 @@ namespace KinKal {
     fittraj_ = std::move(fittraj); // steal the underlying object
     // truncate the domains and fit trajectory to be within the detector range
     auto detrange = detectorRange(hits,exings,true);
-    auto idom = domains.begin();
-    // stop at the 1st domain overlaping the detector range, and erase all elements up to that point
-    while(idom != domains.end() && !(detrange.overlaps((*idom)->range())))++idom;
-    if(idom != domains.begin())domains.erase(domains.begin(),--idom);// leave the overlapping piece
-    auto jdom= domains.rbegin();
-    while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom;
-    domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator
-   // trim the trajectory to this range
-    detrange.combine((*domains.begin())->range());
-    detrange.combine((*domains.rbegin())->range());
+    if(domains.size() > 0){
+      auto idom = domains.begin();
+      // stop at the 1st domain overlaping the detector range, and erase all elements up to that point
+      while(idom != domains.end() && !(detrange.overlaps((*idom)->range())))++idom;
+      if(idom != domains.begin())domains.erase(domains.begin(),--idom);// leave the overlapping piece
+      auto jdom= domains.rbegin();
+      while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom;
+      domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator
+      // trim the trajectory to this range
+      detrange.combine((*domains.begin())->range());
+      detrange.combine((*domains.rbegin())->range());
+    }
     fittraj_->setRange(detrange,true);
     createEffects(hits,exings,domains);
     fit();
   }
 
   // copy constructor
   template<class KTRAJ> Track<KTRAJ>::Track(const Track& rhs, CloneContext& context) :
-      config_(rhs.configs()),
-      bfield_(rhs.bfield()),
-      history_(rhs.history())
+    config_(rhs.configs()),
+    bfield_(rhs.bfield()),
+    history_(rhs.history())
   {
     fittraj_ = std::make_unique<PKTRAJ>(rhs.fitTraj());
     hits_.reserve(rhs.hits().size());
@@ -231,12 +233,12 @@ namespace KinKal {
       exings_.push_back(xng);
     }
     for (const auto& ptr: rhs.domains()){
-        auto dmn = context.get(ptr);
-        domains_.insert(dmn);
+      auto dmn = context.get(ptr);
+      domains_.insert(dmn);
     }
     for (const auto& ptr: rhs.effects()){
-        auto eff = ptr->clone(context);
-        effects_.push_back(std::move(eff));
+      auto eff = ptr->clone(context);
+      effects_.push_back(std::move(eff));
     }
   };
 
@@ -684,8 +686,10 @@ namespace KinKal {
     // sort effects in case ends have migrated
     effects_.sort(KKEFFComp());
     // update domains as needed to cover the end effects
-    TimeRange endrange(effects_.front()->time(),effects_.back()->time());
-    extendDomains(endrange);
+    if(config().bfcorr_){
+      TimeRange endrange(effects_.front()->time(),effects_.back()->time());
+      extendDomains(endrange);
+    }
     KKEFFFWDBND fwdbnds; // bounding sites used for fitting
     KKEFFREVBND revbnds;
     setBounds(fwdbnds,revbnds);

Trajectory/LoopHelix.cc

@@ -72,12 +72,12 @@ namespace KinKal {
 
   void LoopHelix::setBNom(double time, VEC3 const& newbnom) {
     auto db = newbnom - bnom_;
-//    PSMAT dpdpdb =ROOT::Math::SMatrixIdentity();
-//    PSMAT dpdpdb = dPardPardB(time,db);
-//    std::cout << "dpdpdb = " << dpdpdb << std::endl;
-//    pars_.covariance() = ROOT::Math::Similarity(dpdpdb,pars_.covariance());
+    // rotate the covariance to this new basis: this is still work in progress TODO
+    //    PSMAT dpdpdb =ROOT::Math::SMatrixIdentity();
+    //    PSMAT dpdpdb = dPardPardB(time,db);
+    //    std::cout << "dpdpdb = " << dpdpdb << std::endl;
+    //    pars_.covariance() = ROOT::Math::Similarity(dpdpdb,pars_.covariance());
     pars_.parameters() += dPardB(time,db);
-    // rotate covariance: for now, just for the magnitude change.  Rotation is still TODO
     resetBNom(newbnom);
   }
 

Trajectory/PiecewiseTrajectory.hh

@@ -47,6 +47,8 @@ namespace KinKal {
       void append(KTRAJ const& newpiece, bool allowremove=false);
       void prepend(KTRAJ const& newpiece, bool allowremove=false);
       void add(KTRAJ const& newpiece, TimeDir tdir=TimeDir::forwards, bool allowremove=false);
+      // literally append a piece.
+      void add(KTRAJPTR const& pieceptr, TimeDir tdir=TimeDir::forwards);
       // Find the piece associated with a particular time
       KTRAJPTR const& nearestTraj(double time) const { return pieces_[nearestIndex(time)]; }
       KTRAJPTR const& indexTraj(size_t index) const { return pieces_[index]; }
@@ -174,6 +176,24 @@ namespace KinKal {
     }
   }
 
+  template <class KTRAJ> void PiecewiseTrajectory<KTRAJ>::add(KTRAJPTR const& pieceptr, TimeDir tdir) {
+    if(pieces_.empty()){
+      pieces_.push_back(pieceptr);
+    } else if(!(this->range().contains(pieceptr->range()))){
+      if(tdir == TimeDir::forwards){
+        // synchronize phase variables
+        pieceptr->syncPhi0(*pieces_.back());
+        // fine-adjust the time to have 0 gap
+        pieceptr->setRange(TimeRange(pieces_.back()->range().end(),pieceptr->range().end()));
+        pieces_.push_back(pieceptr);
+      } else {
+        pieceptr->syncPhi0(*pieces_.front());
+        pieceptr->setRange(TimeRange(pieceptr->range().begin(),pieces_.front()->range().begin()));
+        pieces_.push_front(pieceptr);
+      }
+    }
+  }
+
   template <class KTRAJ> size_t PiecewiseTrajectory<KTRAJ>::nearestIndex(double time) const {
     size_t retval;
     if(pieces_.empty())throw std::length_error("Empty PiecewiseTrajectory!");