Merge pull request #4447 from VinInn/KUcleanup

KfComponentsHolder cleanup

Alignment/ReferenceTrajectories/src/ReferenceTrajectory.cc

@@ -955,18 +955,17 @@ ReferenceTrajectory::getHitProjectionMatrixT
 {
   // define variables that will be used to setup the KfComponentsHolder
   // (their allocated memory is needed by 'hitPtr->getKfComponents(..)'
-  // ProjectMatrix<double,5,N>  pf; // not needed
-  typename AlgebraicROOTObject<N,5>::Matrix H; 
+
+  ProjectMatrix<double,5,N>  pf; 
   typename AlgebraicROOTObject<N>::Vector r, rMeas; 
   typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas;
   // input for the holder - but dummy is OK here to just get the projection matrix:
   const AlgebraicVector5 dummyPars;
   const AlgebraicSymMatrix55 dummyErr;
-  ProjectMatrix<double,5,N> dummyProjFunc;
 
   // setup the holder with the correct dimensions and get the values
   KfComponentsHolder holder;
-  holder.setup<N>(&r, &V, &H, &dummyProjFunc, &rMeas, &VMeas, dummyPars, dummyErr);
+  holder.setup<N>(&r, &V, &pf, &rMeas, &VMeas, dummyPars, dummyErr);
   hitPtr->getKfComponents(holder);
 
   return asHepMatrix<N,5>(holder.projection<N>());

DataFormats/Math/interface/ProjectMatrix.h

@@ -11,11 +11,30 @@ struct ProjectMatrix{
   typedef  ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> > SMatDD;
   typedef  ROOT::Math::SMatrix<T,N,N > SMatNN;
   typedef  ROOT::Math::SMatrix<T,N,D > SMatND;
+  typedef  ROOT::Math::SMatrix<T,D,N > SMatDN;
 
+
   // no constructor
 
+  // H
+  SMatDN matrix() const {
+    SMatDN r;
+    for (unsigned int i=0; i<D; i++)
+        r(i,index[i]) = T(1.);
+    return r;
+  }
+
+  // constuct given H
+  void fromH(SMatDN const & H) {
+    for (unsigned int i=0; i<D; i++) {
+      index[i]=N;
+      for (unsigned int j=0; j<N; j++)
+        if ( H(i,j)>0 ) { index[i]=j; break;}
+   }
+  }
+
   // H*S
-  SMatND project(SMatDD const & s) {
+  SMatND project(SMatDD const & s) const {
     SMatND r;
     for (unsigned int i=0; i<D; i++)
       for (unsigned int j=0; j<D; j++)
@@ -24,7 +43,7 @@ struct ProjectMatrix{
   }
 
   // K*H
-  SMatNN project(SMatND const & k) {
+  SMatNN project(SMatND const & k) const {
     SMatNN s;
     for (unsigned int i=0; i<N; i++)
       for (unsigned int j=0; j<D; j++)
@@ -33,13 +52,13 @@ struct ProjectMatrix{
   }
 
   // S-K*H
-  void projectAndSubtractFrom(SMatNN & __restrict__ s, SMatND const & __restrict__ k) {
+  void projectAndSubtractFrom(SMatNN & __restrict__ s, SMatND const & __restrict__ k) const {
     for (unsigned int i=0; i<N; i++)
       for (unsigned int j=0; j<D; j++)
 	s(i,index[j]) -= k(i,j);
   }
 
-  // only H(i,index(i))=1.
+  // only H(i,index[i])=1.
   unsigned int index[D];
 
 };

DataFormats/Math/test/ProjectMatrix_t.cpp

@@ -37,6 +37,7 @@ int main() {
   SMatDD S(v,3);
 
   std::cout << S << std::endl;
+  std::cout << std::endl;
 
 
   {
@@ -45,8 +46,10 @@ int main() {
 
     SMatNN V = K*H; 
 
+    std::cout << H << std::endl;
     std::cout << K << std::endl;
     std::cout << V << std::endl;
+    std::cout << std::endl;
 
   }
   {
@@ -55,13 +58,28 @@ int main() {
     SMatND K = H.project(S);  
 
     SMatNN V = H.project(K); 
-
+
+    std::cout << H.matrix() << std::endl;
     std::cout << K << std::endl;
     std::cout << V << std::endl;
+    std::cout << std::endl;
 
   }
 
+  {
+
+    SMatDN HH; HH(0,3)=1; HH(1,4)=1;
+    ProjectMatrix<double,5,2> H; H.fromH(HH);
+    SMatND K = H.project(S);                                                         
+
+    SMatNN V = H.project(K);
 
+    std::cout << H.matrix() << std::endl;
+    std::cout << K << std::endl;
+    std::cout << V << std::endl;
+    std::cout << std::endl;
+
+  }
 
 
   return 0;

DataFormats/TrackerRecHit2D/src/BaseTrackerRecHit.cc

@@ -45,11 +45,11 @@ BaseTrackerRecHit::getKfComponents1D( KfComponentsHolder & holder ) const
   AlgebraicSymMatrix11 & errs = holder.errors<1>();
   errs(0,0) = err_.xx();
 
-  AlgebraicMatrix15 & proj = holder.projection<1>();
-  proj(0,3) = 1;
-  
-   holder.measuredParams<1>() = AlgebraicVector1( holder.tsosLocalParameters().At(3) );
-   holder.measuredErrors<1>() = holder.tsosLocalErrors().Sub<AlgebraicSymMatrix11>( 3, 3 );
+  ProjectMatrix<double,5,1>  & pf = holder.projFunc<1>();
+  pf.index[0] = 3;
+
+  holder.measuredParams<1>() = AlgebraicVector1( holder.tsosLocalParameters().At(3) );
+  holder.measuredErrors<1>() = holder.tsosLocalErrors().Sub<AlgebraicSymMatrix11>( 3, 3 );
 }
 
 void
@@ -65,15 +65,9 @@ BaseTrackerRecHit::getKfComponents2D( KfComponentsHolder & holder ) const
    errs(0,1) = err_.xy();
    errs(1,1) = err_.yy();
 
-
-   AlgebraicMatrix25 & proj = holder.projection<2>();
-   proj(0,3) = 1;
-   proj(1,4) = 1;
-
    ProjectMatrix<double,5,2>  & pf = holder.projFunc<2>();
    pf.index[0] = 3;
    pf.index[1] = 4;
-   holder.doUseProjFunc();
 
    holder.measuredParams<2>() = AlgebraicVector2( & holder.tsosLocalParameters().At(3), 2 );
    holder.measuredErrors<2>() = holder.tsosLocalErrors().Sub<AlgebraicSymMatrix22>( 3, 3 );

DataFormats/TrackingRecHit/interface/KfComponentsHolder.h

@@ -4,26 +4,24 @@
 
 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
 #include "DataFormats/Math/interface/ProjectMatrix.h"
-#include <stdint.h>
+#include <cstdint>
 
 class TrackingRecHit;
 
-#define Debug_KfComponentsHolder
+// #define Debug_KfComponentsHolder
 
 class KfComponentsHolder {
     public:
-  KfComponentsHolder() : params_(0), errors_(0), projection_(0), useProjFunc_(false)
-        {
+  KfComponentsHolder(){
 #ifdef Debug_KfComponentsHolder
-            size_ = 0;
+   size_ = 0;
 #endif
-        }
+}
 
         template <unsigned int D>
         void setup(
             typename AlgebraicROOTObject<D>::Vector       *params,
             typename AlgebraicROOTObject<D,D>::SymMatrix  *errors,
-            typename AlgebraicROOTObject<D,5>::Matrix     *projection,
 	    ProjectMatrix<double,5,D>                     *projFunc, 
             typename AlgebraicROOTObject<D>::Vector       *measuredParams,
             typename AlgebraicROOTObject<D,D>::SymMatrix  *measuredErrors,
@@ -32,19 +30,6 @@ class KfComponentsHolder {
         ) ;
 
 
-  // backward compatible call
-        template <unsigned int D>
-        void setup(
-            typename AlgebraicROOTObject<D>::Vector       *params,
-            typename AlgebraicROOTObject<D,D>::SymMatrix  *errors,
-            typename AlgebraicROOTObject<D,5>::Matrix     *projection,
-            typename AlgebraicROOTObject<D>::Vector       *measuredParams,
-            typename AlgebraicROOTObject<D,D>::SymMatrix  *measuredErrors,
-            const AlgebraicVector5 & tsosLocalParameters, 
-            const AlgebraicSymMatrix55 & tsosLocalErrors 
-        ) ;
-
-
         template <unsigned int D>
         typename AlgebraicROOTObject<D>::Vector & params() { 
 #ifdef Debug_KfComponentsHolder
@@ -61,14 +46,18 @@ class KfComponentsHolder {
             return  * reinterpret_cast<typename AlgebraicROOTObject<D,D>::SymMatrix *>(errors_);
         }
 
+
         template <unsigned int D>
-        typename AlgebraicROOTObject<D,5>::Matrix & projection() { 
+        typename AlgebraicROOTObject<D,5>::Matrix  projection() { 
 #ifdef Debug_KfComponentsHolder
             assert(size_ == D);
 #endif
-            return  * reinterpret_cast<typename AlgebraicROOTObject<D,5>::Matrix *>(projection_);
+            return this->projFunc<D>().matrix();
+            //return  * reinterpret_cast<typename AlgebraicROOTObject<D,5>::Matrix *>(projection_);
         }
 
+
+
         template <unsigned int D>
         ProjectMatrix<double,5,D> & projFunc() { 
 #ifdef Debug_KfComponentsHolder
@@ -99,23 +88,18 @@ class KfComponentsHolder {
         const AlgebraicVector5     & tsosLocalParameters() const { return *tsosLocalParameters_; }
         const AlgebraicSymMatrix55 & tsosLocalErrors()     const { return *tsosLocalErrors_;     }
 
-        bool useProjFunc() const { return useProjFunc_;}
-        void doUseProjFunc() {  useProjFunc_ = true; }
 
         template<unsigned int D> void dump() ;
     private:
 #ifdef Debug_KfComponentsHolder
         uint16_t size_;
 #endif
-  void *params_, *errors_, *projection_, *projFunc_, *measuredParams_, *measuredErrors_;
+  void *params_, *errors_, *projFunc_, *measuredParams_, *measuredErrors_;
   const AlgebraicVector5 * tsosLocalParameters_;
   const AlgebraicSymMatrix55 * tsosLocalErrors_;
 
-  bool useProjFunc_;
-
   template<unsigned int D>
   void genericFill_(const TrackingRecHit &hit);
-
 
 };
 
@@ -124,7 +108,6 @@ template<unsigned int D>
 void KfComponentsHolder::setup(
         typename AlgebraicROOTObject<D>::Vector       *params,
         typename AlgebraicROOTObject<D,D>::SymMatrix  *errors,
-        typename AlgebraicROOTObject<D,5>::Matrix     *projection,
 	ProjectMatrix<double,5,D>                     *projFunc, 
         typename AlgebraicROOTObject<D>::Vector       *measuredParams,
         typename AlgebraicROOTObject<D,D>::SymMatrix  *measuredErrors,
@@ -137,7 +120,6 @@ void KfComponentsHolder::setup(
 #endif
     params_     = params;
     errors_     = errors;
-    projection_ = projection;
     projFunc_ = projFunc;
     measuredParams_ = measuredParams;
     measuredErrors_ = measuredErrors;
@@ -163,3 +145,4 @@ void KfComponentsHolder::dump() {
 #endif 
 
 #endif
+

DataFormats/TrackingRecHit/src/KfComponentsHolder.cc

@@ -11,9 +11,9 @@ void KfComponentsHolder::genericFill_(const TrackingRecHit &hit) {
 
    params<D>()     = asSVector<D>(hit.parameters());
    errors<D>()     = asSMatrix<D>(hit.parametersError());
-   projection<D>() = asSMatrix<D,5>(hit.projectionMatrix());
 
-   const MatD5 & H = projection<D>();
+   MatD5 && H = asSMatrix<D,5>(hit.projectionMatrix());
+   projFunc<D>().fromH(H);
 
    measuredParams<D>() = H * (*tsosLocalParameters_);
    measuredErrors<D>() = ROOT::Math::Similarity(H, (*tsosLocalErrors_));

RecoTracker/SiTrackerMRHTools/src/SiTrackerMultiRecHitUpdator.cc

@@ -192,15 +192,14 @@ double SiTrackerMultiRecHitUpdator::ComputeWeight(const TrajectoryStateOnSurface
 
   // define variables that will be used to setup the KfComponentsHolder
   ProjectMatrix<double,5,N>  pf;
-  typename AlgebraicROOTObject<N,5>::Matrix H;
   typename AlgebraicROOTObject<N>::Vector r, rMeas;
   typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas, W;
   AlgebraicVector5 x = tsos.localParameters().vector();
   const AlgebraicSymMatrix55 &C = (tsos.localError().matrix());
 
   // setup the holder with the correct dimensions and get the values
   KfComponentsHolder holder;
-  holder.template setup<N>(&r, &V, &H, &pf, &rMeas, &VMeas, x, C);
+  holder.template setup<N>(&r, &V,  &pf, &rMeas, &VMeas, x, C);
   aRecHit.getKfComponents(holder);
 
   typename AlgebraicROOTObject<N>::Vector diff;
@@ -289,15 +288,14 @@ std::pair<AlgebraicVector2,AlgebraicSymMatrix22> SiTrackerMultiRecHitUpdator::Co
 
   // define variables that will be used to setup the KfComponentsHolder
   ProjectMatrix<double,5,N>  pf;
-  typename AlgebraicROOTObject<N,5>::Matrix H;
   typename AlgebraicROOTObject<N>::Vector r, rMeas;
   typename AlgebraicROOTObject<N,N>::SymMatrix V, VMeas, Wtemp;
   AlgebraicVector5 x = tsos.localParameters().vector();
   const AlgebraicSymMatrix55 &C = (tsos.localError().matrix());
 
   // setup the holder with the correct dimensions and get the values
   KfComponentsHolder holder;
-  holder.template setup<N>(&r, &V, &H, &pf, &rMeas, &VMeas, x, C);
+  holder.template setup<N>(&r, &V, &pf, &rMeas, &VMeas, x, C);
   aRecHit.getKfComponents(holder);
 
   Wtemp = V;

TrackingTools/GsfTracking/src/PosteriorWeightsCalculator.cc

@@ -23,7 +23,8 @@ std::vector<double> PosteriorWeightsCalculator::weights(const TrackingRecHit& re
   typedef typename AlgebraicROOTObject<5,D>::Matrix Mat5D;
   typedef typename AlgebraicROOTObject<D,D>::SymMatrix SMatDD;
   typedef typename AlgebraicROOTObject<D>::Vector VecD;
-
+  using ROOT::Math::SMatrixNoInit;
+
   std::vector<double> weights;
   if ( predictedComponents.empty() )  {
     edm::LogError("EmptyPredictedComponents")<<"a multi state is empty. cannot compute any weight.";
@@ -36,9 +37,8 @@ std::vector<double> PosteriorWeightsCalculator::weights(const TrackingRecHit& re
   std::vector<double> chi2s;
   chi2s.reserve(predictedComponents.size());
 
-  MatD5 H; 
   VecD r, rMeas; 
-  SMatDD V, R;
+  SMatDD V(SMatrixNoInit{}), R(SMatrixNoInit{});
   AlgebraicVector5 x;
   ProjectMatrix<double,5,D> p;
   //
@@ -47,18 +47,10 @@ std::vector<double> PosteriorWeightsCalculator::weights(const TrackingRecHit& re
   //  
   double chi2Min(DBL_MAX);
   for ( unsigned int i=0; i<predictedComponents.size(); i++ ) {
-//     MeasurementExtractor me(predictedComponents[i]);
-//     // Residuals of aPredictedState w.r.t. aRecHit, 
-//     //!!!     AlgebraicVector r(recHit.parameters(predictedComponents[i]) - me.measuredParameters(recHit));
-//     VecD r = asSVector<D>(recHit.parameters()) - me.measuredParameters<D>(recHit);
-//     // and covariance matrix of residuals
-//     //!!!     AlgebraicSymMatrix V(recHit.parametersError(predictedComponents[i]));
-//     SMatDD V = asSMatrix<D>(recHit.parametersError());
-//     SMatDD R = V + me.measuredError<D>(recHit);
 
     KfComponentsHolder holder; 
     x = predictedComponents[i].localParameters().vector();
-    holder.template setup<D>(&r, &V, &H, &p, &rMeas, &R, 
+    holder.template setup<D>(&r, &V, &p, &rMeas, &R, 
 			     x, predictedComponents[i].localError().matrix());
     recHit.getKfComponents(holder);
 

TrackingTools/KalmanUpdators/BuildFile.xml

@@ -1,3 +1,4 @@
+#<flags   CXXFLAGS="-Ofast -mrecip -funroll-all-loops"/>
 <flags   CXXFLAGS="-O3"/>
 <use   name="boost"/>
 <use   name="clhep"/>

TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h

@@ -25,8 +25,6 @@ class Chi2MeasurementEstimator : public Chi2MeasurementEstimatorBase {
 
   virtual std::pair<bool,double> estimate(const TrajectoryStateOnSurface&,
 				     const TrackingRecHit&) const;
-  template <unsigned int D> std::pair<bool,double> estimate(const TrajectoryStateOnSurface&,
-				     const TrackingRecHit&) const;
 
   virtual Chi2MeasurementEstimator* clone() const {
     return new Chi2MeasurementEstimator(*this);

TrackingTools/KalmanUpdators/interface/KFUpdator.h

@@ -40,8 +40,6 @@ class KFUpdator GCC11_FINAL : public TrajectoryStateUpdator {
   TrajectoryStateOnSurface update(const TrajectoryStateOnSurface&,
                                   const TrackingRecHit&) const;
 
-  template <unsigned int D> TrajectoryStateOnSurface update(const TrajectoryStateOnSurface&,
-                                  const TrackingRecHit&) const;
 
   virtual KFUpdator * clone() const {
     return new KFUpdator(*this);