Belle II Software light-2406-ragdoll
test_gaussian_generators.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 <framework/utilities/MultivariateNormalGenerator.h>
11
12#include <TRandom.h>
13#include <gtest/gtest.h>
14#include <cmath>
15#include <iostream>
16
17using namespace Belle2;
18
19namespace {
20
22 class GaussGeneratorsTests : public ::testing::Test {
23 protected:
24
26 virtual void SetUp() { }
27
29 virtual void TearDown() { }
30
31 };
32
34 TEST_F(GaussGeneratorsTests, Generator6DTest)
35 {
36 Eigen::MatrixXd cov(6, 6);
37
38 // generate random matrix
39 for (int i = 0; i < 6; ++i)
40 for (int j = 0; j < 6; ++j)
41 cov(i, j) = gRandom->Poisson(20) - 20;
42
43 // symmetrise it to obtain covariance
44 cov *= cov.transpose();
45
46 // some random central values
47 Eigen::VectorXd mu(6);
48 mu << 3, 0, -2, 7, -5, -1; // cppcheck-suppress constStatement
49
50 // init the generator
51 MultivariateNormalGenerator gen(mu, cov);
52
53 // calculate the mean and cov of the generated events
54 int N = 50000000;
55 Eigen::VectorXd muS = Eigen::VectorXd::Zero(6);
56 Eigen::MatrixXd covS = Eigen::MatrixXd::Zero(6, 6);
57 for (int i = 0; i < N; ++i) {
58 // generate whole random vector
59 Eigen::VectorXd x = gen.generate();
60 muS += x / N;
61 covS += (x - mu) * (x - mu).transpose() / N;
62 }
63
64 // check that mean vector and calculated cov agree with expectation
65 EXPECT_NEAR((muS - mu).squaredNorm(), 0.0, 1e-2);
66 EXPECT_NEAR((covS - cov).squaredNorm(), 0.0, 5e-2);
67
68 }
69
70
72 TEST_F(GaussGeneratorsTests, ConditionalGenerator2DTest)
73 {
74 //some random cov matrix
75 Eigen::MatrixXd cov(2, 2);
76 cov << 9, -4,
77 -4, 4; // cppcheck-suppress constStatement
78
79 // some mean
80 Eigen::VectorXd mu(2);
81 mu << 5., 9.; // cppcheck-suppress constStatement
82
83 // init the generator
84 ConditionalGaussGenerator gen(mu, cov);
85
86 // generate the events using conditional generator
87 int N = 10000000;
88 Eigen::VectorXd muS = Eigen::VectorXd::Zero(2);
89 Eigen::MatrixXd covS = Eigen::MatrixXd::Zero(2, 2);
90 for (int i = 0; i < N; ++i) {
91 // generate x0 using "external" random generator
92 double x0 = gRandom->Gaus(mu[0], sqrt(cov(0, 0)));
93
94 // generate whole random vector based on x0
95 Eigen::VectorXd x = gen.generate(x0);
96 muS += x / N;
97 covS += (x - mu) * (x - mu).transpose() / N;
98 }
99
100 // check that mean vector and calculated cov agree with expectation
101 EXPECT_NEAR((muS - mu).squaredNorm(), 0.0, 1e-3);
102 EXPECT_NEAR((covS - cov).squaredNorm(), 0.0, 1e-3);
103
104 }
105
107 TEST_F(GaussGeneratorsTests, ConditionalGenerator6DTest)
108 {
109 Eigen::MatrixXd cov(6, 6);
110
111 // generate random matrix
112 for (int i = 0; i < 6; ++i)
113 for (int j = 0; j < 6; ++j)
114 cov(i, j) = gRandom->Poisson(20) - 20;
115
116 // symmetrise it to obtain covariance
117 cov *= cov.transpose();
118
119 // some random central values
120 Eigen::VectorXd mu(6);
121 mu << 3, 0, -2, 7, -5, -1; // cppcheck-suppress constStatement
122
123 // init the generator
124 ConditionalGaussGenerator gen(mu, cov);
125
126 // calculate the mean and cov of the generated events
127 int N = 50000000;
128 Eigen::VectorXd muS = Eigen::VectorXd::Zero(6);
129 Eigen::MatrixXd covS = Eigen::MatrixXd::Zero(6, 6);
130 for (int i = 0; i < N; ++i) {
131 // generate x0 using "external" random generator
132 double x0 = gRandom->Gaus(mu[0], sqrt(cov(0, 0)));
133
134 // generate whole random vector based on x0
135 Eigen::VectorXd x = gen.generate(x0);
136 muS += x / N;
137 covS += (x - mu) * (x - mu).transpose() / N;
138 }
139
140 // check that mean vector and calculated cov agree with expectation
141 EXPECT_NEAR((muS - mu).squaredNorm(), 0.0, 1e-2);
142 EXPECT_NEAR((covS - cov).squaredNorm(), 0.0, 5e-2);
143
144 }
145
146}
Class implementing n-dimensional random number generator from Gaussian distribution where the first c...
Class to generate normal distributed, correlated random numbers given the mean values and the covaria...
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24