BlockSolver Class Reference

solves at best the matrix invertion for calibration More...

#include <BlockSolver.h>

Public Member Functions
int	operator() (const CLHEP::HepMatrix &matrix, const CLHEP::HepVector &vector, CLHEP::HepVector &result)

Private Member Functions
void	pour (CLHEP::HepVector &result, const CLHEP::HepVector &output, const std::vector< int > &where)
	pour results in bigger vector More...
void	shrink (const CLHEP::HepMatrix &matrix, CLHEP::HepMatrix &solution, const CLHEP::HepVector &result, CLHEP::HepVector &input, const std::vector< int > &where)
	eliminate empty columns and rows More...

Detailed Description

solves at best the matrix invertion for calibration

Definition at line 21 of file BlockSolver.h.

Member Function Documentation

int BlockSolver::operator()	(	const CLHEP::HepMatrix &	matrix,
		const CLHEP::HepVector &	vector,
		CLHEP::HepVector &	result
	)

Definition at line 7 of file BlockSolver.cc.

References cuy::col, mps_fire::i, input, convertSQLitetoXML_cfg::output, pour(), shrink(), and electronIdCutBased_cfi::threshold.

{
  double threshold = matrix.trace () / 
                     static_cast<double> (matrix.num_row ()) ;
  std::vector<int> holes ;
  // loop over the matrix rows
  for (int row = 0 ; row < matrix.num_row () ; ++row)
    {
      double sumAbs = 0. ;
      for (int col = 0 ; col < matrix.num_col () ; ++col)
        sumAbs += fabs (matrix[row][col]) ;
      if (sumAbs < threshold) holes.push_back (row) ;  
    } // loop over the matrix rows

  int dim = matrix.num_col () - holes.size () ;

  if (holes.size () == 0) //PG exceptional case!
    {
      for (int i = 0 ; i < result.num_row () ; ++i)
        result[i] = 1. ;
    }    
  else if (dim > 0) 
    {
      CLHEP::HepMatrix solution (dim, dim, 0) ;
      CLHEP::HepVector input (dim, 0) ;
      shrink (matrix, solution, vector, input, holes) ;
      CLHEP::HepVector output = solve (solution,input) ;
      pour (result, output, holes) ;
    } 
  return holes.size () ;
}

void BlockSolver::pour ( CLHEP::HepVector &  result, const CLHEP::HepVector &  output, const std::vector< int > &  where  )

private

pour results in bigger vector

Definition at line 88 of file BlockSolver.cc.

References reco::return().

Referenced by operator()().

{
  std::vector<int>::const_iterator whereCols = where.begin () ;
  int readingIndex = 0 ;
  //PG loop over the output crystals
  for (int xtal = 0 ; xtal < result.num_row () ; ++xtal)
    {
      if (xtal == *whereCols)
        {
          result[xtal] = 1. ;
          ++whereCols ;        
        }
      else
        {
          result[xtal] = output[readingIndex] ;
          ++readingIndex ;
        }    
    } //PG loop over the output crystals

  return ;  
}                

void BlockSolver::shrink ( const CLHEP::HepMatrix &  matrix, CLHEP::HepMatrix &  solution, const CLHEP::HepVector &  result, CLHEP::HepVector &  input, const std::vector< int > &  where  )

private

eliminate empty columns and rows

Definition at line 46 of file BlockSolver.cc.

References cuy::col, and reco::return().

Referenced by operator()().

{
  
  int offsetRow = 0 ;
  std::vector<int>::const_iterator whereRows = where.begin () ;
  // loop over rows
  for (int row = 0 ; row < matrix.num_row () ; ++row)
    {
      if (row == *whereRows) 
        {
//          std::cerr << "        DEBUG shr hole found " << std::endl ;
          ++offsetRow ;
          ++whereRows ;
          continue ;
        } 
      input[row-offsetRow] = result[row] ;        
      int offsetCol = 0 ;
      std::vector<int>::const_iterator whereCols = where.begin () ;
      // loop over columns
      for (int col = 0 ; col < matrix.num_col () ; ++col)
        {
          if (col == *whereCols) 
            {
              ++offsetCol ;
              ++whereCols ;
              continue ;
            }
          solution[row-offsetRow][col-offsetCol] = matrix[row][col] ;    
        }
    } // loop over rows
  return ;  
}