CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

RandomMultiGauss Class Reference

#include <RandomMultiGauss.h>

List of all members.

Public Member Functions

AlgebraicVector fire ()
 RandomMultiGauss (const AlgebraicSymMatrix &aMatrix)
 RandomMultiGauss (const AlgebraicVector &aVector, const AlgebraicSymMatrix &aMatrix)
 ~RandomMultiGauss ()

Private Member Functions

void initialise (const AlgebraicSymMatrix &)

Private Attributes

AlgebraicVector theMeans
int theSize
AlgebraicMatrix theTriangle

Detailed Description

Vector of random numbers according to covariance matrix. Generates vectors of random numbers given a vector of mean values (optional) and a covariance matrix. Will accept empty rows/columns in the input matrix. Uses CLHEP::RandGauss with default engine for generation.

Definition at line 14 of file RandomMultiGauss.h.


Constructor & Destructor Documentation

RandomMultiGauss::RandomMultiGauss ( const AlgebraicVector aVector,
const AlgebraicSymMatrix aMatrix 
)

constructor with explicit vector of mean values

Definition at line 9 of file RandomMultiGauss.cc.

References initialise(), theMeans, and theSize.

                                                                                                     :
  theSize(aMatrix.num_row()),
  theMeans(aVector),
  theTriangle(theSize,theSize,0) {
  //
  // Check consistency
  //
  if ( theMeans.num_row() == theSize ) {
    initialise(aMatrix);
  }
  else {
//    throw DetLogicError("RandomMultiGauss: size of vector and matrix do not match");
    theMeans = AlgebraicVector(theSize,0);
  }
}
RandomMultiGauss::RandomMultiGauss ( const AlgebraicSymMatrix aMatrix)

constructor with covariance matrix only (all means = 0)

Definition at line 27 of file RandomMultiGauss.cc.

References initialise().

                                                                     :
  theSize(aMatrix.num_row()),
  theMeans(theSize,0),
  theTriangle(theSize,theSize,0) {
  //
  initialise(aMatrix);
}
RandomMultiGauss::~RandomMultiGauss ( ) [inline]

Definition at line 23 of file RandomMultiGauss.h.

{}

Member Function Documentation

AlgebraicVector RandomMultiGauss::fire ( )

Generation of a vector of random numbers.

Definition at line 74 of file RandomMultiGauss.cc.

References i, theMeans, theSize, and theTriangle.

                                       {
  AlgebraicVector vRandom(theSize,0);
  for ( int i=0; i<theSize; i++ ) {
    if ( theTriangle[i][i]!=0 ) 
      vRandom[i] = CLHEP::RandGauss::shoot();
  }
  return theTriangle*vRandom+theMeans;
}
void RandomMultiGauss::initialise ( const AlgebraicSymMatrix aMatrix) [private]

generation of the cholesky decomposition

Definition at line 37 of file RandomMultiGauss.cc.

References mathSSE::sqrt(), theSize, and theTriangle.

Referenced by RandomMultiGauss().

                                                                    {
  //
  // Cholesky decomposition with protection against empty rows/columns
  //
  for ( int i1=0; i1<theSize; i1++ ) {
    if ( fabs(aMatrix[i1][i1])<FLT_MIN )  continue;

    for ( int i2=i1; i2<theSize; i2++ ) {
      if ( fabs(aMatrix[i2][i2])<FLT_MIN )  continue;

      double sum = aMatrix[i2][i1];
      for ( int i3=i1-1; i3>=0; i3-- ) {
        if ( fabs(aMatrix[i3][i3])<FLT_MIN )  continue;
        sum -= theTriangle[i1][i3]*theTriangle[i2][i3];
      }
      
      if ( i1==i2 ) {
        //
        // check for positive definite input matrix, but allow for effects
        // due to finite precision
        //
        if ( sum<=0 ) {
//        if ( sum<-FLT_MIN )  throw DetLogicError("RandomMultiGauss: input matrix is not positive definite");
          sum = FLT_MIN;
        }
        theTriangle[i1][i1] = sqrt(sum);
      }
      else {
        theTriangle[i1][i2] = 0.;
        theTriangle[i2][i1] = sum / theTriangle[i1][i1];
      }
    }
  }
}

Member Data Documentation

Definition at line 35 of file RandomMultiGauss.h.

Referenced by fire(), and RandomMultiGauss().

Definition at line 34 of file RandomMultiGauss.h.

Referenced by fire(), initialise(), and RandomMultiGauss().

Definition at line 36 of file RandomMultiGauss.h.

Referenced by fire(), and initialise().