Skip to content

Use trajectory ptrs in ClosestApproach #215

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
brownd1978 merged 4 commits into from
Jan 12, 2026
Merged

Use trajectory ptrs in ClosestApproach #215

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ae4e3eb
Convert PiecewiseClosestApproach to use KTRAJ instead of PKTRAJ to av…
brownd1978 Jan 9, 2026
51c26a9
Make SensorTraj a ptr in ClosestApproach
brownd1978 Jan 10, 2026
37d0b3e
Fix names
brownd1978 Jan 11, 2026
8faca38
Remove constref constructors
brownd1978 Jan 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Detector/Hit.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion Detector/ParameterHit.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
44 changes: 19 additions & 25 deletions Examples/ScintHit.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,18 +17,13 @@ namespace KinKal {
using RESIDHIT = ResidualHit<KTRAJ>;
using HIT = Hit<KTRAJ>;
using KTRAJPTR = std::shared_ptr<KTRAJ>;
using STRAJPTR = std::shared_ptr<SensorLine>;
// copy constructor
ScintHit(ScintHit<KTRAJ> const& rhs):
ResidualHit<KTRAJ>(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()){
/**/
};
Expand All @@ -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
Expand All @@ -70,8 +64,8 @@ namespace KinKal {
};

template <class KTRAJ> ScintHit<KTRAJ>::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())
tvar_(tvar), wvar_(wvar),
tpca_(static_cast<CA const&>(pca))
{}

template <class KTRAJ> Residual const& ScintHit<KTRAJ>::refResidual(unsigned ires) const {
Expand All @@ -81,9 +75,9 @@ namespace KinKal {

template <class KTRAJ> void ScintHit<KTRAJ>::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());
tpca_ = pca.localClosestApproach();
CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(sensorAxis()->measurementTime(), sensorAxis()->measurementTime());
PCA pca(ptraj,sensorAxis(),tphint,precision());
tpca_ = static_cast<CA>(pca);
if(!tpca_.usable())throw std::runtime_error("ScintHit TPOCA failure");
}

Expand All @@ -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);
}
Expand All @@ -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);
}
}

Expand Down
55 changes: 25 additions & 30 deletions Examples/SimpleWireHit.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ namespace KinKal {
using CA = ClosestApproach<KTRAJ,SensorLine>;
using KTRAJPTR = std::shared_ptr<KTRAJ>;
using PTRAJ = ParticleTrajectory<KTRAJ>;
using STRAJPTR = std::shared_ptr<SensorLine>;

enum Dimension { dresid=0, tresid=1}; // residual dimensions

SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca,
Expand All @@ -32,13 +34,7 @@ namespace KinKal {
ResidualHit<KTRAJ>(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()),
dvel_(driftVelocity()),
Expand All @@ -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;
Expand All @@ -73,24 +69,23 @@ namespace KinKal {
double minDOCA() const { return mindoca_; }
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<Residual,2> rresid_; // residuals WRT most recent reference
double mindoca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties
double dvel_; // constant drift speed
Expand All @@ -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
Expand All @@ -113,18 +108,18 @@ namespace KinKal {
template <class KTRAJ> SimpleWireHit<KTRAJ>::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_(pca.localTraj(),wire_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()), // must be explicit to get the right sensor traj reference
whstate_(whstate),
tpca_(static_cast<CA const&>(pca)),
mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) {
}

template <class KTRAJ> void SimpleWireHit<KTRAJ>::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_ = pca.localClosestApproach();
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<CA>(pca);
if(!tpca_.usable())throw std::runtime_error("WireHit TPOCA failure");
}

template <class KTRAJ> void SimpleWireHit<KTRAJ>::updateState(MetaIterConfig const& miconfig, bool first) {
Expand All @@ -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
Expand All @@ -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();
}
}
}
Expand All @@ -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<KTRAJ>(uparams,ca.particleTraj());
return CA(utraj,this->wire(),ca.hint(),ca.precision());
}

Expand Down Expand Up @@ -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;
}
}

Expand Down
22 changes: 8 additions & 14 deletions Examples/StrawXing.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ namespace KinKal {
public:
using PTRAJ = ParticleTrajectory<KTRAJ>;
using KTRAJPTR = std::shared_ptr<KTRAJ>;
using STRAJPTR = std::shared_ptr<SensorLine>;
using EXING = ElementXing<KTRAJ>;
using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
using CA = ClosestApproach<KTRAJ,SensorLine>;
Expand All @@ -24,14 +25,8 @@ namespace KinKal {
// copy constructor
StrawXing(StrawXing const& rhs):
ElementXing<KTRAJ>(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_),
Expand Down Expand Up @@ -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
Expand All @@ -79,17 +74,16 @@ namespace KinKal {
};

template <class KTRAJ> StrawXing<KTRAJ>::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<CA const&>(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)
{}

template <class KTRAJ> void StrawXing<KTRAJ>::updateReference(PTRAJ const& ptraj) {
CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(axis_.timeAtMidpoint(),axis_.timeAtMidpoint());
PCA pca(ptraj,axis_,tphint,precision());
tpca_ = pca.localClosestApproach();
CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(axis()->timeAtMidpoint(),axis()->timeAtMidpoint());
PCA pca(ptraj,axis(),tphint,precision());
tpca_ = static_cast<CA>(pca);
if(!tpca_.usable())throw std::runtime_error("StrawXing TPOCA failure");
}

Expand Down Expand Up @@ -131,7 +125,7 @@ namespace KinKal {
}
if(detail > 1){
ost << " Axis ";
axis_.print(ost,0);
axis()->print(ost,0);
}
ost << std::endl;
}
Expand Down
Loading