diff --git a/Detector/ElementXing.hh b/Detector/ElementXing.hh index 7cc8c906..87a29b08 100644 --- a/Detector/ElementXing.hh +++ b/Detector/ElementXing.hh @@ -33,8 +33,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 c03cec9c..5c5d4b16 100644 --- a/Examples/StrawXing.hh +++ b/Examples/StrawXing.hh @@ -57,6 +57,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/Domain.hh b/Fit/Domain.hh index d9a24adf..9fde2100 100644 --- a/Fit/Domain.hh +++ b/Fit/Domain.hh @@ -11,8 +11,8 @@ 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) {} // clone op for reinstantiation Domain(Domain const&); std::shared_ptr< Domain > clone(CloneContext&) const; @@ -22,21 +22,17 @@ namespace KinKal { 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 }; // 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/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 4f88b5d1..e3500983 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -87,8 +87,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; @@ -98,8 +98,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 ); + // reconstitute from a previous track + Track(Config const& config, BFieldMap const& bfield, PKTRAJPTR& fittraj, HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains); // copy constructor Track(const Track& rhs, CloneContext& context); // extend an existing track with either new configuration, new hits, and/or new material xings @@ -113,8 +115,7 @@ 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_; } + 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_; } @@ -125,77 +126,90 @@ 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, 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); + void fit(HITCOL& hits, EXINGCOL& exings, DOMAINCOL& domains, PKTRAJPTR& fittraj); 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 + 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 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 + 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 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 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 detrange = detectorRange(hits,exings,true); + // 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()); + convertSeed(seedtraj,detrange,domains); + // convert all the primary info to fit effects + createEffects(hits, exings, domains); + // now fit + fit(); + } + + 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 + // trim the trajectory to this range + detrange.combine((*domains.begin())->range()); + detrange.combine((*domains.rbegin())->range()); + fittraj_->setRange(detrange,true); createEffects(hits,exings,domains); - // now fit the track fit(); } @@ -203,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); @@ -229,18 +242,18 @@ namespace KinKal { // 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){ - TimeRange newrange = getRange(hits,exings); + TimeRange newrange = detectorRange(hits,exings); exrange.combine(newrange); } DOMAINCOL domains; @@ -249,21 +262,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()); } } @@ -273,7 +286,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); @@ -301,31 +314,36 @@ namespace KinKal { } } } - // create new traj - auto newtraj = std::make_unique(); + 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); + } + ++olditer; } } - // 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); } @@ -338,34 +356,41 @@ 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"); - fittraj_ = std::make_unique(); + // convert the seedtraj to a PKTRAJ + PKTRAJ 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); + fittraj_ = std::make_unique(firstpiece); } } @@ -451,8 +476,8 @@ namespace KinKal { return; } if(config().bfcorr_){ - // update the limits if new DW effects were added - if(extendDomains(fitrange,config().tol_))setBounds(fwdbnds,revbnds); + // update the limits if new DW effects were added + if(extendDomains(fitrange))setBounds(fwdbnds,revbnds); } FitStateArray states; initFitState(states, fitrange, config().dwt_/miconfig.varianceScale()); @@ -514,7 +539,7 @@ namespace KinKal { } // set the range front.setRange(fitrange); - return std::make_unique(front); + return std::make_unique(front); } // initialize states used before iteration @@ -531,7 +556,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 @@ -615,13 +640,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()); @@ -632,9 +657,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(); } @@ -644,9 +669,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(); } @@ -660,7 +685,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); @@ -705,46 +730,52 @@ 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(PKTRAJ 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)))); + // 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::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((!active) || hit->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((!active) || exing->active() ){ + tmin = std::min(tmin,exing->time()); + tmax = std::max(tmax,exing->time()); + } } return TimeRange(tmin,tmax); } 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 @@ -757,7 +788,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(); } @@ -804,7 +835,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 @@ -836,4 +867,5 @@ namespace KinKal { return TimeRange(tmin,tmax); } } + #endif 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/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 f1c5ecc0..b9e25ccc 100644 --- a/Trajectory/PiecewiseTrajectory.hh +++ b/Trajectory/PiecewiseTrajectory.hh @@ -10,6 +10,7 @@ #include "KinKal/General/MomBasis.hh" #include "KinKal/General/TimeRange.hh" #include +#include #include #include #include @@ -20,6 +21,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); } @@ -55,6 +58,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 @@ -68,13 +73,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("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())); + pieces_.front()->range().restrict(trange); + pieces_.back()->range().restrict(trange); } template PiecewiseTrajectory::PiecewiseTrajectory(KTRAJ const& piece) : pieces_(1,std::make_shared(piece)) @@ -89,13 +98,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 +113,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 +131,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 +139,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 +148,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 +168,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"); } } } @@ -253,6 +262,35 @@ namespace KinKal { return ost; } + template void PiecewiseTrajectory::pieceRange(TimeRange const& range, + std::deque>::const_iterator& first, + std::deque>::const_iterator& last ) const { + 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 = 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(); + } + } + // clone op for reinstantiation template PiecewiseTrajectory::PiecewiseTrajectory(PiecewiseTrajectory const& rhs){