Merge pull request #215 from brownd1978/trajptr

Use trajectory ptrs in ClosestApproach

Author: brownd1978

Detector/Hit.hh

@@ -35,7 +35,7 @@ namespace KinKal {
       virtual void print(std::ostream& ost=std::cout,int detail=0) const = 0;
       // update to a new reference, without changing internal state
       virtual void updateReference(PTRAJ const& ptraj) = 0;
-      virtual KTRAJPTR const& refTrajPtr() const = 0;
+      virtual KTRAJPTR refTrajPtr() const = 0;
       // update the internals of the hit, specific to this meta-iteraion
       virtual void updateState(MetaIterConfig const& config,bool first) = 0;
       // The following provides the constraint/information content of this hit in the trajectory weight space

Detector/ParameterHit.hh

@@ -40,7 +40,7 @@ namespace KinKal {
       // parameter constraints are absolute and can't be updated
       void print(std::ostream& ost=std::cout,int detail=0) const override;
       void updateReference(PTRAJ const& ptraj) override { reftraj_ = ptraj.nearestTraj(time()); }
-      KTRAJPTR const& refTrajPtr() const override { return reftraj_; }
+      KTRAJPTR refTrajPtr() const override { return reftraj_; }
       // ParameterHit-specfic interface
       // construct from constraint values, time, and mask of which parameters to constrain
       ParameterHit(double time, PTRAJ const& ptraj, Parameters const& params, PMASK const& pmask);

Examples/ScintHit.hh

@@ -17,18 +17,13 @@ namespace KinKal {
       using RESIDHIT = ResidualHit<KTRAJ>;
       using HIT = Hit<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
+      using STRAJPTR = std::shared_ptr<SensorLine>;
       // copy constructor
       ScintHit(ScintHit<KTRAJ> const& rhs):
           ResidualHit<KTRAJ>(rhs),
-          saxis_(rhs.sensorAxis()),
           tvar_(rhs.timeVariance()),
           wvar_(rhs.widthVariance()),
-          tpca_(
-                rhs.closestApproach().particleTraj(),
-                saxis_,
-                rhs.closestApproach().hint(),
-                rhs.closestApproach().precision()
-          ),
+          tpca_(rhs.tpca_),
           rresid_(rhs.refResidual()){
         /**/
       };
@@ -46,20 +41,19 @@ namespace KinKal {
       Residual const& refResidual(unsigned ires=0) const override;
       double time() const override { return tpca_.particleToca(); }
       void updateReference(PTRAJ const& ptraj) override;
-      KTRAJPTR const& refTrajPtr() const override { return tpca_.particleTrajPtr(); }
+      KTRAJPTR refTrajPtr() const override { return tpca_.particleTrajPtr(); }
       void updateState(MetaIterConfig const& config,bool first) override;
       void print(std::ostream& ost=std::cout,int detail=0) const override;
      // scintHit explicit interface
       ScintHit(PCA const& pca, double tvar, double wvar);
       virtual ~ScintHit(){}
       // the line encapsulates both the measurement value (through t0), and the light propagation model (through the velocity)
-      auto const& sensorAxis() const { return saxis_; }
+      auto sensorAxis() const { return tpca_.sensorTrajPtr(); }
       auto const& closestApproach() const { return tpca_; }
       double timeVariance() const { return tvar_; }
       double widthVariance() const { return wvar_; }
       auto precision() const { return tpca_.precision(); }
     private:
-      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
@@ -70,8 +64,8 @@ namespace KinKal {
   };
 
   template <class KTRAJ> ScintHit<KTRAJ>::ScintHit(PCA const& pca, double tvar, double wvar) :
-    saxis_(pca.sensorTraj()), tvar_(tvar), wvar_(wvar),
-    tpca_(pca.localTraj(),saxis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP())
+    tvar_(tvar), wvar_(wvar),
+    tpca_(static_cast<CA const&>(pca))
   {}
 
   template <class KTRAJ> Residual const& ScintHit<KTRAJ>::refResidual(unsigned ires) const {
@@ -81,9 +75,9 @@ namespace KinKal {
 
   template <class KTRAJ> void ScintHit<KTRAJ>::updateReference(PTRAJ const& ptraj) {
     // use previous hint, or initialize from the sensor time
-    CAHint tphint = tpca_.usable() ?  tpca_.hint() : CAHint(saxis_.measurementTime(), saxis_.measurementTime());
-    PCA pca(ptraj,saxis_,tphint,precision());
-    tpca_ = pca.localClosestApproach();
+    CAHint tphint = tpca_.usable() ?  tpca_.hint() : CAHint(sensorAxis()->measurementTime(), sensorAxis()->measurementTime());
+    PCA pca(ptraj,sensorAxis(),tphint,precision());
+    tpca_ = static_cast<CA>(pca);
     if(!tpca_.usable())throw std::runtime_error("ScintHit TPOCA failure");
   }
 
@@ -92,26 +86,26 @@ 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 const& sstart = saxis_.start();
-    auto const& send = saxis_.end();
+    auto const& sstart = sensorAxis()->start();
+    auto const& send = sensorAxis()->end();
     // tolerance should come from the config.  Should also test relative to the error. FIXME
-    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) {
+    double tol = sensorAxis()->length()*1.0;
+    double sdist = (ppos - sensorAxis()->middle()).Dot(sensorAxis()->direction());
+    if( (ppos-sstart).Dot(sensorAxis()->direction()) < -tol || (ppos-send).Dot(sensorAxis()->direction()) > tol) {
       // adjust hint to the middle and try agian
-      double sspeed = tpca_.particleTraj().velocity(tpca_.particleToca()).Dot(saxis_.direction());
+      double sspeed = tpca_.particleTraj().velocity(tpca_.particleToca()).Dot(sensorAxis()->direction());
       auto tphint = tpca_.hint();
       tphint.particleToca_ -= sdist/sspeed;
-      tpca_ = CA(tpca_.particleTrajPtr(),saxis_,tphint,precision());
+      tpca_ = CA(tpca_.particleTrajPtr(),sensorAxis(),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());
+      sdist = (tpca_.particlePoca().Vect() - sensorAxis()->middle()).Dot(sensorAxis()->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
     double dd2 = tpca_.dirDot()*tpca_.dirDot();
-    double totvar = tvar_ + wvar_*dd2/(saxis_.speed(sdist)*saxis_.speed(sdist)*(1.0-dd2));
+    double totvar = tvar_ + wvar_*dd2/(sensorAxis()->speed(sdist)*sensorAxis()->speed(sdist)*(1.0-dd2));
     rresid_ = Residual(tpca_.deltaT(),totvar,0.0,tpca_.dTdP());
     this->updateWeight(config);
   }
@@ -124,7 +118,7 @@ namespace KinKal {
     ost << " ScintHit  tvar " << tvar_ << " wvar " << wvar_ << std::endl;
     if(detail > 0){
       ost << "Line ";
-      saxis_.print(ost,detail);
+      sensorAxis()->print(ost,detail);
     }
   }
 

Examples/SimpleWireHit.hh

@@ -22,6 +22,8 @@ namespace KinKal {
       using CA = ClosestApproach<KTRAJ,SensorLine>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       using PTRAJ = ParticleTrajectory<KTRAJ>;
+      using STRAJPTR = std::shared_ptr<SensorLine>;
+
       enum Dimension { dresid=0, tresid=1};  // residual dimensions
 
       SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca,
@@ -32,13 +34,7 @@ namespace KinKal {
           ResidualHit<KTRAJ>(rhs),
           bfield_(rhs.bfield()),
           whstate_(rhs.hitState()),
-          wire_(rhs.wire()),
-          ca_(
-            rhs.closestApproach().particleTraj(),
-            wire_,
-            rhs.closestApproach().hint(),
-            rhs.closestApproach().precision()
-          ),
+          tpca_(rhs.tpca_),
           rresid_(rhs.residuals()),
           mindoca_(rhs.minDOCA()),
           dvel_(driftVelocity()),
@@ -57,11 +53,11 @@ namespace KinKal {
         rv->setClosestApproach(ca);
         return rv;
       };
-      double time() const override { return ca_.particleToca(); }
+      double time() const override { return tpca_.particleToca(); }
       VEC3 dRdX(unsigned ires) const;
       Residual const& refResidual(unsigned ires=dresid) const override;
       void updateReference(PTRAJ const& ptraj) override;
-      KTRAJPTR const& refTrajPtr() const override { return ca_.particleTrajPtr(); }
+      KTRAJPTR refTrajPtr() const override { return tpca_.particleTrajPtr(); }
       void print(std::ostream& ost=std::cout,int detail=0) const override;
       // Use dedicated updater
       void updateState(MetaIterConfig const& config,bool first) override;
@@ -73,24 +69,23 @@ namespace KinKal {
       double minDOCA() const { return mindoca_; }
       int id() const { return id_; }
       CA unbiasedClosestApproach() const;
-      auto const& closestApproach() const { return ca_; }
+      auto const& closestApproach() const { return tpca_; }
       auto const& hitState() const { return whstate_; }
-      auto const& wire() const { return wire_; }
       auto const& bfield() const { return bfield_; }
-      auto precision() const { return ca_.precision(); }
+      auto precision() const { return tpca_.precision(); }
       auto const& residuals() const { return rresid_; }
       double tot() const { return tot_; }
       double totVariance() const { return totvar_; }
+      auto wire() const { return tpca_.sensorTrajPtr(); }
 
     private:
       BFieldMap const& bfield_; // drift calculation requires the BField for ExB effects
       WireHitState whstate_; // current state
-      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
                   // (when multiplied by the propagation velocity).
-      CA ca_; // reference time and position of closest approach to the wire; this is generally biased by the hit
+      CA tpca_; // reference time and position of closest approach to the wire; this is generally biased by the hit
       std::array<Residual,2> rresid_; // residuals WRT most recent reference
       double mindoca_; // effective minimum DOCA used when assigning LR ambiguity, used to define null hit properties
       double dvel_; // constant drift speed
@@ -101,7 +96,7 @@ namespace KinKal {
       void updateResiduals();
 
       // modifiers to support cloning
-      void setClosestApproach(const CA& ca){ ca_ = ca; }
+      void setClosestApproach(const CA& ca){ tpca_ = ca; }
   };
 
   //trivial 'updater' that sets the wire hit state to null
@@ -113,18 +108,18 @@ namespace KinKal {
   template <class KTRAJ> SimpleWireHit<KTRAJ>::SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate,
       double mindoca, double driftspeed, double tvar, double tot, double totvar, double rcell, int id) :
     bfield_(bfield),
-    whstate_(whstate), wire_(pca.sensorTraj()),
-    ca_(pca.localTraj(),wire_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()), // must be explicit to get the right sensor traj reference
+    whstate_(whstate),
+    tpca_(static_cast<CA const&>(pca)),
     mindoca_(mindoca), dvel_(driftspeed), tvar_(tvar), tot_(tot), totvar_(totvar), rcell_(rcell), id_(id) {
     }
 
   template <class KTRAJ> void SimpleWireHit<KTRAJ>::updateReference(PTRAJ const& ptraj) {
     // if we already computed PCA in the previous iteration, use that to set the hint.  This speeds convergence
     // otherwise use the time at the center of the wire
-    CAHint tphint = ca_.usable() ?  ca_.hint() : CAHint(wire_.timeAtMidpoint(),wire_.timeAtMidpoint());
-    PCA pca(ptraj,wire_,tphint,precision());
-    ca_ = pca.localClosestApproach();
-    if(!ca_.usable())throw std::runtime_error("WireHit TPOCA failure");
+    CAHint tphint = tpca_.usable() ?  tpca_.hint() : CAHint(wire()->timeAtMidpoint(),wire()->timeAtMidpoint());
+    PCA pca(ptraj,wire(),tphint,precision());
+    tpca_ = static_cast<CA>(pca);
+    if(!tpca_.usable())throw std::runtime_error("WireHit TPOCA failure");
   }
 
   template <class KTRAJ> void SimpleWireHit<KTRAJ>::updateState(MetaIterConfig const& miconfig, bool first) {
@@ -148,16 +143,16 @@ namespace KinKal {
     }
     rresid_[tresid] = rresid_[dresid] = Residual();
     if(whstate_.active()){
-      rresid_[tresid] = Residual(ca_.deltaT() - tot_, totvar_,0.0,ca_.dTdP()); // always constrain to TOT; this stabilizes the fit
+      rresid_[tresid] = Residual(tpca_.deltaT() - tot_, totvar_,0.0,tpca_.dTdP()); // always constrain to TOT; this stabilizes the fit
       if(whstate_.useDrift()){
         // translate PCA to residual. Use ambiguity assignment to convert drift time to a drift radius
-        double dr = dvel_*whstate_.lrSign()*ca_.deltaT() -ca_.doca();
-        DVEC dRdP = dvel_*whstate_.lrSign()*ca_.dTdP() -ca_.dDdP();
+        double dr = dvel_*whstate_.lrSign()*tpca_.deltaT() -tpca_.doca();
+        DVEC dRdP = dvel_*whstate_.lrSign()*tpca_.dTdP() -tpca_.dDdP();
         rresid_[dresid] = Residual(dr,tvar_*dvel_*dvel_,0.0,dRdP);
       } else {
         // interpret DOCA against the wire directly as a residuals
-        double nulldvar = dvel_*dvel_*(ca_.deltaT()*ca_.deltaT()+0.8);
-        rresid_[dresid] = Residual(ca_.doca(),nulldvar,0.0,ca_.dDdP());
+        double nulldvar = dvel_*dvel_*(tpca_.deltaT()*tpca_.deltaT()+0.8);
+        rresid_[dresid] = Residual(tpca_.doca(),nulldvar,0.0,tpca_.dDdP());
       }
     }
     // now update the weight
@@ -168,9 +163,9 @@ namespace KinKal {
     if (whstate_.active()){
       if (ires == dresid){
         if (whstate_.useDrift()){
-          return ca_.lSign()*ca_.delta().Vect().Unit();
+          return tpca_.lSign()*tpca_.delta().Vect().Unit();
         }else{
-          return -1*ca_.lSign()*ca_.delta().Vect().Unit();
+          return -1*tpca_.lSign()*tpca_.delta().Vect().Unit();
         }
       }
     }
@@ -186,7 +181,7 @@ namespace KinKal {
     // compute the unbiased closest approach; this is brute force, but works
     auto const& ca = this->closestApproach();
     auto uparams = HIT::unbiasedParameters();
-    KTRAJ utraj(uparams,ca.particleTraj());
+    auto utraj = std::make_shared<KTRAJ>(uparams,ca.particleTraj());
     return CA(utraj,this->wire(),ca.hint(),ca.precision());
   }
 
@@ -214,7 +209,7 @@ namespace KinKal {
       ost << std::endl;
     }
     if(detail > 1) {
-      ost << "Approximate Propagation speed " << wire_.speed(100) << " TPOCA " << ca_.tpData() << std::endl;
+      ost << "Approximate Propagation speed " << wire()->speed(100) << " TPOCA " << tpca_.tpData() << std::endl;
     }
   }
 

Examples/StrawXing.hh

@@ -15,6 +15,7 @@ namespace KinKal {
     public:
       using PTRAJ = ParticleTrajectory<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
+      using STRAJPTR = std::shared_ptr<SensorLine>;
       using EXING = ElementXing<KTRAJ>;
       using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
       using CA = ClosestApproach<KTRAJ,SensorLine>;
@@ -24,14 +25,8 @@ namespace KinKal {
       // copy constructor
       StrawXing(StrawXing const& rhs):
           ElementXing<KTRAJ>(rhs),
-          axis_(rhs.axis_),
           smat_(rhs.smat_),
-          tpca_(
-                rhs.closestApproach().particleTraj(),
-                axis_,
-                rhs.closestApproach().hint(),
-                rhs.closestApproach().precision()
-          ),
+          tpca_(rhs.tpca_),
           toff_(rhs.toff_),
           sxconfig_(rhs.sxconfig_),
           varscale_(rhs.varscale_),
@@ -60,12 +55,12 @@ namespace KinKal {
       bool active() const override { return mxings_.size() > 0; }
       // accessors
       auto const& closestApproach() const { return tpca_; }
+      auto axis() const { return tpca_.sensorTrajPtr(); }
       auto const& strawMaterial() const { return smat_; }
       auto const& config() const { return sxconfig_; }
       auto precision() const { return tpca_.precision(); }
 
     private:
-      SensorLine axis_; // straw axis, expressed as a timeline
       StrawMaterial const& smat_;
       CA tpca_; // result of most recent TPOCA
       double toff_; // small time offset
@@ -79,17 +74,16 @@ namespace KinKal {
   };
 
   template <class KTRAJ> StrawXing<KTRAJ>::StrawXing(PCA const& pca, StrawMaterial const& smat) :
-    axis_(pca.sensorTraj()),
     smat_(smat),
-    tpca_(pca.localTraj(),axis_,pca.precision(),pca.tpData(),pca.dDdP(),pca.dTdP()),
+    tpca_(static_cast<CA const&>(pca)),
     toff_(smat.wireRadius()/pca.particleTraj().speed(pca.particleToca())), // locate the effect to 1 side of the wire to avoid overlap with hits
     varscale_(1.0)
   {}
 
   template <class KTRAJ> void StrawXing<KTRAJ>::updateReference(PTRAJ const& ptraj) {
-    CAHint tphint = tpca_.usable() ?  tpca_.hint() : CAHint(axis_.timeAtMidpoint(),axis_.timeAtMidpoint());
-    PCA pca(ptraj,axis_,tphint,precision());
-    tpca_ = pca.localClosestApproach();
+    CAHint tphint = tpca_.usable() ?  tpca_.hint() : CAHint(axis()->timeAtMidpoint(),axis()->timeAtMidpoint());
+    PCA pca(ptraj,axis(),tphint,precision());
+    tpca_ = static_cast<CA>(pca);
     if(!tpca_.usable())throw std::runtime_error("StrawXing TPOCA failure");
   }
 
@@ -131,7 +125,7 @@ namespace KinKal {
     }
     if(detail > 1){
       ost << " Axis ";
-      axis_.print(ost,0);
+      axis()->print(ost,0);
     }
     ost << std::endl;
   }