From ae4e3ebb0bea3a82743c832f4cc0eb708bb8ed32 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 9 Jan 2026 15:15:58 -0800 Subject: [PATCH 1/4] Convert PiecewiseClosestApproach to use KTRAJ instead of PKTRAJ to avoid copy problems --- Examples/ScintHit.hh | 4 +- Examples/SimpleWireHit.hh | 4 +- Examples/StrawXing.hh | 4 +- Trajectory/ClosestApproach.hh | 6 +-- Trajectory/PiecewiseClosestApproach.hh | 42 +++++++-------------- Trajectory/PointClosestApproach.hh | 31 ++++++++------- Trajectory/PointPiecewiseClosestApproach.hh | 36 ++++++------------ 7 files changed, 52 insertions(+), 75 deletions(-) diff --git a/Examples/ScintHit.hh b/Examples/ScintHit.hh index 8a370e8e..84ec4595 100644 --- a/Examples/ScintHit.hh +++ b/Examples/ScintHit.hh @@ -71,7 +71,7 @@ namespace KinKal { template ScintHit::ScintHit(PCA const& pca, double tvar, double wvar) : saxis_(pca.sensorTraj()), tvar_(tvar), wvar_(wvar), - tpca_(pca.localTraj(),saxis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()) + tpca_(static_cast(pca)) {} template Residual const& ScintHit::refResidual(unsigned ires) const { @@ -83,7 +83,7 @@ namespace KinKal { // use previous hint, or initialize from the sensor time CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(saxis_.measurementTime(), saxis_.measurementTime()); PCA pca(ptraj,saxis_,tphint,precision()); - tpca_ = pca.localClosestApproach(); + tpca_ = static_cast(pca); if(!tpca_.usable())throw std::runtime_error("ScintHit TPOCA failure"); } diff --git a/Examples/SimpleWireHit.hh b/Examples/SimpleWireHit.hh index e5657db7..fda05ab6 100644 --- a/Examples/SimpleWireHit.hh +++ b/Examples/SimpleWireHit.hh @@ -114,7 +114,7 @@ namespace KinKal { double mindoca, double driftspeed, double tvar, double tot, double totvar, double rcell, int id) : bfield_(bfield), whstate_(whstate), wire_(pca.sensorTraj()), - ca_(pca.localTraj(),wire_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()), // must be explicit to get the right sensor traj reference + ca_(static_cast(pca)), mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) { } @@ -123,7 +123,7 @@ namespace KinKal { // otherwise use the time at the center of the wire CAHint tphint = ca_.usable() ? ca_.hint() : CAHint(wire_.timeAtMidpoint(),wire_.timeAtMidpoint()); PCA pca(ptraj,wire_,tphint,precision()); - ca_ = pca.localClosestApproach(); + ca_ = static_cast(pca); if(!ca_.usable())throw std::runtime_error("WireHit TPOCA failure"); } diff --git a/Examples/StrawXing.hh b/Examples/StrawXing.hh index 5c5d4b16..3344ab0c 100644 --- a/Examples/StrawXing.hh +++ b/Examples/StrawXing.hh @@ -81,7 +81,7 @@ namespace KinKal { template StrawXing::StrawXing(PCA const& pca, StrawMaterial const& smat) : axis_(pca.sensorTraj()), smat_(smat), - tpca_(pca.localTraj(),axis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()), + tpca_(static_cast(pca)), toff_(smat.wireRadius()/pca.particleTraj().speed(pca.particleToca())), // locate the effect to 1 side of the wire to avoid overlap with hits varscale_(1.0) {} @@ -89,7 +89,7 @@ namespace KinKal { template void StrawXing::updateReference(PTRAJ const& ptraj) { CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(axis_.timeAtMidpoint(),axis_.timeAtMidpoint()); PCA pca(ptraj,axis_,tphint,precision()); - tpca_ = pca.localClosestApproach(); + tpca_ = static_cast(pca); if(!tpca_.usable())throw std::runtime_error("StrawXing TPOCA failure"); } diff --git a/Trajectory/ClosestApproach.hh b/Trajectory/ClosestApproach.hh index 0c95ef22..4401b776 100644 --- a/Trajectory/ClosestApproach.hh +++ b/Trajectory/ClosestApproach.hh @@ -74,14 +74,14 @@ namespace KinKal { void setTrajectory(KTRAJPTR ktrajptr){ ktrajptr_ = ktrajptr; } private: double precision_; // precision used to define convergence + DVEC dDdP_; // derivative of DOCA WRT Parameters + DVEC dTdP_; // derivative of TOCA WRT Parameters + protected: KTRAJPTR ktrajptr_; // kinematic particle trajectory STRAJ const& straj_; // sensor trajectory This should be a shared ptr TODO - protected: // calculate CA given the hint, and fill the state void findTCA(CAHint const& hint); ClosestApproachData tpdata_; // data payload of CA calculation - DVEC dDdP_; // derivative of DOCA WRT Parameters - DVEC dTdP_; // derivative of TOCA WRT Parameters void setParticleTOCA( double ptoca); void setSensorTOCA( double stoca); }; diff --git a/Trajectory/PiecewiseClosestApproach.hh b/Trajectory/PiecewiseClosestApproach.hh index fc0b04be..b98fe3fb 100644 --- a/Trajectory/PiecewiseClosestApproach.hh +++ b/Trajectory/PiecewiseClosestApproach.hh @@ -8,46 +8,32 @@ #include namespace KinKal { - template class PiecewiseClosestApproach : public ClosestApproach,STRAJ> { + template class PiecewiseClosestApproach : public ClosestApproach { public: using PTRAJ = ParticleTrajectory; using KTCA = ClosestApproach; using KTRAJPTR = std::shared_ptr; PiecewiseClosestApproach(PTRAJ const& ptraj, STRAJ const& straj, CAHint const& hint, double precision); - // provide access to the local (non-piecewise) information implicit in this class - size_t particleTrajIndex() const { return pindex_; } - KTRAJ const& localParticleTraj() const { return this->particleTraj().piece(pindex_); } - KTRAJPTR const& localTraj() const { return this->particleTraj().indexTraj(pindex_); } - KTCA localClosestApproach() const { return KTCA(localTraj(),this->sensorTraj(),this->precision(),this->tpData(),this->dDdP(),this->dTdP()); } - private: - size_t pindex_; // indices to the local traj used in TCA calculation }; - template PiecewiseClosestApproach::PiecewiseClosestApproach(ParticleTrajectory const& ptraj, STRAJ const& straj, CAHint const& hint, double prec) : ClosestApproach,STRAJ>(ptraj,straj,prec) { - // iteratively find the nearest piece, and CA for that piece. Start at hints if availalble, otherwise the middle + template PiecewiseClosestApproach::PiecewiseClosestApproach(ParticleTrajectory const& ptraj, STRAJ const& straj, CAHint const& hint, double prec) : KTCA(ptraj.nearestTraj(hint.particleToca_),straj,prec) { + // iteratively find the nearest piece, and CA for that piece. Start at hint static const unsigned maxiter=10; // don't allow infinite iteration. This should be a parameter TODO unsigned niter=0; - size_t oldindex= this->particleTraj().pieces().size(); - this->pindex_ = this->particleTraj().nearestIndex(hint.particleToca_); + size_t oldindex= ptraj.pieces().size(); + auto pindex = ptraj.nearestIndex(hint.particleToca_); // copy over the hint: it needs to evolve CAHint phint = hint; - // iterate until TCA is on the same piece + // iterate until the local TCA is on the same piece do{ - KTCA tpoca(this->particleTraj().piece(pindex_),this->sensorTraj(),phint,prec); - // copy the state - this->tpdata_ = tpoca.tpData(); - this->dDdP_ = tpoca.dDdP(); - this->dTdP_ = tpoca.dTdP(); - // inrange = tpoca.inRange(); - // update the hint - phint.particleToca_ = tpoca.particleToca(); - phint.sensorToca_ = tpoca.sensorToca(); - // update the piece (if needed) - oldindex = pindex_; - pindex_ = this->particleTraj().nearestIndex(tpoca.particlePoca().T()); - } while( pindex_ != oldindex && this->usable() && niter++ < maxiter); - // overwrite the status if we oscillated on the piece - if(this->tpdata_.status() == ClosestApproachData::converged && niter >= maxiter) + this->findTCA(phint); + phint = this->hint(); + oldindex = pindex; + pindex = ptraj.nearestIndex(hint.particleToca_); + this->ktrajptr_ = ptraj.indexTraj(pindex); + } while( pindex != oldindex && this->usable() && niter++ < maxiter); + // overwrite the status if we didn't converge on the piece + if(this->status() == ClosestApproachData::converged && pindex != oldindex) this->tpdata_.status_ = ClosestApproachData::unconverged; // should test explicitly for piece oscillation FIXME! // test if the solution is on a cusp and if so, chose the one with the smallest DOCA TODO diff --git a/Trajectory/PointClosestApproach.hh b/Trajectory/PointClosestApproach.hh index 8e7fb21f..0ce3c3ed 100644 --- a/Trajectory/PointClosestApproach.hh +++ b/Trajectory/PointClosestApproach.hh @@ -20,17 +20,19 @@ namespace KinKal { }; template class PointClosestApproach { public: + using KTRAJPTR = std::shared_ptr; // construct from the particle and sensor trajectories; TCA is computed on construction, given a hint as to where // to start looking, which disambiguates functions with multiple solutions PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, PCAHint const& hint, double precision); // construct without a hint - PointClosestApproach(KTRAJ const& ptraj, VEC4 const& point, double precision); + PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, double precision); // construct with full data - PointClosestApproach(KTRAJ const& ptraj, VEC4 const& point, double precision, + PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, double precision, ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP); // accessors ClosestApproachData const& tpData() const { return tpdata_; } - KTRAJ const& particleTraj() const { return ktraj_; } + KTRAJ const& particleTraj() const { return *ktrajptr_; } + KTRAJPTR const& particleTrajPtr() const { return ktrajptr_; } // derviatives of TOCA and DOCA WRT particle trajectory parameters DVEC const& dDdP() const { return dDdP_; } DVEC const& dTdP() const { return dTdP_; } @@ -58,10 +60,11 @@ namespace KinKal { // calculate CA given the hint, and fill the state private: double precision_; // precision used to define convergence - KTRAJ const& ktraj_; // kinematic particle trajectory ClosestApproachData tpdata_; // data payload of CA calculation DVEC dDdP_; // derivative of DOCA WRT Parameters DVEC dTdP_; // derivative of TOCA WRT Parameters + protected: + KTRAJPTR ktrajptr_; // kinematic particle trajectory }; template PointClosestApproach::PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, double prec) : @@ -69,10 +72,10 @@ namespace KinKal { template PointClosestApproach::PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, double prec, ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP) : - precision_(prec), ktraj_(ktraj), tpdata_(tpdata), dDdP_(dDdP), dTdP_(dTdP) {} + precision_(prec), ktrajptr_(new KTRAJ(ktraj)), tpdata_(tpdata), dDdP_(dDdP), dTdP_(dTdP) {} template PointClosestApproach::PointClosestApproach(KTRAJ const& ktraj, VEC4 const& point, PCAHint const& hint, - double prec) : precision_(prec), ktraj_(ktraj) { + double prec) : precision_(prec), ktrajptr_(new KTRAJ(ktraj)) { // sensor CA is fixed to the point tpdata_.sensCA_ = point; findTCA(hint); @@ -86,13 +89,13 @@ namespace KinKal { static const unsigned maxiter=100; // don't allow infinite iteration. This should be a parameter FIXME! unsigned niter(0); // speed doesn't change - double pspeed = ktraj_.speed(particleToca()); + double pspeed = ktrajptr_->speed(particleToca()); // iterate until change in TOCA is less than precision double dptoca(std::numeric_limits::max()); while(tpdata_.usable() && fabs(dptoca) > precision() && niter++ < maxiter) { // find positions and directions at the current TOCA estimate - tpdata_.partCA_ = ktraj_.position4(tpdata_.particleToca()); - tpdata_.pdir_ = ktraj_.direction(particleToca()); + tpdata_.partCA_ = ktrajptr_->position4(tpdata_.particleToca()); + tpdata_.pdir_ = ktrajptr_->direction(particleToca()); auto dpos = point().Vect()-particlePoca().Vect(); // compute the change in times dptoca = dpos.Dot(tpdata_.pdir_)/pspeed; @@ -104,8 +107,8 @@ namespace KinKal { else tpdata_.status_ = ClosestApproachData::unconverged; // final update - tpdata_.partCA_ = ktraj_.position4(tpdata_.particleToca()); - tpdata_.pdir_ = ktraj_.direction(particleToca()); + tpdata_.partCA_ = ktrajptr_->position4(tpdata_.particleToca()); + tpdata_.pdir_ = ktrajptr_->direction(particleToca()); tpdata_.sdir_ = tpdata_.delta().Vect(); // fill the rest of the state if(usable()){ @@ -116,13 +119,13 @@ namespace KinKal { VEC3 dvechat = dvec.Unit(); // now variances due to the particle trajectory parameter covariance // for DOCA, project the spatial position derivative along the delta-CA direction - DVDP dxdp = ktraj_.dXdPar(particleToca()); + DVDP dxdp = ktrajptr_->dXdPar(particleToca()); SVEC3 dv(dvechat.X(),dvechat.Y(),dvechat.Z()); dDdP_ = -dv*dxdp; dTdP_[KTRAJ::t0Index()] = -1.0; // TOCA is 100% anti-correlated with the (mandatory) t0 component. // project the parameter covariance onto DOCA and TOCA - tpdata_.docavar_ = ROOT::Math::Similarity(dDdP(),ktraj_.params().covariance()); - tpdata_.tocavar_ = ROOT::Math::Similarity(dTdP(),ktraj_.params().covariance()); + tpdata_.docavar_ = ROOT::Math::Similarity(dDdP(),ktrajptr_->params().covariance()); + tpdata_.tocavar_ = ROOT::Math::Similarity(dTdP(),ktrajptr_->params().covariance()); } } diff --git a/Trajectory/PointPiecewiseClosestApproach.hh b/Trajectory/PointPiecewiseClosestApproach.hh index 24c88da2..71a36190 100644 --- a/Trajectory/PointPiecewiseClosestApproach.hh +++ b/Trajectory/PointPiecewiseClosestApproach.hh @@ -8,19 +8,13 @@ #include namespace KinKal { - template class PointPiecewiseClosestApproach : public PointClosestApproach> { + template class PointPiecewiseClosestApproach : public PointClosestApproach { public: using PTRAJ = ParticleTrajectory; using KTCA = PointClosestApproach; // the constructor is the only non-inherited function PointPiecewiseClosestApproach(PTRAJ const& ptraj, VEC4 const& point, PCAHint const& hint, double precision); PointPiecewiseClosestApproach(PTRAJ const& ptraj, VEC4 const& point, double precision); - // provide access to the local (non-piecewise) information implicit in this class - size_t particleTrajIndex() const { return pindex_; } - KTRAJ const& localParticleTraj() const { return this->particleTraj().piece(pindex_); } - KTCA localClosestApproach() const { return KTCA(localParticleTraj(),this->point(),this->precision(),this->tpData(),this->dDdP(),this->dTdP()); } - private: - size_t pindex_; // indices to the local traj used in TCA calculation }; template PointPiecewiseClosestApproach::PointPiecewiseClosestApproach(PTRAJ const& ptraj, VEC4 const& point, double prec) : @@ -28,30 +22,24 @@ namespace KinKal { // iteratively find the nearest piece, and CA for that piece. Start at hints if availalble, otherwise the middle template PointPiecewiseClosestApproach::PointPiecewiseClosestApproach(PTRAJ const& ptraj, VEC4 const& point, PCAHint const& hint, double prec) : - KTCA(ptraj,point,prec) { - // iteratively find the nearest piece, and CA for that piece. Start at hints if availalble, otherwise the middle - static const unsigned maxiter=10; // don't allow infinite iteration. This should be a parameter FIXME! + KTCA(ptraj.nearestTraj(hint.particleToca_),point,prec) { + // iteratively find the nearest piece, and CA for that piece + static const unsigned maxiter=10; // don't allow infinite iteration. unsigned niter=0; - size_t oldindex= this->particleTraj().pieces().size(); - pindex_ = this->particleTraj().nearestIndex(hint.particleToca_); + size_t oldindex= ptraj.pieces().size(); + auto pindex = ptraj.nearestIndex(hint.particleToca_); // copy over the hint: it needs to evolve PCAHint phint = hint; // iterate until TCA is on the same piece do{ - KTCA tpoca(this->particleTraj().piece(pindex_),this->point(),phint,this->precision()); - // copy the state - this->tpdata_ = tpoca.tpData(); - this->dDdP_ = tpoca.dDdP(); - this->dTdP_ = tpoca.dTdP(); - // inrange = tpoca.inRange(); - // update the hint + this->findTCA(phint); phint.particleToca_ = tpoca.particleToca(); - // update the piece (if needed) - oldindex = pindex_; - pindex_ = this->particleTraj().nearestIndex(tpoca.particlePoca().T()); - } while( pindex_ != oldindex && usable() && niter++ < maxiter); + oldindex = pindex; + pindex = ptraj.nearestIndex(tpoca.particlePoca().T()); + this->ktrajptr_ = ptraj.indexTraj(pindex); + } while( pindex != oldindex && usable() && niter++ < maxiter); // overwrite the status if we oscillated on the piece - if(this->tpdata_.status() == ClosestApproachData::converged && niter >= maxiter) + if(this->tpdata_.status() == ClosestApproachData::converged && pindex != oldindex) this->tpdata_.status_ = ClosestApproachData::unconverged; // should test explicitly for piece oscillation FIXME! // test if the solution is on a cusp and if so, chose the one with the smallest DOCA TODO From 51c26a93b1281fa4befbf0d7a5050c6e2a1f4fc6 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sat, 10 Jan 2026 14:10:11 -0800 Subject: [PATCH 2/4] Make SensorTraj a ptr in ClosestApproach --- Detector/Hit.hh | 2 +- Detector/ParameterHit.hh | 2 +- Examples/ScintHit.hh | 40 +++++++-------- Examples/SimpleWireHit.hh | 67 ++++++++++++-------------- Examples/StrawXing.hh | 18 +++---- Tests/ClosestApproachTest.hh | 22 ++++----- Tests/FitTest.hh | 14 +++--- Tests/HitTest.hh | 10 ++-- Tests/ParticleTrajectoryTest.hh | 7 +-- Tests/ToyMC.hh | 9 ++-- Trajectory/ClosestApproach.hh | 36 +++++++------- Trajectory/PiecewiseClosestApproach.hh | 5 +- 12 files changed, 110 insertions(+), 122 deletions(-) diff --git a/Detector/Hit.hh b/Detector/Hit.hh index 56c5aa8b..e4f96ce1 100644 --- a/Detector/Hit.hh +++ b/Detector/Hit.hh @@ -35,7 +35,7 @@ namespace KinKal { virtual void print(std::ostream& ost=std::cout,int detail=0) const = 0; // update to a new reference, without changing internal state virtual void updateReference(PTRAJ const& ptraj) = 0; - virtual KTRAJPTR const& refTrajPtr() const = 0; + virtual KTRAJPTR refTrajPtr() const = 0; // update the internals of the hit, specific to this meta-iteraion virtual void updateState(MetaIterConfig const& config,bool first) = 0; // The following provides the constraint/information content of this hit in the trajectory weight space diff --git a/Detector/ParameterHit.hh b/Detector/ParameterHit.hh index 39831f05..f4125e2b 100644 --- a/Detector/ParameterHit.hh +++ b/Detector/ParameterHit.hh @@ -40,7 +40,7 @@ namespace KinKal { // parameter constraints are absolute and can't be updated void print(std::ostream& ost=std::cout,int detail=0) const override; void updateReference(PTRAJ const& ptraj) override { reftraj_ = ptraj.nearestTraj(time()); } - KTRAJPTR const& refTrajPtr() const override { return reftraj_; } + KTRAJPTR refTrajPtr() const override { return reftraj_; } // ParameterHit-specfic interface // construct from constraint values, time, and mask of which parameters to constrain ParameterHit(double time, PTRAJ const& ptraj, Parameters const& params, PMASK const& pmask); diff --git a/Examples/ScintHit.hh b/Examples/ScintHit.hh index 84ec4595..c845310c 100644 --- a/Examples/ScintHit.hh +++ b/Examples/ScintHit.hh @@ -17,18 +17,13 @@ namespace KinKal { using RESIDHIT = ResidualHit; using HIT = Hit; using KTRAJPTR = std::shared_ptr; + using STRAJPTR = std::shared_ptr; // copy constructor ScintHit(ScintHit const& rhs): ResidualHit(rhs), - saxis_(rhs.sensorAxis()), tvar_(rhs.timeVariance()), wvar_(rhs.widthVariance()), - tpca_( - rhs.closestApproach().particleTraj(), - saxis_, - rhs.closestApproach().hint(), - rhs.closestApproach().precision() - ), + tpca_(rhs.tpca_), rresid_(rhs.refResidual()){ /**/ }; @@ -46,20 +41,19 @@ namespace KinKal { Residual const& refResidual(unsigned ires=0) const override; double time() const override { return tpca_.particleToca(); } void updateReference(PTRAJ const& ptraj) override; - KTRAJPTR const& refTrajPtr() const override { return tpca_.particleTrajPtr(); } + KTRAJPTR refTrajPtr() const override { return tpca_.particleTrajPtr(); } void updateState(MetaIterConfig const& config,bool first) override; void print(std::ostream& ost=std::cout,int detail=0) const override; // scintHit explicit interface ScintHit(PCA const& pca, double tvar, double wvar); virtual ~ScintHit(){} // the line encapsulates both the measurement value (through t0), and the light propagation model (through the velocity) - auto const& sensorAxis() const { return saxis_; } + auto sensorAxis() const { return tpca_.sensorTrajPtr(); } auto const& closestApproach() const { return tpca_; } double timeVariance() const { return tvar_; } double widthVariance() const { return wvar_; } auto precision() const { return tpca_.precision(); } private: - SensorLine saxis_; // symmetry axis of this sensor double tvar_; // variance in the time measurement: assumed independent of propagation distance/time double wvar_; // variance in transverse position of the sensor/measurement in mm. Assumes cylindrical error, could be more general CA tpca_; // reference time and position of closest approach to the axis @@ -70,7 +64,7 @@ namespace KinKal { }; template ScintHit::ScintHit(PCA const& pca, double tvar, double wvar) : - saxis_(pca.sensorTraj()), tvar_(tvar), wvar_(wvar), + tvar_(tvar), wvar_(wvar), tpca_(static_cast(pca)) {} @@ -81,8 +75,8 @@ namespace KinKal { template void ScintHit::updateReference(PTRAJ const& ptraj) { // use previous hint, or initialize from the sensor time - CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(saxis_.measurementTime(), saxis_.measurementTime()); - PCA pca(ptraj,saxis_,tphint,precision()); + CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(sensorAxis()->measurementTime(), sensorAxis()->measurementTime()); + PCA pca(ptraj,sensorAxis(),tphint,precision()); tpca_ = static_cast(pca); if(!tpca_.usable())throw std::runtime_error("ScintHit TPOCA failure"); } @@ -92,26 +86,26 @@ namespace KinKal { // early in the fit when t0 has very large errors. // If it is unphysical try to adjust it back using a better hint. auto ppos = tpca_.particlePoca().Vect(); - auto const& sstart = saxis_.start(); - auto const& send = saxis_.end(); + auto const& sstart = sensorAxis()->start(); + auto const& send = sensorAxis()->end(); // tolerance should come from the config. Should also test relative to the error. FIXME - double tol = saxis_.length()*1.0; - double sdist = (ppos - saxis_.middle()).Dot(saxis_.direction()); - if( (ppos-sstart).Dot(saxis_.direction()) < -tol || (ppos-send).Dot(saxis_.direction()) > tol) { + double tol = sensorAxis()->length()*1.0; + double sdist = (ppos - sensorAxis()->middle()).Dot(sensorAxis()->direction()); + if( (ppos-sstart).Dot(sensorAxis()->direction()) < -tol || (ppos-send).Dot(sensorAxis()->direction()) > tol) { // adjust hint to the middle and try agian - double sspeed = tpca_.particleTraj().velocity(tpca_.particleToca()).Dot(saxis_.direction()); + double sspeed = tpca_.particleTraj().velocity(tpca_.particleToca()).Dot(sensorAxis()->direction()); auto tphint = tpca_.hint(); tphint.particleToca_ -= sdist/sspeed; - tpca_ = CA(tpca_.particleTrajPtr(),saxis_,tphint,precision()); + tpca_ = CA(tpca_.particleTrajPtr(),sensorAxis(),tphint,precision()); // should check if this is still unphysical and disable the hit if so FIXME - sdist = (tpca_.particlePoca().Vect() - saxis_.middle()).Dot(saxis_.direction()); + sdist = (tpca_.particlePoca().Vect() - sensorAxis()->middle()).Dot(sensorAxis()->direction()); } // residual is just delta-T at CA. // the variance includes the measurement variance and the tranvserse size (which couples to the relative direction) // 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)); + double totvar = tvar_ + wvar_*dd2/(sensorAxis()->speed(sdist)*sensorAxis()->speed(sdist)*(1.0-dd2)); rresid_ = Residual(tpca_.deltaT(),totvar,0.0,tpca_.dTdP()); this->updateWeight(config); } @@ -124,7 +118,7 @@ namespace KinKal { ost << " ScintHit tvar " << tvar_ << " wvar " << wvar_ << std::endl; if(detail > 0){ ost << "Line "; - saxis_.print(ost,detail); + sensorAxis()->print(ost,detail); } } diff --git a/Examples/SimpleWireHit.hh b/Examples/SimpleWireHit.hh index fda05ab6..02e9f237 100644 --- a/Examples/SimpleWireHit.hh +++ b/Examples/SimpleWireHit.hh @@ -22,6 +22,8 @@ namespace KinKal { using CA = ClosestApproach; using KTRAJPTR = std::shared_ptr; using PTRAJ = ParticleTrajectory; + using STRAJPTR = std::shared_ptr; + enum Dimension { dresid=0, tresid=1}; // residual dimensions SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca, @@ -32,15 +34,9 @@ namespace KinKal { ResidualHit(rhs), bfield_(rhs.bfield()), whstate_(rhs.hitState()), - wire_(rhs.wire()), - ca_( - rhs.closestApproach().particleTraj(), - wire_, - rhs.closestApproach().hint(), - rhs.closestApproach().precision() - ), + tpca_(rhs.tpca_), rresid_(rhs.residuals()), - mindoca_(rhs.minDOCA()), + mindotpca_(rhs.minDOCA()), dvel_(driftVelocity()), tvar_(timeVariance()), tot_(rhs.tot()), @@ -57,11 +53,11 @@ namespace KinKal { rv->setClosestApproach(ca); return rv; }; - double time() const override { return ca_.particleToca(); } + double time() const override { return tpca_.particleToca(); } VEC3 dRdX(unsigned ires) const; Residual const& refResidual(unsigned ires=dresid) const override; void updateReference(PTRAJ const& ptraj) override; - KTRAJPTR const& refTrajPtr() const override { return ca_.particleTrajPtr(); } + KTRAJPTR refTrajPtr() const override { return tpca_.particleTrajPtr(); } void print(std::ostream& ost=std::cout,int detail=0) const override; // Use dedicated updater void updateState(MetaIterConfig const& config,bool first) override; @@ -70,29 +66,28 @@ namespace KinKal { virtual ~SimpleWireHit(){} double driftVelocity() const { return dvel_; } double timeVariance() const { return tvar_; } - double minDOCA() const { return mindoca_; } + double minDOCA() const { return mindotpca_; } int id() const { return id_; } CA unbiasedClosestApproach() const; - auto const& closestApproach() const { return ca_; } + auto const& closestApproach() const { return tpca_; } auto const& hitState() const { return whstate_; } - auto const& wire() const { return wire_; } auto const& bfield() const { return bfield_; } - auto precision() const { return ca_.precision(); } + auto precision() const { return tpca_.precision(); } auto const& residuals() const { return rresid_; } double tot() const { return tot_; } double totVariance() const { return totvar_; } + auto wire() const { return tpca_.sensorTrajPtr(); } private: BFieldMap const& bfield_; // drift calculation requires the BField for ExB effects WireHitState whstate_; // current state - SensorLine wire_; // local linear approximation to the wire of this hit, encoding all (local) position and time information. // the start time is the measurement time, the direction is from // the physical source of the signal (particle) to the measurement recording location (electronics), the direction magnitude // is the effective signal propagation velocity along the wire, and the time range describes the active wire length // (when multiplied by the propagation velocity). - CA ca_; // reference time and position of closest approach to the wire; this is generally biased by the hit + CA tpca_; // reference time and position of closest approach to the wire; this is generally biased by the hit std::array rresid_; // residuals WRT most recent reference - double mindoca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties + double mindotpca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties double dvel_; // constant drift speed double tvar_; // constant time variance double tot_, totvar_; // TimeOverThreshold and variance @@ -101,7 +96,7 @@ namespace KinKal { void updateResiduals(); // modifiers to support cloning - void setClosestApproach(const CA& ca){ ca_ = ca; } + void setClosestApproach(const CA& ca){ tpca_ = ca; } }; //trivial 'updater' that sets the wire hit state to null @@ -113,18 +108,18 @@ namespace KinKal { template SimpleWireHit::SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca, double driftspeed, double tvar, double tot, double totvar, double rcell, int id) : bfield_(bfield), - whstate_(whstate), wire_(pca.sensorTraj()), - ca_(static_cast(pca)), - mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) { + whstate_(whstate), + tpca_(static_cast(pca)), + mindotpca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) { } template void SimpleWireHit::updateReference(PTRAJ const& ptraj) { // if we already computed PCA in the previous iteration, use that to set the hint. This speeds convergence // otherwise use the time at the center of the wire - CAHint tphint = ca_.usable() ? ca_.hint() : CAHint(wire_.timeAtMidpoint(),wire_.timeAtMidpoint()); - PCA pca(ptraj,wire_,tphint,precision()); - ca_ = static_cast(pca); - if(!ca_.usable())throw std::runtime_error("WireHit TPOCA failure"); + CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(wire()->timeAtMidpoint(),wire()->timeAtMidpoint()); + PCA pca(ptraj,wire(),tphint,precision()); + tpca_ = static_cast(pca); + if(!tpca_.usable())throw std::runtime_error("WireHit TPOCA failure"); } template void SimpleWireHit::updateState(MetaIterConfig const& miconfig, bool first) { @@ -134,12 +129,12 @@ namespace KinKal { auto dwhu = miconfig.findUpdater(); if(nwhu != 0 && dwhu != 0)throw std::invalid_argument(">1 SimpleWireHit updater specified"); if(nwhu != 0){ - mindoca_ = cellRadius(); + mindotpca_ = cellRadius(); whstate_ = nwhu->wireHitState(); // set the residuals based on this state } else if(dwhu != 0){ // update minDoca (for null ambiguity error estimate) - mindoca_ = std::min(dwhu->minDOCA(),cellRadius()); + mindotpca_ = std::min(dwhu->minDOCA(),cellRadius()); // compute the unbiased closest approach. This is brute-force // a more clever solution is to linearly correct the residuals for the change in parameters auto uca = this->unbiasedClosestApproach(); @@ -148,16 +143,16 @@ namespace KinKal { } rresid_[tresid] = rresid_[dresid] = Residual(); if(whstate_.active()){ - rresid_[tresid] = Residual(ca_.deltaT() - tot_, totvar_,0.0,ca_.dTdP()); // always constrain to TOT; this stabilizes the fit + rresid_[tresid] = Residual(tpca_.deltaT() - tot_, totvar_,0.0,tpca_.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(); + double dr = dvel_*whstate_.lrSign()*tpca_.deltaT() -tpca_.doca(); + DVEC dRdP = dvel_*whstate_.lrSign()*tpca_.dTdP() -tpca_.dDdP(); 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,ca_.dDdP()); + double nulldvar = dvel_*dvel_*(tpca_.deltaT()*tpca_.deltaT()+0.8); + rresid_[dresid] = Residual(tpca_.doca(),nulldvar,0.0,tpca_.dDdP()); } } // now update the weight @@ -168,9 +163,9 @@ namespace KinKal { if (whstate_.active()){ if (ires == dresid){ if (whstate_.useDrift()){ - return ca_.lSign()*ca_.delta().Vect().Unit(); + return tpca_.lSign()*tpca_.delta().Vect().Unit(); }else{ - return -1*ca_.lSign()*ca_.delta().Vect().Unit(); + return -1*tpca_.lSign()*tpca_.delta().Vect().Unit(); } } } @@ -186,7 +181,7 @@ namespace KinKal { // compute the unbiased closest approach; this is brute force, but works auto const& ca = this->closestApproach(); auto uparams = HIT::unbiasedParameters(); - KTRAJ utraj(uparams,ca.particleTraj()); + auto utraj = std::make_shared(uparams,ca.particleTraj()); return CA(utraj,this->wire(),ca.hint(),ca.precision()); } @@ -214,7 +209,7 @@ namespace KinKal { ost << std::endl; } if(detail > 1) { - ost << "Approximate Propagation speed " << wire_.speed(100) << " TPOCA " << ca_.tpData() << std::endl; + ost << "Approximate Propagation speed " << wire()->speed(100) << " TPOCA " << tpca_.tpData() << std::endl; } } diff --git a/Examples/StrawXing.hh b/Examples/StrawXing.hh index 3344ab0c..c342fdcc 100644 --- a/Examples/StrawXing.hh +++ b/Examples/StrawXing.hh @@ -15,6 +15,7 @@ namespace KinKal { public: using PTRAJ = ParticleTrajectory; using KTRAJPTR = std::shared_ptr; + using STRAJPTR = std::shared_ptr; using EXING = ElementXing; using PCA = PiecewiseClosestApproach; using CA = ClosestApproach; @@ -24,14 +25,8 @@ namespace KinKal { // copy constructor StrawXing(StrawXing const& rhs): ElementXing(rhs), - axis_(rhs.axis_), smat_(rhs.smat_), - tpca_( - rhs.closestApproach().particleTraj(), - axis_, - rhs.closestApproach().hint(), - rhs.closestApproach().precision() - ), + tpca_(rhs.tpca_), toff_(rhs.toff_), sxconfig_(rhs.sxconfig_), varscale_(rhs.varscale_), @@ -60,12 +55,12 @@ namespace KinKal { bool active() const override { return mxings_.size() > 0; } // accessors auto const& closestApproach() const { return tpca_; } + auto axis() const { return tpca_.sensorTrajPtr(); } auto const& strawMaterial() const { return smat_; } auto const& config() const { return sxconfig_; } auto precision() const { return tpca_.precision(); } private: - SensorLine axis_; // straw axis, expressed as a timeline StrawMaterial const& smat_; CA tpca_; // result of most recent TPOCA double toff_; // small time offset @@ -79,7 +74,6 @@ namespace KinKal { }; template StrawXing::StrawXing(PCA const& pca, StrawMaterial const& smat) : - axis_(pca.sensorTraj()), smat_(smat), tpca_(static_cast(pca)), toff_(smat.wireRadius()/pca.particleTraj().speed(pca.particleToca())), // locate the effect to 1 side of the wire to avoid overlap with hits @@ -87,8 +81,8 @@ namespace KinKal { {} template void StrawXing::updateReference(PTRAJ const& ptraj) { - CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(axis_.timeAtMidpoint(),axis_.timeAtMidpoint()); - PCA pca(ptraj,axis_,tphint,precision()); + CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(axis()->timeAtMidpoint(),axis()->timeAtMidpoint()); + PCA pca(ptraj,axis(),tphint,precision()); tpca_ = static_cast(pca); if(!tpca_.usable())throw std::runtime_error("StrawXing TPOCA failure"); } @@ -131,7 +125,7 @@ namespace KinKal { } if(detail > 1){ ost << " Axis "; - axis_.print(ost,0); + axis()->print(ost,0); } ost << std::endl; } diff --git a/Tests/ClosestApproachTest.hh b/Tests/ClosestApproachTest.hh index 99e796b4..531ae496 100644 --- a/Tests/ClosestApproachTest.hh +++ b/Tests/ClosestApproachTest.hh @@ -118,15 +118,15 @@ int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ if(tr_.Uniform(-1.0,1.0) < 0.0)cost *= -1.0; double sint = sqrt(1.0-cost*cost); MOM4 momv(mom*sint*cos(phi),mom*sint*sin(phi),mom*cost,pmass); - KTRAJ ktraj(origin,momv,icharge,bnom); + auto ktraj = std::make_shared(origin,momv,icharge,bnom); for(unsigned itime=0;itime < ntstep;itime++){ double time = tmin + itime*(tmax-tmin)/(ntstep-1); // create tline perp to trajectory at the specified time, separated by the specified gap VEC3 ppos, pdir; - ppos = ktraj.position3(time); - pdir = ktraj.direction(time); - VEC3 perp1 = ktraj.direction(time,MomBasis::perpdir_); - VEC3 perp2 = ktraj.direction(time,MomBasis::phidir_); + ppos = ktraj->position3(time); + pdir = ktraj->direction(time); + VEC3 perp1 = ktraj->direction(time,MomBasis::perpdir_); + VEC3 perp2 = ktraj->direction(time,MomBasis::phidir_); // choose a specific direction for DOCA // the line traj must be perp. to this and perp to the track double eta = tr_.Uniform(-3.14,3.14); @@ -140,7 +140,7 @@ int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ double gap = tr_.Uniform(0.01,maxgap); VEC3 spos = ppos + gap*docadir; // create the Line - SensorLine tline(spos, time, svel, wlen); + auto tline = std::make_shared(spos, time, svel, wlen); // create ClosestApproach from these CAHint tphint(time,time); TCA tp(ktraj,tline,tphint,1e-8); @@ -160,14 +160,14 @@ int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ } // test PointClosestApproach VEC4 pt(spos.X(),spos.Y(),spos.Z(),time-1.0); - TCAP tpp(ktraj,pt,1e-8); + TCAP tpp(*ktraj,pt,1e-8); if(fabs(fabs(tpp.doca()) - gap) > 1e-8){ cout << "Point DOCA not correct, doca = " << tpp.doca() << " gap " << gap << endl; status = 3; } // test against a piece-traj - PTRAJ ptraj(ktraj); + PTRAJ ptraj(*ktraj); PCA pca(ptraj,tline,tphint,1e-8); if(tp.status() != ClosestApproachData::converged)cout << "ClosestApproach status " << tp.statusName() << " doca " << tp.doca() << " dt " << tp.deltaT() << endl; if(tpp.status() != ClosestApproachData::converged)cout << "PointClosestApproach status " << tpp.statusName() << " doca " << tpp.doca() << " dt " << tpp.deltaT() << endl; @@ -187,11 +187,11 @@ int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ double dstart = -0.5*pchange[ipar]; for(unsigned istep=0;istepparams().parameters(); double dpar = dstart + dstep*istep; dvec[ipar] += dpar; - Parameters pdata(dvec,ktraj.params().covariance()); - KTRAJ dktraj(pdata,ktraj); + Parameters pdata(dvec,ktraj->params().covariance()); + auto dktraj = std::make_shared(pdata,*ktraj); TCA dtp(dktraj,tline,tphint,1e-9); double xd = dtp.doca(); double xt = dtp.deltaT(); diff --git a/Tests/FitTest.hh b/Tests/FitTest.hh index f2762718..c047c3c6 100644 --- a/Tests/FitTest.hh +++ b/Tests/FitTest.hh @@ -517,11 +517,11 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { WIREHITPTR shptr = std::dynamic_pointer_cast (thit); SCINTHITPTR lhptr = std::dynamic_pointer_cast (thit); if(shptr.use_count() > 0){ - auto const& tline = shptr->wire(); - plow = tline.start(); - phigh = tline.end(); + auto tline = shptr->wire(); + plow = tline->start(); + phigh = tline->end(); line->SetLineColor(kRed); - auto hitpos = tline.position3(shptr->closestApproach().sensorToca()); + auto hitpos = tline->position3(shptr->closestApproach().sensorToca()); auto trkpos = fptraj.position3(shptr->closestApproach().particleToca()); hpos->SetPoint(1,hitpos.X(),hitpos.Y(),hitpos.Z()); hpos->SetMarkerColor(kRed); @@ -529,10 +529,10 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { tpos->SetMarkerColor(kGreen); } else if (lhptr.use_count() > 0){ auto const& tline = lhptr->sensorAxis(); - plow = tline.start(); - phigh = tline.end(); + plow = tline->start(); + phigh = tline->end(); line->SetLineColor(kCyan); - auto hitpos = tline.position3(lhptr->closestApproach().sensorToca()); + auto hitpos = tline->position3(lhptr->closestApproach().sensorToca()); auto trkpos = fptraj.position3(lhptr->closestApproach().particleToca()); hpos->SetPoint(1,hitpos.X(),hitpos.Y(),hitpos.Z()); hpos->SetMarkerColor(kCyan); diff --git a/Tests/HitTest.hh b/Tests/HitTest.hh index c5b6d32d..9650cc7e 100644 --- a/Tests/HitTest.hh +++ b/Tests/HitTest.hh @@ -190,14 +190,14 @@ int HitTest(int argc, char **argv, const vector& delpars) { WIREHITPTR shptr = std::dynamic_pointer_cast (thit); SCINTHITPTR lhptr = std::dynamic_pointer_cast (thit); if((bool)shptr){ - auto const& tline = shptr->wire(); - plow = tline.start(); - phigh = tline.end(); + auto tline = shptr->wire(); + plow = tline->start(); + phigh = tline->end(); line->SetLineColor(kRed); } else if ((bool)lhptr){ auto const& tline = lhptr->sensorAxis(); - plow = tline.start(); - phigh = tline.end(); + plow = tline->start(); + phigh = tline->end(); line->SetLineColor(kCyan); } line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); diff --git a/Tests/ParticleTrajectoryTest.hh b/Tests/ParticleTrajectoryTest.hh index 202e3370..85d89567 100644 --- a/Tests/ParticleTrajectoryTest.hh +++ b/Tests/ParticleTrajectoryTest.hh @@ -38,6 +38,7 @@ using namespace KinKal; using namespace std; // avoid confusion with root using KinKal::SensorLine; +using STRAJPTR = std::shared_ptr; void print_usage() { printf("Usage: ParticleTrajectoryTest --changedir i --delta f --tstep f --nsteps i \n"); @@ -225,7 +226,7 @@ int ParticleTrajectoryTest(int argc, char **argv) { // shift the position VEC3 lpos = midpos + gap*rdir; double wlen(1000.0); - SensorLine tline(lpos, ptraj.range().mid(), pvel, wlen); + auto tline = std::make_shared(lpos, ptraj.range().mid(), pvel, wlen); // create ClosestApproach from these CAHint tphint(ptraj.range().mid(),0.0); PCA pca(ptraj,tline, tphint,1e-8); @@ -240,8 +241,8 @@ int ParticleTrajectoryTest(int argc, char **argv) { // draw the line and ClosestApproach TPolyLine3D* line = new TPolyLine3D(2); VEC3 plow, phigh; - plow = tline.start(); - phigh = tline.end(); + plow = tline->start(); + phigh = tline->end(); line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); line->SetPoint(1,phigh.X(),phigh.Y(), phigh.Z()); line->SetLineColor(kOrange); diff --git a/Tests/ToyMC.hh b/Tests/ToyMC.hh index d0891748..b86558b2 100644 --- a/Tests/ToyMC.hh +++ b/Tests/ToyMC.hh @@ -39,6 +39,7 @@ namespace KKTest { using STRAWXING = StrawXing; using STRAWXINGPTR = std::shared_ptr; using PCA = PiecewiseClosestApproach; + using STRAJPTR = std::shared_ptr; // create from aseed ToyMC(BFieldMap const& bfield, double mom, int icharge, double zrange, int iseed, unsigned nhits, bool simmat, bool scinthit, double ambigdoca ,double simmass) : bfield_(bfield), matdb_(sfinder_), @@ -55,7 +56,7 @@ namespace KKTest { } // generate a straw at the given time. direction and drift distance are random - SensorLine generateStraw(PTRAJ const& traj, double htime); + STRAJPTR generateStraw(PTRAJ const& traj, double htime); // create a seed by randomizing the parameters void createSeed(KTRAJ& seed,DVEC const& sigmas, double seedsmear); void extendTraj(PTRAJ& ptraj,double htime); @@ -107,7 +108,7 @@ namespace KKTest { MetaIterConfig miconfig_; // configuration used when calculating initial effect }; - template SensorLine ToyMC::generateStraw(PTRAJ const& traj, double htime) { + template std::shared_ptr ToyMC::generateStraw(PTRAJ const& traj, double htime) { // start with the true helix position at this time auto hpos = traj.position4(htime); auto tdir = traj.direction(htime); @@ -127,7 +128,7 @@ namespace KKTest { // smear measurement time tmeas = tr_.Gaus(tmeas,sigt_); // construct the trajectory for this hit - return SensorLine(mpos,tmeas,vprop,wlen_); + return std::make_shared(mpos,tmeas,vprop,wlen_); } template void ToyMC::simulateParticle(PTRAJ& ptraj,HITCOL& thits, EXINGCOL& dxings, bool addmat) { @@ -224,7 +225,7 @@ namespace KKTest { double tmeas = tr_.Gaus(shmaxtime+(shmeas.Z() - shmax.Z())/cprop_,scitsig_); // create the ttraj for the light propagation VEC3 lvel(0.0,0.0,cprop_); - SensorLine lline(shmeas, tmeas, lvel, clen_); + auto lline = std::make_shared(shmeas, tmeas, lvel, clen_); // then create the hit and add it; the hit has no material CAHint tphint(tmeas,tmeas); PCA pca(ptraj,lline,tphint,tprec_); diff --git a/Trajectory/ClosestApproach.hh b/Trajectory/ClosestApproach.hh index 4401b776..72e5aaa3 100644 --- a/Trajectory/ClosestApproach.hh +++ b/Trajectory/ClosestApproach.hh @@ -25,24 +25,26 @@ namespace KinKal { template class ClosestApproach { public: using KTRAJPTR = std::shared_ptr; + using STRAJPTR = std::shared_ptr; // construct from the particle and sensor trajectories; TCA is computed on construction, given a hint as to where // to start looking, which disambiguates functions with multiple solutions ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, CAHint const& hint, double precision); // same, using Ptrs - ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, CAHint const& hint, double precision); + ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, CAHint const& hint, double precision); // construct without a hint: TCA isn't calculated, state is invalid ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, double precision); - ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, double precision); + ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR straj, double precision); // explicitly construct from all content (no calculation) - ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, double precision, + ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double precision, ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP); // copy constructor ClosestApproach(ClosestApproach const&) = default; // accessors ClosestApproachData const& tpData() const { return tpdata_; } KTRAJ const& particleTraj() const { return *ktrajptr_; } - KTRAJPTR const& particleTrajPtr() const { return ktrajptr_; } - STRAJ const& sensorTraj() const { return straj_; } + KTRAJPTR particleTrajPtr() const { return ktrajptr_; } + STRAJ const& sensorTraj() const { return *strajptr_; } + STRAJPTR sensorTrajPtr() const { return strajptr_; } // derviatives of TOCA and DOCA WRT particle trajectory parameters DVEC const& dDdP() const { return dDdP_; } DVEC const& dTdP() const { return dTdP_; } @@ -78,7 +80,7 @@ namespace KinKal { DVEC dTdP_; // derivative of TOCA WRT Parameters protected: KTRAJPTR ktrajptr_; // kinematic particle trajectory - STRAJ const& straj_; // sensor trajectory This should be a shared ptr TODO + STRAJPTR strajptr_; // sensor trajectory // calculate CA given the hint, and fill the state void findTCA(CAHint const& hint); ClosestApproachData tpdata_; // data payload of CA calculation @@ -87,22 +89,22 @@ namespace KinKal { }; template ClosestApproach::ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, double prec) : - precision_(prec),ktrajptr_(new KTRAJ(ktraj)), straj_(straj) {} + precision_(prec),ktrajptr_(new KTRAJ(ktraj)), strajptr_(new STRAJ(straj)) {} - template ClosestApproach::ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, double prec) : - precision_(prec),ktrajptr_(ktrajptr), straj_(straj) {} + template ClosestApproach::ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double prec) : + precision_(prec),ktrajptr_(ktrajptr), strajptr_(strajptr) {} - template ClosestApproach::ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, double prec, + template ClosestApproach::ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double prec, ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP) : - precision_(prec),ktrajptr_(ktrajptr), straj_(straj), tpdata_(tpdata),dDdP_(dDdP), dTdP_(dTdP) {} + precision_(prec),ktrajptr_(ktrajptr), strajptr_(strajptr), tpdata_(tpdata),dDdP_(dDdP), dTdP_(dTdP) {} template ClosestApproach::ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, CAHint const& hint, double prec) : ClosestApproach(ktraj,straj,prec) { findTCA(hint); } - template ClosestApproach::ClosestApproach(KTRAJPTR const& ktrajptr, STRAJ const& straj, CAHint const& hint, - double prec) : ClosestApproach(ktrajptr,straj,prec) { + template ClosestApproach::ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, CAHint const& hint, + double prec) : ClosestApproach(ktrajptr,strajptr,prec) { findTCA(hint); } @@ -113,7 +115,7 @@ namespace KinKal { dTdP_ = other.dTdP(); ktrajptr_ = other.ktrajptr_; // make sure the sensor traj is the same - if(&straj_ != &other.sensorTraj()) throw std::invalid_argument("Inconsistent ClosestApproach SensorTraj"); + if(strajptr_ != other.sensorTrajPtr()) throw std::invalid_argument("Inconsistent ClosestApproach SensorTraj"); return *this; } @@ -127,7 +129,7 @@ namespace KinKal { unsigned niter(0); // speed doesn't change double pspeed = ktrajptr_->speed(particleToca()); - double sspeed = straj_.speed(sensorToca()); + double sspeed = strajptr_->speed(sensorToca()); // iterate until change in TOCA is less than precision double dptoca(std::numeric_limits::max()), dstoca(std::numeric_limits::max()); while(tpdata_.usable() && (fabs(dptoca) > precision() || fabs(dstoca) > precision()) && niter++ < maxiter) { @@ -193,8 +195,8 @@ namespace KinKal { } template void ClosestApproach::setSensorTOCA(double stoca) { - tpdata_.sensCA_ = straj_.position4(stoca); - tpdata_.sdir_ = straj_.direction(stoca); + tpdata_.sensCA_ = strajptr_->position4(stoca); + tpdata_.sdir_ = strajptr_->direction(stoca); } diff --git a/Trajectory/PiecewiseClosestApproach.hh b/Trajectory/PiecewiseClosestApproach.hh index b98fe3fb..334e33df 100644 --- a/Trajectory/PiecewiseClosestApproach.hh +++ b/Trajectory/PiecewiseClosestApproach.hh @@ -13,10 +13,11 @@ namespace KinKal { using PTRAJ = ParticleTrajectory; using KTCA = ClosestApproach; using KTRAJPTR = std::shared_ptr; - PiecewiseClosestApproach(PTRAJ const& ptraj, STRAJ const& straj, CAHint const& hint, double precision); + using STRAJPTR = std::shared_ptr; + PiecewiseClosestApproach(PTRAJ const& ptraj, STRAJPTR strajptr, CAHint const& hint, double precision); }; - template PiecewiseClosestApproach::PiecewiseClosestApproach(ParticleTrajectory const& ptraj, STRAJ const& straj, CAHint const& hint, double prec) : KTCA(ptraj.nearestTraj(hint.particleToca_),straj,prec) { + template PiecewiseClosestApproach::PiecewiseClosestApproach(ParticleTrajectory const& ptraj, STRAJPTR strajptr, CAHint const& hint, double prec) : KTCA(ptraj.nearestTraj(hint.particleToca_),strajptr,prec) { // iteratively find the nearest piece, and CA for that piece. Start at hint static const unsigned maxiter=10; // don't allow infinite iteration. This should be a parameter TODO unsigned niter=0; From 37d0b3e08b5b6d6979b69aaece879a1937846e3e Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sun, 11 Jan 2026 07:27:03 -0800 Subject: [PATCH 3/4] Fix names --- Examples/SimpleWireHit.hh | 12 ++++++------ Trajectory/ClosestApproach.hh | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Examples/SimpleWireHit.hh b/Examples/SimpleWireHit.hh index 02e9f237..fb152671 100644 --- a/Examples/SimpleWireHit.hh +++ b/Examples/SimpleWireHit.hh @@ -36,7 +36,7 @@ namespace KinKal { whstate_(rhs.hitState()), tpca_(rhs.tpca_), rresid_(rhs.residuals()), - mindotpca_(rhs.minDOCA()), + mindoca_(rhs.minDOCA()), dvel_(driftVelocity()), tvar_(timeVariance()), tot_(rhs.tot()), @@ -66,7 +66,7 @@ namespace KinKal { virtual ~SimpleWireHit(){} double driftVelocity() const { return dvel_; } double timeVariance() const { return tvar_; } - double minDOCA() const { return mindotpca_; } + double minDOCA() const { return mindoca_; } int id() const { return id_; } CA unbiasedClosestApproach() const; auto const& closestApproach() const { return tpca_; } @@ -87,7 +87,7 @@ namespace KinKal { // (when multiplied by the propagation velocity). CA tpca_; // reference time and position of closest approach to the wire; this is generally biased by the hit std::array rresid_; // residuals WRT most recent reference - double mindotpca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties + double mindoca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties double dvel_; // constant drift speed double tvar_; // constant time variance double tot_, totvar_; // TimeOverThreshold and variance @@ -110,7 +110,7 @@ namespace KinKal { bfield_(bfield), whstate_(whstate), tpca_(static_cast(pca)), - mindotpca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) { + mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) { } template void SimpleWireHit::updateReference(PTRAJ const& ptraj) { @@ -129,12 +129,12 @@ namespace KinKal { auto dwhu = miconfig.findUpdater(); if(nwhu != 0 && dwhu != 0)throw std::invalid_argument(">1 SimpleWireHit updater specified"); if(nwhu != 0){ - mindotpca_ = cellRadius(); + mindoca_ = cellRadius(); whstate_ = nwhu->wireHitState(); // set the residuals based on this state } else if(dwhu != 0){ // update minDoca (for null ambiguity error estimate) - mindotpca_ = std::min(dwhu->minDOCA(),cellRadius()); + mindoca_ = std::min(dwhu->minDOCA(),cellRadius()); // compute the unbiased closest approach. This is brute-force // a more clever solution is to linearly correct the residuals for the change in parameters auto uca = this->unbiasedClosestApproach(); diff --git a/Trajectory/ClosestApproach.hh b/Trajectory/ClosestApproach.hh index 72e5aaa3..cb11c4f5 100644 --- a/Trajectory/ClosestApproach.hh +++ b/Trajectory/ClosestApproach.hh @@ -33,7 +33,7 @@ namespace KinKal { ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, CAHint const& hint, double precision); // construct without a hint: TCA isn't calculated, state is invalid ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, double precision); - ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR straj, double precision); + ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double precision); // explicitly construct from all content (no calculation) ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double precision, ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP); From 8faca387884a03228135e13a6f53d57659490317 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Sun, 11 Jan 2026 09:36:03 -0800 Subject: [PATCH 4/4] Remove constref constructors --- Tests/Trajectory.hh | 9 +++++---- Trajectory/ClosestApproach.hh | 11 ----------- 2 files changed, 5 insertions(+), 15 deletions(-) diff --git a/Tests/Trajectory.hh b/Tests/Trajectory.hh index 93077a67..2f66f22a 100644 --- a/Tests/Trajectory.hh +++ b/Tests/Trajectory.hh @@ -276,16 +276,17 @@ int TrajectoryTest(int argc, char **argv,KinKal::DVEC sigmas) { // shift the position VEC3 perpdir(-sin(phi),cos(phi),0.0); VEC3 ppos = pos + gap*perpdir; - SensorLine tline(ppos, ltime, pvel, wlen); + auto tline = std::make_shared(ppos, ltime, pvel, wlen); // find ClosestApproach CAHint hint(ltime,ltime); - ClosestApproach tp(ktraj,tline,hint, 1e-6); + auto ktrajptr = std::make_shared(ktraj); + ClosestApproach tp(ktrajptr,tline,hint, 1e-6); // cout << "ClosestApproach status " << tp.statusName() << " doca " << tp.doca() << " dt " << tp.deltaT() << endl; if(tp.status() == ClosestApproachData::converged) { // draw the line and ClosestApproach TPolyLine3D* line = new TPolyLine3D(2); - auto plow = tline.start(); - auto phigh = tline.end(); + auto plow = tline->start(); + auto phigh = tline->end(); line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); line->SetPoint(1,phigh.X(),phigh.Y(), phigh.Z()); line->SetLineColor(kOrange); diff --git a/Trajectory/ClosestApproach.hh b/Trajectory/ClosestApproach.hh index cb11c4f5..da5ac9df 100644 --- a/Trajectory/ClosestApproach.hh +++ b/Trajectory/ClosestApproach.hh @@ -28,11 +28,8 @@ namespace KinKal { using STRAJPTR = std::shared_ptr; // construct from the particle and sensor trajectories; TCA is computed on construction, given a hint as to where // to start looking, which disambiguates functions with multiple solutions - ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, CAHint const& hint, double precision); - // same, using Ptrs ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, CAHint const& hint, double precision); // construct without a hint: TCA isn't calculated, state is invalid - ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, double precision); ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double precision); // explicitly construct from all content (no calculation) ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double precision, @@ -88,9 +85,6 @@ namespace KinKal { void setSensorTOCA( double stoca); }; - template ClosestApproach::ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, double prec) : - precision_(prec),ktrajptr_(new KTRAJ(ktraj)), strajptr_(new STRAJ(straj)) {} - template ClosestApproach::ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, double prec) : precision_(prec),ktrajptr_(ktrajptr), strajptr_(strajptr) {} @@ -98,11 +92,6 @@ namespace KinKal { ClosestApproachData const& tpdata, DVEC const& dDdP, DVEC const& dTdP) : precision_(prec),ktrajptr_(ktrajptr), strajptr_(strajptr), tpdata_(tpdata),dDdP_(dDdP), dTdP_(dTdP) {} - template ClosestApproach::ClosestApproach(KTRAJ const& ktraj, STRAJ const& straj, CAHint const& hint, - double prec) : ClosestApproach(ktraj,straj,prec) { - findTCA(hint); - } - template ClosestApproach::ClosestApproach(KTRAJPTR ktrajptr, STRAJPTR strajptr, CAHint const& hint, double prec) : ClosestApproach(ktrajptr,strajptr,prec) { findTCA(hint);