From d7e11680b347fe15ab05f7efe6a76abb48b8706b Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Mon, 1 Sep 2025 20:43:25 -0700 Subject: [PATCH 1/8] Simplifications to smooth regrowing seeds --- Fit/Status.hh | 2 +- Fit/Track.hh | 197 ++++++++++++++++-------------- Tests/FitTest.hh | 2 +- Trajectory/PiecewiseTrajectory.hh | 26 ++++ 4 files changed, 130 insertions(+), 97 deletions(-) diff --git a/Fit/Status.hh b/Fit/Status.hh index 321d798d..036e907a 100644 --- a/Fit/Status.hh +++ b/Fit/Status.hh @@ -17,7 +17,7 @@ namespace KinKal { std::string comment_; // further information about the status bool usable() const { return status_ > unfit && status_ < lowNDOF; } bool needsFit() const { return status_ == unfit || status_ == unconverged; } - Status(unsigned miter,unsigned iter=0) : miter_(miter), iter_(iter), status_(unfit), chisq_(NParams()){} + Status(unsigned miter,unsigned iter=0,status stat=unfit,const char* comment="") : miter_(miter), iter_(iter), status_(stat), chisq_(NParams()),comment_(comment){} static std::string statusName(status stat); }; std::ostream& operator <<(std::ostream& os, Status const& fitstatus ); diff --git a/Fit/Track.hh b/Fit/Track.hh index 159bd81b..e41edca0 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -97,8 +97,10 @@ namespace KinKal { using DOMAINCOL = std::set; using CONFIGCOL = std::vector; using FitStateArray = std::array; - // construct from a set of hits and passive material crossings - Track(Config const& config, BFieldMap const& bfield, PTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings ); + // construct from a set of hits and passive material crossings. This will compute the magnetic domains implicitly from the configuration + Track(Config const& config, BFieldMap const& bfield, KTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings ); + // construct from all effects plus a fit trajectory. This is a reconstitution construction + Track(Config const& config, BFieldMap const& bfield, PTRAJ const& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains ); // extend an existing track with either new configuration, new hits, and/or new material xings void extend(Config const& config, HITCOL& hits, EXINGCOL& exings ); // extrapolate the fit through the magnetic field with the given config until the given predicate is satisfied. This function requires @@ -110,7 +112,6 @@ namespace KinKal { // accessors std::vector const& history() const { return history_; } Status const& fitStatus() const { return history_.back(); } // most recent status - PTRAJ const& seedTraj() const { return seedtraj_; } PTRAJ const& fitTraj() const { return *fittraj_; } KKEFFCOL const& effects() const { return effects_; } Config const& config() const { return config_.back(); } @@ -123,12 +124,14 @@ namespace KinKal { TimeRange activeRange() const; // time range of active hits void extendTraj(TimeRange const& newrange); protected: - Track(Config const& cfg, BFieldMap const& bfield, PTRAJ const& seedtraj ); - void fit(HITCOL& hits, EXINGCOL& exings ); + Track(Config const& cfg, BFieldMap const& bfield); + void fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj); private: - // helper functions - TimeRange getRange(HITCOL& hits, EXINGCOL& exings) const; + void createEffects( HITCOL& hits, EXINGCOL& exings, DOMAINCOL const& domains); + void convertSeed(KTRAJ const& seedtraj,TimeRange const& refrange, DOMAINCOL& domains); void fit(); // process the effects and create the trajectory. This executes the current schedule + static TimeRange getRange(HITCOL& hits, EXINGCOL& exings); + bool createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const; bool setBounds(KKEFFFWDBND& fwdbnds,KKEFFREVBND& revbnds); // set the bounds. Returns false if the bounds are empty bool extendDomains(TimeRange const& fitrange,double tol); // extend domains if the fit range changes. Return value says if domains were added void updateDomains(PTRAJ const& ptraj); // Update domains between iterations @@ -137,75 +140,65 @@ namespace KinKal { void initFitState(FitStateArray& states, TimeRange const& fitrange, double dwt=1.0); PTRAJPTR initTraj(FitState& state, TimeRange const& fitrange); bool canIterate() const; - void createEffects( HITCOL& hits, EXINGCOL& exings, DOMAINCOL const& domains); - void createTraj(PTRAJ const& seedtraj,TimeRange const& refrange, DOMAINCOL const& domains); void replaceDomains(DOMAINCOL const& domains); void extendTraj(DOMAINCOL const& domains); void processEnds(); // add a single domain within the tolerance and extend the fit in the specified direction. void addDomain(Domain const& domain,TimeDir const& tdir,bool exact=false); auto& status() { return history_.back(); } // most recent status - // divide the range into magnetic 'domains' within which the BField can be considered constant (parameter change is within tolerance) - bool createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains, double tol) const; + // divide a trajectory into magnetic 'domains' within which the BField can be considered constant (parameter change is within tolerance) // payload CONFIGCOL config_; // configuration BFieldMap const& bfield_; // magnetic field map std::vector history_; // fit status history; records the current iteration - PTRAJ seedtraj_; // seed for the fit PTRAJPTR fittraj_; // result of the current fit KKEFFCOL effects_; // effects used in this fit, sorted by time HITCOL hits_; // hits used in this fit EXINGCOL exings_; // material xings used in this fit DOMAINCOL domains_; // DomainWall domains used in this fit, contiguous and sorted by time - }; - // sub-class constructor, based just on the seed. It requires added hits to create a functional track - template Track::Track(Config const& cfg, BFieldMap const& bfield, PTRAJ const& seedtraj ) : - bfield_(bfield), seedtraj_(seedtraj) + // minimal constructor for subclasses. The resulting object has no fit. + template Track::Track(Config const& cfg, BFieldMap const& bfield) : + config_{cfg}, bfield_(bfield) { - config_.push_back(cfg); if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule"); + history_.push_back(Status(0,0,Status::unfit, "Construction")); } - // construct from configuration, reference (seed) fit, hits,and materials specific to this fit. - template Track::Track(Config const& cfg, BFieldMap const& bfield, PTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings) : Track(cfg,bfield,seedtraj) { - fit(hits,exings); + // construct from configuration, reference (seed) fit, hits,and materials specific to this fit. This will compute the domains according to the configuration before fitting. + template Track::Track(Config const& cfg, BFieldMap const& bfield, KTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings) : Track(cfg,bfield) { + fit(hits,exings, seedtraj); } - template void Track::fit(HITCOL& hits, EXINGCOL& exings) { - // set the seed time based on the min and max time from the inputs - TimeRange refrange = getRange(hits,exings); - seedtraj_.setRange(refrange); - // if correcting for magnetic field effects, define the domains + + template void Track::fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj) { + auto herange = this->getRange(hits,exings); + // convert the seed traj to a piecewaise traj. This creates the domains DOMAINCOL domains; - if(config().bfcorr_ ) { - bool dok = createDomains(seedtraj_, refrange, domains, config().tol_); - if(!dok){ - // failed initialization - history_.push_back(Status(0)); - status().status_ = Status::outsidemap; - status().comment_ = std::string("Initialization error"); - return; - } - } - // Create the initial reference trajectory from the seed trajectory - createTraj(seedtraj_,refrange,domains); - // create all the other effects -// effects_.reserve(hits.size()+exings.size()+domains.size()); - createEffects(hits,exings,domains); - // now fit the track - fit(); + convertSeed(seedtraj,herange,domains); + // convert all the primary info to fit effects + this->createEffects(hits, exings, domains); + // now fit + this->fit(); } + template Track::Track(Config const& cfg, BFieldMap const& bfield, PTRAJ const& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains) : + Track(cfg,bfield) { + fittraj_ = std::make_unique(fittraj); + // check consistency of domains and initial fittraj TODO + createEffects(hits,exings,domains); + fit(); + } + // extend an existing track template void Track::extend(Config const& cfg, HITCOL& hits, EXINGCOL& exings) { + // require the existing fit to be usable + if(!fitStatus().usable())return; // remember the previous config auto const& oldconfig = config_.back(); // update the configuration config_.push_back(cfg); // configuation check if(config().schedule().size() ==0)throw std::invalid_argument("Invalid configuration: no schedule"); - // require the existing fit to be usable - if(!fitStatus().usable())throw std::invalid_argument("Cannot extend unusable fit"); // find the range of the added information, and extend as needed TimeRange exrange = fittraj_->range(); if(hits.size() >0 || exings.size() > 0){ @@ -218,21 +211,21 @@ namespace KinKal { // if the previous configuration didn't have domains, then create them for the full reference range if(!oldconfig.bfcorr_ || oldconfig.tol_ != config().tol_){ // create domains for the whole range - dok &= createDomains(*fittraj_, exrange, domains, config().tol_); - // replace the domains. This also replace the trajectory, as that must reference the new domains + dok &= createDomains(*fittraj_,exrange, domains); + // replace previous domains with these. This replaces the trajectory and bfield-related effects if(dok)replaceDomains(domains); } else { // create domains just for the extensions TimeRange exlow(exrange.begin(),fittraj_->range().begin()); if(exlow.range()>0.0) { DOMAINCOL lowdomains; - dok &= createDomains(*fittraj_, exlow, lowdomains, config().tol_); + dok &= createDomains(*fittraj_,exlow, lowdomains); if(dok)domains.insert(lowdomains.begin(),lowdomains.end()); } TimeRange exhigh(fittraj_->range().end(),exrange.end()); if(exhigh.range()>0.0){ DOMAINCOL highdomains; - dok &= createDomains(*fittraj_, exhigh, highdomains, config().tol_); + dok &= createDomains(*fittraj_,exhigh, highdomains); if(dok)domains.insert(highdomains.begin(),highdomains.end()); } } @@ -242,7 +235,7 @@ namespace KinKal { history_.push_back(Status(0)); status().status_ = Status::outsidemap; status().comment_ = std::string("Extension error"); - return; + return; } // Extend the traj and create the effects for the new info and the new domains extendTraj(domains); @@ -270,31 +263,35 @@ namespace KinKal { } } } - // create new traj auto newtraj = std::make_unique(); // loop over domains for(auto const& domain : domains) { - // loop until we're either out of this domain or the piece is out of this domain - double dtime = domain->begin(); - while(dtime < domain->end()){ - // find the nearest piece of the traj - static double epsilon(1e-10); - auto index = fittraj_->nearestIndex(dtime+epsilon); // make sure step to the next segment at boundaries - auto const& oldpiece = *fittraj_->pieces()[index]; - // create a new piece - KTRAJ newpiece(oldpiece,domain->bnom(),dtime); + // find the range of existing ptraj pieces that overlaps with this domain's range + using KTRAJPTR = std::shared_ptr; + using DKTRAJ = std::deque; + using DKTRAJCITER = DKTRAJ::const_iterator; + DKTRAJCITER first,last; + fittraj_->pieceRange(domain->range(),first,last); + // loop over these pieces + auto olditer = first; + while(olditer != last){ + auto const& oldpiece = **olditer; + // copy this piece, translating bnom to this domain's field + KTRAJ newpiece(oldpiece,domain->bnom(),domain->range().mid()); // set the range for this piece, making sure it is non-zero - double endtime = (index < fittraj_->pieces().size()-1) ? std::max(std::min(domain->end(),oldpiece.range().end()),dtime+epsilon) : domain->end(); - newpiece.range() = TimeRange(dtime,endtime); - newtraj->append(newpiece); - dtime = newpiece.range().end(); + double tstart = std::max(domain->begin(), oldpiece.range().begin()); + double tend = std::min(domain->end(),oldpiece.range().end()); + if(tstart < tend){ + newpiece.range() = TimeRange(tstart,tend); + newtraj->append(newpiece); + } } } - // update all effects to reference this trajectory + // switch over any existing effects to reference this traj (could be none) for (auto& eff : effects_) { eff->updateReference(*newtraj); } - // swap + // swap out the fit fittraj_.swap(newtraj); } @@ -307,31 +304,38 @@ namespace KinKal { } template void Track::extendTraj(TimeRange const& newrange){ - TimeRange temprange(std::min(fittraj_->range().begin(),newrange.begin()), - std::max(fittraj_->range().end(),newrange.end())); - fittraj_->setRange(temprange); + TimeRange temprange(std::min(fittraj_->range().begin(),newrange.begin()), + std::max(fittraj_->range().end(),newrange.end())); + fittraj_->setRange(temprange); } - template void Track::createTraj(PTRAJ const& seedtraj , TimeRange const& range, DOMAINCOL const& domains ) { + template void Track::convertSeed(KTRAJ const& seedtraj,TimeRange const& range, DOMAINCOL& domains) { // if we're making local DomainWall corrections, divide the trajectory into domain pieces. Each will have equivalent parameters, but relative // to the local field if(config().bfcorr_ ) { - if(fittraj_)throw std::invalid_argument("Initial reference trajectory must be empty"); - if(domains.size() == 0)throw std::invalid_argument("Empty domain collection"); + // convert the seedtraj to a PTRAJ + PTRAJ straj(seedtraj); + bool dok = createDomains(straj, range, domains); + if(!dok){ + // failed initialization + history_.emplace_back(0,0,Status::outsidemap, "Domain initialization error"); + return; + } + // create the initial fittraj from the seed. Each piece will have 'the same' physical trajectory as the seed, but reference the local domain BField fittraj_ = std::make_unique(); for(auto const& domain : domains) { // Set the DomainWall to the start of this domain - auto bf = bfield_.fieldVect(seedtraj.position3(domain->begin())); - KTRAJ newpiece(seedtraj.nearestPiece(domain->begin()),bf,domain->begin()); + auto bf = bfield_.fieldVect(seedtraj.position3(domain->mid())); + KTRAJ newpiece(seedtraj,bf,domain->mid()); newpiece.range() = domain->range(); fittraj_->append(newpiece); } } else { - // use the middle of the range as the nominal DomainWall for this fit: + // use the middle of the range as the nominal BField for this fit: double tref = range.mid(); VEC3 bf = bfield_.fieldVect(seedtraj.position3(tref)); // create the first piece. Note this constructor adjusts the parameters according to the local field - KTRAJ firstpiece(seedtraj.nearestPiece(tref),bf,tref); + KTRAJ firstpiece(seedtraj,bf,tref); firstpiece.range() = range; // create the piecewise trajectory from this fittraj_ = std::make_unique(firstpiece); @@ -420,7 +424,7 @@ namespace KinKal { return; } if(config().bfcorr_){ - // update the limits if new DW effects were added + // update the limits if new DW effects were added if(extendDomains(fitrange,config().tol_))setBounds(fwdbnds,revbnds); } FitStateArray states; @@ -674,30 +678,32 @@ namespace KinKal { } } // divide a trajectory into magnetic 'domains' used to apply the DomainWall corrections - template bool Track::createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains, double tol) const { + template bool Track::createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const { bool retval(true); - auto const& ktraj = ptraj.nearestPiece(range.begin()); - // catch exceptions if the fit extends beyond the range of the field map - try { - double trange = bfield_.rangeInTolerance(ktraj,range.begin(),tol); - // define 1st domain to have the 1st effect in the middle. This avoids effects having exactly the same time - double tstart = range.begin() - 0.5*trange; - do { - // see how far we can go on the current traj before the DomainWall change causes the momentum estimate to go out of tolerance - // note this assumes the trajectory is accurate (geometric extrapolation only) - auto const& ktraj = ptraj.nearestPiece(tstart); - trange = bfield_.rangeInTolerance(ktraj,tstart,tol); - domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)),tol)); - // start the next domain at the end of this one - tstart += trange; - } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect - } catch (std::exception const& error) { - retval = false; + if(config().bfcorr_ ) { + auto const& ktraj = ptraj.nearestPiece(range.begin()); + // catch exceptions if the fit extends beyond the range of the field map + try { + double trange = bfield_.rangeInTolerance(ktraj,range.begin(),config().tol_); + // define 1st domain to have the 1st effect in the middle. This avoids effects having exactly the same time + double tstart = range.begin() - 0.5*trange; + do { + // see how far we can go on the current traj before the DomainWall change causes the momentum estimate to go out of tolerance + // note this assumes the trajectory is accurate (geometric extrapolation only) + auto const& ktraj = ptraj.nearestPiece(tstart); + trange = bfield_.rangeInTolerance(ktraj,tstart,config().tol_); + domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)),config().tol_)); + // start the next domain at the end of this one + tstart += trange; + } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect + } catch (std::exception const& error) { + retval = false; + } } return retval; } - template TimeRange Track::getRange(HITCOL& hits, EXINGCOL& exings) const { + template TimeRange Track::getRange(HITCOL& hits, EXINGCOL& exings) { double tmin = std::numeric_limits::max(); double tmax = -std::numeric_limits::max(); // can't assume effects are sorted @@ -805,4 +811,5 @@ namespace KinKal { return TimeRange(tmin,tmax); } } + #endif diff --git a/Tests/FitTest.hh b/Tests/FitTest.hh index 809e5b59..f2762718 100644 --- a/Tests/FitTest.hh +++ b/Tests/FitTest.hh @@ -464,7 +464,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { } } // create and fit the track - KKTRK kktrk(config,*BF,seedtraj,thits,dxings); + KKTRK kktrk(config,*BF,straj,thits,dxings); if(extend && kktrk.fitStatus().usable())kktrk.extend(exconfig,exthits, exdxings); if(!printbad)kktrk.print(cout,detail); TFile fitfile((KTRAJ::trajName() + string("FitTest") + tfname + string(".root")).c_str(),"RECREATE"); diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index 2806accd..234f32d0 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -9,6 +9,7 @@ #include "KinKal/General/MomBasis.hh" #include "KinKal/General/TimeRange.hh" #include +#include #include #include #include @@ -19,6 +20,8 @@ namespace KinKal { public: using KTRAJPTR = std::shared_ptr; using DKTRAJ = std::deque; + using DKTRAJITER = DKTRAJ::iterator; + using DKTRAJCITER = DKTRAJ::const_iterator; // forward calls to the pieces VEC3 position3(double time) const { return nearestPiece(time).position3(time); } VEC3 velocity(double time) const { return nearestPiece(time).velocity(time); } @@ -50,6 +53,8 @@ namespace KinKal { KTRAJ& back() { return *pieces_.back(); } KTRAJPTR const& frontPtr() const { return pieces_.front(); } KTRAJPTR const& backPtr() const { return pieces_.back(); } + void pieceRange(TimeRange const& range, DKTRAJCITER& first, DKTRAJCITER& last ) const; + void pieceRange(TimeRange const& range,DKTRAJITER& first, DKTRAJITER& last ); size_t nearestIndex(double time) const; DKTRAJ const& pieces() const { return pieces_; } // test for spatial gaps @@ -248,6 +253,27 @@ namespace KinKal { return ost; } + template void PiecewiseTrajectory::pieceRange(TimeRange const& range, + std::deque>::const_iterator& first, + std::deque>::const_iterator& last ) const { + // find the first and last pieces which overlap with the range. They can be the same piece. + first = pieces_.cbegin(); + while(first != pieces_.cend() && !((*first)->range().overlaps(range))) ++first; + auto rlast = pieces_.crbegin(); + while(rlast != pieces_.crend() && !((*rlast)->range().overlaps(range))) ++rlast; + last = (rlast+1).base(); // convert back to forwards-iterator + } + + template void PiecewiseTrajectory::pieceRange(TimeRange const& range, + std::deque>::iterator& first, + std::deque>::iterator& last) { + first = pieces_.begin(); + while(first != pieces_.end() && !((*first)->range().overlaps(range))) ++first; + auto rlast = pieces_.rbegin(); + while(rlast != pieces_.rend() && !((*rlast)->range().overlaps(range))) ++rlast; + last= (rlast+1).base(); + } + } #endif From 859da66da47e57607965d55fd8acba4c1f472f4f Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Tue, 2 Sep 2025 22:15:23 -0700 Subject: [PATCH 2/8] Bug fixes after restructuring --- Fit/Status.hh | 2 +- Fit/Track.hh | 1 + Trajectory/PiecewiseTrajectory.hh | 46 ++++++++++++++++++------------- 3 files changed, 29 insertions(+), 20 deletions(-) diff --git a/Fit/Status.hh b/Fit/Status.hh index 036e907a..c81ab76f 100644 --- a/Fit/Status.hh +++ b/Fit/Status.hh @@ -9,7 +9,7 @@ namespace KinKal { // struct to define fit status struct Status { - enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status + enum status {unfit=-1,incomplete,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status unsigned miter_; // meta-iteration number; unsigned iter_; // iteration number; status status_; // current status diff --git a/Fit/Track.hh b/Fit/Track.hh index e41edca0..cd0e0ef2 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -285,6 +285,7 @@ namespace KinKal { newpiece.range() = TimeRange(tstart,tend); newtraj->append(newpiece); } + ++olditer; } } // switch over any existing effects to reference this traj (could be none) diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index 234f32d0..826f0869 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -71,7 +71,7 @@ namespace KinKal { while(pieces_.size() > 1 && trange.begin() > pieces_.front()->range().end() ) pieces_.pop_front(); while(pieces_.size() > 1 && trange.end() < pieces_.back()->range().begin() ) pieces_.pop_back(); } else if(trange.begin() > pieces_.front()->range().end() || trange.end() < pieces_.back()->range().begin()) - throw std::invalid_argument("Invalid Range"); + throw std::invalid_argument("PiecewiseTrajectory::setRange; Invalid Range"); // update piece range pieces_.front()->setRange(TimeRange(trange.begin(),pieces_.front()->range().end())); pieces_.back()->setRange(TimeRange(pieces_.back()->range().begin(),trange.end())); @@ -89,13 +89,13 @@ namespace KinKal { prepend(newpiece,allowremove); break; default: - throw std::invalid_argument("Invalid direction"); + throw std::invalid_argument("PiecewiseTrajectory::add; Invalid direction"); } } template void PiecewiseTrajectory::prepend(KTRAJ const& newpiece, bool allowremove) { // new piece can't have null range - if(newpiece.range().null())throw std::invalid_argument("Can't prepend null range traj"); + if(newpiece.range().null())throw std::invalid_argument("PiecewiseTrajectory::prepend; Can't prepend null range traj"); if(pieces_.empty()){ pieces_.push_back(std::make_shared(newpiece)); } else { @@ -104,7 +104,7 @@ namespace KinKal { if(allowremove) *this = PiecewiseTrajectory(newpiece); else - throw std::invalid_argument("range overlap"); + throw std::invalid_argument("PiecewiseTrajector::prepend; range overlap"); } else { // find the piece that needs to be modified size_t ipiece = nearestIndex(newpiece.range().end()); @@ -122,7 +122,7 @@ namespace KinKal { pieces_.push_front(std::make_shared(newpiece)); pieces_.front()->range() = TimeRange(tmin,pieces_.front()->range().end()); } else { - throw std::invalid_argument("range error"); + throw std::invalid_argument("PiecewiseTrajectory::prepend; range error"); } } } @@ -130,7 +130,7 @@ namespace KinKal { template void PiecewiseTrajectory::append(KTRAJ const& newpiece, bool allowremove) { // new piece can't have null range - if(newpiece.range().null())throw std::invalid_argument("Can't append null range traj"); + if(newpiece.range().null())throw std::invalid_argument("PiecewiseTrajectory::append; Can't append null range traj"); if(pieces_.empty()){ pieces_.push_back(std::make_shared(newpiece)); } else { @@ -139,7 +139,7 @@ namespace KinKal { if(allowremove) *this = PiecewiseTrajectory(newpiece); else - throw std::invalid_argument("range overlap"); + throw std::invalid_argument("PiecewiseTrajectory::append; range overlap"); } else { // find the piece that needs to be modified size_t ipiece = nearestIndex(newpiece.range().begin()); @@ -159,7 +159,7 @@ namespace KinKal { pieces_.push_back(std::make_shared(newpiece)); pieces_.back()->range() = TimeRange(pieces_.back()->range().begin(),tmax); } else { - throw std::invalid_argument("range error"); + throw std::invalid_argument("PiecewiseTrajectory::append; range error"); } } } @@ -256,22 +256,30 @@ namespace KinKal { template void PiecewiseTrajectory::pieceRange(TimeRange const& range, std::deque>::const_iterator& first, std::deque>::const_iterator& last ) const { - // find the first and last pieces which overlap with the range. They can be the same piece. - first = pieces_.cbegin(); - while(first != pieces_.cend() && !((*first)->range().overlaps(range))) ++first; - auto rlast = pieces_.crbegin(); - while(rlast != pieces_.crend() && !((*rlast)->range().overlaps(range))) ++rlast; - last = (rlast+1).base(); // convert back to forwards-iterator + first = last = pieces_.end(); + // check for no overlap + if(this->range().overlaps(range)){ + // find the first and last pieces which overlap with the range. They can be the same piece. + first = pieces_.cbegin(); + while(first != pieces_.cend() && !((*first)->range().overlaps(range))) ++first; + auto rlast = pieces_.crbegin(); + while(rlast != pieces_.crend() && !((*rlast)->range().overlaps(range))) ++rlast; + last = (rlast+1).base(); // convert back to forwards-iterator + } } template void PiecewiseTrajectory::pieceRange(TimeRange const& range, std::deque>::iterator& first, std::deque>::iterator& last) { - first = pieces_.begin(); - while(first != pieces_.end() && !((*first)->range().overlaps(range))) ++first; - auto rlast = pieces_.rbegin(); - while(rlast != pieces_.rend() && !((*rlast)->range().overlaps(range))) ++rlast; - last= (rlast+1).base(); + first = last = pieces_.end(); + // check for no overlap + if(this->range().overlaps(range)){ + first = pieces_.begin(); + while(first != pieces_.end() && !((*first)->range().overlaps(range))) ++first; + auto rlast = pieces_.rbegin(); + while(rlast != pieces_.rend() && !((*rlast)->range().overlaps(range))) ++rlast; + last= (rlast+1).base(); + } } } From d534211c157f7222eb0faa045f0f830182e1b8c3 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Tue, 2 Sep 2025 22:34:40 -0700 Subject: [PATCH 3/8] remove unsed value --- Fit/Status.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Fit/Status.hh b/Fit/Status.hh index c81ab76f..036e907a 100644 --- a/Fit/Status.hh +++ b/Fit/Status.hh @@ -9,7 +9,7 @@ namespace KinKal { // struct to define fit status struct Status { - enum status {unfit=-1,incomplete,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status + enum status {unfit=-1,converged,unconverged,lowNDOF,gapdiverged,paramsdiverged,chisqdiverged,outsidemap,failed}; // fit status unsigned miter_; // meta-iteration number; unsigned iter_; // iteration number; status status_; // current status From 9d39fe3af4a3c17fc2db15e621ee53593ff62154 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Wed, 3 Sep 2025 17:47:13 -0700 Subject: [PATCH 4/8] Simplify domain. Fix regrow constructor --- Fit/Domain.hh | 6 ++-- Fit/Track.hh | 82 ++++++++++++++++++++++++++++----------------------- 2 files changed, 47 insertions(+), 41 deletions(-) diff --git a/Fit/Domain.hh b/Fit/Domain.hh index 1ae624ab..86fa8b47 100644 --- a/Fit/Domain.hh +++ b/Fit/Domain.hh @@ -9,21 +9,19 @@ namespace KinKal { class Domain : public TimeRange { public: - Domain(double lowtime, double range, VEC3 const& bnom, double tol) : range_(lowtime,lowtime+range), bnom_(bnom), tol_(tol) {} - Domain(TimeRange const& range, VEC3 const& bnom, double tol) : range_(range), bnom_(bnom), tol_(tol) {} + Domain(double lowtime, double range, VEC3 const& bnom) : range_(lowtime,lowtime+range), bnom_(bnom) {} + Domain(TimeRange const& range, VEC3 const& bnom) : range_(range), bnom_(bnom) {} bool operator < (Domain const& other) const {return begin() < other.begin(); } auto const& range() const { return range_; } // forward range functions double begin() const { return range_.begin();} double end() const { return range_.end();} double mid() const { return range_.mid();} - double tolerance() const { return tol_; } auto const& bnom() const { return bnom_; } void updateBNom( VEC3 const& bnom) { bnom_ = bnom; }; // used in DomainWall updating private: TimeRange range_; // range of this domain VEC3 bnom_; // nominal BField for this domain - double tol_; // tolerance used to create this domain }; } #endif diff --git a/Fit/Track.hh b/Fit/Track.hh index cd0e0ef2..78ad855c 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -86,8 +86,8 @@ namespace KinKal { using KKMEAS = Measurement; using KKMAT = Material; using KKDW = DomainWall; - using PTRAJ = ParticleTrajectory; - using PTRAJPTR = std::unique_ptr; + using PKTRAJ = ParticleTrajectory; + using PKTRAJPTR = std::unique_ptr; using HIT = Hit; using HITPTR = std::shared_ptr; using HITCOL = std::vector; @@ -99,8 +99,8 @@ namespace KinKal { using FitStateArray = std::array; // construct from a set of hits and passive material crossings. This will compute the magnetic domains implicitly from the configuration Track(Config const& config, BFieldMap const& bfield, KTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings ); - // construct from all effects plus a fit trajectory. This is a reconstitution construction - Track(Config const& config, BFieldMap const& bfield, PTRAJ const& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains ); + // reconstitute from a previous track + Track(Config const& config, BFieldMap const& bfield, PKTRAJPTR& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains); // extend an existing track with either new configuration, new hits, and/or new material xings void extend(Config const& config, HITCOL& hits, EXINGCOL& exings ); // extrapolate the fit through the magnetic field with the given config until the given predicate is satisfied. This function requires @@ -112,7 +112,7 @@ namespace KinKal { // accessors std::vector const& history() const { return history_; } Status const& fitStatus() const { return history_.back(); } // most recent status - PTRAJ const& fitTraj() const { return *fittraj_; } + PKTRAJ const& fitTraj() const { return *fittraj_; } KKEFFCOL const& effects() const { return effects_; } Config const& config() const { return config_.back(); } CONFIGCOL const& configs() const { return config_; } @@ -126,19 +126,20 @@ namespace KinKal { protected: Track(Config const& cfg, BFieldMap const& bfield); void fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj); + void fit(HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains, PKTRAJPTR& fittraj); private: void createEffects( HITCOL& hits, EXINGCOL& exings, DOMAINCOL const& domains); void convertSeed(KTRAJ const& seedtraj,TimeRange const& refrange, DOMAINCOL& domains); void fit(); // process the effects and create the trajectory. This executes the current schedule static TimeRange getRange(HITCOL& hits, EXINGCOL& exings); - bool createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const; + bool createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const; bool setBounds(KKEFFFWDBND& fwdbnds,KKEFFREVBND& revbnds); // set the bounds. Returns false if the bounds are empty - bool extendDomains(TimeRange const& fitrange,double tol); // extend domains if the fit range changes. Return value says if domains were added - void updateDomains(PTRAJ const& ptraj); // Update domains between iterations + bool extendDomains(TimeRange const& fitrange); // extend domains if the fit range changes. Return value says if domains were added + void updateDomains(PKTRAJ const& ptraj); // Update domains between iterations void iterate(MetaIterConfig const& miconfig); - void setStatus(PTRAJPTR& ptrajptr); + void setStatus(PKTRAJPTR& ptrajptr); void initFitState(FitStateArray& states, TimeRange const& fitrange, double dwt=1.0); - PTRAJPTR initTraj(FitState& state, TimeRange const& fitrange); + PKTRAJPTR initTraj(FitState& state, TimeRange const& fitrange); bool canIterate() const; void replaceDomains(DOMAINCOL const& domains); void extendTraj(DOMAINCOL const& domains); @@ -151,7 +152,7 @@ namespace KinKal { CONFIGCOL config_; // configuration BFieldMap const& bfield_; // magnetic field map std::vector history_; // fit status history; records the current iteration - PTRAJPTR fittraj_; // result of the current fit + PKTRAJPTR fittraj_; // result of the current fit KKEFFCOL effects_; // effects used in this fit, sorted by time HITCOL hits_; // hits used in this fit EXINGCOL exings_; // material xings used in this fit @@ -166,10 +167,19 @@ namespace KinKal { } // construct from configuration, reference (seed) fit, hits,and materials specific to this fit. This will compute the domains according to the configuration before fitting. - template Track::Track(Config const& cfg, BFieldMap const& bfield, KTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings) : Track(cfg,bfield) { + // + template Track::Track(Config const& cfg, BFieldMap const& bfield, KTRAJ const& seedtraj, HITCOL& hits, EXINGCOL& exings) + : Track(cfg,bfield) + { fit(hits,exings, seedtraj); } + template Track::Track(Config const& cfg, BFieldMap const& bfield, PKTRAJPTR& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains) + : Track(cfg,bfield) + { + fit(hits,exings,domains,fittraj); + } + template void Track::fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj) { auto herange = this->getRange(hits,exings); // convert the seed traj to a piecewaise traj. This creates the domains @@ -181,13 +191,11 @@ namespace KinKal { this->fit(); } - template Track::Track(Config const& cfg, BFieldMap const& bfield, PTRAJ const& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains) : - Track(cfg,bfield) { - fittraj_ = std::make_unique(fittraj); - // check consistency of domains and initial fittraj TODO - createEffects(hits,exings,domains); - fit(); - } + template void Track::fit(HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains, PKTRAJPTR& fittraj) { + fittraj_ = std::move(fittraj); // steal the underlying object + createEffects(hits,exings,domains); + fit(); + } // extend an existing track template void Track::extend(Config const& cfg, HITCOL& hits, EXINGCOL& exings) { @@ -263,7 +271,7 @@ namespace KinKal { } } } - auto newtraj = std::make_unique(); + auto newtraj = std::make_unique(); // loop over domains for(auto const& domain : domains) { // find the range of existing ptraj pieces that overlaps with this domain's range @@ -314,8 +322,8 @@ namespace KinKal { // if we're making local DomainWall corrections, divide the trajectory into domain pieces. Each will have equivalent parameters, but relative // to the local field if(config().bfcorr_ ) { - // convert the seedtraj to a PTRAJ - PTRAJ straj(seedtraj); + // convert the seedtraj to a PKTRAJ + PKTRAJ straj(seedtraj); bool dok = createDomains(straj, range, domains); if(!dok){ // failed initialization @@ -323,7 +331,7 @@ namespace KinKal { return; } // create the initial fittraj from the seed. Each piece will have 'the same' physical trajectory as the seed, but reference the local domain BField - fittraj_ = std::make_unique(); + fittraj_ = std::make_unique(); for(auto const& domain : domains) { // Set the DomainWall to the start of this domain auto bf = bfield_.fieldVect(seedtraj.position3(domain->mid())); @@ -339,7 +347,7 @@ namespace KinKal { KTRAJ firstpiece(seedtraj,bf,tref); firstpiece.range() = range; // create the piecewise trajectory from this - fittraj_ = std::make_unique(firstpiece); + fittraj_ = std::make_unique(firstpiece); } } @@ -426,7 +434,7 @@ namespace KinKal { } if(config().bfcorr_){ // update the limits if new DW effects were added - if(extendDomains(fitrange,config().tol_))setBounds(fwdbnds,revbnds); + if(extendDomains(fitrange))setBounds(fwdbnds,revbnds); } FitStateArray states; initFitState(states, fitrange, config().dwt_/miconfig.varianceScale()); @@ -488,7 +496,7 @@ namespace KinKal { } // set the range front.setRange(fitrange); - return std::make_unique(front); + return std::make_unique(front); } // initialize states used before iteration @@ -505,7 +513,7 @@ namespace KinKal { } // finalize after iteration - template void Track::setStatus(PTRAJPTR& ptraj) { + template void Track::setStatus(PKTRAJPTR& ptraj) { // compute parameter difference WRT previous iteration. Compare at front and back ends auto const& ffront = ptraj->front(); // test covariance @@ -589,13 +597,13 @@ namespace KinKal { return retval; } - template void Track::updateDomains(PTRAJ const& ptraj) { + template void Track::updateDomains(PKTRAJ const& ptraj) { for(auto& domain : domains_) { domain->updateBNom(bfield_.fieldVect(ptraj.position3(domain->mid()))); } } - template bool Track::extendDomains(TimeRange const& fitrange, double tol) { + template bool Track::extendDomains(TimeRange const& fitrange) { bool retval(false); // then, check if the range has TimeRange drange(domains().begin()->get()->begin(),domains().rbegin()->get()->end()); @@ -606,9 +614,9 @@ namespace KinKal { double time = drange.begin(); while(time > fitrange.begin()){ auto const& ktraj = fittraj_->nearestPiece(time); - double dt = bfield_.rangeInTolerance(ktraj,time,tol); + double dt = bfield_.rangeInTolerance(ktraj,time,config().tol_); TimeRange range(time-dt,time); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),tol); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); addDomain(domain,TimeDir::backwards); time = domain.begin(); } @@ -618,9 +626,9 @@ namespace KinKal { double time = drange.end(); while(time < fitrange.end()){ auto const& ktraj = fittraj_->nearestPiece(time); - double dt = bfield_.rangeInTolerance(ktraj,time,tol); + double dt = bfield_.rangeInTolerance(ktraj,time,config().tol_); TimeRange range(time,time+dt); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),tol); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); addDomain(domain,TimeDir::forwards); time = domain.end(); } @@ -634,7 +642,7 @@ namespace KinKal { effects_.sort(KKEFFComp()); // update domains as needed to cover the end effects TimeRange endrange(effects_.front()->time(),effects_.back()->time()); - extendDomains(endrange,config().tol_); + extendDomains(endrange); KKEFFFWDBND fwdbnds; // bounding sites used for fitting KKEFFREVBND revbnds; setBounds(fwdbnds,revbnds); @@ -679,7 +687,7 @@ namespace KinKal { } } // divide a trajectory into magnetic 'domains' used to apply the DomainWall corrections - template bool Track::createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const { + template bool Track::createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const { bool retval(true); if(config().bfcorr_ ) { auto const& ktraj = ptraj.nearestPiece(range.begin()); @@ -693,7 +701,7 @@ namespace KinKal { // note this assumes the trajectory is accurate (geometric extrapolation only) auto const& ktraj = ptraj.nearestPiece(tstart); trange = bfield_.rangeInTolerance(ktraj,tstart,config().tol_); - domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)),config().tol_)); + domains.emplace(std::make_shared(tstart,trange,bfield_.fieldVect(ktraj.position3(tstart+0.5*trange)))); // start the next domain at the end of this one tstart += trange; } while(tstart < range.end() + 0.5*trange); // ensure the last domain fully covers the last effect @@ -733,7 +741,7 @@ namespace KinKal { auto const& ktraj = fittraj_->nearestPiece(time); double dt = bfield_.rangeInTolerance(ktraj,time,xtest.dpTolerance()); // always positive TimeRange range = tdir == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); - Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),xtest.dpTolerance()); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid()))); addDomain(domain,tdir,true); // use exact transport time = tdir == TimeDir::forwards ? domain.end() : domain.begin(); } From 77311ad32fd1149a9d7e4b696a6c11ef5916a1d6 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 5 Sep 2025 09:47:09 -0700 Subject: [PATCH 5/8] Truncate to remove extrapolation from regrown fit --- Detector/ElementXing.hh | 3 +-- Examples/StrawXing.hh | 1 + Fit/Track.hh | 33 +++++++++++++++++++++++-------- General/TimeRange.hh | 2 +- Trajectory/PiecewiseTrajectory.hh | 12 +++++++---- 5 files changed, 36 insertions(+), 15 deletions(-) diff --git a/Detector/ElementXing.hh b/Detector/ElementXing.hh index 75d5f42f..eff41193 100644 --- a/Detector/ElementXing.hh +++ b/Detector/ElementXing.hh @@ -27,8 +27,7 @@ namespace KinKal { virtual KTRAJ const& referenceTrajectory() const =0; // trajectory WRT which the xing is defined virtual std::vectorconst& matXings() const =0; // Effect of each physical material component of this detector element on this trajectory virtual void print(std::ostream& ost=std::cout,int detail=0) const =0; - // crossings without material are inactive - bool active() const { return matXings().size() > 0; } + virtual bool active() const =0; // momentum change and variance increase associated with crossing this element forwards in time, in spatial basis void momentumChange(SVEC3& dmom, SMAT& dmomvar) const; // parameter change associated with crossing this element forwards in time diff --git a/Examples/StrawXing.hh b/Examples/StrawXing.hh index 6efada09..4a6a744f 100644 --- a/Examples/StrawXing.hh +++ b/Examples/StrawXing.hh @@ -30,6 +30,7 @@ namespace KinKal { KTRAJ const& referenceTrajectory() const override { return tpca_.particleTraj(); } std::vectorconst& matXings() const override { return mxings_; } void print(std::ostream& ost=std::cout,int detail=0) const override; + bool active() const override { return mxings_.size() > 0; } // accessors auto const& closestApproach() const { return tpca_; } auto const& strawMaterial() const { return smat_; } diff --git a/Fit/Track.hh b/Fit/Track.hh index 78ad855c..8f69e3a5 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -123,6 +123,7 @@ namespace KinKal { void print(std::ostream& ost=std::cout,int detail=0) const; TimeRange activeRange() const; // time range of active hits void extendTraj(TimeRange const& newrange); + static TimeRange detectorRange(HITCOL& hits, EXINGCOL& exings,bool active=false); protected: Track(Config const& cfg, BFieldMap const& bfield); void fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj); @@ -131,7 +132,6 @@ namespace KinKal { void createEffects( HITCOL& hits, EXINGCOL& exings, DOMAINCOL const& domains); void convertSeed(KTRAJ const& seedtraj,TimeRange const& refrange, DOMAINCOL& domains); void fit(); // process the effects and create the trajectory. This executes the current schedule - static TimeRange getRange(HITCOL& hits, EXINGCOL& exings); bool createDomains(PKTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains) const; bool setBounds(KKEFFFWDBND& fwdbnds,KKEFFREVBND& revbnds); // set the bounds. Returns false if the bounds are empty bool extendDomains(TimeRange const& fitrange); // extend domains if the fit range changes. Return value says if domains were added @@ -181,7 +181,7 @@ namespace KinKal { } template void Track::fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj) { - auto herange = this->getRange(hits,exings); + auto herange = this->detectorRange(hits,exings); // convert the seed traj to a piecewaise traj. This creates the domains DOMAINCOL domains; convertSeed(seedtraj,herange,domains); @@ -193,6 +193,19 @@ namespace KinKal { template void Track::fit(HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains, PKTRAJPTR& fittraj) { fittraj_ = std::move(fittraj); // steal the underlying object + // truncate the domains and fit trajectory to be within the detector range + auto detrange = detectorRange(hits,exings,true); + auto idom = domains.begin(); + // stop at the 1st domain overlaping the detector range, and erase all elements up to that point + while(idom != domains.end() && !(detrange.overlaps((*idom)->range())))++idom; + if(idom != domains.begin())domains.erase(domains.begin(),--idom);// leave the overlapping piece + auto jdom= domains.rbegin(); + while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom; + domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator + // extend the range to include the first and last domains themselves, and update the traj accordingly + detrange.combine((*domains.begin())->range()); + detrange.combine((*domains.rbegin())->range()); + fittraj_->setRange(detrange,true); createEffects(hits,exings,domains); fit(); } @@ -210,7 +223,7 @@ namespace KinKal { // find the range of the added information, and extend as needed TimeRange exrange = fittraj_->range(); if(hits.size() >0 || exings.size() > 0){ - TimeRange newrange = getRange(hits,exings); + TimeRange newrange = detectorRange(hits,exings); exrange.combine(newrange); } DOMAINCOL domains; @@ -712,17 +725,21 @@ namespace KinKal { return retval; } - template TimeRange Track::getRange(HITCOL& hits, EXINGCOL& exings) { + template TimeRange Track::detectorRange(HITCOL& hits, EXINGCOL& exings,bool active) { double tmin = std::numeric_limits::max(); double tmax = -std::numeric_limits::max(); // can't assume effects are sorted for(auto const& hit : hits){ - tmin = std::min(tmin,hit->time()); - tmax = std::max(tmax,hit->time()); + if(hit->active() || !active){ + tmin = std::min(tmin,hit->time()); + tmax = std::max(tmax,hit->time()); + } } for(auto const& exing : exings){ - tmin = std::min(tmin,exing->time()); - tmax = std::max(tmax,exing->time()); + if(exing->active() || !active){ + tmin = std::min(tmin,exing->time()); + tmax = std::max(tmax,exing->time()); + } } return TimeRange(tmin,tmax); } diff --git a/General/TimeRange.hh b/General/TimeRange.hh index f1606081..00673042 100644 --- a/General/TimeRange.hh +++ b/General/TimeRange.hh @@ -20,7 +20,7 @@ namespace KinKal { bool inRange(double t) const {return t >= begin() && t < end(); } bool null() const { return end() == begin(); } bool overlaps(TimeRange const& other ) const { - return inRange(other.begin()) || inRange(other.end()) || contains(other) || other.contains(*this); } + return inRange(other.begin()) || other.inRange(begin()) || contains(other) || other.contains(*this); } bool contains(TimeRange const& other) const { return (begin() <= other.begin() && end() >= other.end()); } // force time to be in range diff --git a/Trajectory/PiecewiseTrajectory.hh b/Trajectory/PiecewiseTrajectory.hh index 826f0869..14c66654 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -68,13 +68,17 @@ namespace KinKal { template void PiecewiseTrajectory::setRange(TimeRange const& trange, bool trim) { // trim pieces as necessary if(trim){ - while(pieces_.size() > 1 && trange.begin() > pieces_.front()->range().end() ) pieces_.pop_front(); - while(pieces_.size() > 1 && trange.end() < pieces_.back()->range().begin() ) pieces_.pop_back(); + auto ipiece = pieces_.begin(); + while (ipiece != pieces_.end() && !(trange.overlaps((*ipiece)->range())))++ipiece; + if(ipiece != pieces_.begin())pieces_.erase(pieces_.begin(),--ipiece); + auto jpiece=pieces_.rbegin(); + while(jpiece != pieces_.rend() && !(trange.overlaps((*jpiece)->range())))++jpiece; + pieces_.erase(jpiece.base(),pieces_.end()); } else if(trange.begin() > pieces_.front()->range().end() || trange.end() < pieces_.back()->range().begin()) throw std::invalid_argument("PiecewiseTrajectory::setRange; Invalid Range"); // update piece range - pieces_.front()->setRange(TimeRange(trange.begin(),pieces_.front()->range().end())); - pieces_.back()->setRange(TimeRange(pieces_.back()->range().begin(),trange.end())); + pieces_.front()->range().restrict(trange); + pieces_.back()->range().restrict(trange); } template PiecewiseTrajectory::PiecewiseTrajectory(KTRAJ const& piece) : pieces_(1,std::make_shared(piece)) From 380feebdc2166917b03edc2d89879a5ef89c80f0 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 5 Sep 2025 13:05:59 -0700 Subject: [PATCH 6/8] Cleanup --- Fit/Track.hh | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index 8f69e3a5..016bcfd2 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -181,14 +181,14 @@ namespace KinKal { } template void Track::fit(HITCOL& hits, EXINGCOL& exings, KTRAJ const& seedtraj) { - auto herange = this->detectorRange(hits,exings); + auto detrange = detectorRange(hits,exings,true); // convert the seed traj to a piecewaise traj. This creates the domains DOMAINCOL domains; - convertSeed(seedtraj,herange,domains); + convertSeed(seedtraj,detrange,domains); // convert all the primary info to fit effects - this->createEffects(hits, exings, domains); + createEffects(hits, exings, domains); // now fit - this->fit(); + fit(); } template void Track::fit(HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains, PKTRAJPTR& fittraj) { @@ -202,7 +202,7 @@ namespace KinKal { auto jdom= domains.rbegin(); while(jdom != domains.rend() && !(detrange.overlaps((*jdom)->range())))++jdom; domains.erase(jdom.base(),domains.end()); // base points 1 past the reverse iterator - // extend the range to include the first and last domains themselves, and update the traj accordingly + // trim the trajectory to this range detrange.combine((*domains.begin())->range()); detrange.combine((*domains.rbegin())->range()); fittraj_->setRange(detrange,true); @@ -745,7 +745,7 @@ namespace KinKal { } template template bool Track::extrapolate(TimeDir tdir, XTEST const& xtest) { - bool retval = this->fitStatus().usable(); + bool retval = fitStatus().usable(); if(retval){ if(config().bfcorr_){ // test for extrapolation outside the bfield map range @@ -805,7 +805,7 @@ namespace KinKal { } template bool Track::extrapolate(EXINGPTR const& exingptr,TimeDir const& tdir) { - bool retval = this->fitStatus().usable(); + bool retval = fitStatus().usable(); if(retval){ if(tdir == TimeDir::forwards){ // make sure the time is legal From d071e58bdc92a283026467a056dcbb0e8acf7cdf Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 5 Sep 2025 16:06:55 -0700 Subject: [PATCH 7/8] fix merge --- Fit/Domain.hh | 4 +--- Fit/Track.hh | 5 ++--- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/Fit/Domain.hh b/Fit/Domain.hh index 7d22d6f4..9fde2100 100644 --- a/Fit/Domain.hh +++ b/Fit/Domain.hh @@ -32,9 +32,7 @@ namespace KinKal { // clone op for reinstantiation Domain::Domain(Domain const& rhs): range_(rhs.range()), - bnom_(rhs.bnom()), - tol_(rhs.tolerance()){ - /**/ + bnom_(rhs.bnom()){ } std::shared_ptr Domain::clone(CloneContext& context) const{ diff --git a/Fit/Track.hh b/Fit/Track.hh index b7189c3e..315c527c 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -217,10 +217,9 @@ namespace KinKal { template Track::Track(const Track& rhs, CloneContext& context) : config_(rhs.configs()), bfield_(rhs.bfield()), - history_(rhs.history()), - seedtraj_(rhs.seedTraj()) + history_(rhs.history()) { - fittraj_ = std::make_unique(rhs.fitTraj()); + fittraj_ = std::make_unique(rhs.fitTraj()); hits_.reserve(rhs.hits().size()); for (const auto& ptr: rhs.hits()){ auto hit = ptr->clone(context); From 93ea284f9d9b49fa010b36cf8ea8ae3321edfffd Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Fri, 12 Sep 2025 16:19:06 -0700 Subject: [PATCH 8/8] Reorder arguments --- Fit/Track.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Fit/Track.hh b/Fit/Track.hh index 315c527c..e3500983 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -760,13 +760,13 @@ namespace KinKal { double tmax = -std::numeric_limits::max(); // can't assume effects are sorted for(auto const& hit : hits){ - if(hit->active() || !active){ + if((!active) || hit->active()){ tmin = std::min(tmin,hit->time()); tmax = std::max(tmax,hit->time()); } } for(auto const& exing : exings){ - if(exing->active() || !active){ + if((!active) || exing->active() ){ tmin = std::min(tmin,exing->time()); tmax = std::max(tmax,exing->time()); }