9 #include <framework/utilities/ConditionalGaussGenerator.h>
13 using Eigen::MatrixXd;
14 using Eigen::VectorXd;
18 static std::vector<VectorXd> toVectors(
const MatrixXd& mat)
20 std::vector<VectorXd> vList;
21 for (
int j = 0; j < mat.cols(); ++j)
22 vList.push_back(mat.col(j));
28 static MatrixXd toMatrix(
const std::vector<VectorXd>& vList)
30 if (vList.size() == 0)
33 MatrixXd mat(vList[0].rows(), vList.size());
34 for (
unsigned j = 0; j < vList.size(); ++j)
35 mat.col(j) = vList[j];
41 static std::vector<VectorXd> getOrthogonalSpace(VectorXd v0)
44 return toVectors(MatrixXd::Identity(v0.rows(), v0.rows()));
46 VectorXd v0Norm = v0.normalized();
48 std::vector<VectorXd> vList;
49 vList.push_back(v0Norm);
51 while (
int(vList.size()) < v0.rows()) {
53 VectorXd v = VectorXd::Random(v0.rows());
55 for (
const VectorXd& vi : vList)
56 v -= vi.dot(v) * vi / vi.dot(vi);
59 vList.push_back(v.normalized());
64 vList.erase(vList.begin());
74 Eigen::SelfAdjointEigenSolver<MatrixXd> sol(m_covMat);
75 VectorXd vals = sol.eigenvalues();
76 MatrixXd vecs = sol.eigenvectors();
79 std::vector<VectorXd> vBase;
81 for (
int i = 0; i < vals.rows(); ++i)
83 vBase.push_back(sqrt(vals[i]) * vecs.col(i));
87 if (vBase.size() == 0)
91 m_vBaseMat = toMatrix(vBase);
94 VectorXd v0 = m_vBaseMat.row(0);
97 m_cOrt = getOrthogonalSpace(v0);
99 m_v0norm = v0.dot(v0) > 0 ? v0 / v0.dot(v0) : v0;
103 VectorXd ConditionalGaussGenerator::generate(
double x0)
const
105 double dx0 = x0 -
m_mu[0];
113 for (
const VectorXd& co :
m_cOrt) {
114 double rr = gRandom->Gaus();
Eigen::VectorXd m_mu
central values of the distribution
Eigen::MatrixXd m_vBaseMat
transformation matrix between eigen-system of m_covMat and nominal system
std::vector< Eigen::VectorXd > m_cOrt
array of vectors describing free degrees of freedom of random generator
Eigen::VectorXd m_v0norm
vector which scales with dx0
ConditionalGaussGenerator()
dummy constructor without arguments
Abstract base class for different kinds of events.