GsfMatrixTools.h

Go to the documentation of this file.

#ifndef GsfMatrixTools_h_
#define GsfMatrixTools_h_
 
#define SMATRIX_USE_CONSTEXPR
 
#include "Math/SVector.h"
#include "Math/SMatrix.h"
 
namespace GsfMatrixTools {
 
  //   template <unsigned int N>
  //   double similarity (ROOT::Math::SMatrix<double, N, N, ROOT::Math::MatRepSym<double, N> >,
  //             ROOT::Math::SVector<double, N> vector) {
  //     return vector*matrix*vector;
  //   }
 
  template <typename T, unsigned int N>
  T trace(ROOT::Math::SMatrix<T, N, N> const& matrix) {
    T result = 0;
    for (unsigned int i = 0; i < N; ++i)
      result += matrix(i, i);
    return result;
  }
 
  /* compute the trace of a product of two sym matrices
   *   a(i,j)*b(j,i) = a(i,j)*b(i,j) sum over i and j
   */
  template <typename T, unsigned int N>
  T trace(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& a,
          ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& b) {
    auto i1 = a.begin();
    auto e1 = a.end();
    auto i2 = b.begin();
 
    T res = 0;
    // sum of the lower triangle;
    for (; i1 != e1; i1++, i2++)
      res += (*i1) * (*i2);
    res *= T(2.);
    // remove the duplicated diagonal...
    for (unsigned int i = 0; i < N; ++i)
      res -= a(i, i) * b(i, i);
    return res;
  }
 
}  // namespace GsfMatrixTools
#endif