Merge pull request #205 from edcallaghan/track_copy_constructor

Implement Track<> copy-constructor and backing graph-copy infrastructure

Author: brownd1978

Detector/ElementXing.hh

@@ -4,13 +4,16 @@
 //  Describe the material effects of a particle crossing a detector element (piece of the detector)
 //  Used in the kinematic Kalman fit
 //
+#include "KinKal/General/CloneContext.hh"
 #include "KinKal/General/MomBasis.hh"
 #include "KinKal/Detector/MaterialXing.hh"
 #include "KinKal/Trajectory/ParticleTrajectory.hh"
 #include "KinKal/Fit/MetaIterConfig.hh"
 #include <vector>
 #include <array>
 #include <ostream>
+#include <stdexcept>
+#include <string>
 
 namespace KinKal {
   template <class KTRAJ> class ElementXing {
@@ -19,6 +22,9 @@ namespace KinKal {
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       ElementXing() {}
       virtual ~ElementXing() {}
+      // clone op for reinstantiation
+      ElementXing(ElementXing const& rhs) = default;
+      virtual std::shared_ptr< ElementXing<KTRAJ> > clone(CloneContext&) const;
       virtual void updateReference(PTRAJ const& ptraj) = 0; // update the trajectory reference
       virtual void updateState(MetaIterConfig const& config,bool first) =0; // update the state according to this meta-config
       virtual Parameters params() const =0; // parameter change induced by this element crossing WRT the reference parameters going forwards in time
@@ -40,6 +46,14 @@ namespace KinKal {
     private:
   };
 
+  // cloning requires domain knowledge of pointer members of the cloned object,
+  // which must be reassigned explicitly; the default action is thus to throw
+  // an error if a clone routine has not been explicitly provided.
+  template<class KTRAJ> std::shared_ptr< ElementXing<KTRAJ> > ElementXing<KTRAJ>::clone(CloneContext&) const{
+    std::string msg = "Attempt to clone KinKal::Hit subclass with no clone() implementation";
+    throw std::runtime_error(msg);
+  }
+
   template <class KTRAJ> void ElementXing<KTRAJ>::momentumChange(SVEC3& dmom, SMAT& dmomvar) const {
     // compute the parameter effect for forwards time
     double dm, paramomvar, perpmomvar;

Detector/Hit.hh

@@ -4,13 +4,16 @@
 //  Base class to describe a measurement that constrains some aspect of the track fit
 //  The constraint is expressed as a weight WRT a set of reference parameters.
 //
+#include "KinKal/General/CloneContext.hh"
 #include "KinKal/General/Weights.hh"
 #include "KinKal/General/Parameters.hh"
 #include "KinKal/General/Chisq.hh"
 #include "KinKal/Trajectory/ParticleTrajectory.hh"
 #include "KinKal/Fit/MetaIterConfig.hh"
 #include <memory>
 #include <ostream>
+#include <stdexcept>
+#include <string>
 
 namespace KinKal {
   template <class KTRAJ> class Hit {
@@ -22,6 +25,8 @@ namespace KinKal {
       // disallow copy and equivalence
       Hit(Hit const& ) = delete;
       Hit& operator =(Hit const& ) = delete;
+      // clone op for reinstantiation
+      virtual std::shared_ptr< Hit<KTRAJ> > clone(CloneContext&) const;
      // hits may be active (used in the fit) or inactive; this is a pattern recognition feature
       virtual bool active() const =0;
       virtual unsigned nDOF() const=0;
@@ -45,6 +50,14 @@ namespace KinKal {
       Chisq chisquared() const;
   };
 
+  // cloning requires domain knowledge of pointer members of the cloned object,
+  // which must be reassigned explicitly; the default action is thus to throw
+  // an error if a clone routine has not been explicitly provided.
+  template<class KTRAJ> std::shared_ptr< Hit<KTRAJ> > Hit<KTRAJ>::clone(CloneContext& context) const{
+    std::string msg = "Attempt to clone KinKal::Hit subclass with no clone() implementation";
+    throw std::runtime_error(msg);
+  }
+
   template<class KTRAJ> Parameters Hit<KTRAJ>::unbiasedParameters() const {
     if(active()){
       // convert the parameters to a weight, and subtract this hit's weight

Detector/ParameterHit.hh

@@ -15,6 +15,22 @@ namespace KinKal {
       using PTRAJ = ParticleTrajectory<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
 
+      // copy constructor
+      ParameterHit(ParameterHit<KTRAJ> const& rhs):
+          time_(rhs.time()),
+          params_(rhs.constraintParameters()),
+          pweight_(rhs.parameterWeights()),
+          weight_(rhs.weight()),
+          pmask_(rhs.constraintMask()),
+          mask_(rhs.maskMatrix()),
+          ncons_(rhs.numConstrainedParameters()){
+        /**/
+      };
+      // clone op for reinstantiation
+      std::shared_ptr< Hit<KTRAJ> > clone(CloneContext& context) const override{
+        auto rv = std::make_shared< ParameterHit<KTRAJ> >(*this);
+        return rv;
+      };
       // Hit interface overrrides
       bool active() const override { return ncons_ > 0; }
       Chisq chisq(Parameters const& pdata) const override;
@@ -32,6 +48,9 @@ namespace KinKal {
       unsigned nDOF() const override { return ncons_; }
       Parameters const& constraintParameters() const { return params_; }
       PMASK const& constraintMask() const { return pmask_; }
+      Weights parameterWeights() const { return pweight_; };
+      DMAT maskMatrix() const { return mask_; }
+      unsigned numConstrainedParameters() const { return ncons_; };
     private:
       double time_; // time of this constraint: must be supplied on construction and does not change
       KTRAJPTR reftraj_; // reference WRT this hits weight was calculated

Detector/ResidualHit.hh

@@ -10,6 +10,10 @@ namespace KinKal {
 
   template <class KTRAJ> class ResidualHit : public Hit<KTRAJ> {
     public:
+      // copy constructor
+      ResidualHit(ResidualHit<KTRAJ> const& rhs): weight_(rhs.weight()){
+        /**/
+      };
       // override of some Hit interface.  Subclasses must still implement update and material methods
       using HIT = Hit<KTRAJ>;
       using PTRAJ = ParticleTrajectory<KTRAJ>;

Examples/ScintHit.hh

@@ -17,6 +17,30 @@ namespace KinKal {
       using RESIDHIT = ResidualHit<KTRAJ>;
       using HIT = Hit<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
+      // 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()
+          ),
+          rresid_(rhs.refResidual()){
+        /**/
+      };
+      // clone op for reinstantiation
+      std::shared_ptr< Hit<KTRAJ> > clone(CloneContext& context) const override{
+        auto rv = std::make_shared< ScintHit<KTRAJ> >(*this);
+        auto ca = rv->closestApproach();
+        auto trajectory = std::make_shared<KTRAJ>(ca.particleTraj());
+        ca.setTrajectory(trajectory);
+        rv->setClosestApproach(ca);
+        return rv;
+      };
       // ResidualHit interface implementation
       unsigned nResid() const override { return 1; } // 1 time residual
       Residual const& refResidual(unsigned ires=0) const override;
@@ -40,6 +64,9 @@ namespace KinKal {
       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
       Residual rresid_; // residual WRT most recent reference parameters
+
+      // modifiers to support cloning
+      void setClosestApproach(const CA& ca){ tpca_ = ca; }
   };
 
   template <class KTRAJ> ScintHit<KTRAJ>::ScintHit(PCA const& pca, double tvar, double wvar) :

Examples/SimpleWireHit.hh

@@ -27,6 +27,36 @@ namespace KinKal {
       SimpleWireHit(BFieldMap const& bfield, PCA const& pca, WireHitState const& whstate, double mindoca,
           double driftspeed, double tvar, double tot, double totvar, double rcell,int id);
       unsigned nResid() const override { return 2; } // 2 residuals
+      // clone op for reinstantiation
+      SimpleWireHit(SimpleWireHit<KTRAJ> const& rhs):
+          ResidualHit<KTRAJ>(rhs),
+          bfield_(rhs.bfield()),
+          whstate_(rhs.hitState()),
+          wire_(rhs.wire()),
+          ca_(
+            rhs.closestApproach().particleTraj(),
+            wire_,
+            rhs.closestApproach().hint(),
+            rhs.closestApproach().precision()
+          ),
+          rresid_(rhs.residuals()),
+          mindoca_(rhs.minDOCA()),
+          dvel_(driftVelocity()),
+          tvar_(timeVariance()),
+          tot_(rhs.tot()),
+          totvar_(rhs.totVariance()),
+          rcell_(rhs.cellRadius()),
+          id_(rhs.id()){
+        /**/
+      };
+      std::shared_ptr< Hit<KTRAJ> > clone(CloneContext& context) const override{
+        auto rv = std::make_shared< SimpleWireHit<KTRAJ> >(*this);
+        auto ca = rv->closestApproach();
+        auto trajectory = std::make_shared<KTRAJ>(ca.particleTraj());
+        ca.setTrajectory(trajectory);
+        rv->setClosestApproach(ca);
+        return rv;
+      };
       double time() const override { return ca_.particleToca(); }
       VEC3 dRdX(unsigned ires) const;
       Residual const& refResidual(unsigned ires=dresid) const override;
@@ -48,6 +78,10 @@ namespace KinKal {
       auto const& wire() const { return wire_; }
       auto const& bfield() const { return bfield_; }
       auto precision() const { return ca_.precision(); }
+      auto const& residuals() const { return rresid_; }
+      double tot() const { return tot_; }
+      double totVariance() const { return totvar_; }
+
     private:
       BFieldMap const& bfield_; // drift calculation requires the BField for ExB effects
       WireHitState whstate_; // current state
@@ -65,6 +99,9 @@ namespace KinKal {
       double rcell_; // straw radius
       int id_; // id
       void updateResiduals();
+
+      // modifiers to support cloning
+      void setClosestApproach(const CA& ca){ ca_ = ca; }
   };
 
   //trivial 'updater' that sets the wire hit state to null

Examples/StrawXing.hh

@@ -21,6 +21,33 @@ namespace KinKal {
       // construct from PCA and material
       StrawXing(PCA const& pca, StrawMaterial const& smat);
       virtual ~StrawXing() {}
+      // 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()
+          ),
+          toff_(rhs.toff_),
+          sxconfig_(rhs.sxconfig_),
+          varscale_(rhs.varscale_),
+          mxings_(rhs.mxings_),
+          fparams_(rhs.fparams_){
+        /**/
+      }
+      // clone op for reinstantiation
+      std::shared_ptr< ElementXing<KTRAJ> > clone(CloneContext& context) const override{
+        auto rv = std::make_shared< StrawXing<KTRAJ> >(*this);
+        auto ca = rv->closestApproach();
+        auto trajectory = std::make_shared<KTRAJ>(ca.particleTraj());
+        ca.setTrajectory(trajectory);
+        rv->setClosestApproach(ca);
+        return rv;
+      }
       // ElementXing interface
       void updateReference(PTRAJ const& ptraj) override;
       void updateState(MetaIterConfig const& config,bool first) override;
@@ -35,6 +62,7 @@ namespace KinKal {
       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_;
@@ -44,6 +72,9 @@ namespace KinKal {
       double varscale_; // variance scale
       std::vector<MaterialXing> mxings_;
       Parameters fparams_; // parameter change for forwards time
+
+      // modifiers to support cloning
+      void setClosestApproach(const CA& ca){ tpca_ = ca; }
   };
 
   template <class KTRAJ> StrawXing<KTRAJ>::StrawXing(PCA const& pca, StrawMaterial const& smat) :

Fit/Domain.hh

@@ -4,13 +4,18 @@
 //  domain used to compute magnetic field corrections. Magnetic bending not described by the intrinsic parameterization is assumed
 //  to be negligible over the domain
 //
+#include <memory>
+#include "KinKal/General/CloneContext.hh"
 #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) {}
+      // clone op for reinstantiation
+      Domain(Domain const&);
+      std::shared_ptr< Domain > clone(CloneContext&) const;
       bool operator < (Domain const& other) const {return begin() < other.begin(); }
       auto const& range() const { return range_; }
       // forward range functions
@@ -25,6 +30,19 @@ namespace KinKal {
       VEC3 bnom_; // nominal BField for this domain
       double tol_; // tolerance used to create this domain
   };
+
+  // clone op for reinstantiation
+  Domain::Domain(Domain const& rhs):
+      range_(rhs.range()),
+      bnom_(rhs.bnom()),
+      tol_(rhs.tolerance()){
+    /**/
+  }
+
+  std::shared_ptr<Domain> Domain::clone(CloneContext& context) const{
+    auto rv = std::make_shared<Domain>(*this);
+    return rv;
+  }
 }
 #endif
 

Fit/DomainWall.hh

@@ -34,20 +34,33 @@ namespace KinKal {
       auto const& parameterChange() const { return dpfwd_; }
       virtual ~DomainWall(){}
       // disallow copy and equivalence
-      DomainWall(DomainWall const& ) = delete;
+//    DomainWall(DomainWall const& ) = delete;
       DomainWall& operator =(DomainWall const& ) = delete;
+      // clone op for reinstantiation
+      DomainWall(DomainWall const&);
+      std::unique_ptr< Effect<KTRAJ> > clone(CloneContext&) const override;
       // specific DomainWall interface
       DomainWall(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_; }
+      // other accessors
+      auto const& prevPtr() const { return prev_; }
+      auto const& nextPtr() const { return next_; }
+      auto const& fwdChange() const { return dpfwd_; }
+      auto const& prevWeight() const { return prevwt_; }
+      auto const& nextWeight() const { return nextwt_; }
+      auto const& fwdCovarianceRotation() const { return dpdpdb_; }
 
     private:
       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
 
+      // modifiers to support cloning
+      void setPrevPtr(DOMAINPTR const& ptr){ prev_ = ptr; }
+      void setNextPtr(DOMAINPTR const& ptr){ next_ = ptr; }
   };
 
   template<class KTRAJ> DomainWall<KTRAJ>::DomainWall(
@@ -132,5 +145,26 @@ namespace KinKal {
     return ost;
   }
 
+  // clone op for reinstantiation
+  template <class KTRAJ>
+  DomainWall<KTRAJ>::DomainWall(DomainWall const& rhs):
+      dpfwd_(rhs.fwdChange()),
+      prevwt_(rhs.prevWeight()),
+      nextwt_(rhs.nextWeight()),
+      dpdpdb_(rhs.fwdCovarianceRotation()){
+    /**/
+  }
+
+  template <class KTRAJ>
+  std::unique_ptr< Effect<KTRAJ> > DomainWall<KTRAJ>::clone(CloneContext& context) const{
+    auto casted = std::make_unique< DomainWall<KTRAJ> >(*this);
+    DOMAINPTR prev = context.get(prev_);
+    casted->setPrevPtr(prev);
+    DOMAINPTR next = context.get(next_);
+    casted->setNextPtr(next);
+    //auto rv = std::make_unique< Effect<KTRAJ> >(casted);
+    auto rv = std::move(casted);
+    return rv;
+  }
 }
 #endif