JEigenValues2D.hh

Go to the documentation of this file.

#ifndef __JEIGENVALUES2D__
#define __JEIGENVALUES2D__
 
#include <cmath>
 
#include "TMatrixDSymEigen.h"
 
#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JGeometry2DToolkit.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JGEOMETRY2D {}
namespace JPP { using namespace JGEOMETRY2D; }
 
namespace JGEOMETRY2D {
 
  /**
   * Eigen values in 2D.
   */
  class JEigenValues2D :
    public JVector2D
  {
  public:
    /**
     * Constructor.
     *
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class T>
    JEigenValues2D(T __begin,
                   T __end) :
      JVector2D(),
      N(0)
    {
      const JCenter2D center(__begin, __end);
        
      // RMS matrix
        
      TMatrixDSym A(2);
        
      for (int i = 0; i != 2; ++i) {
        for (int j = 0; j != 2; ++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();
          
        A(0,0) += (dx * dx);
        A(0,1) += (dx * dy);
        
        A(1,1) += (dy * dy);
      }
      
      A(1,0) = A(0,1);
 
      const TVectorD V = TMatrixDSymEigen(A).GetEigenValues();   // copy!
 
      N = V.GetNoElements();
 
      switch (N) {
 
      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 (n > N)
        return N;
      else
        return n;
    }
 
 
  protected:
    int N;
  };
}
 
#endif