diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 994e98bf..afb9b54f 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -16,7 +16,7 @@ jobs: matrix: os: ["ubuntu-latest", "macos-latest"] python-version: ["3.9"] - root-version: ["6.26.10"] + root-version: ["6.30.2"] build-type: ["Debug", "Release"] steps: diff --git a/.gitignore b/.gitignore index e3d77230..6a22363e 100644 --- a/.gitignore +++ b/.gitignore @@ -5,5 +5,5 @@ KinKal/UnitTests/Dict.cc debug/* *.root .DS_Store -.vscode +.vscode/* .*.swp diff --git a/Detector/ElementXing.hh b/Detector/ElementXing.hh index ec2d8844..316cba7a 100644 --- a/Detector/ElementXing.hh +++ b/Detector/ElementXing.hh @@ -7,7 +7,6 @@ #include "KinKal/General/MomBasis.hh" #include "KinKal/Detector/MaterialXing.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/General/TimeDir.hh" #include "KinKal/Detector/Hit.hh" #include "KinKal/Fit/MetaIterConfig.hh" #include @@ -25,7 +24,7 @@ namespace KinKal { virtual ~ElementXing() {} virtual void updateReference(KTRAJPTR const& ktrajptr) = 0; // update the trajectory reference virtual void updateState(MetaIterConfig const& config,bool first) =0; // update the state according to this meta-config - virtual Parameters parameters(TimeDir tdir) const =0; // parameter change induced by this element crossing WRT the reference + virtual Parameters params() const =0; // parameter change induced by this element crossing WRT the reference parameters going forwards in time virtual double time() const=0; // time the particle crosses thie element virtual double transitTime() const=0; // time to cross this element virtual KTRAJ const& referenceTrajectory() const =0; // trajectory WRT which the xing is defined @@ -33,19 +32,18 @@ namespace KinKal { 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; } - // calculate the cumulative material effect from these crossings - void materialEffects(TimeDir tdir, std::array& dmom, std::array& momvar) const; + // calculate the cumulative material effect from all the materials in this element crossing going forwards in time + void materialEffects(std::array& dmom, std::array& momvar) const; // sum radiation fraction double radiationFraction() const; - static double elossFactor(TimeDir const& tdir) { return tdir == TimeDir::forwards ? 1.0 : -1.0; } private: }; - template void ElementXing::materialEffects(TimeDir tdir, std::array& dmom, std::array& momvar) const { + template void ElementXing::materialEffects(std::array& dmom, std::array& momvar) const { // compute the derivative of momentum to energy, at the reference trajectory double mom = referenceTrajectory().momentum(time()); double mass = referenceTrajectory().mass(); - double dmFdE = elossFactor(tdir)*sqrt(mom*mom+mass*mass)/(mom*mom); // dimension of 1/E + double dmFdE = sqrt(mom*mom+mass*mass)/(mom*mom); // dimension of 1/E // loop over crossings for this detector piece for(auto const& mxing : matXings()){ // compute FRACTIONAL momentum change and variance on that in the given direction diff --git a/Detector/StrawXing.hh b/Detector/StrawXing.hh index fa8dc460..0d860868 100644 --- a/Detector/StrawXing.hh +++ b/Detector/StrawXing.hh @@ -7,7 +7,7 @@ #include "KinKal/Detector/ElementXing.hh" #include "KinKal/Detector/StrawMaterial.hh" #include "KinKal/Detector/StrawXingConfig.hh" -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" namespace KinKal { @@ -16,15 +16,15 @@ namespace KinKal { using PTRAJ = ParticleTrajectory; using KTRAJPTR = std::shared_ptr; using EXING = ElementXing; - using PCA = PiecewiseClosestApproach; - using CA = ClosestApproach; + using PCA = PiecewiseClosestApproach; + using CA = ClosestApproach; // construct from PCA and material StrawXing(PCA const& pca, StrawMaterial const& smat); virtual ~StrawXing() {} // ElementXing interface void updateReference(KTRAJPTR const& ktrajptr) override; void updateState(MetaIterConfig const& config,bool first) override; - Parameters parameters(TimeDir tdir) const override; + Parameters params() const override; double time() const override { return tpca_.particleToca() + toff_; } // offset time WRT TOCA to avoid exact overlapp with the wire hit double transitTime() const override; // time to cross this element KTRAJ const& referenceTrajectory() const override { return tpca_.particleTraj(); } @@ -36,7 +36,7 @@ namespace KinKal { auto const& config() const { return sxconfig_; } auto precision() const { return tpca_.precision(); } private: - Line axis_; // straw axis, expressed as a timeline + SensorLine axis_; // straw axis, expressed as a timeline StrawMaterial const& smat_; CA tpca_; // result of most recent TPOCA double toff_; // small time offset @@ -78,7 +78,7 @@ namespace KinKal { if(mxings_.size() > 0){ // compute the parameter effect for forwards time std::array dmom = {0.0,0.0,0.0}, momvar = {0.0,0.0,0.0}; - this->materialEffects(TimeDir::forwards, dmom, momvar); + this->materialEffects(dmom, momvar); // get the parameter derivative WRT momentum DPDV dPdM = referenceTrajectory().dPardM(time()); double mommag = referenceTrajectory().momentum(time()); @@ -101,11 +101,8 @@ namespace KinKal { } } - template Parameters StrawXing::parameters(TimeDir tdir) const { - if(tdir == TimeDir::forwards) - return fparams_; - else - return Parameters(-fparams_.parameters(),fparams_.covariance()); + template Parameters StrawXing::params() const { + return fparams_; } template double StrawXing::transitTime() const { diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index c3d827d2..c4a30bfa 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -5,7 +5,7 @@ ROOT_GENERATE_DICTIONARY(ExamplesDict # headers to include in ROOT dictionary HitInfo.hh - BFieldInfo.hh + DomainWallInfo.hh ParticleTrajectoryInfo.hh MaterialInfo.hh LINKDEF LinkDef.h diff --git a/Examples/BFieldInfo.hh b/Examples/DomainWallInfo.hh similarity index 53% rename from Examples/BFieldInfo.hh rename to Examples/DomainWallInfo.hh index 6f10693c..6dca403b 100644 --- a/Examples/BFieldInfo.hh +++ b/Examples/DomainWallInfo.hh @@ -1,14 +1,14 @@ -#ifndef KinKal_BFieldInfo_hh -#define KinKal_BFieldInfo_hh +#ifndef KinKal_DomainWallInfo_hh +#define KinKal_DomainWallInfo_hh #include namespace KinKal { - struct BFieldInfo { - BFieldInfo(){}; + struct DomainWallInfo { + DomainWallInfo(){}; Int_t active_; Float_t time_, range_; static std::string leafnames() { return std::string("active/i:time/f:range/f"); } }; - typedef std::vector KKBFIV; + typedef std::vector KKBFIV; } #endif diff --git a/Examples/LinkDef.h b/Examples/LinkDef.h index 6f4ac83f..54116fb8 100644 --- a/Examples/LinkDef.h +++ b/Examples/LinkDef.h @@ -2,8 +2,8 @@ #pragma link C++ class KinKal::HitInfo+; #pragma link C++ class vector+; -#pragma link C++ class KinKal::BFieldInfo+; -#pragma link C++ class vector+; +#pragma link C++ class KinKal::DomainWallInfo+; +#pragma link C++ class vector+; #pragma link C++ class KinKal::MaterialInfo+; #pragma link C++ class vector+; #pragma link C++ class KinKal::ParticleTrajectoryInfo+; diff --git a/Examples/ScintHit.hh b/Examples/ScintHit.hh index c4d8e256..10fbaf3b 100644 --- a/Examples/ScintHit.hh +++ b/Examples/ScintHit.hh @@ -4,7 +4,7 @@ // simple hit subclass representing a time measurement using scintillator light from a crystal or plastic scintillator // #include "KinKal/Detector/ResidualHit.hh" -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" #include namespace KinKal { @@ -12,8 +12,8 @@ namespace KinKal { template class ScintHit : public ResidualHit { public: using PTRAJ = ParticleTrajectory; - using PCA = PiecewiseClosestApproach; - using CA = ClosestApproach; + using PCA = PiecewiseClosestApproach; + using CA = ClosestApproach; using RESIDHIT = ResidualHit; using HIT = Hit; using KTRAJPTR = std::shared_ptr; @@ -35,7 +35,7 @@ namespace KinKal { double widthVariance() const { return wvar_; } auto precision() const { return tpca_.precision(); } private: - Line saxis_; // symmetry axis of this sensor + SensorLine saxis_; // symmetry axis of this sensor double tvar_; // variance in the time measurement: assumed independent of propagation distance/time double wvar_; // variance in transverse position of the sensor/measurement in mm. Assumes cylindrical error, could be more general CA tpca_; // reference time and position of closest approach to the axis @@ -54,7 +54,7 @@ namespace KinKal { template void ScintHit::updateReference(KTRAJPTR const& ktrajptr) { // use previous hint, or initialize from the sensor time - CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(saxis_.t0(), saxis_.t0()); + CAHint tphint = tpca_.usable() ? tpca_.hint() : CAHint(saxis_.measurementTime(), saxis_.measurementTime()); tpca_ = CA(ktrajptr,saxis_,tphint,precision()); if(!tpca_.usable())throw std::runtime_error("ScintHit TPOCA failure"); } @@ -64,22 +64,21 @@ namespace KinKal { // early in the fit when t0 has very large errors. // If it is unphysical try to adjust it back using a better hint. auto ppos = tpca_.particlePoca().Vect(); - auto sstart = saxis_.startPosition(); - auto send = saxis_.endPosition(); - double slen = (send-sstart).R(); + auto const& sstart = saxis_.start(); + auto const& send = saxis_.end(); // tolerance should come from the config. Should also test relative to the error. FIXME - double tol = slen*1.0; - double sdist = (ppos - saxis_.position3(saxis_.timeAtMidpoint())).Dot(saxis_.direction()); - if( (ppos-sstart).Dot(saxis_.direction()) < -tol || - (ppos-send).Dot(saxis_.direction()) > tol) { + double tol = saxis_.length()*1.0; + double sdist = (ppos - saxis_.middle()).Dot(saxis_.direction()); + if( (ppos-sstart).Dot(saxis_.direction()) < -tol || (ppos-send).Dot(saxis_.direction()) > tol) { // adjust hint to the middle and try agian double sspeed = tpca_.particleTraj().velocity(tpca_.particleToca()).Dot(saxis_.direction()); - auto tphint = tpca_.hint(); tphint.particleToca_ -= sdist/sspeed; tpca_ = CA(tpca_.particleTrajPtr(),saxis_,tphint,precision()); // should check if this is still unphysical and disable the hit if so FIXME + sdist = (tpca_.particlePoca().Vect() - saxis_.middle()).Dot(saxis_.direction()); } + // residual is just delta-T at CA. // the variance includes the measurement variance and the tranvserse size (which couples to the relative direction) // Might want to do more updating (set activity) based on DOCA in future: TODO diff --git a/Examples/SimpleWireHit.hh b/Examples/SimpleWireHit.hh index bd787a53..23a7f6f2 100644 --- a/Examples/SimpleWireHit.hh +++ b/Examples/SimpleWireHit.hh @@ -7,7 +7,7 @@ #include "KinKal/Examples/DOCAWireHitUpdater.hh" #include "KinKal/Examples/WireHitStructs.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" #include "KinKal/Trajectory/ClosestApproach.hh" #include "KinKal/General/BFieldMap.hh" @@ -18,8 +18,8 @@ namespace KinKal { template class SimpleWireHit : public ResidualHit { public: using HIT = Hit; - using PCA = PiecewiseClosestApproach; - using CA = ClosestApproach; + using PCA = PiecewiseClosestApproach; + using CA = ClosestApproach; using KTRAJPTR = std::shared_ptr; enum Dimension { dresid=0, tresid=1}; // residual dimensions @@ -49,7 +49,7 @@ namespace KinKal { private: BFieldMap const& bfield_; // drift calculation requires the BField for ExB effects WireHitState whstate_; // current state - Line wire_; // local linear approximation to the wire of this hit, encoding all (local) position and time information. + SensorLine wire_; // local linear approximation to the wire of this hit, encoding all (local) position and time information. // the start time is the measurement time, the direction is from // the physical source of the signal (particle) to the measurement recording location (electronics), the direction magnitude // is the effective signal propagation velocity along the wire, and the time range describes the active wire length @@ -129,7 +129,7 @@ namespace KinKal { return rresid_[ires]; } - template ClosestApproach SimpleWireHit::unbiasedClosestApproach() const { + template ClosestApproach SimpleWireHit::unbiasedClosestApproach() const { // compute the unbiased closest approach; this is brute force, but works auto const& ca = this->closestApproach(); auto uparams = HIT::unbiasedParameters(); diff --git a/Fit/BField.hh b/Fit/BField.hh deleted file mode 100644 index 2ec7c24f..00000000 --- a/Fit/BField.hh +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef KinKal_BField_hh -#define KinKal_BField_hh -// -// Effect to correct the fit parameters for the change in BField along a small piece of the trajectory. -// This effect adds no information content or noise (presently), just transports the parameters -// -#include "KinKal/Fit/Effect.hh" -#include "KinKal/General/TimeDir.hh" -#include "KinKal/General/BFieldMap.hh" -#include "KinKal/Fit/Config.hh" -#include -#include -#include -#include - -namespace KinKal { - template class BField : public Effect { - public: - using KKEFF = Effect; - using PTRAJ = ParticleTrajectory; - using KTRAJPTR = std::shared_ptr; - double time() const override { return drange_.mid(); } // apply the correction at the middle of the range - bool active() const override { return bfcorr_; } - void process(FitState& kkdata,TimeDir tdir) override; - void updateState(MetaIterConfig const& miconfig,bool first) override {} - void updateConfig(Config const& config) override { bfcorr_ = config.bfcorr_; } - void updateReference(KTRAJPTR const& ltrajptr) override {} // nothing explicit here - void print(std::ostream& ost=std::cout,int detail=0) const override; - void append(PTRAJ& fit,TimeDir tdir) override; - Chisq chisq(Parameters const& pdata) const override { return Chisq();} - auto const& parameterChange() const { return dpfwd_; } - virtual ~BField(){} - // disallow copy and equivalence - BField(BField const& ) = delete; - BField& operator =(BField const& ) = delete; - // create from the domain range, the effect, and the - BField(Config const& config, BFieldMap const& bfield,TimeRange const& drange) : - bfield_(bfield), drange_(drange), bfcorr_(config.bfcorr_) {} - TimeRange const& range() const { return drange_; } - - private: - BFieldMap const& bfield_; // bfield - TimeRange drange_; // extent of this effect. The middle is at the transition point between 2 bfield domains (domain transition) - DVEC dpfwd_; // aggregate effect in parameter space of BFieldMap change over this domain in the forwards time direction - bool bfcorr_; // apply correction or not - }; - - template void BField::process(FitState& kkdata,TimeDir tdir) { - if(bfcorr_){ - kkdata.append(dpfwd_,tdir); - // rotate the covariance matrix for the change in BField. This requires 2nd derivatives TODO - } - } - - template void BField::append(PTRAJ& ptraj,TimeDir tdir) { - if(bfcorr_){ - double etime = time(); - // make sure the piece is appendable - if((tdir == TimeDir::forwards && ptraj.back().range().begin() > etime) || - (tdir == TimeDir::backwards && ptraj.front().range().end() < etime) ) - throw std::invalid_argument("BField: Can't append piece"); - // assume the next domain has ~about the same range - TimeRange newrange = (tdir == TimeDir::forwards) ? TimeRange(etime,std::max(ptraj.range().end(),drange_.end())) : - TimeRange(std::min(ptraj.range().begin(),drange_.begin()),etime); - // update the parameters according to the change in bnom across this domain - // This corresponds to keeping the physical position and momentum constant, but referring to the BField - // at the end vs the begining of the domain - // Use the 1st order approximation: the exact method tried below doesn't do any better (slightly worse) - VEC3 bend = (tdir == TimeDir::forwards) ? bfield_.fieldVect(ptraj.position3(drange_.end())) : bfield_.fieldVect(ptraj.position3(drange_.begin())); - // update the parameter change due to the BField change. Note this assumes the traj piece - // at the begining of the domain has the same bnom as the BField at that point in space - KTRAJ newpiece = (tdir == TimeDir::forwards) ? ptraj.back() : ptraj.front(); - newpiece.setBNom(etime,bend); - newpiece.range() = newrange; - // extract the parameter change for the next processing BEFORE appending - // This should really be done in updateState, but it's easier here and has no knock-on effects - dpfwd_ = (tdir == TimeDir::forwards) ? newpiece.params().parameters() - ptraj.back().params().parameters() : ptraj.back().params().parameters() - newpiece.params().parameters(); - if( tdir == TimeDir::forwards){ - ptraj.append(newpiece); - } else { - ptraj.prepend(newpiece); - } - // exact calculation (for reference) - // extract the particle state at this transition - // auto pstate = ptraj.back().stateEstimate(etime); - // re-compute the trajectory at the domain end using this state - // KTRAJ newpiece(pstate,bend,newrange); - // set the parameter change for the next processing BEFORE appending - // dpfwd_ = newpiece.params().parameters()-ptraj.back().params().parameters(); - // ptraj.append(newpiece); - } - } - - template void BField::print(std::ostream& ost,int detail) const { - ost << "BField " << static_castconst&>(*this); - ost << " effect " << dpfwd_ << " domain range " << drange_ << std::endl; - } - - template std::ostream& operator <<(std::ostream& ost, BField const& kkmat) { - kkmat.print(ost,0); - return ost; - } - -} -#endif diff --git a/Fit/Config.hh b/Fit/Config.hh index fd19afd8..a73ff62a 100644 --- a/Fit/Config.hh +++ b/Fit/Config.hh @@ -17,25 +17,22 @@ namespace KinKal { struct Config { enum printLevel{none=0,minimal, basic, complete, detailed, extreme}; - using Schedule = std::vector; + using Schedule = std::vector; + explicit Config() {} explicit Config(Schedule const& schedule) : Config() { schedule_ = schedule; } - Config() : maxniter_(10), dwt_(1.0e6), convdchisq_(0.01), divdchisq_(10.0), pdchisq_(1.0e6), divgap_(100.0), - tol_(1.0e-4), minndof_(5), bfcorr_(true), ends_(true), plevel_(none) {} - Schedule& schedule() { return schedule_; } Schedule const& schedule() const { return schedule_; } - // algebraic iteration parameters - unsigned maxniter_; // maximum number of algebraic iterations for this config - double dwt_; // dweighting of initial seed covariance - double convdchisq_; // maximum change in chisquared/dof for convergence - double divdchisq_; // minimum change in chisquared/dof for divergence - double pdchisq_; // maximum allowed parameter change (units of chisqred) WRT previous reference - double divgap_; // maximum average gap of trajectory before calling it diverged (mm) - double tol_; // tolerance on fractional momentum accuracy due to BField domain steps - unsigned minndof_; // minimum number of DOFs to continue fit - bool bfcorr_; // whether to make BFieldMap corrections in the fit - bool ends_; // process the passive effects at each end of the track after schedule completion - printLevel plevel_; // print level + unsigned maxniter_ = 10; // maximum number of algebraic iterations for this config + double dwt_ = 1.0e6; // dweighting of initial seed covariance + double convdchisq_ = 1.0e-2; // maximum change in chisquared/dof for convergence + double divdchisq_ = 1.0e1; // minimum change in chisquared/dof for divergence + double pdchisq_ = 1.0e6; // maximum allowed parameter change (units of chisqred) WRT previous reference + double divgap_ = 1.0e2; // maximum average gap of trajectory before calling it diverged (mm) + double tol_ = 1.0e-4; // tolerance on fractional momentum accuracy due to BField domain steps + unsigned minndof_ = 5; // minimum number of DOFs to continue fit + bool bfcorr_ = true; // whether to make BFieldMap corrections in the fit + bool ends_ = true; // process the passive effects at each end of the track after schedule completion + printLevel plevel_ = none; // print level // schedule of meta-iterations. These will be executed sequentially until completion or failure Schedule schedule_; }; diff --git a/Fit/Domain.hh b/Fit/Domain.hh new file mode 100644 index 00000000..1e0c3d26 --- /dev/null +++ b/Fit/Domain.hh @@ -0,0 +1,29 @@ +#ifndef KinKal_Domain_hh +#define KinKal_Domain_hh +// +// domain used in computing magnetic field corrections: just a time range and the tolerance used to define it +// +#include "KinKal/General/TimeRange.hh" +#include "KinKal/General/Vectors.hh" +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) {} + 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/DomainWall.hh b/Fit/DomainWall.hh new file mode 100644 index 00000000..8d73a51a --- /dev/null +++ b/Fit/DomainWall.hh @@ -0,0 +1,120 @@ +#ifndef KinKal_DomainWall_hh +#define KinKal_DomainWall_hh +// +// Effect describing the change in fit parameters for the change in BField crossing between 2 domains +// This effect adds no information content or noise (presently), just transports the parameters +// +#include "KinKal/Fit/Effect.hh" +#include "KinKal/Fit/Domain.hh" +#include "KinKal/Fit/Config.hh" +#include "KinKal/Fit/FitState.hh" +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/BFieldMap.hh" +#include +#include +#include +#include + +namespace KinKal { + template class DomainWall : public Effect { + public: + using KKEFF = Effect; + using PTRAJ = ParticleTrajectory; + using DOMAINPTR = std::shared_ptr; + double time() const override { return prev_->end(); } + bool active() const override { return true; } // always active + void process(FitState& kkdata,TimeDir tdir) override; + void updateState(MetaIterConfig const& miconfig,bool first) override; + void updateConfig(Config const& config) override {} + void updateReference(PTRAJ const& ptraj) override; + void print(std::ostream& ost=std::cout,int detail=0) const override; + void append(PTRAJ& fit,TimeDir tdir) override; + Chisq chisq(Parameters const& pdata) const override { return Chisq();} // no information added + auto const& parameterChange() const { return dpfwd_; } + virtual ~DomainWall(){} + // disallow copy and equivalence + DomainWall(DomainWall const& ) = delete; + DomainWall& operator =(DomainWall const& ) = delete; + // specific DomainWall interface + // create from the domain and BField + DomainWall(BFieldMap const& bfield,DOMAINPTR const& prevdomain,DOMAINPTR const& nextdomain, PTRAJ const& ptraj); + // previous and next domains + auto const& prevDomain() const { return *prev_; } + auto const& nextDomain() const { return *next_; } + + private: + BFieldMap const& bfield_; // bfield + DOMAINPTR prev_, next_; // pointers to previous and next domains + DVEC dpfwd_; // parameter change across this domain wall in the forwards time direction + Weights prevwt_, nextwt_; // cache of weights + PSMAT dpdpdb_; // forward rotation of covariance matrix going in the forwards direction + + }; + + template DomainWall::DomainWall(BFieldMap const& bfield, + DOMAINPTR const& prevdomain, DOMAINPTR const& nextdomain, PTRAJ const& ptraj) : + bfield_(bfield), prev_(prevdomain), next_(nextdomain) { + updateReference(ptraj); + } + + template void DomainWall::process(FitState& fstate,TimeDir tdir) { + if(tdir == TimeDir::forwards) { + prevwt_ += fstate.wData(); + fstate.append(dpfwd_,tdir); + // fstate.pData().covariance() = ROOT::Math::Similarity(dpdpdb_,fstate.pData().covariance()); Not tested TODO + nextwt_ += fstate.wData(); + } else { + nextwt_ += fstate.wData(); + fstate.append(dpfwd_,tdir); + // fstate.pData().covariance() = ROOT::Math::SimilarityT(dpdpdb_,fstate.pData().covariance()); + prevwt_ += fstate.wData(); + } + } + + template void DomainWall::updateState(MetaIterConfig const& miconfig,bool first) { + // reset the cached weights + prevwt_ = nextwt_ = Weights(); + } + + template void DomainWall::updateReference(PTRAJ const& ptraj) { + // update change in parameters for passage through this DW + auto const& refpiece = ptraj.nearestPiece(time()-1e-5); // disambiguate derivativates + auto db = next_->bnom() - prev_->bnom(); + dpfwd_ = refpiece.dPardB(this->time(),db); + dpdpdb_ = refpiece.dPardPardB(this->time(),db); + } + + template void DomainWall::append(PTRAJ& ptraj,TimeDir tdir) { + double etime = time(); + // make sure the piece is appendable + if((tdir == TimeDir::forwards && ptraj.back().range().begin() > etime) || + (tdir == TimeDir::backwards && ptraj.front().range().end() < etime) ) + throw std::invalid_argument("DomainWall: Can't append piece"); + auto const& oldpiece = (tdir == TimeDir::forwards) ? ptraj.back() : ptraj.front(); + KTRAJ newpiece(oldpiece); + if( tdir == TimeDir::forwards){ + newpiece.range() = TimeRange(next_->begin(),std::max(ptraj.range().end(),next_->end())); + newpiece.params() = Parameters(nextwt_); + newpiece.resetBNom(next_->bnom()); + ptraj.append(newpiece); + } else { + newpiece.range() = TimeRange(std::min(ptraj.range().begin(),prev_->begin()),prev_->end()); + newpiece.params() = Parameters(prevwt_); + newpiece.resetBNom(prev_->bnom()); + ptraj.prepend(newpiece); + } + } + + template void DomainWall::print(std::ostream& ost,int detail) const { + ost << "DomainWall " << static_castconst&>(*this); + ost << " previous domain " << *prev_ << " next domain " << *next_; + ost << " effect " << dpfwd_ << std::endl; + } + + template std::ostream& operator <<(std::ostream& ost, DomainWall const& kkmat) { + kkmat.print(ost,0); + return ost; + } + +} +#endif diff --git a/Fit/Effect.hh b/Fit/Effect.hh index 209125b7..547f2f77 100644 --- a/Fit/Effect.hh +++ b/Fit/Effect.hh @@ -20,7 +20,6 @@ namespace KinKal { public: // type of the data payload used for processing the fit using PTRAJ = ParticleTrajectory; - using KTRAJPTR = std::shared_ptr; Effect() {} virtual ~Effect(){} // Effect interface @@ -35,7 +34,7 @@ namespace KinKal { // add this effect to a trajectory in the given direction virtual void append(PTRAJ& fit,TimeDir tdir) =0; // update the reference trajectory for this effect - virtual void updateReference(KTRAJPTR const& ltraj) =0; + virtual void updateReference(PTRAJ const& ptraj) =0; // chisquared WRT a given local parameter set, assumed uncorrelatedd This is used for convergence testing virtual Chisq chisq(Parameters const& pdata) const = 0; // diagnostic printout diff --git a/Fit/ExtraConfig.hh b/Fit/ExtraConfig.hh new file mode 100644 index 00000000..4f913230 --- /dev/null +++ b/Fit/ExtraConfig.hh @@ -0,0 +1,15 @@ +#ifndef KinKal_ExtraConfig_hh +#define KinKal_ExtraConfig_hh +// +// This class defines the configuration for extrapolating a fit beyond the measurement region +// +#include "KinKal/General/TimeDir.hh" + +namespace KinKal { + struct ExtraConfig { + TimeDir xdir_ = TimeDir::forwards; // time direction to extend + double tol_ = 1.0e-3; // tolerance on fractional momentum accuracy due to BField domain steps + double maxdt_ = 1.0e2; // maximum time to extend the fit + }; +} +#endif diff --git a/Fit/FitState.hh b/Fit/FitState.hh index 70646b23..747ad10a 100644 --- a/Fit/FitState.hh +++ b/Fit/FitState.hh @@ -17,8 +17,8 @@ namespace KinKal { // accessors bool hasParameters() const { return hasParameters_; } bool hasWeights() const { return hasWeights_; } - // add to either parameters or weights. Parameters can have a direction - void append(Parameters const& pdata,TimeDir tdir=TimeDir::forwards) { + // add to either parameters or weights. Parameters can have a direction. Covariance is always additivie + void append(Parameters const& pdata,TimeDir tdir) { if(tdir==TimeDir::forwards) pData().parameters() += pdata.parameters(); else diff --git a/Fit/Material.hh b/Fit/Material.hh index b80144f8..e12b4b42 100644 --- a/Fit/Material.hh +++ b/Fit/Material.hh @@ -17,7 +17,6 @@ namespace KinKal { public: using KKEFF = Effect; using PTRAJ = ParticleTrajectory; - using KTRAJPTR = std::shared_ptr; using EXING = ElementXing; using EXINGPTR = std::shared_ptr; double time() const override { return exing_->time();} @@ -26,34 +25,34 @@ namespace KinKal { void updateState(MetaIterConfig const& miconfig,bool first) override; void updateConfig(Config const& config) override {} void append(PTRAJ& fit,TimeDir tdir) override; - void updateReference(KTRAJPTR const& ltrajptr) override; + void updateReference(PTRAJ const& ptraj) override; Chisq chisq(Parameters const& pdata) const override { return Chisq();} void print(std::ostream& ost=std::cout,int detail=0) const override; virtual ~Material(){} // create from the material and a trajectory Material(EXINGPTR const& dxing, PTRAJ const& ptraj); // accessors - auto const& cache() const { return cache_; } auto const& elementXing() const { return *exing_; } auto const& elementXingPtr() const { return exing_; } auto const& referenceTrajectory() const { return exing_->referenceTrajectory(); } private: EXINGPTR exing_; // element crossing for this effect - Weights cache_; // cache of weight processing in opposite directions, used to build the fit trajectory + Weights nextwt_; // cache of weight forwards of this effect. }; - template Material::Material(EXINGPTR const& dxing, PTRAJ const& ptraj) : exing_(dxing) {} + template Material::Material(EXINGPTR const& dxing, PTRAJ const& ptraj) : exing_(dxing) { + updateReference(ptraj); + } template void Material::process(FitState& kkdata,TimeDir tdir) { + // The assymetry in the cache processing WRT fit update implements the convention that the cache is forwards in time of this effect, and insures the effect is only included once if(exing_->active()){ - // forwards, set the cache AFTER processing this effect if(tdir == TimeDir::forwards) { - kkdata.append(exing_->parameters(tdir)); - cache_ += kkdata.wData(); + kkdata.append(exing_->params(),tdir); + nextwt_ += kkdata.wData(); } else { - // backwards, set the cache BEFORE processing this effect, to avoid double-counting it - cache_ += kkdata.wData(); - kkdata.append(exing_->parameters(tdir)); + nextwt_ += kkdata.wData(); + kkdata.append(exing_->params(),tdir); } } } @@ -62,7 +61,7 @@ namespace KinKal { // update the ElementXing exing_->updateState(miconfig,first); // reset the cached weights - cache_ = Weights(); + nextwt_ = Weights(); } template void Material::append(PTRAJ& ptraj,TimeDir tdir) { @@ -74,24 +73,25 @@ namespace KinKal { (tdir == TimeDir::backwards && etime > ptraj.front().range().end()) ) throw std::invalid_argument("New piece overlaps existing"); KTRAJ newpiece = (tdir == TimeDir::forwards) ? ptraj.back() : ptraj.front(); - newpiece.params() = Parameters(cache_); // make sure the range includes the transit time newpiece.range() = (tdir == TimeDir::forwards) ? TimeRange(etime,std::max(ptraj.range().end(),etime+exing_->transitTime())) : TimeRange(std::min(ptraj.range().begin(),etime-exing_->transitTime()),etime); - if( tdir == TimeDir::forwards) + // invert the cached weight to define the parameters + newpiece.params() = Parameters(nextwt_); + if( tdir == TimeDir::forwards){ ptraj.append(newpiece); - else + } else { + // Since the cache was forwards of this site, we have to apply the effect of this material xing to the parameters. + auto temp = exing_->params(); + temp.parameters() *= -1; // reverse the sign of the parameter change + newpiece.params() += temp; ptraj.prepend(newpiece); + } } - // update the xing - if( tdir == TimeDir::forwards) - exing_->updateReference(ptraj.backPtr()); - else - exing_->updateReference(ptraj.frontPtr()); } - template void Material::updateReference(KTRAJPTR const& ltrajptr) { - exing_->updateReference(ltrajptr); + template void Material::updateReference(PTRAJ const& ptraj) { + exing_->updateReference(ptraj.nearestTraj(exing_->time())); } template void Material::print(std::ostream& ost,int detail) const { @@ -99,8 +99,8 @@ namespace KinKal { ost << " ElementXing "; exing_->print(ost,detail); if(detail >3){ - ost << " cache "; - cache().print(ost,detail); + ost << " forward cache "; + nextwt_.print(ost,detail); } } diff --git a/Fit/Measurement.hh b/Fit/Measurement.hh index 43bdb132..589b971e 100644 --- a/Fit/Measurement.hh +++ b/Fit/Measurement.hh @@ -15,7 +15,6 @@ namespace KinKal { public: using KKEFF = Effect; using PTRAJ = ParticleTrajectory; - using KTRAJPTR = std::shared_ptr; using HIT = Hit; using HITPTR = std::shared_ptr; // Effect Interface @@ -24,21 +23,23 @@ namespace KinKal { void process(FitState& kkdata,TimeDir tdir) override; void updateState(MetaIterConfig const& miconfig,bool first) override; void updateConfig(Config const& config) override {} - void updateReference(KTRAJPTR const& ltrajptr) override; + void updateReference(PTRAJ const& ptraj) override; void append(PTRAJ& fit,TimeDir tdir) override; Chisq chisq(Parameters const& pdata) const override; void print(std::ostream& ost=std::cout,int detail=0) const override; virtual ~Measurement(){} // local functions // construct from a hit and reference trajectory - Measurement(HITPTR const& hit); + Measurement(HITPTR const& hit,PTRAJ const& ptraj); // access the underlying hit HITPTR const& hit() const { return hit_; } private: - HITPTR hit_ ; // hit used for this constraint + HITPTR hit_ ; // hit used for this measurement }; - template Measurement::Measurement(HITPTR const& hit) : hit_(hit) {} + template Measurement::Measurement(HITPTR const& hit,PTRAJ const& ptraj) : hit_(hit) { + this->updateReference(ptraj); + } template void Measurement::process(FitState& kkdata,TimeDir tdir) { // add this effect's information. direction is irrelevant for processing hits @@ -51,22 +52,17 @@ namespace KinKal { } template void Measurement::append(PTRAJ& ptraj,TimeDir tdir) { - // update the hit to reference this trajectory. Use the end piece - if(tdir == TimeDir::forwards) - hit_->updateReference(ptraj.backPtr()); - else - hit_->updateReference(ptraj.frontPtr()); } - template void Measurement::updateReference(KTRAJPTR const& ltrajptr) { - hit_->updateReference(ltrajptr); + template void Measurement::updateReference(PTRAJ const& ptraj) { + hit_->updateReference(ptraj.nearestTraj(hit_->time())); } template Chisq Measurement::chisq(Parameters const& pdata) const { return hit_->chisq(pdata); } - template void Measurement::print(std::ostream& ost, int detail) const { + template void Measurement::print(std::ostream& ost, int detail) const { ost << "Measurement " << static_cast const&>(*this) << std::endl; if(detail > 0){ hit_->print(ost,detail); diff --git a/Fit/Status.hh b/Fit/Status.hh index 5bb6e093..d514174b 100644 --- a/Fit/Status.hh +++ b/Fit/Status.hh @@ -14,7 +14,7 @@ namespace KinKal { status status_; // current status Chisq chisq_; // current chisquared std::string comment_; // further information about the status - bool usable() const { return status_ < lowNDOF; } + 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){} static std::string statusName(status stat); diff --git a/Fit/Track.hh b/Fit/Track.hh index 732f1923..03d7c900 100644 --- a/Fit/Track.hh +++ b/Fit/Track.hh @@ -2,7 +2,7 @@ #define KinKal_Track_hh // // Primary class of the Kinematic Kalman fit. This class owns the state describing -// the fit inputs (measurements, material interactions, BField corrections, etc), the result of the fit, +// the fit inputs (measurements, material interactions, BField changes, etc), the result of the fit, // and the methods for computing it. The fit result is expressed as a piecewise kinematic covariant // particle trajectory, providing position, momentum etc information about the particle with covariance // as a function of physical time. @@ -40,14 +40,19 @@ #include "KinKal/Fit/Effect.hh" #include "KinKal/Fit/Measurement.hh" #include "KinKal/Fit/Material.hh" -#include "KinKal/Fit/BField.hh" +#include "KinKal/Fit/DomainWall.hh" #include "KinKal/Fit/Config.hh" +#include "KinKal/Fit/ExtraConfig.hh" #include "KinKal/Fit/Status.hh" +#include "KinKal/Fit/Domain.hh" #include "KinKal/General/BFieldMap.hh" +#include "KinKal/General/TimeRange.hh" #include "KinKal/General/TimeDir.hh" #include "TMath.h" #include #include +#include +#include #include #include #include @@ -57,6 +62,11 @@ #include namespace KinKal { + using DOMAINPTR = std::shared_ptr; + bool operator < (DOMAINPTR const& a, DOMAINPTR const& b) { return a->begin() < b->begin(); } + bool operator < (DOMAINPTR const& a, double val) { return a->begin() < val; } + bool operator < (double val, DOMAINPTR const& b) { return val < b->begin(); } + template class Track { public: using KKEFF = Effect; @@ -68,14 +78,15 @@ namespace KinKal { return false; } }; - using KKEFFCOL = std::vector>; // container type for effects - using KKEFFFWD = typename std::vector>::iterator; - using KKEFFREV = typename std::vector>::reverse_iterator; + + using KKEFFCOL = std::list>; + using KKEFFFWD = typename KKEFFCOL::iterator; + using KKEFFREV = typename KKEFFCOL::reverse_iterator; using KKEFFFWDBND = std::array; using KKEFFREVBND = std::array; using KKMEAS = Measurement; using KKMAT = Material; - using KKBFIELD = BField; + using KKDW = DomainWall; using PTRAJ = ParticleTrajectory; using PTRAJPTR = std::unique_ptr; using HIT = Hit; @@ -84,13 +95,17 @@ namespace KinKal { using EXING = ElementXing; using EXINGPTR = std::shared_ptr; using EXINGCOL = std::vector; - using DOMAINCOL = std::vector; + 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 ); // 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 with the given config until the given predicate is satisfied. This function requires + // the fit be valid, otherwise the return code is false. If successful the status, domains, and trajectory of the fit are updated + // Note that the actual fit itself is unchanged + template bool extrapolate(ExtraConfig const& xconfig, XTEST const& XTest); // accessors std::vector const& history() const { return history_; } Status const& fitStatus() const { return history_.back(); } // most recent status @@ -111,19 +126,24 @@ namespace KinKal { // helper functions TimeRange getRange(HITCOL& hits, EXINGCOL& exings) const; void fit(); // process the effects and create the trajectory. This executes the current schedule - void setBounds(KKEFFFWDBND& fwdbnds, KKEFFREVBND& revbnds); + 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 void iterate(MetaIterConfig const& miconfig); void setStatus(PTRAJPTR& ptrajptr); - void initFitState(FitStateArray& states, double dwt=1.0); + 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 replaceTraj(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); auto& status() { return history_.back(); } // most recent status - // divide a kinematic trajectory range into magnetic 'domains' within which the BField inhomogeneity effects are within tolerance - void createDomains(PTRAJ const& ptraj, TimeRange const& range, std::vector& ranges, TimeDir tdir=TimeDir::forwards) const; + // divide the range into magnetic 'domains' within which the BField can be considered constant (parameter change is within tolerance) + void createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains, double tol) const; // payload CONFIGCOL config_; // configuration BFieldMap const& bfield_; // magnetic field map @@ -133,7 +153,8 @@ namespace KinKal { 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_; // BField domains 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 ) : @@ -151,13 +172,13 @@ namespace KinKal { // 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 BField effects, define the domains + // if correcting for DomainWall effects, define the domains DOMAINCOL domains; - if(config().bfcorr_ ) createDomains(seedtraj_, refrange, domains); + if(config().bfcorr_ ) createDomains(seedtraj_, refrange, domains, config().tol_); // 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()); +// effects_.reserve(hits.size()+exings.size()+domains.size()); createEffects(hits,exings,domains); // now fit the track fit(); @@ -165,6 +186,8 @@ namespace KinKal { // extend an existing track template void Track::extend(Config const& cfg, HITCOL& hits, EXINGCOL& exings) { + // remember the previous config + auto const& oldconfig = config_.back(); // update the configuration config_.push_back(cfg); // configuation check @@ -180,31 +203,29 @@ namespace KinKal { DOMAINCOL domains; if(config().bfcorr_ ) { // if the previous configuration didn't have domains, then create them for the full reference range - auto oldconfig = config_.end(); --oldconfig; --oldconfig; - if(!oldconfig->bfcorr_){ + if(!oldconfig.bfcorr_ || oldconfig.tol_ != config().tol_){ // create domains for the whole range - createDomains(*fittraj_, exrange,domains); - // replace the reftraj with one with BField rotations - replaceTraj(domains); + createDomains(*fittraj_, exrange, domains, config().tol_); + // replace the domains. This also replace the trajectory, as that must reference the new domains + replaceDomains(domains); } else { - // create domains just for the extensions, and extend the reftraj as needed + // create domains just for the extensions TimeRange exlow(exrange.begin(),fittraj_->range().begin()); if(exlow.range()>0.0) { DOMAINCOL lowdomains; - createDomains(*fittraj_, exlow, lowdomains, TimeDir::backwards); - extendTraj(domains); - domains.insert(domains.begin(),lowdomains.begin(),lowdomains.end()); + createDomains(*fittraj_, exlow, lowdomains, config().tol_); + domains.insert(lowdomains.begin(),lowdomains.end()); } TimeRange exhigh(fittraj_->range().end(),exrange.end()); if(exhigh.range()>0.0){ DOMAINCOL highdomains; - createDomains(*fittraj_, exhigh, highdomains,TimeDir::forwards); - extendTraj(highdomains); - domains.insert(domains.end(),highdomains.begin(),highdomains.end()); + createDomains(*fittraj_, exhigh, highdomains, config().tol_); + domains.insert(highdomains.begin(),highdomains.end()); } } } - // create the effects for the new info and the new domains + // Extend the traj and create the effects for the new info and the new domains + extendTraj(domains); createEffects(hits,exings,domains); // update all the effects for this new configuration for(auto& ieff : effects_ ) ieff->updateConfig(config()); @@ -212,24 +233,39 @@ namespace KinKal { fit(); } - // replace the traj with one describing the 'same' trajectory in space, but using the local BField as reference - template void Track::replaceTraj(DOMAINCOL const& domains) { + // replace domains when DomainWall correction is added or changed. the traj must also be replaced, so that + // the pieces correspond to the new domains + template void Track::replaceDomains(DOMAINCOL const& domains) { + // if domains exist, clear them and remove all DomainWall effects + if(domains_.size() > 0){ + domains_.clear(); + // remove all existing DomainWall effects + auto ieff = effects_.begin(); + while(ieff != effects_.end()){ + const KKDW* kkbf = dynamic_cast(ieff->get()); + if(kkbf != 0){ + ieff = effects_.erase(ieff); + } else { + ++ieff; + } + } + } // create new traj auto newtraj = std::make_unique(); // loop over domains for(auto const& domain : domains) { - double dtime = domain.begin(); - // Set the BField to the start of this domain + double dtime = domain->begin(); + // Set the DomainWall to the start of this domain auto bf = bfield_.fieldVect(fittraj_->position3(dtime)); // loop until we're either out of this domain or the piece is out of this domain - while(dtime < domain.end()){ - // find the nearest piece of the current reftraj + while(dtime < domain->end()){ + // find the nearest piece of the traj auto index = fittraj_->nearestIndex(dtime); auto const& oldpiece = *fittraj_->pieces()[index]; // create a new piece KTRAJ newpiece(oldpiece,bf,dtime); // set the range as needed - double endtime = (index < fittraj_->pieces().size()-1) ? std::min(domain.end(),oldpiece.range().end()) : domain.end(); + double endtime = (index < fittraj_->pieces().size()-1) ? std::min(domain->end(),oldpiece.range().end()) : domain->end(); newpiece.range() = TimeRange(dtime,endtime); newtraj->append(newpiece); // update the time @@ -237,42 +273,38 @@ namespace KinKal { dtime = newpiece.range().end()+epsilon; // to avoid boundary } } - // update existing effects to reference this trajectory + // update all effects to reference this trajectory for (auto& eff : effects_) { - auto const& ltrajptr = newtraj->nearestTraj(eff->time()); - eff->updateReference(ltrajptr); + eff->updateReference(*newtraj); } - // swap fittraj_.swap(newtraj); } template void Track::extendTraj(DOMAINCOL const& domains ) { - // dummy implementation: need to put in parameter rotation at each domain boundary FIXME! if(domains.size() > 0){ - // extend the reftraj range - TimeRange newrange(std::min(fittraj_->range().begin(),domains.front().begin()), - std::max(fittraj_->range().end(),domains.back().end())); + TimeRange newrange(std::min(fittraj_->range().begin(),domains.begin()->get()->begin()), + std::max(fittraj_->range().end(),domains.rbegin()->get()->end())); fittraj_->setRange(newrange); } } template void Track::createTraj(PTRAJ const& seedtraj , TimeRange const& range, DOMAINCOL const& domains ) { - // if we're making local BField corrections, divide the trajectory into domain pieces. Each will have equivalent parameters, but relative + // 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(); for(auto const& domain : domains) { - // Set the BField to the start of this domain - auto bf = bfield_.fieldVect(seedtraj.position3(domain.begin())); - KTRAJ newpiece(seedtraj.nearestPiece(domain.begin()),bf,domain.begin()); - newpiece.range() = domain; + // 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()); + newpiece.range() = domain->range(); fittraj_->append(newpiece); } } else { - // use the middle of the range as the nominal BField for this fit: + // use the middle of the range as the nominal DomainWall 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 @@ -284,32 +316,31 @@ namespace KinKal { } template void Track::createEffects( HITCOL& hits, EXINGCOL& exings,DOMAINCOL const& domains ) { - // pre-allocate as needed - effects_.reserve(effects_.size()+hits.size()+exings.size()+domains.size()); - // append the effects. First, loop over the hits + // create and append the effects. First, loop over the hits for(auto& hit : hits ) { - // create the hit effects and insert them in the collection - effects_.emplace_back(std::make_unique(hit)); - // update hit reference; this should be done on construction FIXME - hit->updateReference(fittraj_->nearestTraj(hit->time())); + effects_.emplace_back(std::make_unique(hit,*fittraj_)); } //add material effects for(auto& exing : exings) { effects_.emplace_back(std::make_unique(exing,*fittraj_)); - // update xing reference; should be done on construction FIXME - exing->updateReference(fittraj_->nearestTraj(exing->time())); } - // add BField effects - for( auto const& domain : domains) { - // create the BField effect for integrated differences over this range - effects_.emplace_back(std::make_unique(config(),bfield_,domain)); + // add DomainWall effects + if(config().bfcorr_ && domains.size() > 1){ + auto nextdom = domains.cbegin(); + auto prevdom = nextdom; + ++nextdom; + while( nextdom != domains.cend() ){ + if(fabs(prevdom->get()->end()-nextdom->get()->begin())>1e-10)throw std::invalid_argument("Invalid domains"); + + effects_.emplace_back(std::make_unique(bfield_,*prevdom,*nextdom ,*fittraj_)); + prevdom = nextdom; + ++nextdom; + } } - // sort - std::sort(effects_.begin(),effects_.end(),KKEFFComp ()); - // store the inputs; these are just for convenience + // store the inputs; these are just for and may not be in time order hits_.insert(hits_.end(),hits.begin(),hits.end()); exings_.insert(exings_.end(),exings.begin(),exings.end()); - domains_.insert(domains_.end(),domains.begin(),domains.end()); + domains_.insert(domains.cbegin(),domains.cend()); } // fit the track @@ -332,7 +363,7 @@ namespace KinKal { } while(canIterate()); if(!status().usable())break; } - // if the fit is usable, process the passive effects on either end + // if the fit is usable, extend the track as necessary if(config().ends_ && status().usable()){ history_.push_back(fitStatus()); status().comment_ = "EndProcessing"; @@ -342,90 +373,97 @@ namespace KinKal { status().status_ = Status::failed; status().comment_ = status().comment_ + error.what(); } - } + } if(config().plevel_ > Config::none)print(std::cout, config().plevel_); } - // single algebraic iteration + // single algebraic iteration of the fit template void Track::iterate(MetaIterConfig const& miconfig) { if(config().plevel_ >= Config::basic)std::cout << "Processing fit iteration " << fitStatus().iter_ << std::endl; // update the effects for this configuration; this will sort the effects and find the iteration bounds bool first = status().iter_ == 0; // 1st iteration of a meta-iteration: update the effect internals for(auto& ieff : effects_ ) ieff->updateState(miconfig,first); // sort the sites, and set the iteration bounds - std::sort(effects_.begin(),effects_.end(),KKEFFComp ()); - KKEFFFWDBND fwdbnds; + effects_.sort(KKEFFComp()); + KKEFFFWDBND fwdbnds; // bounding sites used for fitting KKEFFREVBND revbnds; - setBounds(fwdbnds,revbnds ); - int ndof(0); ndof -= NParams(); + if(!setBounds(fwdbnds,revbnds)){ + status().status_ = Status::lowNDOF; + return; + } + // update the domains + TimeRange fitrange(fwdbnds[0]->get()->time(),revbnds[0]->get()->time()); + if(config().bfcorr_){ + // update the limits if new DW effects were added + if(extendDomains(fitrange,config().tol_))setBounds(fwdbnds,revbnds); + } + FitStateArray states; + initFitState(states, fitrange, config().dwt_/miconfig.varianceScale()); + // process the relevant effects forwards in time, adding their info to the fit state. Also compute chisquared for(auto feff=fwdbnds[0];feff!=fwdbnds[1];++feff){ - auto const* kkmeas = dynamic_cast(feff->get()); -// if(kkmeas && kkmeas->active())ndof += kkmeas->hit()->nDOF(); - if(kkmeas && kkmeas->active())ndof++; // this is more conservative than the above, but still not a complete test, since some measurements - // have redundant DOFs.FIXME - } - if(ndof >= (int)config().minndof_) { - // initialize the fit state to be used in this iteration, deweighting as specified. Not sure if using variance scale is right TODO - // To be consistent with hit errors I should scale by the ratio of current to previous temperature. Or maybe skip this?? FIXME - FitStateArray states; - initFitState(states, config().dwt_/miconfig.varianceScale()); - // loop over relevant effects, adding their info to the fit state. Also compute chisquared - for(auto feff=fwdbnds[0];feff!=fwdbnds[1];++feff){ - auto effptr = feff->get(); - // update chisquared increment WRT the current state: only needed once - Chisq dchisq = effptr->chisq(states[0].pData()); - status().chisq_ += dchisq; - // process - effptr->process(states[0],TimeDir::forwards); - if(config().plevel_ >= Config::detailed && dchisq.nDOF() > 0){ - std::cout << "Chisq increment " << dchisq << " "; - effptr->print(std::cout,config().plevel_-Config::detailed); - } - } - double mintime(std::numeric_limits::max()); - double maxtime(-std::numeric_limits::max()); - for(auto beff = revbnds[0]; beff!=revbnds[1]; ++beff){ - auto effptr = beff->get(); - effptr->process(states[1],TimeDir::backwards); - mintime = std::min(mintime,effptr->time()); - maxtime = std::max(maxtime,effptr->time()); + auto effptr = feff->get(); + // update chisquared increment WRT the current state: only needed once + Chisq dchisq = effptr->chisq(states[0].pData()); + status().chisq_ += dchisq; + // process + effptr->process(states[0],TimeDir::forwards); + if(config().plevel_ >= Config::detailed && dchisq.nDOF() > 0){ + std::cout << "Chisq increment " << dchisq << " "; + effptr->print(std::cout,config().plevel_-Config::detailed); } - // convert the fit result into a new trajectory - // initialize the parameters to the backward processing end - auto front = fittraj_->front(); - front.params() = states[1].pData(); - // extend range if needed -// TimeRange maxrange(std::min(fittraj_->range().begin(),fwdbnds[0]->get()->time()), -// std::max(fittraj_->range().end(),revbnds[0]->get()->time())); - TimeRange maxrange(mintime-0.1,maxtime+0.1); // FIXME - front.setRange(maxrange); - auto ptraj = std::make_unique(front); - // process forwards, adding pieces as necessary. This also sets the effects to reference the new trajectory - for(auto& ieff=fwdbnds[0]; ieff != fwdbnds[1]; ++ieff) { - ieff->get()->append(*ptraj,TimeDir::forwards); - } - setStatus(ptraj); // set the status for this iteration - // prepare for the next iteration: first, update the references for effects outside the fit range - // (the ones inside the range were updated above in 'append') - if(status().usable()){ - for(auto feff=fwdbnds[1]; feff != effects_.end(); ++feff) - feff->get()->updateReference(ptraj->nearestTraj(feff->get()->time())); - for(auto beff=revbnds[1]; beff != effects_.rend(); ++beff) - beff->get()->updateReference(ptraj->nearestTraj(beff->get()->time())); - } - // now all effects reference the new traj: we can swap it with the old + } + if(status().chisq_.nDOF() < (int)config().minndof_) { // I need a better way to define coverage as this test doesn't guarantee all parameters are constrained TODO + status().status_ = Status::lowNDOF; + return; + } + // Now process effects backwards in time + for(auto beff = revbnds[0]; beff!=revbnds[1]; ++beff){ + auto effptr = beff->get(); + effptr->process(states[1],TimeDir::backwards); + } + // convert the innermost state into a new trajectory + auto ptraj = initTraj(states[1],fitrange); + // append forwards in time + for(auto& ieff=fwdbnds[0]; ieff != fwdbnds[1]; ++ieff) { + ieff->get()->append(*ptraj,TimeDir::forwards); + } + setStatus(ptraj); // set the status for this iteration. This compares the trajs, so it must be done before swapping out the old fit + if(status().usable()){ + // prepare for the next iteration. Update the domains + updateDomains(*ptraj); + // update the effects to reference the new trajectory + for(auto& effptr : effects_) effptr->updateReference(*ptraj); fittraj_.swap(ptraj); + // now all effects reference the new traj: we can swap the fit to that. if(config().plevel_ >= Config::complete)fittraj_->print(std::cout,1); - } else { - status().chisq_ = Chisq(-1.0,ndof); - status().status_ = Status::lowNDOF; } } - // initialize statess used before iteration - template void Track::initFitState(FitStateArray& states, double dwt) { - auto fwdtraj = fittraj_->front(); - auto revtraj = fittraj_->back(); + template std::unique_ptr> Track::initTraj(FitState& state, TimeRange const& fitrange) { + // use the existing fit to define the parameterization + auto front = fittraj_->front(); + // initialize the parameters to the earliest backward processing end + front.params() = state.pData(); + // set bnom for the front piece parameters to the domain used + if(config().bfcorr_){ + // find the relevant domain + double ftime = fitrange.begin(); + for(auto const& domain : domains_) { + if(domain->range().inRange(ftime)){ + front.resetBNom(domain->bnom()); + break; + } + } + } + // set the range + front.setRange(fitrange); + return std::make_unique(front); + } + + // initialize states used before iteration + template void Track::initFitState(FitStateArray& states, TimeRange const& fitrange, double dwt) { + auto fwdtraj = fittraj_->nearestPiece(fitrange.begin()); + auto revtraj = fittraj_->nearestPiece(fitrange.end()); // dweight the covariance, scaled by the temperature. fwdtraj.params().covariance() *= dwt; revtraj.params().covariance() *= dwt; @@ -437,9 +475,16 @@ namespace KinKal { // finalize after iteration template void Track::setStatus(PTRAJPTR& ptraj) { - // to test for compute parameter difference WRT previous iteration. Compare at front and back ends - // to test for compute parameter difference WRT previous iteration. Compare at front and back ends + // compute parameter difference WRT previous iteration. Compare at front and back ends auto const& ffront = ptraj->front(); + // test covariance + auto fdiag = ffront.params().covariance().Diagonal(); + for(size_t ipar = 0; ipar < NParams(); ++ipar){ + if(fdiag[ipar] < 0.0 || std::isnan(fdiag[ipar])){ + status().status_ = Status::failed; + return; + } + } auto const& sfront = fittraj_->nearestPiece(ffront.range().mid()); DVEC dpfront = ffront.params().parameters() - sfront.params().parameters(); DMAT frontwt = sfront.params().covariance(); @@ -447,12 +492,19 @@ namespace KinKal { double dpchisqfront = ROOT::Math::Similarity(dpfront,frontwt); // back auto const& fback = ptraj->back(); + auto bdiag = fback.params().covariance().Diagonal(); + for(size_t ipar = 0; ipar < NParams(); ++ipar){ + if(bdiag[ipar] < 0.0 || std::isnan(bdiag[ipar])){ + status().status_ = Status::failed; + return; + } + } auto const& sback = fittraj_->nearestPiece(fback.range().mid()); DVEC dpback = fback.params().parameters() - sback.params().parameters(); DMAT backwt = sback.params().covariance(); if(! backwt.Invert())throw std::runtime_error("Reference covariance uninvertible"); double dpchisqback = ROOT::Math::Similarity(dpback,backwt); - // fit chisquared chang3 + // fit chisquared change double dchisq = config().convdchisq_ + 1e-4; // initialize to insure 0th iteration doesn't converge if(fitStatus().iter_ > 0){ auto prevstat = history_.rbegin(); @@ -478,9 +530,11 @@ namespace KinKal { status().status_ = Status::unconverged; } - // update between iterations - template void Track::setBounds(KKEFFFWDBND& fwdbnds, KKEFFREVBND& revbnds) { -// find bounds between first and last measurement + template bool Track::setBounds(KKEFFFWDBND& fwdbnds,KKEFFREVBND& revbnds) { + // find the bounding sites for algebraic processing this iteraiton. This excludes effects which can't change the fit + // find first measurement + bool retval(false); + fwdbnds[0] = effects_.end(); for(auto ieff=effects_.begin();ieff!=effects_.end();++ieff){ auto const* kkmeas = dynamic_cast(ieff->get()); if(kkmeas != 0 && kkmeas->active()){ @@ -489,30 +543,75 @@ namespace KinKal { break; } } - for(auto ieff=effects_.rbegin();ieff!=effects_.rend();++ieff){ - auto const* kkmeas = dynamic_cast(ieff->get()); - if(kkmeas != 0 && kkmeas->active()){ - revbnds[0] = ieff; - fwdbnds[1] = ieff.base(); - break; + if(fwdbnds[0] != effects_.end()){ + retval = true; + // now the last measurement + for(auto ieff=effects_.rbegin();ieff!=effects_.rend();++ieff){ + auto const* kkmeas = dynamic_cast(ieff->get()); + if(kkmeas != 0 && kkmeas->active()){ + revbnds[0] = ieff; + fwdbnds[1] = ieff.base(); + break; + } + } + } + return retval; + } + + template void Track::updateDomains(PTRAJ const& ptraj) { + for(auto& domain : domains_) { + domain->updateBNom(bfield_.fieldVect(ptraj.position3(domain->mid()))); + } + } + + template bool Track::extendDomains(TimeRange const& fitrange, double tol) { + bool retval(false); + // then, check if the range has + TimeRange drange(domains().begin()->get()->begin(),domains().rbegin()->get()->end()); + if(!drange.contains(fitrange)){ + retval = true; + // we need to extend the domains. First backwards + if(drange.begin() > fitrange.begin()){ + double time = drange.begin(); + while(time > fitrange.begin()){ + auto const& ktraj = fittraj_->nearestPiece(time); + double dt = bfield_.rangeInTolerance(ktraj,time,tol); + TimeRange range(time-dt,time); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),tol); + addDomain(domain,TimeDir::backwards); + time = domain.begin(); + } + } + // then forwards + if(drange.end() < fitrange.end()){ + double time = drange.end(); + while(time < fitrange.end()){ + auto const& ktraj = fittraj_->nearestPiece(time); + double dt = bfield_.rangeInTolerance(ktraj,time,tol); + TimeRange range(time,time+dt); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),tol); + addDomain(domain,TimeDir::forwards); + time = domain.end(); + } } } + return retval; } template void Track::processEnds() { - // sort the end sites - KKEFFFWDBND fwdbnds; // bounds for iterating + // sort effects in case ends have migrated + effects_.sort(KKEFFComp()); + // update domains as needed to cover the end effects + TimeRange endrange(effects_.front()->time(),effects_.back()->time()); + extendDomains(endrange,config().tol_); + KKEFFFWDBND fwdbnds; // bounding sites used for fitting KKEFFREVBND revbnds; setBounds(fwdbnds,revbnds); - // update the end effects: this makes their internal content consistent with the others - // use the final meta-iteration - for(auto feff=fwdbnds[1]; feff != effects_.end(); ++feff) - feff->get()->updateState(config().schedule().back(),false); - for(auto reff=revbnds[1]; reff != effects_.rend(); ++reff) - reff->get()->updateState(config().schedule().back(),false); - // then process these sites. Start with the state at the apropriate end, but without any deweighting + // initialize the fit state where we left off processing FitStateArray states; - initFitState(states, 1.0); // no deweighting + TimeRange fitrange(fwdbnds[0]->get()->time(),revbnds[0]->get()->time()); + initFitState(states, fitrange, 1.0); // no deweighting + // process forwards and backwards for(auto feff=fwdbnds[1]; feff != effects_.end(); ++feff) feff->get()->process(states[1],TimeDir::forwards); for(auto reff=revbnds[1]; reff != effects_.rend(); ++reff) @@ -548,19 +647,21 @@ namespace KinKal { for(auto const& eff : effects()) eff.get()->print(ost,detail-3); } } - // divide a trajectory into magnetic 'domains' used to apply the BField corrections - template void Track::createDomains(PTRAJ const& ptraj, TimeRange const& range, std::vector& ranges, - TimeDir tdir) const { - double tstart; - tstart = range.begin(); + // divide a trajectory into magnetic 'domains' used to apply the DomainWall corrections + template void Track::createDomains(PTRAJ const& ptraj, TimeRange const& range, DOMAINCOL& domains, double tol) const { + auto const& ktraj = ptraj.nearestPiece(range.begin()); + 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 BField change causes the momentum estimate to go out of tolerance + // 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); - double tend = bfield_.rangeInTolerance(ktraj,tstart,config().tol_); - ranges.emplace_back(tstart,tend); - // start the next domain and the end of this one - tstart = tend; - } while(tstart < range.end()); + 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 } template TimeRange Track::getRange(HITCOL& hits, EXINGCOL& exings) const { @@ -576,5 +677,49 @@ namespace KinKal { } return TimeRange(tmin,tmax); } + + template template bool Track::extrapolate(ExtraConfig const& xconfig, XTEST const& xtest) { + bool retval(false); + if(this->fitStatus().usable()){ + if(config().bfcorr_){ + // iterate until the extrapolation condition is met + double time = xconfig.xdir_ == TimeDir::forwards ? domains_.crbegin()->get()->end() : domains_.cbegin()->get()->begin(); + double tstart = time; + while(fabs(time-tstart) < xconfig.maxdt_ && xtest.needsExtrapolation(time) ){ + // create a domain for this extrapolation + auto const& ktraj = fittraj_->nearestPiece(time); + double dt = bfield_.rangeInTolerance(ktraj,time,xconfig.tol_); // always positive + TimeRange range = xconfig.xdir_ == TimeDir::forwards ? TimeRange(time,time+dt) : TimeRange(time-dt,time); + Domain domain(range,bfield_.fieldVect(ktraj.position3(range.mid())),xconfig.tol_); + addDomain(domain,xconfig.xdir_); + time = xconfig.xdir_ == TimeDir::forwards ? domain.end() : domain.begin(); + } + retval = true; + } + } else { + retval = true; + } + return retval; + } + + template void Track::addDomain(Domain const& domain,TimeDir const& tdir) { + auto dptr = std::make_shared(domain); + if(tdir == TimeDir::forwards){ + auto const& prevdom = *domains_.crbegin(); + auto const& ktraj = fittraj_->nearestPiece(prevdom->end()); + FitState fstate(ktraj.params()); + effects_.emplace_back(std::make_unique(bfield_,prevdom,dptr,ktraj)); + effects_.back()->process(fstate,tdir); + effects_.back()->append(*fittraj_,tdir); + } else { + auto const& nextdom = *domains_.cbegin(); + auto const& ktraj = fittraj_->nearestPiece(nextdom->begin()); + FitState fstate(ktraj.params()); + effects_.emplace_front(std::make_unique(bfield_,dptr,nextdom,ktraj)); + effects_.front()->process(fstate,tdir); + effects_.front()->append(*fittraj_,tdir); + } + domains_.insert(dptr); + } } #endif diff --git a/General/BFieldMap.hh b/General/BFieldMap.hh index 6e2cd7a9..8ed493db 100644 --- a/General/BFieldMap.hh +++ b/General/BFieldMap.hh @@ -56,7 +56,7 @@ namespace KinKal { } // estimate how long the momentum vector from the given trajectory will stay within the given (fractional) tolerance given the field spatial variation - // ie mag(P_true(tstart+dt) - P_traj(tstart+dt)) < tol. This is good to 1st order (ignores trajectory curvature) + // ie mag(P_true(tstart+dt) - P_traj(tstart+dt)) < tol. This is a 2nd order local calculation template double BFieldMap::rangeInTolerance(KTRAJ const& ktraj, double tstart, double tol) const { auto tpos = ktraj.position3(tstart); // starting position double dp = ktraj.momentum(tstart)*tol; // fractional tolerance on momentum @@ -64,9 +64,9 @@ namespace KinKal { auto dBdt = fieldDeriv(tpos,vel); // change in field WRT time along this velocity double d2pdt2 = (dBdt.Cross(vel)).R()*cbar()*fabs(ktraj.charge()); // 2nd derivative of momentum due to B change along the path if(d2pdt2 > 1e-10) - return tstart + sqrt(dp/d2pdt2); + return sqrt(dp/d2pdt2); else - return ktraj.range().end(); + return ktraj.range().range(); } // trivial instance of the above, used for testing diff --git a/General/CMakeLists.txt b/General/CMakeLists.txt index 8de26380..8e392d4e 100644 --- a/General/CMakeLists.txt +++ b/General/CMakeLists.txt @@ -26,6 +26,8 @@ ROOT_GENERATE_DICTIONARY(GeneralDict FitData.hh TimeRange.hh Vectors.hh + PhysicalConstants.h + SystemOfUnits.h LINKDEF LinkDef.h MODULE General ) diff --git a/Geometry/CMakeLists.txt b/Geometry/CMakeLists.txt index b7e5b838..b67fa591 100644 --- a/Geometry/CMakeLists.txt +++ b/Geometry/CMakeLists.txt @@ -13,6 +13,7 @@ add_library(Geometry SHARED Intersection.cc Plane.cc Ray.cc + LineSegment.cc Rectangle.cc ThinSolid.cc ) diff --git a/Geometry/Cylinder.cc b/Geometry/Cylinder.cc index 510b55d8..a816a026 100644 --- a/Geometry/Cylinder.cc +++ b/Geometry/Cylinder.cc @@ -55,7 +55,11 @@ namespace KinKal { } Disk Cylinder::frontDisk() const { - return Disk(axis_,uDirection(),center_-halflen_*axis_,radius_); + return Disk(axis_,uDirection(),center_-halflen_*axis_,radius_); + } + + Disk Cylinder::midDisk() const { + return Disk(axis_,uDirection(),center_,radius_); } Disk Cylinder::backDisk() const { @@ -80,13 +84,13 @@ namespace KinKal { IntersectFlag Cylinder::intersect(Ray const& ray,double& dist, bool forwards, double tol) const { IntersectFlag retval; - double ddot = ray.dir_.Dot(axis_); + double ddot = ray.direction().Dot(axis_); double alpha = (1.0 - ddot*ddot); // always positive // make sure the ray isn't co-linear on the relevant scale if(alpha > tol/std::max(radius_,halflen_)){ - auto rvec = ray.start_ - center_; + auto rvec = ray.start() - center_; double sdot = rvec.Dot(axis_); - double beta = sdot*ddot - rvec.Dot(ray.dir_); + double beta = sdot*ddot - rvec.Dot(ray.direction()); double gamma = rvec.Mag2() - sdot*sdot - radius2_; double beta2 = beta*beta; double ag = alpha*gamma; diff --git a/Geometry/Cylinder.hh b/Geometry/Cylinder.hh index 6c089855..b8dfc295 100644 --- a/Geometry/Cylinder.hh +++ b/Geometry/Cylinder.hh @@ -29,6 +29,7 @@ namespace KinKal { double halfLength() const { return halflen_; } // front and rear boundaries of this cylinder as disks. U direction is arbitrary Disk frontDisk() const; + Disk midDisk() const; Disk backDisk() const; // inscribed rectangle with given normal direction. U direction is long the axis Rectangle inscribedRectangle(VEC3 const& norm) const; diff --git a/Geometry/Frustrum.cc b/Geometry/Frustrum.cc index 0edd3829..25a921eb 100644 --- a/Geometry/Frustrum.cc +++ b/Geometry/Frustrum.cc @@ -79,6 +79,10 @@ namespace KinKal { return Disk(axis_,uDirection(),center_-halflen_*axis_,radius(-halflen_)); } + Disk Frustrum::midDisk() const { + return Disk(axis_,uDirection(),center_,radius(0.0)); + } + Disk Frustrum::backDisk() const { return Disk(axis_,uDirection(),center_+halflen_*axis_,radius(halflen_)); } @@ -100,10 +104,10 @@ namespace KinKal { // exact solution IntersectFlag retval; // displace the ray start to be WRT the base of the cone - auto spos = ray.start_ - center_ + halflen_*axis_; - double z1 = ray.dir_.Dot(axis_); + auto spos = ray.start() - center_ + halflen_*axis_; + double z1 = ray.direction().Dot(axis_); double z0 = spos.Dot(axis_); - auto pvec = ray.dir_ - z1*axis_; + auto pvec = ray.direction() - z1*axis_; auto qvec = spos - z0*axis_; double b2 = pvec.Mag2() - z1*z1*drda_*drda_; double b1 = z1*drda_*(drda_*z0 + r0_) - pvec.Dot(qvec); diff --git a/Geometry/Frustrum.hh b/Geometry/Frustrum.hh index e0cb89c3..f394b054 100644 --- a/Geometry/Frustrum.hh +++ b/Geometry/Frustrum.hh @@ -34,6 +34,7 @@ namespace KinKal { double halfAngle() const { return atan2(sint_,cost_); } // front and rear boundaries of this cylinder as disks. U direction is arbitrary Disk frontDisk() const; + Disk midDisk() const; Disk backDisk() const; private: VEC3 axis_; // symmetry axis of the cylinder diff --git a/Geometry/Intersect.hh b/Geometry/Intersect.hh index 22e3ff7b..8247a49a 100644 --- a/Geometry/Intersect.hh +++ b/Geometry/Intersect.hh @@ -34,7 +34,7 @@ namespace KinKal { stepinside = surf.isInside(pos); } while (startinside == stepinside && trange.inRange(ttest)); if(startinside != stepinside){ - // we crossed the cylinder: backup and do a linear search. + // we crossed the surface: backup and do a linear search. ttest -= tstep; double dist; IntersectFlag rayinter; @@ -142,8 +142,9 @@ namespace KinKal { auto axis = helix.axis(trange.begin()); // test for the helix being circular or tangent to the plane double vz = helix.axisSpeed(); // speed along the helix axis - double ddot = fabs(axis.dir_.Dot(plane.normal())); - if(fabs(vz*trange.range()) > tol && ddot > tol ){ + double ddot = fabs(axis.direction().Dot(plane.normal())); + double zrange = fabs(vz*trange.range()); + if(zrange > tol && ddot > tol/zrange ){ // Find the intersection time of the helix axis (along bnom) with the plane double dist(0.0); auto pinter = plane.intersect(axis,dist,true,tol); @@ -271,12 +272,12 @@ namespace KinKal { // use pointers to cast to avoid avoid a throw const Surface* surfp = &surf; // go through the possibilities: I don't know of anything more elegant + auto plane = dynamic_cast(surfp); + if(plane)return intersect(ktraj,*plane,trange,tol); auto cyl = dynamic_cast(surfp); if(cyl)return intersect(ktraj,*cyl,trange,tol); auto fru = dynamic_cast(surfp); if(fru)return intersect(ktraj,*fru,trange,tol); - auto plane = dynamic_cast(surfp); - if(plane)return intersect(ktraj,*plane,trange,tol); // unknown surface subclass; return failure return Intersection(); } diff --git a/Geometry/LineSegment.cc b/Geometry/LineSegment.cc new file mode 100644 index 00000000..1098e0cf --- /dev/null +++ b/Geometry/LineSegment.cc @@ -0,0 +1,19 @@ +#include "KinKal/Geometry/LineSegment.hh" +#include +#include +#include +#include + +namespace KinKal { + + void LineSegment::print(std::ostream& ost, int detail) const { + ost << *this << std::endl; + } + + std::ostream& operator <<(std::ostream& ost, LineSegment const& lineseg) { + ost << " LineSegment starting " << lineseg.start() + << " direction " << lineseg.direction() + << " length " << lineseg.length(); + return ost; + } +} diff --git a/Geometry/LineSegment.hh b/Geometry/LineSegment.hh new file mode 100644 index 00000000..682d4062 --- /dev/null +++ b/Geometry/LineSegment.hh @@ -0,0 +1,31 @@ +#ifndef KinKal_LineSegment_hh +#define KinKal_LineSegment_hh +// +// directed line segment +// +#include "KinKal/General/Vectors.hh" +#include "KinKal/Geometry/Ray.hh" + +namespace KinKal { + class LineSegment : public Ray { + public: + // construct from a point and direction + LineSegment(VEC3 const& dir, VEC3 const& start, double length) : Ray(dir,start), length_(length) { + end_ = position(length_); + } + // construct from 2 points + LineSegment(VEC3 const& start, VEC3 const& end) : Ray(end-start,start), length_((end-start).R()), end_(end) {} + // accessors + VEC3 const& end() const { return end_; } + VEC3 middle() const { return 0.5*(start() + end()); } + double length() const { return length_; } + void print(std::ostream& ost, int detail) const; + // position from the end + VEC3 endPosition(double dist) const { return end_ - dist*direction(); } + private: + double length_; // segment length + VEC3 end_; // end position + }; + std::ostream& operator <<(std::ostream& ost, LineSegment const& tLineSegment); +} +#endif diff --git a/Geometry/ParticleTrajectoryIntersect.hh b/Geometry/ParticleTrajectoryIntersect.hh index 50c93316..2d1c8f18 100644 --- a/Geometry/ParticleTrajectoryIntersect.hh +++ b/Geometry/ParticleTrajectoryIntersect.hh @@ -9,11 +9,11 @@ #include "KinKal/Geometry/Intersect.hh" namespace KinKal { // Find first intersection of a particle trajectory in the specified range. This generic implementation tests - template Intersection pIntersect(ParticleTrajectory const& ptraj, SURF const& surf, TimeRange trange, double tol) { + template Intersection pIntersect(ParticleTrajectory const& ptraj, SURF const& surf, TimeRange trange, double tstart, double tol) { Intersection retval; // loop over pieces, and test the ones in range VEC3 spos, epos; - auto curr = ptraj.nearestTraj(trange.begin()); + auto curr = ptraj.nearestTraj(tstart); auto prev = curr; // loop until we find the best piece unsigned ntries(0); @@ -35,27 +35,54 @@ namespace KinKal { return retval; } // KinematicLine-based particle trajectory intersect implementation can always use the generic function - Intersection intersect(ParticleTrajectory const& kklptraj, KinKal::Surface const& surf, TimeRange trange,double tol) { - return pIntersect(kklptraj,surf,trange,tol); + Intersection intersect(ParticleTrajectory const& kklptraj, KinKal::Surface const& surf, TimeRange trange, double tol) { + return pIntersect(kklptraj,surf,trange,trange.begin(),tol); } // Helix-based particle trajectory intersect implementation with a plane template Intersection phpIntersect(ParticleTrajectory const& phelix, KinKal::Plane const& plane, TimeRange trange ,double tol) { - // for now, call generic function. In future, we can do a smarter binary search for the correct piece using the 'constant' - // z velocity - return pIntersect(phelix,plane,trange,tol); + // use the middle range time to estimate the start time + auto const& midhelix = phelix.nearestPiece(trange.mid()); + auto axis = midhelix.axis(trange.mid()); + double dist; // distance from ray start to the plane + auto ainter = plane.intersect(axis,dist,false,tol); + double tstart = trange.begin(); // backup if the following fails due to pathological geometry + if(ainter.onsurface_){ + double vz = midhelix.axisSpeed(); // speed along the helix axis + tstart = trange.mid() + dist/vz; + } + return pIntersect(phelix,plane,trange,tstart,tol); } template < class HELIX> Intersection phcIntersect( ParticleTrajectory const& phelix, KinKal::Cylinder const& cyl, TimeRange trange ,double tol) { - // for now, call generic function. In future, we can call the above intersection on the end disks to find the correct range more efficiently - return pIntersect(phelix,cyl,trange,tol); + // use the middle range time to estimate the start time + auto const& midhelix = phelix.nearestPiece(trange.mid()); + auto axis = midhelix.axis(trange.mid()); + double dist; // distance to the midplane + auto mdisk = cyl.midDisk(); + auto ainter = mdisk.intersect(axis,dist,false,tol); + double tstart = trange.begin(); // backup if the following fails due to pathological geometry + if(ainter.onsurface_ ){ + double vz = midhelix.axisSpeed(); // speed along the helix axis + tstart = trange.mid() + dist/vz; // time the axis reaches the midplane + } + return pIntersect(phelix,cyl,trange,tstart,tol); } template < class HELIX> Intersection phfIntersect( ParticleTrajectory const& phelix, KinKal::Frustrum const& fru, TimeRange trange ,double tol) { - // for now, call generic function. In future, we can call the above intersection on the end disks to find the correct range more efficiently - return pIntersect(phelix,fru,trange,tol); + // use the middle range time to estimate the start time + auto const& midhelix = phelix.nearestPiece(trange.mid()); + auto axis = midhelix.axis(trange.mid()); + double dist; // distance to the midplane + auto mdisk = fru.midDisk(); + auto ainter = mdisk.intersect(axis,dist,false,tol); + double tstart = trange.begin(); // backup if the following fails due to pathological geometry + if(ainter.onsurface_ ){ + double vz = midhelix.axisSpeed(); // speed along the helix axis + tstart = trange.mid() + dist/vz; // time the axis reaches the midplane + } + return pIntersect(phelix,fru,trange,tstart,tol); } - // explicit 'specializations' for the different helix types Intersection intersect( ParticleTrajectory const& ploophelix, KinKal::Cylinder const& cyl, TimeRange trange ,double tol) { @@ -82,12 +109,12 @@ namespace KinKal { // use pointers to cast to avoid avoid a throw const Surface* surfp = &surf; // go through the possibilities: I don't know of anything more elegant + auto plane = dynamic_cast(surfp); + if(plane)return intersect(pktraj,*plane,trange,tol); auto cyl = dynamic_cast(surfp); if(cyl)return intersect(pktraj,*cyl,trange,tol); auto fru = dynamic_cast(surfp); if(fru)return intersect(pktraj,*fru,trange,tol); - auto plane = dynamic_cast(surfp); - if(plane)return intersect(pktraj,*plane,trange,tol); // unknown surface subclass; return failure return Intersection(); } diff --git a/Geometry/Plane.cc b/Geometry/Plane.cc index 4a231f69..2b6b2b66 100644 --- a/Geometry/Plane.cc +++ b/Geometry/Plane.cc @@ -22,9 +22,9 @@ namespace KinKal { IntersectFlag Plane::intersect(Ray const& ray,double& dist, bool forwards, double tol) const { IntersectFlag retval; - double ddir = norm_.Dot(ray.dir_); + double ddir = norm_.Dot(ray.direction()); if(fabs(ddir)>0.0) { - double pdist = norm_.Dot(center_ - ray.start_); + double pdist = norm_.Dot(center_ - ray.start()); dist = pdist/ddir; if(dist > 0.0 || !forwards){ retval.onsurface_ = true; diff --git a/Geometry/Plane.hh b/Geometry/Plane.hh index 4a44ce24..19da7b7c 100644 --- a/Geometry/Plane.hh +++ b/Geometry/Plane.hh @@ -24,7 +24,7 @@ namespace KinKal { auto const& uDirection() const { return udir_; } auto const& vDirection() const { return vdir_; } auto const& normal() const { return norm_; } - Plane tangentPlane(VEC3 const& point) const override { return Plane(norm_,udir_,point); } + Plane tangentPlane(VEC3 const& ) const override { return *this; } // plane interfac3 auto const& center() const { return center_; } private: diff --git a/Geometry/Ray.cc b/Geometry/Ray.cc index 517ea164..defa34aa 100644 --- a/Geometry/Ray.cc +++ b/Geometry/Ray.cc @@ -1,5 +1,5 @@ #include "KinKal/Geometry/Ray.hh" std::ostream& operator <<(std::ostream& ost, KinKal::Ray const& ray){ - ost << "Ray starting " << ray.start_ << " direction " << ray.dir_; + ost << "Ray starting " << ray.start() << " direction " << ray.direction(); return ost; } diff --git a/Geometry/Ray.hh b/Geometry/Ray.hh index e25044ce..2261b0bd 100644 --- a/Geometry/Ray.hh +++ b/Geometry/Ray.hh @@ -7,12 +7,18 @@ #include "KinKal/General/Vectors.hh" #include namespace KinKal { - struct Ray { - // construct, making sure direciton is unit - Ray(VEC3 const& dir, VEC3 const& start) : dir_(dir.Unit()), start_(start){} - VEC3 dir_; // direction - VEC3 start_; // starting position - VEC3 position(double distance) const { return start_ + distance*dir_; } + class Ray { + public: + // construct, making sure direciton is unit + Ray(VEC3 const& dir, VEC3 const& start) : dir_(dir.Unit()), start_(start){} + VEC3 position(double distance) const { return start_ + distance*dir_; } + VEC3 const& direction() const { return dir_; } + VEC3 const& start() const { return start_; } + // Distance along the the ray to the point of Closest Approach to a point + double distanceTo(VEC3 const& point) const { return (point - start_).Dot(dir_); } + private: + VEC3 dir_; // direction + VEC3 start_; // starting position }; } std::ostream& operator <<(std::ostream& ost, KinKal::Ray const& ray); diff --git a/Tests/BFieldMapTest.hh b/Tests/BFieldMapTest.hh index 29306dd5..530a811b 100644 --- a/Tests/BFieldMapTest.hh +++ b/Tests/BFieldMapTest.hh @@ -2,7 +2,6 @@ // test basic functions of BFieldMap class // #include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/Trajectory/Line.hh" #include "KinKal/Trajectory/ClosestApproach.hh" #include "KinKal/General/BFieldMap.hh" #include "KinKal/General/Vectors.hh" @@ -52,7 +51,7 @@ int BFieldMapTest(int argc, char **argv) { double pmass(0.511); BFieldMap *BF(0); double Bgrad(0.0), dBx(0.0), dBy(0.0), dBz(0.0); - double tol(0.1); + double tol(1e-4), simtol(1e-5); double zrange(3000.0); // tracker dimension static struct option long_options[] = { @@ -63,6 +62,7 @@ int BFieldMapTest(int argc, char **argv) { {"dBz", required_argument, 0, 'Z' }, {"Bgrad", required_argument, 0, 'g' }, {"Tol", required_argument, 0, 't' }, + {"SimTol", required_argument, 0, 's' }, {NULL, 0,0,0} }; int opt; @@ -82,6 +82,8 @@ int BFieldMapTest(int argc, char **argv) { break; case 't' : tol = atof(optarg); break; + case 's' : simtol = atof(optarg); + break; case 'q' : icharge = atoi(optarg); break; default: print_usage(); @@ -100,11 +102,16 @@ int BFieldMapTest(int argc, char **argv) { } // first, create a traj based on the actual field at this point KKTest::ToyMC toy(*BF, mom, icharge, zrange, iseed, 0, false, false, -1.0, pmass ); - toy.setTolerance(tol); + toy.setTolerance(simtol); PTRAJ tptraj; - HITCOL thits; - EXINGCOL dxings; - toy.simulateParticle(tptraj, thits, dxings); +// HITCOL thits; +// EXINGCOL dxings; +// toy.simulateParticle(tptraj, thits, dxings); + toy.createTraj(tptraj); + double tbeg = tptraj.range().begin(); + auto vel = tptraj.velocity(tbeg); + double tend = tbeg + zrange/vel.Z(); + toy.extendTraj(tptraj,tend); // then, create a piecetraj around the nominal field with corrections, auto pos = tptraj.position4(tptraj.range().begin()); auto momv = tptraj.momentum4(pos.T()); @@ -115,7 +122,7 @@ int BFieldMapTest(int argc, char **argv) { TimeRange prange = start.range(); do { auto const& piece = xptraj.back(); - prange = TimeRange(prange.begin(),BF->rangeInTolerance(piece,prange.begin(), tol)); + prange = TimeRange(prange.begin(),prange.begin() + BF->rangeInTolerance(piece,prange.begin(), tol)); // integrate the momentum change over this range VEC3 dp = BF->integrate(piece,prange); // approximate change in position @@ -154,9 +161,9 @@ int BFieldMapTest(int argc, char **argv) { xdp = BF->integrate( xptraj, xptraj.range()); ldp = BF->integrate( lptraj, lptraj.range()); ndp = BF->integrate( start, start.range()); - cout << "TTraj " << tptraj << " integral " << tdp << endl; - cout << "XTraj " << xptraj << " integral " << xdp << endl; - cout << "LTraj " << lptraj << " integral " << ldp << endl; + cout << "Simulated Traj " << tptraj << " integral " << tdp << endl; + cout << "Extract Traj " << xptraj << " integral " << xdp << endl; + cout << "Linear Traj " << lptraj << " integral " << ldp << endl; cout << "Nominal " << start << " integral " << ndp << endl; // setup histograms diff --git a/Tests/CentralHelixFit_unit.cc b/Tests/CentralHelixFit_unit.cc index 2606bec6..f1e451e3 100644 --- a/Tests/CentralHelixFit_unit.cc +++ b/Tests/CentralHelixFit_unit.cc @@ -6,10 +6,15 @@ int main(int argc, char **argv) { cout << "Testing gradient field, correction 2" << endl; std::vector arguments; arguments.push_back(argv[0]); - arguments.push_back("--Bgrad"); - arguments.push_back("-0.036"); // mu2e-like field gradient +// arguments.push_back("--Bgrad"); +// arguments.push_back("-0.036"); // mu2e-like field gradient arguments.push_back("--Schedule"); - arguments.push_back("driftfit.txt"); + arguments.push_back("MCAmbigFixedField.txt"); +// arguments.push_back("seedfit.txt"); +// arguments.push_back("--extend"); +// arguments.push_back("driftextend.txt"); +// arguments.push_back("--ttree"); +// arguments.push_back("1"); std::vector myargv; for (const auto& arg : arguments) myargv.push_back((char*)arg.data()); diff --git a/Tests/ClosestApproachTest.hh b/Tests/ClosestApproachTest.hh index 79b868c8..99e796b4 100644 --- a/Tests/ClosestApproachTest.hh +++ b/Tests/ClosestApproachTest.hh @@ -1,7 +1,7 @@ // // test basic functions of ClosestApproach using KTraj and Line // -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/ClosestApproach.hh" #include "KinKal/Trajectory/PointClosestApproach.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" @@ -40,8 +40,6 @@ using namespace KinKal; using namespace std; -// avoid confusion with root -using KinKal::Line; void print_usage() { printf("Usage: ClosestApproachTest --charge i--gap f --tmin f --tmax f --vprop f --delta f \n"); @@ -51,9 +49,9 @@ template int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ gROOT->SetBatch(kTRUE); gStyle->SetOptFit(1); - using TCA = ClosestApproach; + using TCA = ClosestApproach; using TCAP = PointClosestApproach; - using PCA = PiecewiseClosestApproach; + using PCA = PiecewiseClosestApproach; using PTRAJ = ParticleTrajectory; int opt; int status(0); @@ -142,7 +140,7 @@ int ClosestApproachTest(int argc, char **argv, KinKal::DVEC pchange ){ double gap = tr_.Uniform(0.01,maxgap); VEC3 spos = ppos + gap*docadir; // create the Line - Line tline(spos, time, svel, wlen); + SensorLine tline(spos, time, svel, wlen); // create ClosestApproach from these CAHint tphint(time,time); TCA tp(ktraj,tline,tphint,1e-8); diff --git a/Tests/Derivatives.hh b/Tests/Derivatives.hh index d58ca713..d1197276 100644 --- a/Tests/Derivatives.hh +++ b/Tests/Derivatives.hh @@ -431,16 +431,17 @@ int test(int argc, char **argv) { << dPdM << endl; } - // test changes due to BFieldMap + // test changing bnom TCanvas* dbcan[3]; // 3 directions std::vector bgraphs[3]; std::array basis; basis[0] = reftraj.bnom().Unit(); basis[1] = reftraj.direction(MomBasis::phidir_); basis[2] = VEC3(basis[1].Y(),-basis[1].X(), 0.0); //perp - std::array bnames{"BDirection", "PhiDirection", "MomPerpendicular"}; + std::array bnames{"BDir", "PhiDir", "PerpDir"}; // gaps TGraph* bgapgraph[3]; + auto state = reftraj.stateEstimate(ttest); for(size_t idir =0; idir<3;idir++){ bgraphs[idir] = std::vector(NParams(),0); for(size_t ipar = 0; ipar < NParams(); ipar++){ @@ -457,24 +458,33 @@ int test(int argc, char **argv) { for(int id=0;id1.0e-6) cout << "State vector " << ipar << " doesn't match: original " - << state.state()[ipar] << " rotated " << newstate.state()[ipar] << endl; - } - dpx = newbfhel.params().parameters() - reftraj.params().parameters(); - // 1st order change trajectory - KTRAJ dbtraj(reftraj,bf,ttest); + KTRAJ xdbtraj(state,bf); + dpx = xdbtraj.params().parameters() - reftraj.params().parameters(); + // test 1st order change + auto dbtraj = reftraj; + dbtraj.setBNom(ttest, bf); dpdb = dbtraj.params().parameters() - reftraj.params().parameters(); for(size_t ipar = 0; ipar < NParams(); ipar++){ bgraphs[idir][ipar]->SetPoint(id,dpx[ipar], dpdb[ipar]); } - bgapgraph[idir]->SetPoint(id,dval,(dbtraj.position3(ttest)-newbfhel.position3(ttest)).R()); - // BField derivatives + bgapgraph[idir]->SetPoint(id,dval,(dbtraj.position3(ttest)-reftraj.position3(ttest)).R()); + // test state + if(id==0){ + // exact state test + auto xdbstate = xdbtraj.stateEstimate(ttest); + for(size_t ipar=0;ipar < ParticleState::dimension(); ipar++){ + if(fabs(state.state()[ipar] - xdbstate.state()[ipar])>1.0e-6) cout << "Exact State vector " << ipar << " doesn't match: dir " + << idir << " original " << state.state()[ipar] << " changed " << xdbstate.state()[ipar] << endl; + } + // 1st order state test + auto dbstate = dbtraj.stateEstimate(ttest); + for(size_t ipar=0;ipar < ParticleState::dimension(); ipar++){ + if(fabs(state.state()[ipar] - dbstate.state()[ipar])>1.0e-6) cout << "1st order State vector " << ipar << " doesn't match: dir " + << idir << " original " << state.state()[ipar] << " changed " << dbstate.state()[ipar] << endl; + } + } } char gtitle[80]; @@ -504,6 +514,33 @@ int test(int argc, char **argv) { dbcan[idir]->Draw(); dbcan[idir]->Write(); } + +// test covariance rotation +// double db = 0.001; +// // PSMAT dpdbdp = traj.dPardBdPar(traj.t0()); +// KTRAJ dbtraj(traj); +//// DVEC dpdb = traj.dPardB(traj.t0()); +// auto bnomp = (1.0+db)*traj.bnom(); +// dbtraj.setBNom(traj.t0(),bnomp); +//// dbtraj.resetBNom(bnomp); +//// dbtraj.params().parameters() += dpdb*db/(1.0+db); +// auto dbpstate = dbtraj.stateEstimate(traj.t0()); +// for(size_t ipar=0; ipar < NParams(); ipar++){ +// auto dp = fabs(pstate.state()[ipar]-dbpstate.state()[ipar])/db; +// if(dp > 1.0e-8){ +// std::cout << "State mismatch, par " << ipar << " diff " << dp << " pars " << pstate.state()[ipar] << " " << dbpstate.state()[ipar] << std::endl; +// retval = -1; +// } +// for(size_t jpar=0; jpar < NParams(); jpar++){ +// auto dc = fabs(pstate.stateCovariance()[ipar][jpar]-dbpstate.stateCovariance()[ipar][jpar])/sqrt(pstate.stateCovariance()[ipar][ipar]*pstate.stateCovariance()[jpar][jpar]); +// if( dc > 1.0e2) { // temporarily disable FIXME! +// std::cout << "State Covariance mismatch par " << ipar << " , " << jpar << " diff " << dc << " covs " << pstate.stateCovariance()[ipar][jpar] << " " << dbpstate.stateCovariance()[ipar][jpar] << std::endl; +// retval = -1; +// } +// } +// } +// } + TCanvas* dbgcan = new TCanvas("dbgcan","DB Gap",800,800); dbgcan->Divide(2,2); for(int idir=0;idir<3;++idir){ diff --git a/Tests/FitTest.hh b/Tests/FitTest.hh index 16ba0b58..03743127 100644 --- a/Tests/FitTest.hh +++ b/Tests/FitTest.hh @@ -3,7 +3,6 @@ // #include "KinKal/General/Vectors.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/Trajectory/Line.hh" #include "KinKal/Examples/SimpleWireHit.hh" #include "KinKal/Examples/ScintHit.hh" #include "KinKal/Detector/StrawXing.hh" @@ -14,12 +13,12 @@ #include "KinKal/Fit/Config.hh" #include "KinKal/Fit/Measurement.hh" #include "KinKal/Fit/Material.hh" -#include "KinKal/Fit/BField.hh" +#include "KinKal/Fit/DomainWall.hh" #include "KinKal/Fit/Track.hh" #include "KinKal/Tests/ToyMC.hh" #include "KinKal/Examples/HitInfo.hh" #include "KinKal/Examples/MaterialInfo.hh" -#include "KinKal/Examples/BFieldInfo.hh" +#include "KinKal/Examples/DomainWallInfo.hh" #include "KinKal/Examples/ParticleTrajectoryInfo.hh" #include "KinKal/Examples/DOCAWireHitUpdater.hh" #include "KinKal/General/PhysicalConstants.h" @@ -66,7 +65,6 @@ using namespace MatEnv; using namespace KinKal; using namespace std; // avoid confusion with root -using KinKal::Line; void print_usage() { printf("Usage: FitTest --momentum f --simparticle i --fitparticle i--charge i --nhits i --hres f --seed i -ambigdoca f --nevents i --simmat i--fitmat i --ttree i --Bz f --dBx f --dBy f --dBz f--Bgrad f --tolerance f --TFilesuffix c --PrintBad i --PrintDetail i --ScintHit i --invert i --Schedule a --ssmear i --constrainpar i --inefficiency f --extend s --TimeBuffer f --matvarscale i\n"); } @@ -95,6 +93,60 @@ double dTraj(KTRAJ const& kt1, KTRAJ const& kt2, double t1, double& t2) { return ((pos2-pos1).Cross(dir1)).R(); } +template +int testState(KinKal::Track const& kktrk,DVEC sigmas) { + // test parameterstate + int retval(0); + auto const& traj = kktrk.fitTraj().nearestPiece(kktrk.fitTraj().range().mid()); + auto pstate = traj.stateEstimate(traj.t0()); + auto pos1 = traj.position3(traj.t0()); + auto pos2 = pstate.position3(); + if(fabs(pos1.Dot(pos2) - pos1.Mag2()) >1e-10){ + std::cout << "Position difference " << pos1 << " " << pos2 << std::endl; + retval = -1; + } + auto mom1 = traj.momentum3(traj.t0()); + auto mom2 = pstate.momentum3(); + if(fabs(mom1.Dot(mom2)-mom1.Mag2())>1e-10){ + std::cout << "Momentum difference " << mom1 << " " << mom2 << std::endl; + retval = -1; + } + double momvar1 = traj.momentumVariance(traj.t0()); + double momvar2 = pstate.momentumVariance(); + if(fabs(momvar1-momvar2)>1e-10){ + std::cout << "Momentum variance differencer " << momvar1 << " " << momvar2 << std::endl; + retval = -1; + } + // test reversibility + KTRAJ testtraj(pstate,traj.bnom(),traj.range()); + testtraj.syncPhi0(traj); + + for(size_t ipar=0; ipar < NParams(); ipar++){ + auto dp = fabs(traj.paramVal(ipar)-testtraj.paramVal(ipar))/sigmas(ipar); + if( dp > 1.0e-10){ + if(ipar == KTRAJ::phi0Index()){ + if(fabs(traj.paramVal(ipar)-testtraj.paramVal(ipar))-2*M_PI > 1.0e-8){ + std::cout << "Parameter mismatch, par " << ipar << " diff " << traj.paramVal(ipar) << " " << testtraj.paramVal(ipar) << std::endl; + retval = -1; + } + // allow phi0 to be off by 2pi + } else { + std::cout << "Parameter mismatch, par " << ipar << " diff " << traj.paramVal(ipar) << " " << testtraj.paramVal(ipar) << std::endl; + retval = -1; + } + } + for(size_t jpar=0; jpar < NParams(); jpar++){ + auto dc = fabs(traj.params().covariance()(ipar,jpar)-testtraj.params().covariance()(ipar,jpar))/sigmas(ipar)*sigmas(jpar); + if(dc > 1.0e-8){ + std::cout << "Covariance mismatch par " << ipar << " , " << jpar << " diff " << traj.params().covariance()(ipar,jpar) << " " << testtraj.params().covariance()(ipar,jpar) << std::endl; + retval = -1; + } + } + } + + return retval; +} + int makeConfig(string const& cfile, KinKal::Config& config,bool mvarscale=true) { string fullfile; if(strncmp(cfile.c_str(),"/",1) == 0) { @@ -167,7 +219,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { using KKEFF = Effect; using KKMEAS = Measurement; using KKMAT = Material; - using KKBFIELD = BField; + using KKDW = DomainWall; using PTRAJ = ParticleTrajectory; using MEAS = Hit; using MEASPTR = std::shared_ptr; @@ -213,7 +265,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { unsigned nhits(40); unsigned nsteps(200); // steps for traj comparison double seedsmear(10.0); - double momsigma(0.2); + double momsigma(0.3); double ineff(0.05); bool simmat(true), scinthit(true); int retval(EXIT_SUCCESS); @@ -420,35 +472,14 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { KTRAJPars ftpars_, mtpars_, btpars_, spars_, ffitpars_, ffiterrs_, mfitpars_, mfiterrs_, bfitpars_, bfiterrs_; float chisq_, btmom_, mtmom_, ftmom_, ffmom_, mfmom_, bfmom_, ffmomerr_, mfmomerr_, bfmomerr_, chiprob_; float fft_,mft_, bft_; - int ndof_, niter_, status_, igap_, nmeta_, nkkbf_, nkkhit_, nkkmat_; + int ndof_, niter_, status_, igap_, nmeta_, nkkdw_, nkkhit_, nkkmat_; int nactivehit_, nstrawhit_, nscinthit_, nnull_; float sbeg_, send_, fbeg_, fend_; float maxgap_, avgap_; + int ts = testState(kktrk,sigmas); +// if(ts != 0)return ts; + if(ts != 0)std::cout << "state test failed" << std::endl; - // test parameterstate - auto const& traj = kktrk.fitTraj().front(); - auto pstate = traj.stateEstimate(traj.t0()); - double momvar1 = traj.momentumVariance(traj.t0()); - double momvar2 = pstate.momentumVariance(); - if(fabs(momvar1-momvar2)>1e-10){ - std::cout << "Momentum variance error " << momvar1 << " " << momvar2 << std::endl; - return -3; - } - // full reversibility - KTRAJ testtraj(pstate,traj.bnom(),traj.range()); - for(size_t ipar=0; ipar < NParams(); ipar++){ - if(fabs(traj.paramVal(ipar)-testtraj.paramVal(ipar)) > 1.0e-10){ - std::cout << "Parameter error " << traj.paramVal(ipar) << " " << testtraj.paramVal(ipar) << std::endl; - return -3; - } - for(size_t jpar=0; jpar < NParams(); jpar++){ - if(fabs(traj.params().covariance()(ipar,jpar)-testtraj.params().covariance()(ipar,jpar)) > 1.0e-6){ - std::cout << "Covariance error " << traj.paramVal(ipar) << " " << testtraj.paramVal(ipar) << std::endl; - return -3; - } - } - } - std::cout << "Passed ParameterState tests" << std::endl; if(nevents ==0 ){ // draw the fit result TCanvas* pttcan = new TCanvas("pttcan","PieceKTRAJ",1000,1000); @@ -487,8 +518,8 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { SCINTHITPTR lhptr = std::dynamic_pointer_cast (thit); if(shptr.use_count() > 0){ auto const& tline = shptr->wire(); - plow = tline.startPosition(); - phigh = tline.endPosition(); + plow = tline.start(); + phigh = tline.end(); line->SetLineColor(kRed); auto hitpos = tline.position3(shptr->closestApproach().sensorToca()); auto trkpos = fptraj.position3(shptr->closestApproach().particleToca()); @@ -498,8 +529,8 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { tpos->SetMarkerColor(kGreen); } else if (lhptr.use_count() > 0){ auto const& tline = lhptr->sensorAxis(); - plow = tline.startPosition(); - phigh = tline.endPosition(); + plow = tline.start(); + phigh = tline.end(); line->SetLineColor(kCyan); auto hitpos = tline.position3(lhptr->closestApproach().sensorToca()); auto trkpos = fptraj.position3(lhptr->closestApproach().particleToca()); @@ -547,7 +578,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { ftree->Branch("bferrs.", &bfiterrs_,KTRAJPars::leafnames().c_str()); ftree->Branch("chisq", &chisq_,"chisq/F"); ftree->Branch("ndof", &ndof_,"ndof/I"); - ftree->Branch("nkkbf", &nkkbf_,"nkkbf/I"); + ftree->Branch("nkkdw", &nkkdw_,"nkkdw/I"); ftree->Branch("nkkmat", &nkkmat_,"nkkmat/I"); ftree->Branch("nkkhit", &nkkhit_,"nkkhit/I"); ftree->Branch("nactivehit", &nactivehit_,"nactivehit/I"); @@ -714,7 +745,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { maxgap_ = avgap_ = -1; igap_ = -1; // fill effect information - nkkbf_ = 0; nkkhit_ = 0; nkkmat_ = 0; + nkkdw_ = 0; nkkhit_ = 0; nkkmat_ = 0; // accumulate chisquared info chisq_ = fstat.chisq_.chisq(); ndof_ = fstat.chisq_.nDOF(); @@ -756,8 +787,17 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { ffmomerr_ = -1.0; mfmomerr_ = -1.0; bfmomerr_ = -1.0; + // gaps + double maxgap, avgap; + size_t igap; + kktrk.fitTraj().gaps(maxgap, igap, avgap); + maxgap_ = maxgap; + avgap_ = avgap; + igap_ = igap; if(fstat.usable()){ + int ts = testState(kktrk,sigmas); + if(ts != 0)std::cout << "state test failed" << std::endl; // basic info auto const& fptraj = kktrk.fitTraj(); // compare parameters at the first traj of both true and fit @@ -819,11 +859,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { fmompull->Fill((ffmom_-ftmom_)/ffmomerr_); mmompull->Fill((mfmom_-mtmom_)/mfmomerr_); bmompull->Fill((bfmom_-btmom_)/bfmomerr_); - // state space parameter difference and errors - // ParticleStateEstimate tslow = tptraj.state(tlow); - // ParticleStateEstimate tshigh = tptraj.state(thigh); - // ParticleStateEstimate slow = fptraj.stateEstimate(tlow); - // ParticleStateEstimate shigh = fptraj.stateEstimate(thigh); + // if(ttree && fstat.usable()){ fbeg_ = fptraj.range().begin(); fend_ = fptraj.range().end(); @@ -831,7 +867,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { nactivehit_ = nstrawhit_ = nnull_ = nscinthit_ = 0; for(auto const& eff: kktrk.effects()) { const KKMEAS* kkhit = dynamic_cast(eff.get()); - const KKBFIELD* kkbf = dynamic_cast(eff.get()); + const KKDW* kkdw = dynamic_cast(eff.get()); const KKMAT* kkmat = dynamic_cast(eff.get()); if(kkhit != 0){ nkkhit_++; @@ -914,7 +950,7 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { minfo.active_ = kkmat->active(); minfo.nxing_ = kkmat->elementXing().matXings().size(); std::array dmom = {0.0,0.0,0.0}, momvar = {0.0,0.0,0.0}; - kkmat->elementXing().materialEffects(TimeDir::forwards, dmom, momvar); + kkmat->elementXing().materialEffects(dmom, momvar); minfo.dmomf_ = dmom[MomBasis::momdir_]; minfo.momvar_ = momvar[MomBasis::momdir_]; minfo.perpvar_ = momvar[MomBasis::perpdir_]; @@ -926,12 +962,12 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { } minfovec.push_back(minfo); } - if(kkbf != 0){ - nkkbf_++; - BFieldInfo bfinfo; - bfinfo.active_ = kkbf->active(); - bfinfo.time_ = kkbf->time(); - bfinfo.range_ = kkbf->range().range(); + if(kkdw != 0){ + nkkdw_++; + DomainWallInfo bfinfo; + bfinfo.active_ = kkdw->active(); + bfinfo.time_ = kkdw->time(); + bfinfo.range_ = kkdw->nextDomain().mid()-kkdw->prevDomain().mid(); bfinfovec.push_back(bfinfo); } } @@ -961,12 +997,6 @@ int FitTest(int argc, char *argv[],KinKal::DVEC const& sigmas) { ktinfo.dt_= tstep-ttrue; tinfovec.push_back(ktinfo); } - double maxgap, avgap; - size_t igap; - fptraj.gaps(maxgap, igap, avgap); - maxgap_ = maxgap; - avgap_ = avgap; - igap_ = igap; } } else if(printbad){ cout << "Bad Fit event " << ievent << " status " << kktrk.fitStatus() << endl; diff --git a/Tests/HitTest.hh b/Tests/HitTest.hh index 7d8d97a0..1fd67b56 100644 --- a/Tests/HitTest.hh +++ b/Tests/HitTest.hh @@ -3,7 +3,6 @@ // #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/Trajectory/LoopHelix.hh" -#include "KinKal/Trajectory/Line.hh" #include "KinKal/Trajectory/ClosestApproach.hh" #include "KinKal/Examples/SimpleWireHit.hh" #include "KinKal/Examples/ScintHit.hh" @@ -41,11 +40,9 @@ using namespace MatEnv; using namespace KinKal; using namespace std; -// avoid confusion with root -using KinKal::Line; void print_usage() { - printf("Usage: HitTest --momentum f --particle i --charge i --strawhit i --scinthit i --zrange f --nhits i --hres f --seed i --ambigdoca f --By f --Bgrad f --simmat_ i --prec f\n"); + printf("Usage: HitTest --momentum f --particle i --charge i --strawhit i --scinthit i --zrange f --nhits i --hres f --seed i --ambigdoca f --By f --Bgrad f --simmat_ i --prec f --maxdr f \n"); } template @@ -77,6 +74,7 @@ int HitTest(int argc, char **argv, const vector& delpars) { double Bgrad(0.0), By(0.0); bool simmat_(true), scinthit_(true), strawhit_(true); double zrange(3000.0); // tracker dimension + double maxdr(1.0); static struct option long_options[] = { {"momentum", required_argument, 0, 'm' }, @@ -93,6 +91,7 @@ int HitTest(int argc, char **argv, const vector& delpars) { {"By", required_argument, 0, 'y' }, {"Bgrad", required_argument, 0, 'g' }, {"prec", required_argument, 0, 'P' }, + {"maxdr", required_argument, 0, 'M' }, {NULL, 0,0,0} }; @@ -104,6 +103,8 @@ int HitTest(int argc, char **argv, const vector& delpars) { break; case 'p' : imass = atoi(optarg); break; + case 'M' : maxdr = atof(optarg); + break; case 'q' : icharge = atoi(optarg); break; case 'z' : zrange = atof(optarg); @@ -190,13 +191,13 @@ int HitTest(int argc, char **argv, const vector& delpars) { SCINTHITPTR lhptr = std::dynamic_pointer_cast (thit); if((bool)shptr){ auto const& tline = shptr->wire(); - plow = tline.startPosition(); - phigh = tline.endPosition(); + plow = tline.start(); + phigh = tline.end(); line->SetLineColor(kRed); } else if ((bool)lhptr){ auto const& tline = lhptr->sensorAxis(); - plow = tline.startPosition(); - phigh = tline.endPosition(); + plow = tline.start(); + phigh = tline.end(); line->SetLineColor(kCyan); } line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); @@ -285,7 +286,7 @@ int HitTest(int argc, char **argv, const vector& delpars) { // compare the change with the expected from the derivatives double ddr = ROOT::Math::Dot(pder,dpvec); hderivg[ipar]->SetPoint(ipt++,dr,ddr); - if(fabs(dr - ddr) > 1.0 ){ + if(fabs(dr - ddr) > maxdr ){ cout << "Large ddiff " << KTRAJ::paramName(tpar) << " " << *thit << " delta " << dpar << " doca " << tpdata.doca() << " DirDot " << tpdata.dirDot() <<" Exact change " << dr << " deriv " << ddr << endl; status = 2; diff --git a/Tests/LoopHelixBField_unit.cc b/Tests/LoopHelixBField_unit.cc index dde979d5..b0d5d387 100644 --- a/Tests/LoopHelixBField_unit.cc +++ b/Tests/LoopHelixBField_unit.cc @@ -1,5 +1,20 @@ #include "KinKal/Trajectory/LoopHelix.hh" #include "KinKal/Tests/BFieldMapTest.hh" int main(int argc, char **argv) { - return BFieldMapTest(argc,argv); + if(argc == 1){ + cout << "Testing gradient field, correction 2" << endl; + std::vector arguments; + arguments.push_back(argv[0]); + arguments.push_back("--Bgrad"); + arguments.push_back("-0.036"); // mu2e-like field gradient + arguments.push_back("--Tol"); + arguments.push_back("1.0e-4"); + arguments.push_back("1"); + std::vector myargv; + for (const auto& arg : arguments) + myargv.push_back((char*)arg.data()); + myargv.push_back(nullptr); + return BFieldMapTest(myargv.size()-1,myargv.data()); + } else + return BFieldMapTest(argc,argv); } diff --git a/Tests/LoopHelixFit_unit.cc b/Tests/LoopHelixFit_unit.cc index ff00ba6d..a63acd24 100644 --- a/Tests/LoopHelixFit_unit.cc +++ b/Tests/LoopHelixFit_unit.cc @@ -6,17 +6,20 @@ int main(int argc, char **argv) { cout << "Testing gradient field, correction 2" << endl; std::vector arguments; arguments.push_back(argv[0]); - arguments.push_back("--Bgrad"); - arguments.push_back("-0.036"); // mu2e-like field gradient +// arguments.push_back("--Bgrad"); +// arguments.push_back("-0.036"); // mu2e-like field gradient arguments.push_back("--Schedule"); - arguments.push_back("driftfit.txt"); - arguments.push_back("--ttree"); - arguments.push_back("1"); + arguments.push_back("MCAmbigFixedField.txt"); +// arguments.push_back("seedfit.txt"); +// arguments.push_back("--extend"); +// arguments.push_back("driftextend.txt"); +// arguments.push_back("--ttree"); +// arguments.push_back("1"); std::vector myargv; for (const auto& arg : arguments) myargv.push_back((char*)arg.data()); myargv.push_back(nullptr); return FitTest(myargv.size()-1,myargv.data(),sigmas); } else - return FitTest(argc,argv,sigmas); + return FitTest(argc,argv,sigmas); } diff --git a/Tests/MCAmbig.txt b/Tests/MCAmbig.txt index 5211fa32..e64917c5 100644 --- a/Tests/MCAmbig.txt +++ b/Tests/MCAmbig.txt @@ -1,7 +1,7 @@ # # Test Configuration file for iteration schedule -# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tbuff tol minndof bfcor plevel -10 1.0e6 1.0 1.0e8 1.0e8 100.0 1.0e-4 5 0 0 0 +# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tol minndof bfcor ends plevel +10 1.0e6 1.0 50.0 1.0e6 100.0 1.0e-4 5 1 0 0 # Then, meta-iteration specific parameters: 2.0 1.0 diff --git a/Tests/MCAmbigFixedField.txt b/Tests/MCAmbigFixedField.txt new file mode 100644 index 00000000..5211fa32 --- /dev/null +++ b/Tests/MCAmbigFixedField.txt @@ -0,0 +1,8 @@ +# +# Test Configuration file for iteration schedule +# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tbuff tol minndof bfcor plevel +10 1.0e6 1.0 1.0e8 1.0e8 100.0 1.0e-4 5 0 0 0 +# Then, meta-iteration specific parameters: +2.0 +1.0 +0.0 diff --git a/Tests/ParticleTrajectoryTest.hh b/Tests/ParticleTrajectoryTest.hh index 8496da50..202e3370 100644 --- a/Tests/ParticleTrajectoryTest.hh +++ b/Tests/ParticleTrajectoryTest.hh @@ -2,7 +2,7 @@ // test basic functions of the ParticleTrajectory, using the LoopHelix class // #include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" #include "KinKal/General/BFieldMap.hh" #include "KinKal/General/PhysicalConstants.h" @@ -37,7 +37,7 @@ using namespace KinKal; using namespace std; // avoid confusion with root -using KinKal::Line; +using KinKal::SensorLine; void print_usage() { printf("Usage: ParticleTrajectoryTest --changedir i --delta f --tstep f --nsteps i \n"); @@ -46,7 +46,7 @@ void print_usage() { template int ParticleTrajectoryTest(int argc, char **argv) { using PTRAJ = ParticleTrajectory; - using PCA = PiecewiseClosestApproach; + using PCA = PiecewiseClosestApproach; double mom(105.0), cost(0.7), phi(0.5); unsigned npts(50); int icharge(-1); @@ -225,7 +225,7 @@ int ParticleTrajectoryTest(int argc, char **argv) { // shift the position VEC3 lpos = midpos + gap*rdir; double wlen(1000.0); - Line tline(lpos, ptraj.range().mid(), pvel, wlen); + SensorLine tline(lpos, ptraj.range().mid(), pvel, wlen); // create ClosestApproach from these CAHint tphint(ptraj.range().mid(),0.0); PCA pca(ptraj,tline, tphint,1e-8); @@ -240,8 +240,8 @@ int ParticleTrajectoryTest(int argc, char **argv) { // draw the line and ClosestApproach TPolyLine3D* line = new TPolyLine3D(2); VEC3 plow, phigh; - plow = tline.startPosition(); - phigh = tline.endPosition(); + plow = tline.start(); + phigh = tline.end(); line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); line->SetPoint(1,phigh.X(),phigh.Y(), phigh.Z()); line->SetLineColor(kOrange); diff --git a/Tests/ToyMC.hh b/Tests/ToyMC.hh index ab7f1aa9..98afb659 100644 --- a/Tests/ToyMC.hh +++ b/Tests/ToyMC.hh @@ -9,7 +9,7 @@ #include "KinKal/MatEnv/DetMaterial.hh" #include "KinKal/MatEnv/SimpleFileFinder.hh" #include "KinKal/Examples/CaloDistanceToTime.hh" -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/Trajectory/PiecewiseClosestApproach.hh" #include "KinKal/Examples/SimpleWireHit.hh" @@ -37,7 +37,7 @@ namespace KKTest { using SCINTHITPTR = std::shared_ptr; using STRAWXING = StrawXing; using STRAWXINGPTR = std::shared_ptr; - using PCA = PiecewiseClosestApproach; + using PCA = PiecewiseClosestApproach; // create from aseed ToyMC(BFieldMap const& bfield, double mom, int icharge, double zrange, int iseed, unsigned nhits, bool simmat, bool scinthit, double ambigdoca ,double simmass) : bfield_(bfield), matdb_(sfinder_,MatEnv::DetMaterial::moyalmean), // use the moyal based eloss model @@ -52,7 +52,7 @@ namespace KKTest { } // generate a straw at the given time. direction and drift distance are random - Line generateStraw(PTRAJ const& traj, double htime); + SensorLine generateStraw(PTRAJ const& traj, double htime); // create a seed by randomizing the parameters void createSeed(KTRAJ& seed,DVEC const& sigmas, double seedsmear); void extendTraj(PTRAJ& ptraj,double htime); @@ -99,19 +99,19 @@ namespace KKTest { }; - template Line ToyMC::generateStraw(PTRAJ const& traj, double htime) { + template SensorLine ToyMC::generateStraw(PTRAJ const& traj, double htime) { // start with the true helix position at this time auto hpos = traj.position4(htime); - auto hdir = traj.direction(htime); - // generate a random direction for the straw - double eta = tr_.Uniform(-M_PI,M_PI); - VEC3 sdir(cos(eta),sin(eta),0.0); + auto tdir = traj.direction(htime); + // generate a random azimuth direction for the straw + double azimuth = tr_.Uniform(-M_PI,M_PI); + VEC3 sdir(cos(azimuth),sin(azimuth),0.0); // generate a random drift perp to this and the trajectory double rdrift = tr_.Uniform(-rstraw_,rstraw_); - VEC3 drift = (sdir.Cross(hdir)).Unit(); - VEC3 dpos = hpos.Vect() + rdrift*drift; + VEC3 driftdir = (sdir.Cross(tdir)).Unit(); + VEC3 dpos = hpos.Vect() + rdrift*driftdir; // cout << "Generating hit at position " << dpos << endl; - double dprop = tr_.Uniform(0.0,0.5*wlen_); + double dprop = tr_.Uniform(0.0,wlen_); VEC3 mpos = dpos + sdir*dprop; VEC3 vprop = sdir*sprop_; // measured time is after propagation and drift @@ -119,7 +119,7 @@ namespace KKTest { // smear measurement time tmeas = tr_.Gaus(tmeas,sigt_); // construct the trajectory for this hit - return Line(mpos,tmeas,vprop,wlen_); + return SensorLine(mpos,tmeas,vprop,wlen_); } template void ToyMC::simulateParticle(PTRAJ& ptraj,HITCOL& thits, EXINGCOL& dxings, bool addmat) { @@ -175,7 +175,7 @@ namespace KKTest { auto endmom = endpiece.momentum4(tstraw); auto endpos = endpiece.position4(tstraw); std::array dmom {0.0,0.0,0.0}, momvar {0.0,0.0,0.0}; - sxing->materialEffects(TimeDir::forwards, dmom, momvar); + sxing->materialEffects(dmom, momvar); for(int idir=0;idir<=MomBasis::phidir_; idir++) { auto mdir = static_cast(idir); double momsig = sqrt(momvar[idir]); @@ -206,8 +206,8 @@ namespace KKTest { template void ToyMC::createScintHit(PTRAJ& ptraj, HITCOL& thits) { // create a ScintHit at the end, axis parallel to z - // first, find the position at showermax_. - VEC3 shmaxTrue,shmaxMeas; + // first, find the position at showermax. + VEC3 shmax,shmeas; double tend = thits.back()->time(); // extend to the calorimeter z VEC3 pvel = ptraj.velocity(tend); @@ -218,27 +218,17 @@ namespace KKTest { pvel = ptraj.velocity(shstart); // compute time at showermax double shmaxtime = shstart + shmax_/pvel.R(); - auto endpos = ptraj.position4(shstart); - shmaxTrue = ptraj.position3(shmaxtime); // true position at shower-max - // smear the x-y position by the transverse variance. - shmaxMeas.SetX(tr_.Gaus(shmaxTrue.X(),shPosSig_)); - shmaxMeas.SetY(tr_.Gaus(shmaxTrue.Y(),shPosSig_)); - // set the z position to the sensor plane (end of the crystal) - shmaxMeas.SetZ(endpos.Z()+clen_); - // set the measurement time to correspond to the light propagation from showermax_, smeared by the resolution - //double tmeas = tr_.Gaus(shmaxtime+(shmaxMeas.Z()-shmaxTrue.Z())/cprop_,scitsig_); + shmax = ptraj.position3(shmaxtime); // true position at shower-max + // Compute measurement position: smear the x-y position by the transverse variance. + shmeas.SetX(tr_.Gaus(shmax.X(),shPosSig_)); + shmeas.SetY(tr_.Gaus(shmax.Y(),shPosSig_)); + // set the z position to the sensor plane (forward end of the crystal) + shmeas.SetZ(cz_+clen_); + // set the measurement time to correspond to the light propagation from showermax_, smeared by the time resolution + double tmeas = tr_.Gaus(shmaxtime+(shmeas.Z() - shmax.Z())/cprop_,scitsig_); // create the ttraj for the light propagation VEC3 lvel(0.0,0.0,cprop_); - - // put in manual values - //std::shared_ptr calod2t = std::make_shared(85.76, clen_-27.47); - std::shared_ptr calod2t = std::make_shared(cprop_, clen_-27.47, 0.001); - //CaloDistanceToTime calod2t(tmeas, 85.76, 27.47); - double tmeas = shmaxtime + calod2t->time(shmaxMeas.Z() - shmaxTrue.Z()); - Line lline(shmaxMeas, tmeas, lvel, clen_, calod2t); - - // original - //Line lline(shmaxMeas,tmeas,lvel,clen_); + SensorLine lline(shmeas, tmeas, lvel, clen_); // then create the hit and add it; the hit has no material CAHint tphint(tmeas,tmeas); PCA pca(ptraj,lline,tphint,tprec_); @@ -257,7 +247,6 @@ namespace KKTest { template void ToyMC::extendTraj(PTRAJ& ptraj,double htime) { if(htime > ptraj.range().end()){ - ROOT::Math::SMatrix bgrad; VEC3 pos,vel, dBdt; pos = ptraj.position3(htime); vel = ptraj.velocity(htime); @@ -266,13 +255,15 @@ namespace KKTest { if(dBdt.R() != 0.0){ double tbeg = ptraj.range().end(); while(tbeg < htime) { - double tend = bfield_.rangeInTolerance(ptraj.back(),tbeg,tol_); + double tend = std::min(tbeg + bfield_.rangeInTolerance(ptraj.back(),tbeg,tol_),htime); double tmid = 0.5*(tbeg+tend); auto bf = bfield_.fieldVect(ptraj.position3(tmid)); - auto pos = ptraj.position4(tmid); - auto mom = ptraj.momentum4(tmid); + auto pos = ptraj.position4(tend); + auto mom = ptraj.momentum4(tend); TimeRange prange(tbeg,tend); KTRAJ newend(pos,mom,ptraj.charge(),bf,prange); + // make sure phi0 stays continuous + newend.syncPhi0(ptraj.back()); ptraj.append(newend); tbeg = tend; } @@ -287,10 +278,14 @@ namespace KKTest { double tsint = sqrt(1.0-tcost*tcost); MOM4 tmomv(mom_*tsint*cos(tphi),mom_*tsint*sin(tphi),mom_*tcost,simmass_); double tmax = fabs(zrange_/(CLHEP::c_light*tcost)); - VEC4 torigin(tr_.Gaus(0.0,osig_), tr_.Gaus(0.0,osig_), tr_.Gaus(-0.5*zrange_,osig_),tr_.Uniform(t0off_-tmax,t0off_+tmax)); + double tbeg = tr_.Uniform(t0off_-tmax,t0off_+tmax); + VEC4 torigin(tr_.Gaus(0.0,osig_), tr_.Gaus(0.0,osig_), tr_.Gaus(-0.5*zrange_,osig_),tbeg); VEC3 bsim = bfield_.fieldVect(torigin.Vect()); - KTRAJ ktraj(torigin,tmomv,icharge_,bsim,TimeRange(torigin.T(),torigin.T()+tmax)); - ptraj = PTRAJ(ktraj); + KTRAJ ktraj(torigin,tmomv,icharge_,bsim,TimeRange(tbeg,tbeg+tmax)); + // set initial range + double tend = tbeg + bfield_.rangeInTolerance(ktraj,tbeg,tol_); + ktraj.setRange(TimeRange(tbeg,tend)); + ptraj.append(ktraj); } } #endif diff --git a/Tests/Trajectory.hh b/Tests/Trajectory.hh index 4172b162..62a6f3b3 100644 --- a/Tests/Trajectory.hh +++ b/Tests/Trajectory.hh @@ -1,7 +1,7 @@ // // test basic functions of kinematic trajectory class // -#include "KinKal/Trajectory/Line.hh" +#include "KinKal/Trajectory/SensorLine.hh" #include "KinKal/Trajectory/ClosestApproach.hh" #include "KinKal/General/ParticleStateEstimate.hh" #include "KinKal/General/PhysicalConstants.h" @@ -30,8 +30,6 @@ using namespace KinKal; using namespace std; -// avoid confusion with root -using KinKal::Line; void print_usage() { printf("Usage: Trajectory --momentum f --costheta f --azimuth f --particle i --charge i --xorigin f -- yorigin f --zorigin f --torigin --tmin f--tmax f --ltime f --By f --invert i\n"); @@ -134,6 +132,17 @@ int TrajectoryTest(int argc, char **argv,KinKal::DVEC sigmas) { double sint = sqrt(1.0-cost*cost); MOM4 momv(mom*sint*cos(phi),mom*sint*sin(phi),mom*cost,pmass); KTRAJ ktraj(origin,momv,icharge,bnom,TimeRange(-10,10)); + +// test state + auto state = ktraj.state(3.0); + KTRAJ straj(state,bnom,ktraj.range()); + for(size_t ipar=0; ipar < NParams(); ipar++){ + if(fabs(ktraj.paramVal(ipar)-straj.paramVal(ipar)) > 1.0e-10){ + std::cout << "Parameter mismatch, par " << ipar << " diff " << ktraj.paramVal(ipar) << " " << straj.paramVal(ipar) << std::endl; + return -1; + } + } + if(invert)ktraj.invertCT(); auto testmom = ktraj.momentum4(ot); // cout << "KTRAJ with momentum " << momv.Vect() << " position " << origin << " has parameters: " << ktraj << endl; @@ -267,16 +276,16 @@ int TrajectoryTest(int argc, char **argv,KinKal::DVEC sigmas) { // shift the position VEC3 perpdir(-sin(phi),cos(phi),0.0); VEC3 ppos = pos + gap*perpdir; - Line tline(ppos, ltime, pvel, wlen); + SensorLine tline(ppos, ltime, pvel, wlen); // find ClosestApproach CAHint hint(ltime,ltime); - ClosestApproach tp(ktraj,tline,hint, 1e-6); + ClosestApproach tp(ktraj,tline,hint, 1e-6); // cout << "ClosestApproach status " << tp.statusName() << " doca " << tp.doca() << " dt " << tp.deltaT() << endl; if(tp.status() == ClosestApproachData::converged) { // draw the line and ClosestApproach TPolyLine3D* line = new TPolyLine3D(2); - auto plow = tline.startPosition(); - auto phigh = tline.endPosition(); + auto plow = tline.start(); + auto phigh = tline.end(); line->SetPoint(0,plow.X(),plow.Y(), plow.Z()); line->SetPoint(1,phigh.X(),phigh.Y(), phigh.Z()); line->SetLineColor(kOrange); @@ -339,9 +348,9 @@ int TrajectoryTest(int argc, char **argv,KinKal::DVEC sigmas) { // test axis auto axis = ktraj.axis(ltime); auto bdir = ktraj.bnom().Unit(); - auto rtest = (axis.start_-ktraj.position3(ltime)).R(); - if( fabs(axis.dir_.Dot(acc)) > 1e-9 || fabs(rtest-ktraj.bendRadius()) > 1e-9 || - (acc.R() != 0 && fabs(fabs(axis.dir_.Dot(bdir))-1.0)>1e-9) ){ + auto rtest = (axis.start()-ktraj.position3(ltime)).R(); + if( fabs(axis.direction().Dot(acc)) > 1e-9 || fabs(rtest-ktraj.bendRadius()) > 1e-9 || + (acc.R() != 0 && fabs(fabs(axis.direction().Dot(bdir))-1.0)>1e-9) ){ cout << "Axis check failed " << endl; return -2; } diff --git a/Tests/debug.txt b/Tests/debug.txt new file mode 100644 index 00000000..529abef1 --- /dev/null +++ b/Tests/debug.txt @@ -0,0 +1,5 @@ +# Test Configuration iteration schedule for a null (no drift) fit +# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tol minndof bfcor ends plevel +10 1.0e6 1.0 50.0 1.0e6 100.0 1e-2 5 1 0 0 +# Then, meta-iteration specific parameters: temperature and update type (null hit) +2.0 0 diff --git a/Tests/driftextend.txt b/Tests/driftextend.txt index db02a8f0..cd39dbb5 100644 --- a/Tests/driftextend.txt +++ b/Tests/driftextend.txt @@ -1,7 +1,7 @@ # # Test Configuration iteration schedule for a drift extension -# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge tol minndof bfcor ends plevel -10 1.0e6 1.0 50.0 1.0e6 1e-4 5 1 0 0 +# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tol minndof bfcor ends plevel +10 1.0e6 1.0 50.0 1.0e6 100.0 1e-4 5 1 1 0 # Order: # temperature updater mindoca maxdoca 2.0 1 1.5 5 diff --git a/Tests/driftfit.txt b/Tests/driftfit.txt index 720c0e61..a6af9ff5 100644 --- a/Tests/driftfit.txt +++ b/Tests/driftfit.txt @@ -1,12 +1,12 @@ # # Test Configuration iteration schedule for a drift fit # first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge divgap tol minndof bfcor ends plevel -10 1.0e7 1.0 50.0 1.0e6 100.0 1e-4 5 1 1 0 +10 1.0e7 1.0 50.0 1.0e6 100.0 1e-4 5 1 0 0 # Order: # temperature updater (mindoca maxdoca) 2.0 0 1.0 0 -0.5 0 +0.0 0 2.0 1 1.5 5 1.0 1 1.0 3.5 0.5 1 0.5 2.8 diff --git a/Tests/fixedBfit.txt b/Tests/fixedBfit.txt new file mode 100644 index 00000000..6864cd33 --- /dev/null +++ b/Tests/fixedBfit.txt @@ -0,0 +1,15 @@ + +# +# Test Configuration iteration schedule for a drift fit +# first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge divgap tol minndof bfcor ends plevel +10 1.0e7 1.0 50.0 1.0e6 100.0 1e-4 5 0 0 0 +# Order: +# temperature updater (mindoca maxdoca) +2.0 0 +1.0 0 +0.5 0 +2.0 1 1.5 5 +1.0 1 1.0 3.5 +0.5 1 0.5 2.8 +0.2 1 0.5 2.8 +0.0 1 0.5 2.8 diff --git a/Tests/seedfit.txt b/Tests/seedfit.txt index 2c2df8a0..8795db05 100644 --- a/Tests/seedfit.txt +++ b/Tests/seedfit.txt @@ -1,6 +1,6 @@ # Test Configuration iteration schedule for a null (no drift) fit # first global parameters: maxniter dewight dchisquared_converge dchisquared_diverge dchisq_paramdiverge gapdiv tol minndof bfcor ends plevel -10 1.0e6 1.0 50.0 1.0e6 100.0 1e-4 5 1 0 0 +10 1.0e6 1.0 50.0 1.0e6 100.0 1e-2 5 1 0 0 # Then, meta-iteration specific parameters: temperature and update type (null hit) 2.0 0 1.0 0 diff --git a/Trajectory/CMakeLists.txt b/Trajectory/CMakeLists.txt index e719b296..d61f45ac 100644 --- a/Trajectory/CMakeLists.txt +++ b/Trajectory/CMakeLists.txt @@ -9,7 +9,6 @@ add_library(Trajectory SHARED ClosestApproachData.cc SensorLine.cc KinematicLine.cc - Line.cc LoopHelix.cc POCAUtil.cc ConstantDistanceToTime.cc diff --git a/Trajectory/CentralHelix.cc b/Trajectory/CentralHelix.cc index ff72c3af..92250991 100644 --- a/Trajectory/CentralHelix.cc +++ b/Trajectory/CentralHelix.cc @@ -1,6 +1,7 @@ #include "KinKal/Trajectory/CentralHelix.hh" #include "Math/AxisAngle.h" #include +#include #include using namespace std; @@ -34,18 +35,14 @@ namespace KinKal { // Transform into the system where Z is along the Bfield. This is a pure rotation about the origin VEC4 pos(pos0); MOM4 mom(mom0); - g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); - if(fabs(g2l_(bnom_).Theta()) > 1.0e-6)throw invalid_argument("Rotation Error"); + setTransforms(); pos = g2l_(pos); mom = g2l_(mom); - // create inverse rotation; this moves back into the original coordinate system - l2g_ = g2l_.Inverse(); // kinematic to geometric conversion double radToMom = BFieldMap::cbar()*charge*bnom_.R(); double momToRad = 1.0/radToMom; abscharge_ = abs(charge); double mbar = -mass_ * momToRad; - absmbar_ = fabs(mbar); // caches double pt = sqrt(mom.perp2()); double radius = fabs(pt*momToRad); @@ -77,51 +74,27 @@ namespace KinKal { param(z0_) = z0 - nwind*deltaz; // t0, also correcting for winding param(t0_) = pos.T() -(dphi + 2*M_PI*nwind)/Omega(); - // test -// auto testpos = position3(pos0.T()); -// auto testmom = momentum3(pos0.T()); -// auto dp = testpos - pos0.Vect(); -// auto dm = testmom - mom0.Vect(); -// if(dp.R() > 1.0e-5 || dm.R() > 1.0e-5)throw invalid_argument("Construction Test Failure"); -// // check -// auto lmom = localMomentum(pos0.T()); -// auto tcent = center(); -// if(fabs(lcent.phi()-tcent.phi())>1e-5 || fabs(lcent.perp2()-tcent.perp2()) > 1e-5){ -// cout << "center " << lcent << " test center " << tcent << endl; -// } -// if(fabs(tan(phi0()) +1.0/tan(lcent.phi())) > 1e-5){ -// cout << "phi0 " << phi0() << " test phi0 " << -1.0/tan(lcent.phi()) << endl; -// } -// double d0t = sign()*sqrt(lcent.perp2())-sqrt(lmom.perp2())/Q(); -// if(fabs(d0t - d0()) > 1e-5){ -// cout << " d0 " << d0() << " d0 test " << d0t << endl; -// } } void CentralHelix::setBNom(double time, VEC3 const& bnom) { - // adjust the parameters for the change in bnom - absmbar_ *= bnom_.R()/bnom.R(); - pars_.parameters() += dPardB(time,bnom); + // adjust the parameters for the change in bnom holding the state constant + VEC3 db = bnom-bnom_; + pars_.parameters() += dPardB(time,db); + resetBNom(bnom); + // rotate covariance TODO + } + + void CentralHelix::resetBNom(VEC3 const& bnom) { bnom_ = bnom; - // adjust rotations to global space - g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); - l2g_ = g2l_.Inverse(); + setTransforms(); } CentralHelix::CentralHelix(CentralHelix const& other, VEC3 const& bnom, double trot) : CentralHelix(other) { - absmbar_ *= bnom_.R()/bnom.R(); - bnom_ = bnom; - pars_.parameters() += other.dPardB(trot,bnom); - g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); - l2g_ = g2l_.Inverse(); + setBNom(trot,bnom); } CentralHelix::CentralHelix(Parameters const &pdata, double mass, int charge, double bnom, TimeRange const& range) : trange_(range), pars_(pdata), mass_(mass), bnom_(VEC3(0.0,0.0,bnom)){ //FIXME for now just ignore sign - // if(abscharge < 0) throw invalid_argument("Central helix charge sign should be defined by omega"); - // compute kinematic cache - double momToRad = 1.0/(BFieldMap::cbar()*charge*bnom); - absmbar_ = fabs(-mass_ * momToRad); abscharge_ = abs(charge); } @@ -140,6 +113,19 @@ namespace KinKal { pars_.covariance() = ROOT::Math::Similarity(dpds,pstate.stateCovariance()); } + void CentralHelix::syncPhi0(CentralHelix const& other) { + // adjust the phi0 of this traj to agree with the reference, keeping its value (mod 2pi) the same. + static double twopi = 2*M_PI; + int nloop = static_cast(std::round( (other.phi0() - phi0())/twopi)); + if(nloop != 0) pars_.parameters()[phi0_] += nloop*twopi; + } + + void CentralHelix::setTransforms() { + // adjust rotations to global space + g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); + l2g_ = g2l_.Inverse(); + } + double CentralHelix::momentumVariance(double time) const { DVEC dMomdP(0.0, 0.0, -1.0/omega() , 0.0 , sinDip()*cosDip() , 0.0); dMomdP *= momentum(time); @@ -437,14 +423,21 @@ namespace KinKal { retval[phi0_] = -sign()*rad*cent.Dot(mperp)/(rcval*rcval); retval[z0_] = lpos.Z()-z0() + tanDip()*retval[phi0_]/omega(); retval[t0_] = time-t0() + retval[phi0_]/Omega(); - return (1.0/bnom_.R())*retval; + return retval; } - DVEC CentralHelix::dPardB(double time, VEC3 const& BPrime) const { + PSMAT CentralHelix::dPardPardB(double time,VEC3 const& db) const { + PSMAT dpdpdb = ROOT::Math::SMatrixIdentity(); + return dpdpdb; + } + + + DVEC CentralHelix::dPardB(double time, VEC3 const& dB) const { // rotate Bfield difference into local coordinate system - VEC3 dB = g2l_(BPrime-bnom_); + VEC3 dBloc = g2l_(dB); // find the parameter change due to BField magnitude change using component parallel to the local nominal Bfield (always along z) - DVEC retval = dPardB(time)*dB.Z(); + double bfrac = dBloc.Z()/bnom_.R(); + DVEC retval = dPardB(time)*bfrac/(1.0 + bfrac); // find the change in (local) position and momentum due to the rotation implied by the B direction change // work in local coordinate system to avoid additional matrix mulitplications auto xvec = localPosition(time); diff --git a/Trajectory/CentralHelix.hh b/Trajectory/CentralHelix.hh index 1d811d1c..403aacb6 100644 --- a/Trajectory/CentralHelix.hh +++ b/Trajectory/CentralHelix.hh @@ -29,6 +29,7 @@ namespace KinKal { // classes implementing the Kalman fit // define the indices and names of the parameters enum ParamIndex {d0_=0,phi0_=1,omega_=2,z0_=3,tanDip_=4,t0_=5,npars_=6}; + constexpr static ParamIndex phi0Index() { return phi0_; } constexpr static ParamIndex t0Index() { return t0_; } static std::vector const ¶mNames(); @@ -54,6 +55,7 @@ namespace KinKal { explicit CentralHelix(ParticleState const& pstate, VEC3 const& bnom, TimeRange const& range=TimeRange()); // same, including covariance information explicit CentralHelix(ParticleStateEstimate const& pstate, VEC3 const& bnom, TimeRange const& range=TimeRange()); + void syncPhi0(CentralHelix const& other); // particle position and momentum as a function of time VEC4 position4(double time) const; VEC3 position3(double time) const; @@ -65,17 +67,19 @@ namespace KinKal { VEC3 acceleration(double time) const; VEC3 direction(double time, MomBasis::Direction mdir= MomBasis::momdir_) const; // scalar momentum and energy in MeV/c units - double momentum(double time=0) const { return fabs(mass_ * pbar() / absmbar_); } + double momentum(double time=0) const { return fabs(mass_ * pbar() / mbar()); } double momentumVariance(double time=0) const; double positionVariance(double time,MomBasis::Direction dir) const; PMAT planeCovariance(double time,Plane const& plane) const; - double energy(double time=0) const { return fabs(mass_ * ebar() / absmbar_); } + double energy(double time=0) const { return fabs(mass_ * ebar() / mbar()); } // local momentum direction basis void print(std::ostream& ost, int detail) const; TimeRange const& range() const { return trange_; } TimeRange& range() { return trange_; } void setRange(TimeRange const& trange) { trange_ = trange; } void setBNom(double time, VEC3 const& bnom); + // change the BField. This also resets the transforms + void resetBNom(VEC3 const& bnom); bool inRange(double time) const { return trange_.inRange(time); } // momentum change derivatives; this is required to instantiate a KalTrk using this KTraj @@ -103,20 +107,20 @@ namespace KinKal { // helicity is defined as the sign of the projection of the angular momentum vector onto the linear momentum vector double helicity() const { return copysign(1.0,tanDip()); } // needs to be checked TODO double pbar() const { return 1./ (omega() * cosDip() ); } // momentum in mm - double ebar2() const { return pbar()*pbar() + absmbar_ * absmbar_; } + double ebar2() const { return pbar()*pbar() + mbar() * mbar(); } double ebar() const { return sqrt(ebar2()); } // energy in mm + double mbar() const { return fabs(mass_/Q()); } // mass in mm + double Q() const { return -BFieldMap::cbar()*charge()*bnom_.R(); } // reduced charge double cosDip() const { return 1./sqrt(1.+ tanDip() * tanDip() ); } double sinDip() const { return tanDip()*cosDip(); } - double mbar() const { return copysign(absmbar_,omega()); } // mass in mm; includes charge information! - double Q() const { return mass_/copysign(absmbar_,omega()); } // reduced charge double omegaZ() const { return omega()/(CLHEP::c_light*beta()*tanDip()); } // dPhi/dz double beta() const { return fabs(pbar()/ebar()); } // relativistic beta - double gamma() const { return fabs(ebar()/absmbar_); } // relativistic gamma - double betaGamma() const { return fabs(pbar()/absmbar_); } // relativistic betagamma + double gamma() const { return fabs(ebar()/mbar()); } // relativistic gamma + double betaGamma() const { return fabs(pbar()/mbar()); } // relativistic betagamma double Omega() const { return Q()*CLHEP::c_light/energy(); } // true angular velocity double dphi(double t) const { return Omega()*(t - t0()); } // rotation WRT 0 at a given time double phi(double t) const { return dphi(t) + phi0(); } // absolute azimuth at a given time - VEC3 const &bnom(double time=0.0) const { return bnom_; } + VEC3 const& bnom(double time=0.0) const { return bnom_; } double bnomR() const { return bnom_.R(); } DPDV dPardX(double time) const; DPDV dPardM(double time) const; @@ -128,6 +132,7 @@ namespace KinKal { // Parameter derivatives given a change in BFieldMap DVEC dPardB(double time) const; DVEC dPardB(double time, VEC3 const& BPrime) const; + PSMAT dPardPardB(double time,VEC3 const& db) const; // Parameter covariance rotation for a change in BField // flip the helix in time and charge; it remains unchanged geometrically void invertCT() { @@ -151,11 +156,11 @@ namespace KinKal { PSMAT dPardStateLoc(double time) const; // derivative of parameters WRT local state double rc() const { return -1.0/omega() - d0(); } VEC3 center() const { return VEC3(rc()*sin(phi0()), -rc()*cos(phi0()), 0.0); } // local circle center + void setTransforms(); // define global to local and local to global given BNom TimeRange trange_; Parameters pars_; // parameters double mass_; // in units of MeV/c^2 - double absmbar_; // reduced mass in units of mm, computed from the mass and nominal field - int abscharge_; // absolute value of charge in units of proton charge + int abscharge_; // absolute value of charge in units of proton charge; we need to take the sign of the charge from omega VEC3 bnom_; // nominal BField vector, from the map ROOT::Math::Rotation3D l2g_, g2l_; // rotations between local and global coordinates const static std::vector paramTitles_; diff --git a/Trajectory/ClosestApproach.hh b/Trajectory/ClosestApproach.hh index 6284c92a..0c0e3ea8 100644 --- a/Trajectory/ClosestApproach.hh +++ b/Trajectory/ClosestApproach.hh @@ -10,6 +10,7 @@ #include #include #include +#include namespace KinKal { // Hint class for TCA calculation. TCA search will start at these TOCA values. This allows to @@ -130,9 +131,13 @@ namespace KinKal { VEC3 delta = sensorPoca().Vect()-particlePoca().Vect(); double ddot = sensorDirection().Dot(particleDirection()); double denom = 1.0 - ddot*ddot; - // check for parallel) + // check for error conditions if(denom<1.0e-5){ - tpdata_.status_ = ClosestApproachData::pocafailed; + tpdata_.status_ = ClosestApproachData::parallel; + break; + } + if(std::isnan(denom)){ + tpdata_.status_ = ClosestApproachData::failed; break; } double pdd = delta.Dot(particleDirection()); @@ -143,17 +148,12 @@ namespace KinKal { // update the TOCA estimates setParticleTOCA(particleToca()+dptoca); setSensorTOCA(sensorToca()+dstoca); - } - if(tpdata_.status_ != ClosestApproachData::pocafailed){ - if(niter < maxiter) - tpdata_.status_ = ClosestApproachData::converged; - else - tpdata_.status_ = ClosestApproachData::unconverged; // need to add divergence and oscillation tests TODO } // fill the rest of the state if(usable()){ - // sign doca by angular momentum projected onto difference vector. This is just a convention + if(fabs(dptoca) < precision() && fabs(dstoca) < precision() )tpdata_.status_ = ClosestApproachData::converged; + // sign doca by angular momentum projected onto difference vector. This is just a convention VEC3 dvec = delta().Vect(); VEC3 pdir = particleDirection(); tpdata_.lsign_ = copysign(1.0,sensorDirection().Cross(pdir).Dot(dvec)); diff --git a/Trajectory/ClosestApproachData.cc b/Trajectory/ClosestApproachData.cc index 0a8a8f72..f4db21ed 100644 --- a/Trajectory/ClosestApproachData.cc +++ b/Trajectory/ClosestApproachData.cc @@ -1,6 +1,6 @@ #include "KinKal/Trajectory/ClosestApproachData.hh" namespace KinKal { - const std::vector ClosestApproachData::statusNames_ = { "converged", "unconverged", "oscillating", "diverged","pocafailed", "invalid"}; + const std::vector ClosestApproachData::statusNames_ = { "converged", "unconverged", "oscillating", "diverged","parallel","failed", "invalid"}; std::string const& ClosestApproachData::statusName(TPStat status) { return statusNames_[static_cast(status)];} std::ostream& operator << (std::ostream& ost, ClosestApproachData const& cadata) { diff --git a/Trajectory/ClosestApproachData.hh b/Trajectory/ClosestApproachData.hh index 6500f93c..d271541e 100644 --- a/Trajectory/ClosestApproachData.hh +++ b/Trajectory/ClosestApproachData.hh @@ -10,7 +10,7 @@ namespace KinKal { struct ClosestApproachData { - enum TPStat{converged=0,unconverged,oscillating,diverged,pocafailed,invalid}; + enum TPStat{converged=0,unconverged,oscillating,diverged,parallel,failed,invalid}; static std::string const& statusName(TPStat status); //accessors VEC4 const& particlePoca() const { return partCA_; } diff --git a/Trajectory/KinematicLine.cc b/Trajectory/KinematicLine.cc index 76198588..f806bfe2 100644 --- a/Trajectory/KinematicLine.cc +++ b/Trajectory/KinematicLine.cc @@ -8,6 +8,7 @@ #include "Math/AxisAngle.h" #include "Math/VectorUtil.h" #include +#include #include using namespace std; @@ -96,6 +97,13 @@ namespace KinKal { KinematicLine::KinematicLine(KinematicLine const& other, VEC3 const& bnom, double trot) : KinematicLine(other) { } + void KinematicLine::syncPhi0(KinematicLine const& other) { +// adjust the phi0 of this traj to agree with the reference, keeping its value (mod 2pi) the same. + static double twopi = 2*M_PI; + int nloop = static_cast(std::round( (other.phi0() - phi0())/twopi)); + if(nloop != 0) pars_.parameters()[phi0_] += nloop*twopi; + } + VEC3 KinematicLine::position3(double time) const { return (pos0() + flightLength(time) * direction()); } diff --git a/Trajectory/KinematicLine.hh b/Trajectory/KinematicLine.hh index 2c7d9183..28429fd8 100644 --- a/Trajectory/KinematicLine.hh +++ b/Trajectory/KinematicLine.hh @@ -39,6 +39,7 @@ class KinematicLine { static std::string const ¶mUnit(ParamIndex index); constexpr static ParamIndex t0Index() { return t0_; } + constexpr static ParamIndex phi0Index() { return phi0_; } static std::string const &trajName(); // This also requires the nominal BField, which can be a vector (3d) or a @@ -54,6 +55,7 @@ class KinematicLine { KinematicLine(Parameters const& pars, double mass, int charge, VEC3 const& bnom, TimeRange const& trange=TimeRange() ); KinematicLine(Parameters const& pars); + void syncPhi0(KinematicLine const& other); explicit KinematicLine(ParticleState const& pstate, VEC3 const& bnom, TimeRange const& range=TimeRange()); // same, including covariance information @@ -100,7 +102,8 @@ class KinematicLine { TimeRange const &range() const { return trange_; } TimeRange &range() {return trange_; } virtual void setRange(TimeRange const &trange) { trange_ = trange; } - void setBNom(double time, VEC3 const& bnom) { bnom_ = bnom; } + void setBNom(double time, VEC3 const& bnom) { resetBNom(bnom); } + void resetBNom(VEC3 const& bnom) { bnom_ = bnom; } bool inRange(double time) const { return trange_.inRange(time); } double speed() const { return ( mom()/ energy()) * CLHEP::c_light; } @@ -139,6 +142,7 @@ class KinematicLine { // Parameter derivatives given a change in BField. These return null for KinematicLine DVEC dPardB(double time) const { return DVEC(); } DVEC dPardB(double time, VEC3 const& BPrime) const { return DVEC(); } + PSMAT dPardPardB(double time,VEC3 const& db) const { return ROOT::Math::SMatrixIdentity(); } // implement 'helix' interface. This has a physically valid interpretion even for a line Ray axis(double time) const; // helix axis in global coordinates double bendRadius() const { return 0.0; } diff --git a/Trajectory/Line.cc b/Trajectory/Line.cc deleted file mode 100644 index 5acebab4..00000000 --- a/Trajectory/Line.cc +++ /dev/null @@ -1,49 +0,0 @@ -#include "KinKal/Trajectory/Line.hh" -#include "KinKal/Trajectory/SensorLine.hh" -#include "KinKal/Trajectory/ConstantDistanceToTime.hh" -#include -#include -#include -#include -#include - -using namespace std; -using namespace ROOT::Math; - -namespace KinKal { - Line::Line(VEC4 const& pos0, VEC3 const& svel, double length ) : Line(pos0.Vect(), pos0.T(), svel, length) {} - Line::Line(VEC3 const& pos0, double tmeas , VEC3 const& svel, double length ) : - t0_(tmeas), d2t_(new ConstantDistanceToTime(sqrt(svel.Mag2()))), sline_(pos0, svel, length) {} - Line::Line(VEC3 const& p0, VEC3 const& p1, double t0, double speed ) : t0_(t0), d2t_(new ConstantDistanceToTime(speed)), sline_(p0, p1) {} - Line::Line(VEC3 const& p0, double t0, VEC3 const& svel, double length, std::shared_ptr d) : t0_(t0), d2t_(d), sline_(p0, svel, length) {} - - VEC3 Line::position3(double time) const { - return sline_.position3(d2t_->distance(time - t0_)); - } - - VEC4 Line::position4(double time) const { - VEC3 pos3 = position3(time); - return VEC4(pos3.X(),pos3.Y(),pos3.Z(),time); - } - - VEC3 Line::velocity(double time) const { - return direction(time)*speed(time); - } - - double Line::TOCA(VEC3 const& point) const { - double s = sline_.DOCA(point); - return s/speed(s) - t0(); - } - - void Line::print(std::ostream& ost, int detail) const { - ost << " Line, initial position " << startPosition() - << " t0 " << t0() - << " direction " << direction() << endl; - } - - std::ostream& operator <<(std::ostream& ost, Line const& tline) { - tline.print(ost,0); - return ost; - } - -} diff --git a/Trajectory/Line.hh b/Trajectory/Line.hh deleted file mode 100644 index 38ec3e71..00000000 --- a/Trajectory/Line.hh +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef KinKal_Line_hh -#define KinKal_Line_hh - -// -// Linear time-based trajectory with a constant velocity. -// Used as part of the kinematic Kalman fit -// -#include -#include "KinKal/General/Vectors.hh" -#include "KinKal/General/TimeRange.hh" -#include "KinKal/Trajectory/DistanceToTime.hh" -#include "KinKal/Trajectory/SensorLine.hh" - -namespace KinKal { - class Line { - public: - // construct from a spacetime point (typically the measurement position and time) and propagation velocity (mm/ns). - Line(VEC4 const& p0, VEC3 const& svel, double length); - Line(VEC3 const& p0, double t0, VEC3 const& svel, double length); - // construct from 2 points plus timing information. P0 is the measurement (near) end, p1 the far end. Signals propagate from far to near - Line(VEC3 const& p0, VEC3 const& p1, double t0, double speed ); - Line(VEC3 const& p0, double t0, VEC3 const& svel, double length, std::shared_ptr d2t); - // accessors - double t0() const { return t0_; } - double& t0() { return t0_; } // detector updates need to refine t0 - - // signal ends at pos0 - VEC3 startPosition() const { return sline_.startPosition(); } - VEC3 const& endPosition() const { return sline_.endPosition(); } - double speed(double time) const { return d2t_->speed(d2t_->distance(time-t0_)); } - double length() const { return sline_.length(); } - VEC3 const& direction() const { return sline_.direction(); } - // TOCA to a point - double TOCA(VEC3 const& point) const; - // geometric accessors - VEC3 position3(double time) const; - VEC4 position4(double time) const; - VEC3 velocity(double time) const; - VEC3 const& direction(double time) const { return sline_.direction(); } - void print(std::ostream& ost, int detail) const; - double timeAtMidpoint() const { return t0_ + d2t_->time(0.5*length()); } - - private: - double t0_; // intial time (at pos0) - std::shared_ptr d2t_; // represents the possibly nonlinear distance to time relationship of the line - SensorLine sline_; // geometic representation of the line - }; - std::ostream& operator <<(std::ostream& ost, Line const& tline); -} -#endif diff --git a/Trajectory/LoopHelix.cc b/Trajectory/LoopHelix.cc index b164fab3..ceba6af6 100644 --- a/Trajectory/LoopHelix.cc +++ b/Trajectory/LoopHelix.cc @@ -28,57 +28,64 @@ namespace KinKal { string const& LoopHelix::trajName() { return trajName_; } LoopHelix::LoopHelix() : mass_(0.0), charge_(0) {} - LoopHelix::LoopHelix( VEC4 const& pos0, MOM4 const& mom0, int charge, double bnom, TimeRange const& range) : LoopHelix(pos0,mom0,charge,VEC3(0.0,0.0,bnom),range) {} - LoopHelix::LoopHelix( VEC4 const& pos0, MOM4 const& mom0, int charge, VEC3 const& bnom, TimeRange const& trange) : trange_(trange), mass_(mom0.M()), charge_(charge), bnom_(bnom) { + LoopHelix::LoopHelix( VEC4 const& gpos, MOM4 const& gmom, int charge, double bnom, TimeRange const& range) : LoopHelix(gpos,gmom,charge,VEC3(0.0,0.0,bnom),range) {} + LoopHelix::LoopHelix( VEC4 const& gpos, MOM4 const& gmom, int charge, VEC3 const& bnom, TimeRange const& trange) : trange_(trange), mass_(gmom.M()), charge_(charge), bnom_(bnom) { static double twopi = 2*M_PI; - // Transform into the system where Z is along the Bfield, which is the implicit coordinate system of the parameterization. - // The transform is a pure rotation about the origin - VEC4 pos(pos0); - MOM4 mom(mom0); - g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); - if(fabs(g2l_(bnom_).Theta()) > 1.0e-6)throw invalid_argument("Rotation Error"); - // to convert global vectors into parameters they must first be rotated into the local system. - pos = g2l_(pos); - mom = g2l_(mom); - // create inverse rotation; this moves back into the global coordinate system - l2g_ = g2l_.Inverse(); + // Transform position and momentum into the system where Z is along the Bfield, which is the implicit coordinate system of the parameterization. + // This is a pure rotation about the origin + setTransforms(); + auto lpos = g2l_(gpos); + auto lmom = g2l_(gmom); // compute some simple useful parameters - double pt = mom.Pt(); - double phibar = mom.Phi(); + double pt = lmom.Pt(); + double lmomphi = lmom.Phi(); // translation factor from MeV/c to curvature radius in mm, B in Tesla; signed by the charge!!! - double momToRad = 1.0/Q(); + double invq = 1.0/Q(); // transverse radius of the helix - param(rad_) = pt*momToRad; + param(rad_) = pt*invq; // longitudinal wavelength - param(lam_) = mom.Z()*momToRad; - // time at z=0 - double om = omega(); - param(t0_) = pos.T() - pos.Z()/(om*lam()); + param(lam_) = lmom.Z()*invq; + // time when particle reaches local z=0 + double zbar = lpos.Z()/lam(); + param(t0_) = lpos.T() - zbar/omega(); // compute winding that puts phi0 in the range -pi,pi - double nwind = rint((pos.Z()/lam() - phibar)/twopi); - // cout << "winding number = " << nwind << endl; - // azimuth at z=0 - param(phi0_) = phibar - om*(pos.T()-t0()) + twopi*nwind; + double nwind = round((zbar - lmomphi)/twopi); + // particle momentum azimuth at z=0 + param(phi0_) = lmomphi - zbar + twopi*nwind; // circle center - param(cx_) = pos.X() - mom.Y()*momToRad; - param(cy_) = pos.Y() + mom.X()*momToRad; - // test -// auto testpos = position3(pos0.T()); -// auto testmom = momentum3(pos0.T()); -// auto dp = testpos - pos0.Vect(); -// auto dm = testmom - mom0.Vect(); -// if(dp.R() > 1.0e-5 || dm.R() > 1.0e-5)throw invalid_argument("Construction Test Failure"); - } - - void LoopHelix::setBNom(double time, VEC3 const& bnom) { - // adjust the parameters for the change in bnom - pars_.parameters() += dPardB(time,bnom); - bnom_ = bnom; + param(cx_) = lpos.X() - lmom.Y()*invq; + param(cy_) = lpos.Y() + lmom.X()*invq; + } + + void LoopHelix::syncPhi0(LoopHelix const& other) { +// adjust the phi0 of this traj to agree with the reference, keeping its value (mod 2pi) the same. + static double twopi = 2*M_PI; + int nloop = static_cast(round( (other.phi0() - phi0())/twopi)); + if(nloop != 0) pars_.parameters()[phi0_] += nloop*twopi; + } + + void LoopHelix::setTransforms() { // adjust rotations to global space g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); l2g_ = g2l_.Inverse(); } + void LoopHelix::setBNom(double time, VEC3 const& newbnom) { + auto db = newbnom - bnom_; +// PSMAT dpdpdb =ROOT::Math::SMatrixIdentity(); +// PSMAT dpdpdb = dPardPardB(time,db); +// std::cout << "dpdpdb = " << dpdpdb << std::endl; +// pars_.covariance() = ROOT::Math::Similarity(dpdpdb,pars_.covariance()); + pars_.parameters() += dPardB(time,db); + // rotate covariance: for now, just for the magnitude change. Rotation is still TODO + resetBNom(newbnom); + } + + void LoopHelix::resetBNom(VEC3 const& newbnom) { + bnom_ = newbnom; + setTransforms(); + } + LoopHelix::LoopHelix(LoopHelix const& other, VEC3 const& bnom, double tref) : LoopHelix(other) { setBNom(tref,bnom); } @@ -89,10 +96,8 @@ namespace KinKal { LoopHelix::LoopHelix( Parameters const& pars, double mass, int charge, VEC3 const& bnom, TimeRange const& trange ) : trange_(trange), pars_(pars), mass_(mass), charge_(charge), bnom_(bnom) { - // set the transforms - g2l_ = Rotation3D(AxisAngle(VEC3(sin(bnom_.Phi()),-cos(bnom_.Phi()),0.0),bnom_.Theta())); - l2g_ = g2l_.Inverse(); - } + setTransforms(); + } LoopHelix::LoopHelix(ParticleState const& pstate, VEC3 const& bnom, TimeRange const& range) : LoopHelix(pstate.position4(),pstate.momentum4(),pstate.charge(),bnom,range) @@ -257,33 +262,60 @@ namespace KinKal { return dPardMLoc(time)*g2lmat; } - DVEC LoopHelix::dPardB(double time) const { + PSMAT LoopHelix::dPardPardB(double time,VEC3 const& db) const { + auto newbnom = bnom_ + db; + double bfrac = bnom_.R()/newbnom.R(); + auto xvec = localPosition(time); double phival = phi(time); + double cphi = cos(phival); + double sphi = sin(phival); + // compute the rows + DVEC dpdpdb_rad; dpdpdb_rad[rad_] = bfrac; + DVEC dpdpdb_lam; dpdpdb_lam[lam_] = bfrac; + DVEC dpdpdb_cx; dpdpdb_cx[cx_] = 1.0; dpdpdb_cx[rad_] = sphi*(1-bfrac); + DVEC dpdpdb_cy; dpdpdb_cy[cy_] = 1.0; dpdpdb_cy[rad_] = -cphi*(1-bfrac); + DVEC dpdpdb_phi0; dpdpdb_phi0[phi0_] = 1.0; dpdpdb_phi0[lam_] = -xvec.Z()*(1-1.0/bfrac)/(lam()*lam()); + DVEC dpdpdb_t0; dpdpdb_t0[t0_] = 1.0; + // build the matrix from these rows + PSMAT dpdpdb; + dpdpdb.Place_in_row(dpdpdb_rad,rad_,0); + dpdpdb.Place_in_row(dpdpdb_lam,lam_,0); + dpdpdb.Place_in_row(dpdpdb_cx,cx_,0); + dpdpdb.Place_in_row(dpdpdb_cy,cy_,0); + dpdpdb.Place_in_row(dpdpdb_phi0,phi0_,0); + dpdpdb.Place_in_row(dpdpdb_t0,t0_,0); + return dpdpdb; + } + + DVEC LoopHelix::dPardB(double time, VEC3 const& db) const { + // record position and momentum in local coordinates; these are constant + auto xvec = localPosition(time); + auto mvec = localMomentum(time); + // translate newbnom to local coordinates + VEC3 dbloc = g2l_(db); + // find the parameter change due to bnom magnitude change. These are exact + double br = bnom_.R(); + auto zhat = VEC3(0.0,0.0,1.0); + auto bloc = br*zhat; + auto newbloc = dbloc + bloc; + double nbr = newbloc.R(); DVEC retval; - retval[rad_] = -rad(); - retval[lam_] = -lam(); - retval[cx_] = rad()*sin(phival); - retval[cy_] = -rad()*cos(phival); - retval[phi0_] = -dphi(time); + double dbf = (br - nbr)/nbr; + retval[rad_] = rad()*dbf; + retval[lam_] = lam()*dbf; + retval[phi0_] = xvec.Z()*(br -nbr)/br/lam(); retval[t0_] = 0.0; - return (1.0/bnom_.R())*retval; - } + retval[cx_] = -sin(mvec.Phi())*retval[rad_]; + retval[cy_] = cos(mvec.Phi())*retval[rad_]; - DVEC LoopHelix::dPardB(double time, VEC3 const& BPrime) const { - // rotate new B field difference into local coordinate system - VEC3 dB = g2l_(BPrime-bnom_); - // find the parameter change due to BField magnitude change usng component parallel to the local nominal Bfield (always along z) - DVEC retval = dPardB(time)*dB.Z(); - // find the change in (local) position and momentum due to the rotation implied by the B direction change + // find the change in local position and momentum due to the rotation implied by the B direction change // work in local coordinate system to avoid additional matrix mulitplications - auto xvec = localPosition(time); - auto mvec = localMomentum(time); - VEC3 BxdB =VEC3(0.0,0.0,1.0).Cross(dB)/bnomR(); + VEC3 BxdB = zhat.Cross(dbloc)/br; VEC3 dx = xvec.Cross(BxdB); VEC3 dm = mvec.Cross(BxdB); - // convert these to a full state vector change + // convert these to the state vector change ParticleState dstate(dx,dm,time,mass(),charge()); - // convert the change in (local) state due to rotation to parameter space + // convert the (local) state due to rotation to parameter space retval += dPardStateLoc(time)*dstate.state(); return retval; } diff --git a/Trajectory/LoopHelix.hh b/Trajectory/LoopHelix.hh index ce5a5766..db832f12 100644 --- a/Trajectory/LoopHelix.hh +++ b/Trajectory/LoopHelix.hh @@ -29,6 +29,7 @@ namespace KinKal { // classes implementing the Kalman fit // define the indices and names of the parameters enum ParamIndex {rad_=0,lam_=1,cx_=2,cy_=3,phi0_=4,t0_=5,npars_=6}; + constexpr static ParamIndex phi0Index() { return phi0_; } constexpr static ParamIndex t0Index() { return t0_; } static std::vector const& paramNames(); @@ -55,6 +56,8 @@ namespace KinKal { LoopHelix( Parameters const& pars, double mass, int charge, VEC3 const& bnom, TimeRange const& trange=TimeRange() ); // copy payload and override the parameters; Is this used? LoopHelix(Parameters const& pdata, LoopHelix const& other); + // synchronize phi0, which has a 2pi wrapping + void syncPhi0(LoopHelix const& other); VEC4 position4(double time) const; VEC3 position3(double time) const; VEC3 velocity(double time) const; @@ -67,6 +70,8 @@ namespace KinKal { void setRange(TimeRange const& trange) { trange_ = trange; } // allow resetting the BField. Note this is time-dependent void setBNom(double time, VEC3 const& bnom); + // change the BField. This also resets the transforms + void resetBNom(VEC3 const& bnom); bool inRange(double time) const { return trange_.inRange(time); } VEC3 momentum3(double time) const; MOM4 momentum4(double time) const; @@ -127,9 +132,8 @@ namespace KinKal { PSMAT dStatedPar(double time) const; // derivative of global state WRT parameters DVEC momDeriv(double time, MomBasis::Direction mdir) const; // projection of M derivatives onto direction basis // package the above for full (global) state - // Parameter derivatives given a change in BField - DVEC dPardB(double time) const; // parameter derivative WRT change in BField magnitude - DVEC dPardB(double time, VEC3 const& BPrime) const; // parameter change given a new BField vector + DVEC dPardB(double time, VEC3 const& db) const; // parameter change given a change in BField vector; this includes the magnitude and direction changes + PSMAT dPardPardB(double time,VEC3 const& db) const; // Parameter covariance rotation for a change in BField // helix interface VEC3 center(double time) const; // helix center in global coordinates Ray axis(double time) const; // helix axis in global coordinates @@ -143,6 +147,7 @@ namespace KinKal { DPDV dPardXLoc(double time) const; // return the derivative of the parameters WRT the local (unrotated) position vector DPDV dPardMLoc(double time) const; // return the derivative of the parameters WRT the local (unrotated) momentum vector PSMAT dPardStateLoc(double time) const; // derivative of parameters WRT local state + void setTransforms(); // define global to local and local to global given BNom TimeRange trange_; Parameters pars_; // parameters diff --git a/Trajectory/POCAUtil.cc b/Trajectory/POCAUtil.cc index 6ef3e350..017671e2 100644 --- a/Trajectory/POCAUtil.cc +++ b/Trajectory/POCAUtil.cc @@ -1,15 +1,3 @@ -/* - - code based on MuteUtilities/src/TwoLinePCA: - - Given two lines in 3D, compute the distance of closest - approach between the two lines. The lines are - specified in point-slope form. - - Original author Rob Kutschke - - */ - #include "KinKal/Trajectory/POCAUtil.hh" #include #include diff --git a/Trajectory/POCAUtil.hh b/Trajectory/POCAUtil.hh index 7ddc9975..e20930f6 100644 --- a/Trajectory/POCAUtil.hh +++ b/Trajectory/POCAUtil.hh @@ -1,5 +1,16 @@ #ifndef KinKal_POCAUtil_hh #define KinKal_POCAUtil_hh +/* + + code based on MuteUtilities/src/TwoLinePCA: + + Given two lines in 3D, compute the distance of closest + approach between the two lines. The lines are + specified in point-slope form. + + Original author Rob Kutschke + + */ #include "KinKal/General/Vectors.hh" #include "Math/Rotation3D.h" diff --git a/Trajectory/SensorLine.cc b/Trajectory/SensorLine.cc index 848c2fe5..f5fab7ff 100644 --- a/Trajectory/SensorLine.cc +++ b/Trajectory/SensorLine.cc @@ -1,34 +1,52 @@ #include "KinKal/Trajectory/SensorLine.hh" +#include "KinKal/Trajectory/ConstantDistanceToTime.hh" #include #include #include #include +#include using namespace std; using namespace ROOT::Math; namespace KinKal { - SensorLine::SensorLine(VEC4 const& p0, VEC3 const& svel, double length) : SensorLine(p0.Vect(), svel, length) {} - SensorLine::SensorLine(VEC3 const& p0, VEC3 const& svel, double length ) : - pos0_(p0), dir_(svel.Unit()), length_(length) {} - SensorLine::SensorLine(VEC3 const& p0, VEC3 const& p1) : pos0_(p0), dir_((p0-p1).Unit()), length_((p1-p0).R()) {} - VEC3 SensorLine::position3(double distance) const { - return pos0_ + distance*dir_; + SensorLine::SensorLine(VEC4 const& mpos, VEC3 const& svel, double length ) : SensorLine(mpos.Vect(), mpos.T(), svel, length) {} + +// note the measurement position is at the END of the directed line segment, since that points in the direction of signal propagation + SensorLine::SensorLine(VEC3 const& mpos, double mtime , VEC3 const& svel, double length ) : + mtime_(mtime), d2t_(new ConstantDistanceToTime(svel.R())), lineseg_(mpos-svel.unit()*length, mpos) {} + + SensorLine::SensorLine(VEC3 const& mpos, VEC3 const& endpos, double mtime, double speed ) : mtime_(mtime), + d2t_(new ConstantDistanceToTime(speed)), lineseg_(endpos,mpos) {} + + SensorLine::SensorLine(VEC3 const& mpos, double mtime, VEC3 const& svel, double length, std::shared_ptr d) : mtime_(mtime), + d2t_(d), lineseg_(mpos-svel.unit()*length, mpos) {} + + VEC3 SensorLine::position3(double time) const { + return lineseg_.endPosition(d2t_->distance(mtime_ - time)); + } + + VEC4 SensorLine::position4(double time) const { + VEC3 pos3 = position3(time); + return VEC4(pos3.X(),pos3.Y(),pos3.Z(),time); + } + + VEC3 SensorLine::velocity(double time) const { + return direction(time)*speed(time); } - double SensorLine::DOCA(VEC3 const& point) const { - return (point - pos0_).Dot(dir_); + double SensorLine::timeTo(VEC3 const& point) const { + double s = lineseg_.distanceTo(point); + return s/speed(s) - measurementTime(); } void SensorLine::print(std::ostream& ost, int detail) const { - ost << " SensorLine, initial position " << pos0_ - << " direction " << dir_ - << " length " << length_ << endl; + ost << *this << endl; } - std::ostream& operator <<(std::ostream& ost, SensorLine const& tSensorLine) { - tSensorLine.print(ost,0); + std::ostream& operator <<(std::ostream& ost, SensorLine const& sline) { + ost << "SensorLine " << sline.line() << " measurement time " << sline.measurementTime(); return ost; } diff --git a/Trajectory/SensorLine.hh b/Trajectory/SensorLine.hh index 684456a0..9e8eadc6 100644 --- a/Trajectory/SensorLine.hh +++ b/Trajectory/SensorLine.hh @@ -1,34 +1,49 @@ #ifndef KinKal_SensorLine_hh #define KinKal_SensorLine_hh - -// Used as part of the kinematic Kalman fit // +// Class describing linear signal propagation in a sensor with an arbitrary distance-to-time relationship +// +#include #include "KinKal/General/Vectors.hh" +#include "KinKal/Trajectory/DistanceToTime.hh" +#include "KinKal/Geometry/LineSegment.hh" namespace KinKal { class SensorLine { public: - // construct from a spacetime point (typically the measurement position and time) and propagation velocity (mm/ns). - SensorLine(VEC4 const& p0, VEC3 const& svel, double length); - SensorLine(VEC3 const& p0, VEC3 const& svel, double length); - // construct from 2 points. P0 is the signal origin position (start position), p1 the measurement position(end position). Signals propagate from start to end. - SensorLine(VEC3 const& p0, VEC3 const& p1); + // construct from the measurement position and time and signal propagation velocity (mm/ns). Note the velocity points TOWARDS the sensor + SensorLine(VEC4 const& mpos, VEC3 const& svel, double length); + SensorLine(VEC3 const& mpos, double mtime, VEC3 const& svel, double length); + // construct from measurement position and time, and opposite (far) end of the sensor. Signals propagate from the far end to the measurement position with the given speed + SensorLine(VEC3 const& mpos, VEC3 const& endpos, double mtime, double speed ); + SensorLine(VEC3 const& mpos, double mtime, VEC3 const& svel, double length, std::shared_ptr d2t); // accessors - // signal ends at pos0 - VEC3 startPosition() const { return pos0_ - length_*dir_; } - VEC3 const& endPosition() const { return pos0_; } - double length() const { return length_; } - VEC3 const& direction() const { return dir_; } - // Distance of Closest Approach to a point - double DOCA(VEC3 const& point) const; + auto const& measurementPosition() const { return lineseg_.end(); } // measurement is at the end (latest time) + double measurementTime() const { return mtime_; } + double& measurementTime() { return mtime_; } // Fit updates need to refine this + // LineSegment positions + auto const& start() const { return lineseg_.start(); } + auto const& end() const { return lineseg_.end(); } + auto middle() const { return lineseg_.middle(); } + // + auto const& line() const { return lineseg_; } + double speed(double time) const { return d2t_->speed(d2t_->distance(mtime_ - time)); } // D2T is relative to measurement end + double length() const { return lineseg_.length(); } + VEC3 const& direction() const { return lineseg_.direction(); } + // time to a point + double timeTo(VEC3 const& point) const; // geometric accessors - VEC3 position3(double distance) const; + VEC3 position3(double time) const; + VEC4 position4(double time) const; + VEC3 velocity(double time) const; // signal velocity + VEC3 const& direction(double time) const { return lineseg_.direction(); } void print(std::ostream& ost, int detail) const; - + double timeAtMidpoint() const { return mtime_ - d2t_->time(0.5*length()); } private: - VEC3 pos0_, dir_; // signal physical origin position and signal propagation direction - double length_; // distance from signal origin to measurement (electronics) point + double mtime_; // measurement time, at the far end of the directed line segment + std::shared_ptr d2t_; // distance to time relationship along the sensor + LineSegment lineseg_; // linear representation of the sensor }; - std::ostream& operator <<(std::ostream& ost, SensorLine const& tSensorLine); + std::ostream& operator <<(std::ostream& ost, SensorLine const& sline); } #endif