From bc9a12d7bec460c8638c88e33034a78592061b6f Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 17 Jun 2024 18:06:22 -0700 Subject: [PATCH] Improve segment accessors --- RecoDataProducts/inc/KalSeed.hh | 4 +- RecoDataProducts/inc/KalSegment.hh | 11 ++- RecoDataProducts/src/KalSeed.cc | 98 +++++---------------- RecoDataProducts/src/KalSegment.cc | 32 ++++--- TEveEventDisplay/src/TEveMu2eCustomHelix.cc | 2 +- 5 files changed, 49 insertions(+), 98 deletions(-) diff --git a/RecoDataProducts/inc/KalSeed.hh b/RecoDataProducts/inc/KalSeed.hh index 8f3160c53c..0e6d364290 100644 --- a/RecoDataProducts/inc/KalSeed.hh +++ b/RecoDataProducts/inc/KalSeed.hh @@ -53,7 +53,8 @@ namespace mu2e { KinKal::TimeRange timeRange() const { return KinKal::TimeRange(_segments.front()._tmin,_segments.back()._tmax); } bool hasCaloCluster() const { return _chit.caloCluster().isNonnull(); } art::Ptr const& caloCluster() const { return _chit.caloCluster(); } - std::vector::const_iterator nearestSeg(double time) const; + std::vector::const_iterator nearestSegment(double time) const; + std::vector::const_iterator t0Segment(double& t0) const; // return the segment associated with the t0 value, or the closest to it. Also return the t0 value std::vector::const_iterator intersection(SurfaceId const& id) const; bool loopHelixFit() const { return _status.hasAllProperties(TrkFitFlag::KKLoopHelix); } bool centralHelixFit() const { return _status.hasAllProperties(TrkFitFlag::KKCentralHelix); } @@ -81,7 +82,6 @@ namespace mu2e { std::vector _hits; // hit seeds for all the hits used in this fit std::vector _straws; // straws interesected by this fit TrkCaloHitSeed _chit; // CaloCluster-based hit. If it has no CaloCluster, this has no content - std::vector::const_iterator nearestSegment(float time) const; // // deprecated BTrk legacy content, DO NOT write any new code which depends on these functions // find the nearest segment to a given the time diff --git a/RecoDataProducts/inc/KalSegment.hh b/RecoDataProducts/inc/KalSegment.hh index ed0bcb5713..714fd3d3be 100644 --- a/RecoDataProducts/inc/KalSegment.hh +++ b/RecoDataProducts/inc/KalSegment.hh @@ -11,6 +11,7 @@ #include "Offline/RecoDataProducts/inc/HelixVal.hh" #include "Offline/DataProducts/inc/GenVector.hh" #include "Offline/RecoDataProducts/inc/HitT0.hh" +#include "Offline/RecoDataProducts/inc/TrkFitFlag.hh" #include "KinKal/General/ParticleState.hh" #include "KinKal/Trajectory/KinematicLine.hh" #include "KinKal/Trajectory/CentralHelix.hh" @@ -43,21 +44,19 @@ namespace mu2e { Float_t tmin() const { return _tmin; } Float_t tmax() const { return _tmax; } KinKal::TimeRange timeRange() const; - auto tref() const { return _pstate.time(); } + auto tref() const { return _pstate.time(); }// note this is NOT t0; it is the reference time at which the trajectory was sampled + double t0Val(TrkFitFlag const& flag = TrkFitFlag(TrkFitFlag::KKLoopHelix)) const;// this will give the local segment's t0, interpreted for the given trajectory type, which is not necessarily the same as the track t0 XYZVectorF const& bnom() const { return _bnom; } KinKal::VEC3 KKbnom() const { return KinKal::VEC3(_bnom); } -// convenience function - double t0Val() const; - Float_t _tmin, _tmax; // time range // main payload is the particle state estimate. this includes all the kinematic information to // interpret as anything else. BField is needed to interpret geometrically XYZVectorF _bnom; // Bfield associated with this segment, needed to reconstitute helix KinKal::ParticleStateEstimate _pstate; // particle state at this sample - // the following are deprecated legacy functions specific to the BTrk fit and will go away once we migrate to KinKal + // the following are deprecated legacy functions specific to the BTrk, these should be removed FIXME + HitT0 t0() const; HelixVal helix() const; HelixCov covar() const; - HitT0 t0() const; void mom(double flt, XYZVectorF& momvec) const; double fmin() const { return timeToFlt(_tmin); } // local 3D flight range double fmax() const { return timeToFlt(_tmax); } diff --git a/RecoDataProducts/src/KalSeed.cc b/RecoDataProducts/src/KalSeed.cc index f5cc80ff8f..9f91d3360b 100644 --- a/RecoDataProducts/src/KalSeed.cc +++ b/RecoDataProducts/src/KalSeed.cc @@ -63,7 +63,7 @@ namespace mu2e { return retval; } - std::vector::const_iterator KalSeed::nearestSeg(double time) const { + std::vector::const_iterator KalSeed::nearestSegment(double time) const { auto retval = segments().end(); float dmin(std::numeric_limits::max()); for(auto ikseg = segments().begin(); ikseg !=segments().end(); ikseg++) { @@ -98,20 +98,6 @@ namespace mu2e { } return retval; } - - std::vector::const_iterator KalSeed::nearestSegment(float time) const { - auto retval = segments().end(); - double tmin(std::numeric_limits::max()); - for(auto ikseg = segments().begin(); ikseg !=segments().end(); ikseg++) { - double dt = fabs(time - 0.5*(ikseg->tmin()+ikseg->tmax())); - if(dt < tmin){ - tmin = dt; - retval = ikseg; - } - } - return retval; - } - std::vector::const_iterator KalSeed::nearestSegment(const XYZVectorF& pos) const { auto retval = segments().end(); double dmin(std::numeric_limits::max()); @@ -125,74 +111,32 @@ namespace mu2e { return retval; } - double KalSeed::t0Val() const { + + std::vector::const_iterator KalSeed::t0Segment(double& t0) const { if(segments().size() <= 0) throw cet::exception("RECO")<<"mu2e::KalSeed: no segments" << std::endl; -// find the segment nearest z=0. If there's one segment, we are done +// start with the first segment auto iseg = segments().begin(); - if(_segments.size() > 1){ - auto pvel = iseg->state().velocity(); - double vz = pvel.Z(); // to a good approximation, B is along Z - auto pref = iseg->position3(); - double zmin = pref.Z() + vz*(iseg->tmin()-iseg->tref()); - double zmax = pref.Z() + vz*(iseg->tmax()-iseg->tref()); - if(zmin > 0.0 || zmax < 0.0){ - double mindz = std::min(fabs(zmin), fabs(zmax)); - // find the segment closest to z=0 - auto jseg = ++iseg; - while(jseg != segments().end()) { - pvel = jseg->state().velocity(); - vz = pvel.Z(); - pref = jseg->position3(); - zmin = pref.Z() + vz*(jseg->tmin()-jseg->tref()); - zmax = pref.Z() + vz*(jseg->tmax()-jseg->tref()); - double dz = std::min(fabs(zmin),fabs(zmax)); - if(zmin < 0.0 && zmax > 0.0){ - iseg = jseg; - break; - } else if(dz < mindz){ - iseg = jseg; - mindz = dz; - } - ++jseg; - } - } + auto oldseg = segments().end(); + t0 = iseg->t0Val(_status); + while(iseg != oldseg && iseg->timeRange().inRange(t0)){ + oldseg = iseg; + iseg = nearestSegment(t0); } - return iseg->t0Val(); + return iseg; + } + + double KalSeed::t0Val() const { + double t0; + t0Segment(t0); + return t0; } + // deprecated legacy functions + HitT0 KalSeed::t0() const { - if(segments().size() <= 0) throw cet::exception("RECO")<<"mu2e::KalSeed: no segments" << std::endl; -// find the segment nearest z=0. If there's one segment, we are done - auto iseg = segments().begin(); - if(_segments.size() > 1){ - auto pvel = iseg->state().velocity(); - double vz = pvel.Z(); // to a good approximation, B is along Z - auto pref = iseg->position3(); - double zmin = pref.Z() + vz*(iseg->tmin()-iseg->tref()); - double zmax = pref.Z() + vz*(iseg->tmax()-iseg->tref()); - if(zmin > 0.0 || zmax < 0.0){ - double mindz = std::min(fabs(zmin), fabs(zmax)); - // find the segment closest to z=0 - auto jseg = ++iseg; - while(jseg != segments().end()) { - pvel = jseg->state().velocity(); - vz = pvel.Z(); - pref = jseg->position3(); - zmin = pref.Z() + vz*(jseg->tmin()-jseg->tref()); - zmax = pref.Z() + vz*(jseg->tmax()-jseg->tref()); - double dz = std::min(fabs(zmin),fabs(zmax)); - if(zmin < 0.0 && zmax > 0.0){ - iseg = jseg; - break; - } else if(dz < mindz){ - iseg = jseg; - mindz = dz; - } - ++jseg; - } - } - } - return HitT0(iseg->t0Val(),1.0); //FIXME + double t0; + t0Segment(t0); + return HitT0(t0,-1.0); } } diff --git a/RecoDataProducts/src/KalSegment.cc b/RecoDataProducts/src/KalSegment.cc index ccf79b4b80..7bae6f4b27 100644 --- a/RecoDataProducts/src/KalSegment.cc +++ b/RecoDataProducts/src/KalSegment.cc @@ -12,6 +12,26 @@ namespace mu2e { } } + double KalSegment::t0Val(TrkFitFlag const& flag) const { + // convert to the appropriate parameterization. + if(flag.hasAllProperties(TrkFitFlag::KKLoopHelix)) { + auto lh = loopHelix(); + return lh.paramVal(KinKal::LoopHelix::t0Index()); + } else if(flag.hasAllProperties(TrkFitFlag::KKCentralHelix)) { + auto ch = loopHelix(); + return ch.paramVal(KinKal::CentralHelix::t0Index()); + } else if (flag.hasAllProperties(TrkFitFlag::KKLine)) { + auto kl = kinematicLine(); + return kl.paramVal(KinKal::KinematicLine::t0Index()); + } + // throw cet::exception("RECO")<<"mu2e::KalSegment: no trajectory specified in flag" << std::endl; + // for now, revert to the legacy implementation. Once BTrk is fully removed this should be removed + auto vel = _pstate.velocity(); + return _pstate.time() - _pstate.position3().Z()/vel.Z(); + } + + // deprecated legacy functions, these should be removed FIXME + HelixVal KalSegment::helix() const { // CentralHelix uses the same parameter convention as BTrk. First, // convert the state estimate into a helix. @@ -56,16 +76,4 @@ namespace mu2e { auto vel = _pstate.velocity(); return (time - t0Val())*vel.R(); } - - double KalSegment::t0Val() const { - auto vel = _pstate.velocity(); - return _pstate.time() - _pstate.position3().Z()/vel.Z(); - } - - HitT0 KalSegment::t0() const { - // convert to LoopHelix. In that parameterization, t0 is defined WRT z=0 - auto lhelix = loopHelix(); - return HitT0(lhelix.params().parameters()(KinKal::LoopHelix::t0_), - sqrt(lhelix.params().covariance()(KinKal::LoopHelix::t0_,KinKal::LoopHelix::t0_))); - } } diff --git a/TEveEventDisplay/src/TEveMu2eCustomHelix.cc b/TEveEventDisplay/src/TEveMu2eCustomHelix.cc index aa6e2cd743..2c0dc45ba3 100644 --- a/TEveEventDisplay/src/TEveMu2eCustomHelix.cc +++ b/TEveEventDisplay/src/TEveMu2eCustomHelix.cc @@ -26,7 +26,7 @@ namespace mu2e{ auto pos = fseg.position3(); auto vel = fseg.velocity(); double tz = fseg.tref() + (zpos-pos.Z())/vel.Z(); - auto zseg = fKalSeed_.nearestSeg(tz); + auto zseg = fKalSeed_.nearestSegment(tz); pos = zseg->position3(); vel = zseg->velocity(); tz = zseg->tref() + (zpos-pos.Z())/vel.Z();