JEigenValues3D.hh

Go to the documentation of this file.

#ifndef __JEIGENVALUES3D__
#define __JEIGENVALUES3D__

#include <cmath>

#include "TMatrixDSymEigen.h"

#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"


/**
 * \author mdejong
 */

namespace JGEOMETRY3D {}
namespace JPP { using namespace JGEOMETRY3D; }

namespace JGEOMETRY3D {

  /**
   * Eigen values in 3D.
   */
  class JEigenValues3D :
    public JVector3D
  {
  public:
    /**
     * Constructor.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class T>
    JEigenValues3D(T __begin,
                   T __end) :
      JVector3D(),
      N(0)
    {
      const JCenter3D center(__begin, __end);
        
      // RMS matrix
        
      TMatrixDSym A(3);
        
      for (int i = 0; i != 3; ++i) {
        for (int j = 0; j != 3; ++j) {
          A(i,j) = 0.0;
        }
      }
        
      for (T i = __begin; i != __end; ++i) {
        
        const double dx = center.getX() - i->getX();
        const double dy = center.getY() - i->getY();
        const double dz = center.getZ() - i->getZ();
          
        A(0,0) += (dx * dx);
        A(0,1) += (dx * dy);
        A(0,2) += (dx * dz);
        
        A(1,1) += (dy * dy);
        A(1,2) += (dy * dz);
        
        A(2,2) += (dz * dz);
      }
      
      A(1,0) = A(0,1);
      A(2,0) = A(0,2);
      A(2,1) = A(1,2);

      const TVectorD V = TMatrixDSymEigen(A).GetEigenValues();   // copy!

      N = V.GetNoElements();

      switch (N) {

      case 3:
        __z = V[2];

      case 2:
        __y = V[1];

      case 1:
        __x = V[0];
        break;

      default:
        break;
      }
    }


    /**
     * Get number of eigen values.
     *
     * \return                number of eigen values
     */
    int getN() const
    {
      return N;
    }


    /**
     * Get number of eigen values with minimal value.
     *
     * \param  precision      minimal value
     * \return                number of eigen values
     */
    int getN(double precision) const
    {
      int n = 0;

      if (fabs(getX()) >= precision) ++n;
      if (fabs(getY()) >= precision) ++n;
      if (fabs(getZ()) >= precision) ++n;

      if (n > N)
        return N;
      else
        return n;
    }

  protected:
    int N;
  };
}
 
#endif