JFIT::JMatrixNZ Class Reference

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z). More...

#include <JMatrixNZ.hh>

Inheritance diagram for JFIT::JMatrixNZ:

Classes
struct	variance
	Auxiliary data structure for co-variance calculation. More...

Public Types
typedef std::vector < std::vector< double > >	matrix_type
typedef matrix_type::const_iterator	const_row_type
typedef std::vector< double > ::const_iterator	const_col_type
typedef matrix_type::iterator	row_type
typedef std::vector< double > ::iterator	col_type

Public Member Functions
	JMatrixNZ ()
	Default contructor. More...
template<class T >
	JMatrixNZ (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
	Constructor. More...
template<class T >
void	set (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
	Set co-variance matrix. More...
void	invert ()
	Invert matrix. More...
void	update (const unsigned int k, const double extra)
	Update inverted matrix at given diagonal element. More...
void	resize (const size_t size)
	Resize matrix. More...
double	at (size_t row, size_t col) const
	Get matrix element. More...
double &	at (size_t row, size_t col)
	Get matrix element. More...
double	operator() (size_t row, size_t col) const
	Get matrix element. More...
double &	operator() (size_t row, size_t col)
	Get matrix element. More...
JMatrixND &	reset ()
	Set matrix to the null matrix. More...
JMatrixND &	setIdentity ()
	Set to identity matrix. More...
JMatrixND &	transpose ()
	Transpose. More...
JMatrixND &	negate ()
	Negate matrix. More...
JMatrixND &	add (const JMatrixND &A)
	Matrix addition. More...
JMatrixND &	sub (const JMatrixND &A)
	Matrix subtraction. More...
JMatrixND &	mul (const double factor)
	Scale matrix. More...
JMatrixND &	mul (const JMatrixND &A, const JMatrixND &B)
	Matrix multiplication. More...
JMatrixND &	mul (const JNullType &object)
	Multiply with object. More...
JMatrixND &	div (const double factor)
	Scale matrix. More...
bool	equals (const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
	Equality. More...
bool	isIdentity (const double eps=std::numeric_limits< double >::min()) const
	Test identity. More...
double	getDeterminant () const
	Get determinant of matrix. More...
bool	isPositiveDefinite () const
	Test positive definiteness. More...
const bool	isPositiveSemiDefinite () const
	Test positive semi-definiteness. More...

Static Private Member Functions
static double	getDot (const variance &first, const variance &second)
	Get dot product. More...

Private Attributes
std::vector< variance >	buffer

Detailed Description

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).

In this, the given angular and time resolution are taken into account.

Definition at line 28 of file JMatrixNZ.hh.

Member Typedef Documentation

typedef std::vector< std::vector<double> > JMATH::JMatrixND_t::matrix_type

inherited

Definition at line 44 of file JMatrixND.hh.

typedef matrix_type::const_iterator JMATH::JMatrixND_t::const_row_type

inherited

Definition at line 46 of file JMatrixND.hh.

typedef std::vector<double>::const_iterator JMATH::JMatrixND_t::const_col_type

inherited

Definition at line 47 of file JMatrixND.hh.

typedef matrix_type::iterator JMATH::JMatrixND_t::row_type

inherited

Definition at line 49 of file JMatrixND.hh.

typedef std::vector<double>::iterator JMATH::JMatrixND_t::col_type

inherited

Definition at line 50 of file JMatrixND.hh.

Constructor & Destructor Documentation

JFIT::JMatrixNZ::JMatrixNZ ( )

inline

Default contructor.

Definition at line 35 of file JMatrixNZ.hh.

                :
      JMatrixNS()
    {}

template<class T >

JFIT::JMatrixNZ::JMatrixNZ ( const JVector3D &  pos, T  __begin, T  __end, const double  alpha, const double  sigma  )

inline

Constructor.

The template argument T refers to an iterator of a data structure which should have the following member methods:

• double getX(); // [m]
• double getY(); // [m]
• double getZ(); // [m]

Parameters

pos	reference position [m]
__begin	begin of data
__end	end of data
alpha	angular resolution [deg]
sigma	time resolution [ns]

Definition at line 55 of file JMatrixNZ.hh.

                                      :
      JMatrixNS()
    {
      set(pos, __begin, __end, alpha, sigma);
    }

Member Function Documentation

template<class T >

void JFIT::JMatrixNZ::set ( const JVector3D &  pos, T  __begin, T  __end, const double  alpha, const double  sigma  )

inline

Set co-variance matrix.

The template argument T refers to an iterator of a data structure which should have the following member methods:

• double getX(); // [m]
• double getY(); // [m]
• double getZ(); // [m]

Parameters

pos	reference position [m]
__begin	begin of data
__end	end of data
alpha	angular resolution [deg]
sigma	time resolution [ns]

Definition at line 81 of file JMatrixNZ.hh.

    {
      using namespace std;
      using namespace JTOOLS;


      const int N = distance(__begin, __end);

      this ->resize(N);
      buffer.resize(N);


      const double ta = alpha * PI / 180.0;
      const double ct = cos(ta);
      const double st = sin(ta);


      // angular resolution

      for (T i = __begin; i != __end; ++i) {
        
        const double dx = i->getX() - pos.getX(); 
        const double dy = i->getY() - pos.getY(); 
        const double dz = i->getZ() - pos.getZ(); 

        const double R  = sqrt(dx*dx + dy*dy);
        
        double x  = ta * getKappaC() * getInverseSpeedOfLight();
        double y  = ta * getKappaC() * getInverseSpeedOfLight();
        double v  = ta * getInverseSpeedOfLight();
        double w  = ta * getInverseSpeedOfLight();

        if (R != 0.0) {
          x *= dx / R;
          y *= dy / R;
        }

        x *=  (dz*ct - dx*st);
        y *=  (dz*ct - dy*st);
        v *= -(dx*ct + dz*st);
        w *= -(dy*ct + dz*st);
        
        buffer[distance(__begin,i)] = variance(x, y, v, w);
      }

      for (int i = 0; i != N; ++i) {

        for (int j = 0; j != i; ++j) {
          (*this)(i, j) = getDot(buffer[i], buffer[j]);
          (*this)(j, i) = (*this)(i, j);
        }

        (*this)(i, i) = getDot(buffer[i], buffer[i]) +  sigma * sigma;
      }
    }

static double JFIT::JMatrixNZ::getDot ( const variance &  first, const variance &  second  )

inlinestaticprivate

Get dot product.

Parameters

first	first variance
second	second variance

Returns

dot product

Definition at line 193 of file JMatrixNZ.hh.

    {
      return (first.x * second.x  +
              first.y * second.y  +
              first.v * second.v  +
              first.w * second.w);
              
    }

void JMATH::JMatrixNS::invert ( )

inlineinherited

Invert matrix.

Definition at line 61 of file JMatrixNS.hh.

    {
      // LDU decomposition

      unsigned int i, j, k;
      double       val, del;
      row_type     row, root;
      col_type     col, cursor;

      JPermutationMatrix P;

      for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {

        row_type pivot;

        val   = (*root)[k];
        pivot = root;

        for (row = root; ++row != this->end(); ) {
          if (fabs((*row)[k]) > fabs(val)) {
            val   = (*row)[k];
            pivot = row;
          }
        }

        if (val == 0) {
          throw JDivisionByZero("LDU decomposition zero pivot");
        }

        if (pivot != root) {

          std::iter_swap(root, pivot);

          P.push_back(std::make_pair(std::distance(this->begin(),root),
                                     std::distance(this->begin(),pivot)));
        }

        // decompose in LDU

        for (row = root; ++row != this->end(); ) {
          
          del = (*row)[k] /= val;
          
          for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); ) {
            (*col) -= del * (*(++cursor));
          }
        }
      }


      // invert D

      for (k = 0; k != this->size(); ++k) {

        double& a = this->at(k,k);
        
        if (a == 0) {
          throw JDivisionByZero("LDU decomposition zero pivot");
        }

        a = 1.0/a;
      }


      // invert U

      for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {

        col_type pivot = root->begin() + k;
        
        for (col = pivot, j = k + 1; ++col != root->end(); ++j) {
          
          (*col) *= -(*pivot);
          
          for (cursor = pivot, row = root; ++cursor != col; ) {
            (*col) -= (*cursor) * (*(++row))[j];
          }
          
          (*col) *= this->at(j,j);
        }
      }


      // invert L

      for (i = size(); i != 0; ) {

        root = this->begin() + --i;

        for (j = i; j != 0; ) {

          double& a = (*root)[--j];

          a = -a;

          for (k = j + 1, cursor = root->begin() + k, row = this->begin() + k; k != i; ++cursor, ++row, ++k) {
            a -= (*cursor) * (*row)[j];
          }
        }
      }


      //  U^-1 x D^-1 x L^-1

      for (root = this->begin(), i = 0; root != this->end(); ++root, ++i) {

        val = (*root)[i];
        
        for (col = root->begin(), j = 0; j != i; ++col, ++j) {
          
          (*col) *= val;
          
          for (row = root, cursor = root->begin() + i; ++cursor != root->end(); ) {
            (*col) += (*cursor) * (*(++row))[j];
          }
        }
        
        for (j = i, col = root->begin() + j; j != this->size() - 1; ++col, ++j) {
          for (k = j + 1, row = this->begin() + k, cursor = root->begin() + k; k != this->size(); ++k, ++cursor, ++row) {
            (*col) += (*(cursor)) * (*(row))[j];
          }
        }
      }


      for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
        for (row = this->begin(); row != this->end(); ++row) {
          std::swap((*row)[p->first], (*row)[p->second]);
        }
      }
    }

void JMATH::JMatrixNS::update ( const unsigned int  k, const double  extra  )

inlineinherited

Update inverted matrix at given diagonal element.

If A is the original matrix and A' is such that for each (i,j) != (k,k):

      A'[i][j] = A[i][j]          
      A'[k][k] = A[k][k] + extra

then JMatrixNS::update(k, extra) will modify the already inverted matrix of A such that it will be the inverse of A'.

Parameters

k	index of diagonal element
extra	value to add

Definition at line 208 of file JMatrixNS.hh.

    {
      if (k >= this->size()) {
        throw JIndexOutOfRange("JMatrixNS::invert(): index out of range.");
      }

      unsigned int i, j;
      row_type     row, root;
      col_type     col;

      const double del  =  this->at(k,k);
      const double din  =  1.0 / (1.0 + extra * del);

      root = this->begin() + k;

      for (i = 0, row = this->begin(); row != this->end(); ++i, ++row) {

        if (i != k) {

          const double val = (*row)[k];

          for (j = 0, col = row->begin(); col != row->end(); ++j, ++col) {

            if (j != k) {
              (*col) -= val * (*root)[j] * extra * din;
            }
          }

          (*row)[k] *= din;
        }
      }

      for (j = 0, col = root->begin(); col != root->end(); ++j, ++col) {

        if (j != k) {
          (*col) *= din;
        }
      }

      (*root)[k] = del * din;
    }

void JMATH::JMatrixND::resize ( const size_t  size )

inlineinherited

Resize matrix.

Note that this method resets the matrix.

Parameters

size

dimension

Definition at line 98 of file JMatrixND.hh.

    {
      matrix_type::resize(size);

      for (row_type row = this->begin(); row != this->end(); ++row) {

        row->resize(size);

        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = 0.0;
        }
      }
    }

double JMATH::JMatrixND::at ( size_t  row, size_t  col  ) const

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns

matrix element at (row,col)

Definition at line 120 of file JMatrixND.hh.

    {
      return (*this)[row][col];
    }

double& JMATH::JMatrixND::at ( size_t  row, size_t  col  )

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns

matrix element at (row,col)

Definition at line 133 of file JMatrixND.hh.

    {
      return (*this)[row][col];
    }

double JMATH::JMatrixND::operator() ( size_t  row, size_t  col  ) const

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns

matrix element at (row,col)

Definition at line 146 of file JMatrixND.hh.

    {
      return at(row, col);
    }

double& JMATH::JMatrixND::operator() ( size_t  row, size_t  col  )

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns

matrix element at (row,col)

Definition at line 159 of file JMatrixND.hh.

    {
      return at(row, col);
    }

JMatrixND& JMATH::JMatrixND::reset ( )

inlineinherited

Set matrix to the null matrix.

Returns

this matrix

Definition at line 170 of file JMatrixND.hh.

    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = 0.0;
        }
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::setIdentity ( )

inlineinherited

Set to identity matrix.

Returns

this matrix

Definition at line 187 of file JMatrixND.hh.

    {
      for (size_t i = 0; i != this->size(); ++i) {

        at(i,i) = 1.0;

        for (size_t j = 0; j != i; ++j) {
          at(i,j) = at(j,i) = 0.0;
        }
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::transpose ( )

inlineinherited

Transpose.

Returns

this matrix

Definition at line 207 of file JMatrixND.hh.

    {
      for (size_t i = 0; i != this->size(); ++i) {
        for (size_t j = 0; j != i; ++j) {
          std::swap(at(i,j), at(j,i));
        }
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::negate ( )

inlineinherited

Negate matrix.

Returns

this matrix

Definition at line 224 of file JMatrixND.hh.

    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = -*col;
        }
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::add ( const JMatrixND &  A )

inlineinherited

Matrix addition.

Parameters

A

matrix

Returns

this matrix

Definition at line 242 of file JMatrixND.hh.

    {
      if (this->size() == A.size()) {

        for (size_t i = 0; i != this->size(); ++i) {
          for (size_t j = 0; j != this->size(); ++j) {
            at(i,j) += A.at(i,j);
          }
        }

      } else {
        THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND &  A )

inlineinherited

Matrix subtraction.

Parameters

A

matrix

Returns

this matrix

Definition at line 266 of file JMatrixND.hh.

    {
      if (this->size() == A.size()) {

        for (size_t i = 0; i != this->size(); ++i) {
          for (size_t j = 0; j != this->size(); ++j) {
            at(i,j) -= A.at(i,j);
          }
        }

      } else {
        THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::mul ( const double  factor )

inlineinherited

Scale matrix.

Parameters

factor

factor

Returns

this matrix

Definition at line 290 of file JMatrixND.hh.

    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col *= factor;
        }
      }

      return *this;
    }

JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND &  A, const JMatrixND &  B  )

inlineinherited

Matrix multiplication.

Parameters

A	matrix
B	matrix

Returns

this matrix

Definition at line 327 of file JMatrixND.hh.

    {
      if (A.size() == B.size()) {
        
        this->resize(A.size());
        
        for (size_t i = 0; i != A.size(); ++i) {
          for (size_t j = 0; j != A.size(); ++j) {
            
            this->at(i,j) = 0.0;
            
            for (size_t k = 0; k != A.size(); ++k) {
              this->at(i,j) += A[i][k] * B[k][j];
            }
          }
        }
        
      } else {

        THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
      }
      
      return *this;
    }

JMatrixND & JMATH::JMath< JMatrixND , JNullType >::mul ( const JNullType &  object )

inlineinherited

Multiply with object.

Parameters

object

object

Returns

result object

Definition at line 273 of file JMath.hh.

    {
      return static_cast<JFirst_t&>(*this) = JCalculator<JFirst_t>::calculator.mul(static_cast<const JFirst_t&>(*this), object);
    }

JMatrixND& JMATH::JMatrixND::div ( const double  factor )

inlineinherited

Scale matrix.

Parameters

factor

factor

Returns

this matrix

Definition at line 308 of file JMatrixND.hh.

    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col /= factor;
        }
      }

      return *this;
    }

bool JMATH::JMatrixND::equals ( const JMatrixND &  A, const double  eps = std::numeric_limits<double>::min()  ) const

inlineinherited

Equality.

Parameters

A	matrix
eps	numerical precision

Returns

true if matrices identical; else false

Definition at line 361 of file JMatrixND.hh.

    {
      if (this->size() == A.size()) {

        for (const_row_type row1 = this->begin(), row2 = A.begin(); row1 != this->end(); ++row1, ++row2) {
          for (const_col_type col1 = row1->begin(), col2 = row2->begin(); col1 != row1->end(); ++col1, ++col2) {
            if (fabs(*col1 - *col2) > eps) {
              return false;
            }
          }
        }

        return true;

      } else {
        THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }
    }

bool JMATH::JMatrixND::isIdentity ( const double  eps = std::numeric_limits<double>::min() ) const

inlineinherited

Test identity.

Parameters

eps	numerical precision

Returns

true if identity matrix; else false

Definition at line 388 of file JMatrixND.hh.

    {
      for (size_t i = 0; i != this->size(); ++i) {

        if (fabs(1.0 - at(i,i)) > eps) {
          return false;
        };
        
        for (size_t j = 0; j != i; ++j) {
          if (fabs(at(i,j)) > eps || fabs(at(j,i)) > eps) {
            return false;
          };
        }
      }

      return true;
    }

double JMATH::JMatrixND::getDeterminant ( ) const

inlineinherited

Get determinant of matrix.

Returns

determinant of matrix

Definition at line 412 of file JMatrixND.hh.

    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);
    
      double det = 1.0;
      
      for (size_t i = 0; i != A.size(); ++i) {
        det *= A[i][i];
      }
      
      return det * A.getSign();
    }

bool JMATH::JMatrixND::isPositiveDefinite ( ) const

inlineinherited

Test positive definiteness.

Returns

true if positive definite; else false

Definition at line 433 of file JMatrixND.hh.

    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);
    
      for (size_t i = 0; i != A.size(); ++i) {
        if (A[i][i] <= 0.0) {
          return false;
        }
      }

      return true;
    }

const bool JMATH::JMatrixND::isPositiveSemiDefinite ( ) const

inlineinherited

Test positive semi-definiteness.

Returns

true if positive semi-definite; else false

Definition at line 454 of file JMatrixND.hh.

    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);    

      for (size_t i = 0; i != A.size(); ++i) {
        if (A[i][i] < 0.0) {
          return false;
        }
      }

      return true;
    }

Member Data Documentation

std::vector<variance> JFIT::JMatrixNZ::buffer

private

Definition at line 183 of file JMatrixNZ.hh.

---

The documentation for this class was generated from the following file:

• JMatrixNZ.hh