Belle II Software  light-2212-foldex
ConditionalGaussGenerator.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <framework/utilities/ConditionalGaussGenerator.h>
10 #include <TRandom.h>
11 
12 using namespace Belle2;
13 using Eigen::MatrixXd;
14 using Eigen::VectorXd;
15 
16 
18 static std::vector<VectorXd> toVectors(const MatrixXd& mat)
19 {
20  std::vector<VectorXd> vList;
21  for (int j = 0; j < mat.cols(); ++j)
22  vList.push_back(mat.col(j));
23 
24  return vList;
25 }
26 
28 static MatrixXd toMatrix(const std::vector<VectorXd>& vList)
29 {
30  if (vList.size() == 0)
31  return MatrixXd();
32 
33  MatrixXd mat(vList[0].rows(), vList.size());
34  for (unsigned j = 0; j < vList.size(); ++j)
35  mat.col(j) = vList[j];
36 
37  return mat;
38 }
39 
41 static std::vector<VectorXd> getOrthogonalSpace(VectorXd v0)
42 {
43  if (v0.dot(v0) == 0)
44  return toVectors(MatrixXd::Identity(v0.rows(), v0.rows()));
45 
46  VectorXd v0Norm = v0.normalized();
47 
48  std::vector<VectorXd> vList;
49  vList.push_back(v0Norm);
50 
51  while (int(vList.size()) < v0.rows()) {
52 
53  VectorXd v = VectorXd::Random(v0.rows());
54 
55  for (const VectorXd& vi : vList)
56  v -= vi.dot(v) * vi / vi.dot(vi);
57 
58  if (v.dot(v) > 0) {
59  vList.push_back(v.normalized());
60  }
61  }
62 
63 
64  vList.erase(vList.begin()); // note len(vList) >= 1
65 
66  return vList;
67 
68 }
69 
70 
71 ConditionalGaussGenerator::ConditionalGaussGenerator(const VectorXd& mu, const MatrixXd& covMat) : m_mu(mu), m_covMat(covMat)
72 {
73 
74  Eigen::SelfAdjointEigenSolver<MatrixXd> sol(m_covMat);
75  VectorXd vals = sol.eigenvalues();
76  MatrixXd vecs = sol.eigenvectors();
77 
78 
79  std::vector<VectorXd> vBase;
80  //keep only vectors with positive eigenvalue
81  for (int i = 0; i < vals.rows(); ++i)
82  if (vals[i] > 0) {
83  vBase.push_back(sqrt(vals[i]) * vecs.col(i));
84  }
85 
86 
87  if (vBase.size() == 0) // matrix is zero, no smearing
88  return;
89 
90 
91  m_vBaseMat = toMatrix(vBase);
92 
93  // get the longitudinal space
94  VectorXd v0 = m_vBaseMat.row(0);
95 
96  // get the orthogonal space
97  m_cOrt = getOrthogonalSpace(v0);
98 
99  m_v0norm = v0.dot(v0) > 0 ? v0 / v0.dot(v0) : v0; //normalize, or keep it zero
100 
101 }
102 
103 VectorXd ConditionalGaussGenerator::generate(double x0) const
104 {
105  double dx0 = x0 - m_mu[0];
106 
107  // in case of zero cov matrix
108  if (m_vBaseMat.cols() == 0)
109  return m_mu;
110 
111  VectorXd x = m_v0norm * dx0;
112 
113  for (const VectorXd& co : m_cOrt) {
114  double rr = gRandom->Gaus();
115  x += rr * co;
116  }
117 
118  VectorXd vec = m_vBaseMat * x;
119 
120  return (m_mu + vec);
121 }
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.
Definition: ClusterUtils.h:23