GsfTools/interface/GsfMatrixTools.h

0001 #ifndef GsfMatrixTools_h_
0002 #define GsfMatrixTools_h_
0003
0004 #define SMATRIX_USE_CONSTEXPR
0005
0006 #include "Math/SVector.h"
0007 #include "Math/SMatrix.h"
0008
0009 namespace GsfMatrixTools {
0010
0011   //   template <unsigned int N>
0012   //   double similarity (ROOT::Math::SMatrix<double, N, N, ROOT::Math::MatRepSym<double, N> >,
0013   //             ROOT::Math::SVector<double, N> vector) {
0014   //     return vector*matrix*vector;
0015   //   }
0016
0017   template <typename T, unsigned int N>
0018   T trace(ROOT::Math::SMatrix<T, N, N> const& matrix) {
0019     T result = 0;
0020     for (unsigned int i = 0; i < N; ++i)
0021       result += matrix(i, i);
0022     return result;
0023   }
0024
0025   /* compute the trace of a product of two sym matrices
0026    *   a(i,j)*b(j,i) = a(i,j)*b(i,j) sum over i and j
0027    */
0028   template <typename T, unsigned int N>
0029   T trace(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& a,
0030           ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& b) {
0031     auto i1 = a.begin();
0032     auto e1 = a.end();
0033     auto i2 = b.begin();
0034
0035     T res = 0;
0036     // sum of the lower triangle;
0037     for (; i1 != e1; i1++, i2++)
0038       res += (*i1) * (*i2);
0039     res *= T(2.);
0040     // remove the duplicated diagonal...
0041     for (unsigned int i = 0; i < N; ++i)
0042       res -= a(i, i) * b(i, i);
0043     return res;
0044   }
0045
0046 }  // namespace GsfMatrixTools
0047 #endif