TEvtGen/EvtGenBase/EvtMatrix.hh

   1 /*****************************************************************************
   2  * Project: BaBar detector at the SLAC PEP-II B-factory
   3  * Package: EvtGenBase
   4  *    File: $Id: EvtMatrix.hh,v 1.5 2008/11/27 16:27:24 jordix Exp $
   5  *
   6  * Description:
   7  *   Class to make simple computations with matrices: assignment, product,
   8  *      determinant, inverse... Still many functions could be implemented.
   9  *
  10  * Modification history:
  11  *   Jordi Garra Ticó     2008/07/03         File created
  12  *****************************************************************************/
  13
  14
  15 #ifndef __EVT_MATRIX_HH__
  16 #define __EVT_MATRIX_HH__
  17
  18 #include <vector>
  19 #include <sstream>
  20
  21
  22
  23 template <class T> class EvtMatrix
  24 {
  25 private:
  26   T** _mat;
  27   int _range;
  28 public:
  29   EvtMatrix() : _range( 0 ) {};
  30   ~EvtMatrix();
  31   inline void setRange( int range );
  32
  33   T& operator()( int row, int col ) { return _mat[ row ][ col ]; }
  34   T* operator[]( int row )          { return _mat[ row ];        }
  35   T det();
  36   EvtMatrix* min( int row, int col );
  37   EvtMatrix* inverse();
  38   std::string dump();
  39
  40   template <class M> friend EvtMatrix< M >* operator*( const EvtMatrix< M >& left, const EvtMatrix< M >& right );
  41 };
  42
  43
  44
  45
  46 template <class T> inline void EvtMatrix< T >::setRange( int range )
  47 {
  48   // If the range is changed, delete any previous matrix stored
  49   //    and allocate elements with the newly specified range.
  50   if ( _range != range )
  51     {
  52       if ( _range )
  53         {
  54           for ( int row = 0; row < _range; row++ )
  55             delete[] _mat[ row ];
  56           delete[] _mat;
  57         }
  58
  59       _mat = new T*[ range ];
  60       for ( int row = 0; row < range; row++ )
  61         _mat[ row ] = new T[ range ];
  62
  63       // Set the new range.
  64       _range = range;
  65     }
  66
  67   // Since user is willing to change the range, reset the matrix elements.
  68   for ( int row = 0; row < _range; row++ )
  69     for ( int col = 0; col < _range; col++ )
  70       _mat[ row ][ col ] = 0.;
  71 }
  72
  73
  74 template <class T> EvtMatrix< T >::~EvtMatrix()
  75 {
  76   for( int row = 0; row < _range; row++ )
  77     delete[] _mat[ row ];
  78   delete[] _mat;
  79 }
  80
  81
  82 template <class T> std::string EvtMatrix< T >::dump()
  83 {
  84   std::ostringstream str;
  85
  86   for ( int row = 0; row < _range; row++ )
  87     {
  88       str << "|";
  89       for ( int col = 0; col < _range; col++ )
  90         str << "\t" << _mat[ row ][ col ];
  91       str << "\t|" << std::endl;
  92     }
  93
  94   return str.str();
  95 }
  96
  97
  98 template <class T> T EvtMatrix< T >::det()
  99 {
 100   if ( _range == 1 )
 101     return _mat[ 0 ][ 0 ];
 102
 103   // There's no need to define the range 2 determinant manually, but it may
 104   //    speed up the calculation.
 105   if ( _range == 2 )
 106     return _mat[ 0 ][ 0 ] * _mat[ 1 ][ 1 ] - _mat[ 0 ][ 1 ] * _mat[ 1 ][ 0 ];
 107
 108   T sum = 0.;
 109
 110   for ( int col = 0; col < _range; col++ )
 111     {
 112       EvtMatrix< T >* minor = min( 0, col );
 113       sum += pow( -1., col ) * _mat[ 0 ][ col ] * minor->det();
 114       delete minor;
 115     }
 116
 117   return sum;
 118 }
 119
 120
 121 // Returns the minor at (i, j).
 122 template <class T> EvtMatrix< T >* EvtMatrix< T >::min( int row, int col )
 123 {
 124   EvtMatrix< T >* minor = new EvtMatrix< T >();
 125   minor->setRange( _range - 1 );
 126
 127   int minIndex = 0;
 128
 129   for ( int r = 0; r < _range; r++ )
 130     for ( int c = 0; c < _range; c++ )
 131       if ( ( r != row ) && ( c != col ) )
 132         {
 133           (*minor)( minIndex / ( _range - 1 ), minIndex % ( _range - 1 ) ) = _mat[ r ][ c ];
 134           minIndex++;
 135         }
 136
 137   return minor;
 138 }
 139
 140
 141 template <class T> EvtMatrix< T >* EvtMatrix< T >::inverse()
 142 {
 143   EvtMatrix< T >* inv = new EvtMatrix< T >();
 144   inv->setRange( _range );
 145
 146   if ( det() == 0 )
 147     {
 148       std::cerr << "This matrix has a null determinant and cannot be inverted. Returning zero matrix." << std::endl;
 149       for ( int row = 0; row < _range; row++ )
 150         for ( int col = 0; col < _range; col++ )
 151           (*inv)( row, col ) = 0.;
 152       return inv;
 153     }
 154
 155   T determinant = det();
 156
 157   for ( int row = 0; row < _range; row++ )
 158     for ( int col = 0; col < _range; col++ )
 159       {
 160         EvtMatrix< T >* minor = min( row, col );
 161         inv->_mat[col][row] = pow( -1., row + col ) * minor->det() / determinant;
 162         delete minor;
 163       }
 164
 165   return inv;
 166 }
 167
 168
 169 template <class T>
 170 EvtMatrix< T >* operator*( const EvtMatrix< T >& left, const EvtMatrix< T >& right )
 171 {
 172   // Chech that the matrices have the correct range.
 173   if ( left._range != right._range )
 174     {
 175       std::cerr << "These matrices cannot be multiplied." << std::endl;
 176       return new EvtMatrix< T >();
 177     }
 178
 179   EvtMatrix< T >* mat = new EvtMatrix< T >();
 180   mat->setRange( left._range );
 181
 182   // Initialize the elements of the matrix.
 183   for ( int row = 0; row < left._range; row++ )
 184     for ( int col = 0; col < right._range; col++ )
 185       (*mat)[ row ][ col ] = 0;
 186
 187   for ( int row = 0; row < left._range; row++ )
 188     for ( int col = 0; col < right._range; col++ )
 189       for ( int line = 0; line < right._range; line++ )
 190         (*mat)[ row ][ col ] += left._mat[ row ][ line ] * right._mat[ line ][ col ];
 191
 192   return mat;
 193 }
 194
 195
 196 #endif