File indexing completed on 2023-10-25 09:36:10
0001
0002 #include "CommonTools/Statistics/interface/RandomMultiGauss.h"
0003 #include "CLHEP/Random/RandGauss.h"
0004
0005 #include <cfloat>
0006
0007
0008
0009 RandomMultiGauss::RandomMultiGauss(const AlgebraicVector& aVector, const AlgebraicSymMatrix& aMatrix)
0010 : theSize(aMatrix.num_row()), theMeans(aVector), theTriangle(theSize, theSize, 0) {
0011
0012
0013
0014 if (theMeans.num_row() == theSize) {
0015 initialise(aMatrix);
0016 } else {
0017
0018 theMeans = AlgebraicVector(theSize, 0);
0019 }
0020 }
0021
0022
0023
0024 RandomMultiGauss::RandomMultiGauss(const AlgebraicSymMatrix& aMatrix)
0025 : theSize(aMatrix.num_row()), theMeans(theSize, 0), theTriangle(theSize, theSize, 0) {
0026
0027 initialise(aMatrix);
0028 }
0029
0030
0031
0032 void RandomMultiGauss::initialise(const AlgebraicSymMatrix& aMatrix) {
0033
0034
0035
0036 for (int i1 = 0; i1 < theSize; i1++) {
0037 if (fabs(aMatrix[i1][i1]) < FLT_MIN)
0038 continue;
0039
0040 for (int i2 = i1; i2 < theSize; i2++) {
0041 if (fabs(aMatrix[i2][i2]) < FLT_MIN)
0042 continue;
0043
0044 double sum = aMatrix[i2][i1];
0045 for (int i3 = i1 - 1; i3 >= 0; i3--) {
0046 if (fabs(aMatrix[i3][i3]) < FLT_MIN)
0047 continue;
0048 sum -= theTriangle[i1][i3] * theTriangle[i2][i3];
0049 }
0050
0051 if (i1 == i2) {
0052
0053
0054
0055
0056 if (sum <= 0) {
0057
0058 sum = FLT_MIN;
0059 }
0060 theTriangle[i1][i1] = sqrt(sum);
0061 } else {
0062 theTriangle[i1][i2] = 0.;
0063 theTriangle[i2][i1] = sum / theTriangle[i1][i1];
0064 }
0065 }
0066 }
0067 }
0068
0069
0070
0071 AlgebraicVector RandomMultiGauss::fire() {
0072 AlgebraicVector vRandom(theSize, 0);
0073 for (int i = 0; i < theSize; i++) {
0074 if (theTriangle[i][i] != 0)
0075 vRandom[i] = CLHEP::RandGauss::shoot();
0076 }
0077 return theTriangle * vRandom + theMeans;
0078 }