From d05ff85198c1f5a0a38bebe1ea8f6d51b9a4079f Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Wed, 24 Sep 2025 17:27:56 -0700 Subject: [PATCH 1/7] Fix phi0 phase issue --- Trajectory/LoopHelix.cc | 10 +++++----- Trajectory/PiecewiseTrajectory.hh | 23 +++++++++++++++++++++++ 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/Trajectory/LoopHelix.cc b/Trajectory/LoopHelix.cc index 1bf2350b..8152ec43 100644 --- a/Trajectory/LoopHelix.cc +++ b/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); } diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index b9e25ccc..f0d461c5 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/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. This will throw if there's a time gap or overlap outside the given tolerance + void add(KTRAJPTR const& pieceptr, double tol, 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,27 @@ namespace KinKal { } } + template void PiecewiseTrajectory::add(KTRAJPTR const& pieceptr,double tol, TimeDir tdir) { + if(pieces_.empty()){ + pieces_.push_back(pieceptr); + } else { + if(pieceptr->range().range() < tol)throw std::invalid_argument("PiecewiseTrajectory::add; piece range too small"); + double dt = tdir == TimeDir::forwards ? fabs(this->range().end() -pieceptr->range().begin()) : fabs(this->range().begin() -pieceptr->range().end()); + if(dt > tol)throw std::invalid_argument("PiecewiseTrajectory::add; range error"); + 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 size_t PiecewiseTrajectory::nearestIndex(double time) const { size_t retval; if(pieces_.empty())throw std::length_error("Empty PiecewiseTrajectory!"); From cb2accdf30440c7af89b9d52b8c712e7fab2cef8 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 29 Sep 2025 12:14:16 -0700 Subject: [PATCH 2/7] Improve comments --- Detector/ElementXing.hh | 16 +++++++++------- Fit/Material.hh | 2 +- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/Detector/ElementXing.hh b/Detector/ElementXing.hh index 87a29b08..afcb4cb0 100644 --- a/Detector/ElementXing.hh +++ b/Detector/ElementXing.hh @@ -54,16 +54,17 @@ namespace KinKal { } template void ElementXing::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 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 Parameters ElementXing::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); } diff --git a/Fit/Material.hh b/Fit/Material.hh index 6dac9534..a4aa205b 100644 --- a/Fit/Material.hh +++ b/Fit/Material.hh @@ -99,7 +99,7 @@ namespace KinKal { } template void Material::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"); From 68708f88165ecc6cb9d42e18a2d8d79038955c0a Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Thu, 2 Oct 2025 17:14:00 -0700 Subject: [PATCH 3/7] Simplify piece adding --- Trajectory/PiecewiseTrajectory.hh | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index f0d461c5..4ee3284b 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -48,7 +48,7 @@ namespace KinKal { void prepend(KTRAJ const& newpiece, bool allowremove=false); void add(KTRAJ const& newpiece, TimeDir tdir=TimeDir::forwards, bool allowremove=false); // literally append a piece. This will throw if there's a time gap or overlap outside the given tolerance - void add(KTRAJPTR const& pieceptr, double tol, TimeDir tdir=TimeDir::forwards); + 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]; } @@ -176,13 +176,10 @@ namespace KinKal { } } - template void PiecewiseTrajectory::add(KTRAJPTR const& pieceptr,double tol, TimeDir tdir) { + template void PiecewiseTrajectory::add(KTRAJPTR const& pieceptr, TimeDir tdir) { if(pieces_.empty()){ pieces_.push_back(pieceptr); - } else { - if(pieceptr->range().range() < tol)throw std::invalid_argument("PiecewiseTrajectory::add; piece range too small"); - double dt = tdir == TimeDir::forwards ? fabs(this->range().end() -pieceptr->range().begin()) : fabs(this->range().begin() -pieceptr->range().end()); - if(dt > tol)throw std::invalid_argument("PiecewiseTrajectory::add; range error"); + } else if(!(this->range().contains(pieceptr->range()))){ if(tdir == TimeDir::forwards){ // synchronize phase variables pieceptr->syncPhi0(*pieces_.back()); From ea6b765f7576e854960d23d656441d5739ae2847 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Thu, 2 Oct 2025 17:14:58 -0700 Subject: [PATCH 4/7] Record unbiased residual failure as a flag --- Detector/Residual.hh | 13 +++++++------ Detector/ResidualHit.hh | 8 ++++++-- Examples/ScintHit.hh | 2 +- Examples/SimpleWireHit.hh | 6 +++--- 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/Detector/Residual.hh b/Detector/Residual.hh index 9514dfc0..f21fe60f 100644 --- a/Detector/Residual.hh +++ b/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); } diff --git a/Detector/ResidualHit.hh b/Detector/ResidualHit.hh index c3b6374f..68236746 100644 --- a/Detector/ResidualHit.hh +++ b/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); } diff --git a/Examples/ScintHit.hh b/Examples/ScintHit.hh index 0c5846d0..8a370e8e 100644 --- a/Examples/ScintHit.hh +++ b/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); } diff --git a/Examples/SimpleWireHit.hh b/Examples/SimpleWireHit.hh index 8f95b9da..e5657db7 100644 --- a/Examples/SimpleWireHit.hh +++ b/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 From 71bd40835eb862e4984bdf62d37aaeb2a3a88ad1 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 3 Oct 2025 08:49:23 -0700 Subject: [PATCH 5/7] Fix comment --- Trajectory/PiecewiseTrajectory.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index 4ee3284b..e7c6aa30 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -47,7 +47,7 @@ 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. This will throw if there's a time gap or overlap outside the given tolerance + // 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)]; } From 0b5a94a603873ac56835059ecb39744567fe1cb7 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 3 Oct 2025 09:53:03 -0700 Subject: [PATCH 6/7] Protect against regrowing fits without domains. Fix indent --- Fit/Track.hh | 44 +++++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index e3500983..62b41547 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -164,10 +164,10 @@ namespace KinKal { // minimal constructor for subclasses. The resulting object has no fit. template Track::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,16 +198,18 @@ 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(); @@ -215,9 +217,9 @@ namespace KinKal { // copy constructor template Track::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(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)); } }; From c72d8cd01193354e04d5fc382b20ff873c5da4e5 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 3 Oct 2025 16:04:43 -0700 Subject: [PATCH 7/7] Add protection --- Fit/Track.hh | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index 62b41547..1f91ebd7 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -686,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);