Skip to content

Commit

Permalink
Improve segment accessors
Browse files Browse the repository at this point in the history
  • Loading branch information
@brownd1978
brownd1978 committed Jun 18, 2024
1 parent 0433a1e commit bc9a12d
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 98 deletions.
4 changes: 2 additions & 2 deletions RecoDataProducts/inc/KalSeed.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<CaloCluster> const& caloCluster() const { return _chit.caloCluster(); }
std::vector<KalSegment>::const_iterator nearestSeg(double time) const;
std::vector<KalSegment>::const_iterator nearestSegment(double time) const;
std::vector<KalSegment>::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<KalIntersection>::const_iterator intersection(SurfaceId const& id) const;
bool loopHelixFit() const { return _status.hasAllProperties(TrkFitFlag::KKLoopHelix); }
bool centralHelixFit() const { return _status.hasAllProperties(TrkFitFlag::KKCentralHelix); }
Expand Down Expand Up @@ -81,7 +82,6 @@ namespace mu2e {
std::vector<TrkStrawHitSeed> _hits; // hit seeds for all the hits used in this fit
std::vector<TrkStraw> _straws; // straws interesected by this fit
TrkCaloHitSeed _chit; // CaloCluster-based hit. If it has no CaloCluster, this has no content
std::vector<KalSegment>::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
Expand Down
11 changes: 5 additions & 6 deletions RecoDataProducts/inc/KalSegment.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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); }
Expand Down
98 changes: 21 additions & 77 deletions RecoDataProducts/src/KalSeed.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ namespace mu2e {
return retval;
}

std::vector<KalSegment>::const_iterator KalSeed::nearestSeg(double time) const {
std::vector<KalSegment>::const_iterator KalSeed::nearestSegment(double time) const {
auto retval = segments().end();
float dmin(std::numeric_limits<float>::max());
for(auto ikseg = segments().begin(); ikseg !=segments().end(); ikseg++) {
Expand Down Expand Up @@ -98,20 +98,6 @@ namespace mu2e {
}
return retval;
}

std::vector<KalSegment>::const_iterator KalSeed::nearestSegment(float time) const {
auto retval = segments().end();
double tmin(std::numeric_limits<float>::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<KalSegment>::const_iterator KalSeed::nearestSegment(const XYZVectorF& pos) const {
auto retval = segments().end();
double dmin(std::numeric_limits<float>::max());
Expand All @@ -125,74 +111,32 @@ namespace mu2e {
return retval;
}

double KalSeed::t0Val() const {

std::vector<KalSegment>::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);
}
}

32 changes: 20 additions & 12 deletions RecoDataProducts/src/KalSegment.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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_)));
}
}
2 changes: 1 addition & 1 deletion TEveEventDisplay/src/TEveMu2eCustomHelix.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down

0 comments on commit bc9a12d

Please sign in to comment.