ProjectMatrix.h

Go to the documentation of this file.

#ifndef DataFormat_Math_ProjectMatrix_H
#define DataFormat_Math_ProjectMatrix_H

#define SMATRIX_USE_CONSTEXPR
#include "Math/SMatrix.h"

template <typename T, unsigned int N, unsigned int D>
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) const {
    SMatND r;
    for (unsigned int i = 0; i < D; i++)
      for (unsigned int j = 0; j < D; j++)
        r(index[i], j) = s(i, j);
    return r;
  }

  // K*H
  SMatNN project(SMatND const& k) const {
    SMatNN s;
    for (unsigned int i = 0; i < N; i++)
      for (unsigned int j = 0; j < D; j++)
        s(i, index[j]) = k(i, j);
    return s;
  }

  // S-K*H
  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.
  unsigned int index[D];
};

#endif