Merge pull request #187 from brownd1978/extrap

Cleanups and improvements to support extrapolation

.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:

.gitignore

@@ -5,5 +5,5 @@ KinKal/UnitTests/Dict.cc
 debug/*
 *.root
 .DS_Store
-.vscode
+.vscode/*
 .*.swp

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 <vector>
@@ -25,27 +24,26 @@ 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
       virtual std::vector<MaterialXing>const&  matXings() const =0; // Effect of each physical material component of this detector element on this trajectory
       virtual void print(std::ostream& ost=std::cout,int detail=0) const =0;
       // crossings  without material are inactive
       bool active() const { return matXings().size() > 0; }
-      // calculate the cumulative material effect from these crossings
-      void materialEffects(TimeDir tdir, std::array<double,3>& dmom, std::array<double,3>& momvar) const;
+      // calculate the cumulative material effect from all the materials in this element crossing going forwards in time
+      void materialEffects(std::array<double,3>& dmom, std::array<double,3>& momvar) const;
       // sum radiation fraction
       double radiationFraction() const;
-      static double elossFactor(TimeDir const& tdir) { return tdir == TimeDir::forwards ? 1.0 : -1.0; }
     private:
   };
 
-  template <class KTRAJ> void ElementXing<KTRAJ>::materialEffects(TimeDir tdir, std::array<double,3>& dmom, std::array<double,3>& momvar) const {
+  template <class KTRAJ> void ElementXing<KTRAJ>::materialEffects(std::array<double,3>& dmom, std::array<double,3>& 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

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<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       using EXING = ElementXing<KTRAJ>;
-      using PCA = PiecewiseClosestApproach<KTRAJ,Line>;
-      using CA = ClosestApproach<KTRAJ,Line>;
+      using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
+      using CA = ClosestApproach<KTRAJ,SensorLine>;
       // 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<double,3> 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 <class KTRAJ> Parameters StrawXing<KTRAJ>::parameters(TimeDir tdir) const {
-    if(tdir == TimeDir::forwards)
-      return fparams_;
-    else
-      return Parameters(-fparams_.parameters(),fparams_.covariance());
+  template <class KTRAJ> Parameters StrawXing<KTRAJ>::params() const {
+    return fparams_;
   }
 
   template <class KTRAJ> double StrawXing<KTRAJ>::transitTime() const {

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

Examples/BFieldInfo.hh → Examples/DomainWallInfo.hh

@@ -1,14 +1,14 @@
 
-#ifndef KinKal_BFieldInfo_hh
-#define KinKal_BFieldInfo_hh
+#ifndef KinKal_DomainWallInfo_hh
+#define KinKal_DomainWallInfo_hh
 #include <vector>
 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<BFieldInfo> KKBFIV;
+  typedef std::vector<DomainWallInfo> KKBFIV;
 }
 #endif

Examples/LinkDef.h

@@ -2,8 +2,8 @@
 
 #pragma link C++ class KinKal::HitInfo+;
 #pragma link C++ class vector<KinKal::HitInfo>+;
-#pragma link C++ class KinKal::BFieldInfo+;
-#pragma link C++ class vector<KinKal::BFieldInfo>+;
+#pragma link C++ class KinKal::DomainWallInfo+;
+#pragma link C++ class vector<KinKal::DomainWallInfo>+;
 #pragma link C++ class KinKal::MaterialInfo+;
 #pragma link C++ class vector<KinKal::MaterialInfo>+;
 #pragma link C++ class KinKal::ParticleTrajectoryInfo+;

Examples/ScintHit.hh

@@ -4,16 +4,16 @@
 //  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 <stdexcept>
 namespace KinKal {
 
   template <class KTRAJ> class ScintHit : public ResidualHit<KTRAJ> {
     public:
       using PTRAJ = ParticleTrajectory<KTRAJ>;
-      using PCA = PiecewiseClosestApproach<KTRAJ,Line>;
-      using CA = ClosestApproach<KTRAJ,Line>;
+      using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
+      using CA = ClosestApproach<KTRAJ,SensorLine>;
       using RESIDHIT = ResidualHit<KTRAJ>;
       using HIT = Hit<KTRAJ>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
@@ -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 <class KTRAJ> void ScintHit<KTRAJ>::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

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 KTRAJ> class SimpleWireHit : public ResidualHit<KTRAJ> {
     public:
       using HIT = Hit<KTRAJ>;
-      using PCA = PiecewiseClosestApproach<KTRAJ,Line>;
-      using CA = ClosestApproach<KTRAJ,Line>;
+      using PCA = PiecewiseClosestApproach<KTRAJ,SensorLine>;
+      using CA = ClosestApproach<KTRAJ,SensorLine>;
       using KTRAJPTR = std::shared_ptr<KTRAJ>;
       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 <class KTRAJ> ClosestApproach<KTRAJ,Line> SimpleWireHit<KTRAJ>::unbiasedClosestApproach() const {
+  template <class KTRAJ> ClosestApproach<KTRAJ,SensorLine> SimpleWireHit<KTRAJ>::unbiasedClosestApproach() const {
     // compute the unbiased closest approach; this is brute force, but works
     auto const& ca = this->closestApproach();
     auto uparams = HIT::unbiasedParameters();