Belle II Software  release-05-01-25
numerics.test.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2014 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost <oliver.frost@desy.de> *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/trackFindingCDC/numerics/SinEqLine.h>
11 
12 #include <tracking/trackFindingCDC/geometry/Line2D.h>
13 
14 #include <tracking/trackFindingCDC/numerics/CovarianceMatrixUtil.h>
15 
16 #include <tracking/trackFindingCDC/numerics/Median.h>
17 #include <tracking/trackFindingCDC/numerics/WithWeight.h>
18 #include <tracking/trackFindingCDC/numerics/ESign.h>
19 
20 #include <framework/gearbox/Unit.h>
21 
22 #include <gtest/gtest.h>
23 
24 
25 using namespace Belle2;
26 using namespace TrackFindingCDC;
27 
28 
29 TEST(TrackFindingCDCTest, numerics_sign)
30 {
31  EXPECT_EQ(ESign::c_Plus, sign(INFINITY));
32  EXPECT_EQ(ESign::c_Plus, sign(2));
33  EXPECT_EQ(ESign::c_Plus, sign(0.0));
34  EXPECT_EQ(ESign::c_Minus, sign(-0.0));
35  EXPECT_EQ(ESign::c_Minus, sign(-2));
36  EXPECT_EQ(ESign::c_Minus, sign(-INFINITY));
37  EXPECT_EQ(ESign::c_Invalid, sign(NAN));
38 
39  EXPECT_TRUE(ESignUtil::isValid(ESign::c_Plus));
40  EXPECT_TRUE(ESignUtil::isValid(ESign::c_Minus));
41  EXPECT_TRUE(ESignUtil::isValid(ESign::c_Zero));
42 
43  EXPECT_FALSE(ESignUtil::isValid(ESign::c_Invalid));
44  EXPECT_FALSE(ESignUtil::isValid(static_cast<ESign>(7)));
45 
46  EXPECT_EQ(ESign::c_Minus, ESignUtil::opposite(ESign::c_Plus));
47  EXPECT_EQ(ESign::c_Plus, ESignUtil::opposite(ESign::c_Minus));
48  EXPECT_EQ(ESign::c_Zero, ESignUtil::opposite(ESign::c_Zero));
49  EXPECT_EQ(ESign::c_Invalid, ESignUtil::opposite(ESign::c_Invalid));
50 }
51 
52 
53 TEST(TrackFindingCDCTest, numerics_SinEqLine_isIncreasing)
54 {
55 
56  Vector2D lower(0.0, 1.0);
57  Vector2D upper(1.0, 2.0);
58 
59  EXPECT_EQ(EIncDec::c_Increasing, SinEqLine::getEIncDec(lower, upper));
60 
61 }
62 
63 
64 TEST(TrackFindingCDCTest, numerics_SinEqLine_getIPeriodFromIHalfPeriod)
65 {
66  int iHalfPeriod = -1;
67  int iPeriod = SinEqLine::getIPeriodFromIHalfPeriod(iHalfPeriod);
68  EXPECT_EQ(-1, iPeriod);
69 }
70 
71 
72 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootInInterval_simple)
73 {
74 
75  // Simple sin.
76  SinEqLine sinEqLine(0.0, 0.0);
77 
78  double rootX = sinEqLine.computeRootInInterval(M_PI / 2, 3 * M_PI / 2);
79  double rootY = sinEqLine.map(rootX);
80 
81  EXPECT_NEAR(M_PI, rootX, 10e-7);
82  EXPECT_NEAR(0.0, rootY, 10e-7);
83 
84 }
85 
86 
87 
88 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootInInterval_const)
89 {
90 
91  // Constant sin.
92  SinEqLine sinEqLine(0.0, 1.0 / 2.0);
93 
94  double rootX = sinEqLine.computeRootInInterval(M_PI / 2, 3 * M_PI / 2);
95  double rootY = sinEqLine.map(rootX);
96 
97  EXPECT_NEAR(150.0 * Unit::deg, rootX, 10e-7);
98  EXPECT_NEAR(0.0, rootY, 10e-7);
99 
100 }
101 
102 
103 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootInInterval_complex)
104 {
105 
106  // Setup a line that is a
107  double rootX = 150.0 * Unit::deg;
108  Line2D line = Line2D::fromSlopeIntercept(1.0 / 2.0, 1.0 / 2.0);
109  line.moveAlongFirst(rootX);
110 
111  SinEqLine sinEqLine(line);
112 
113  double solvedRootX = sinEqLine.computeRootInInterval(M_PI / 2, 3 * M_PI / 2);
114  double solvedRootY = sinEqLine.map(rootX);
115 
116  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
117  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
118 
119 }
120 
121 
122 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootLargerThanExtemumInHalfPeriod_simple)
123 {
124 
125  // Setup a line that is a
126  double rootX = 0.0 * Unit::deg;
127  Line2D line = Line2D::fromSlopeIntercept(-1.0 / 2.0, 0);
128 
129  SinEqLine sinEqLine(line);
130 
131  double solvedRootX = sinEqLine.computeRootLargerThanExtemumInHalfPeriod(-1);
132  double solvedRootY = sinEqLine.map(rootX);
133 
134  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
135  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
136 
137 }
138 
139 
140 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootLargerThanExtemumInHalfPeriod)
141 {
142 
143  // Setup a line that is a
144  double rootX = 150.0 * Unit::deg;
145  Line2D line = Line2D::fromSlopeIntercept(1.0 / 2.0, 1.0 / 2.0);
146  line.moveAlongFirst(rootX);
147 
148  SinEqLine sinEqLine(line);
149 
150  double solvedRootX = sinEqLine.computeRootLargerThanExtemumInHalfPeriod(0);
151  double solvedRootY = sinEqLine.map(rootX);
152 
153  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
154  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
155 
156 }
157 
158 
159 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeSmallestPositiveRoot)
160 {
161 
162  // Setup a line that is a
163  double rootX = 150.0 * Unit::deg;
164  Line2D line = Line2D::fromSlopeIntercept(1.0 / 2.0, 1.0 / 2.0);
165  line.moveAlongFirst(rootX);
166 
167  SinEqLine sinEqLine(line);
168 
169  double solvedRootX = sinEqLine.computeSmallestPositiveRoot();
170  double solvedRootY = sinEqLine.map(rootX);
171 
172  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
173  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
174 
175 }
176 
177 
178 
179 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeSmallestPositiveRoot_largeSlope)
180 {
181 
182  // Setup a line that is a
183  double rootX = 150.0 * Unit::deg;
184  Line2D line = Line2D::fromSlopeIntercept(2.0, 1.0 / 2.0);
185  line.moveAlongFirst(rootX);
186 
187  SinEqLine sinEqLine(line);
188 
189  double solvedRootX = sinEqLine.computeSmallestPositiveRoot();
190  double solvedRootY = sinEqLine.map(rootX);
191 
192  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
193  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
194 
195 }
196 
197 
198 
199 TEST(TrackFindingCDCTest, numerics_SinEqLine_computeRootForLargeSlope)
200 {
201  // Setup a line that is a
202  double rootX = 150.0 * Unit::deg;
203  Line2D line = Line2D::fromSlopeIntercept(2.0, 1.0 / 2.0);
204  line.moveAlongFirst(rootX);
205 
206  SinEqLine sinEqLine(line);
207 
208  double solvedRootX = sinEqLine.computeRootForLargeSlope();
209  double solvedRootY = sinEqLine.map(rootX);
210 
211  EXPECT_NEAR(rootX, solvedRootX, 10e-7);
212  EXPECT_NEAR(0.0, solvedRootY, 10e-7);
213 }
214 
215 TEST(TrackFindingCDCTest, numerics_median)
216 {
217  std::vector<double> fourValues{0.0, 3.0, 5.0, 100};
218  EXPECT_EQ(4.0, median(std::move(fourValues)));
219 
220  std::vector<double> threeValues{0.0, 3.0, 5.0, 100.0, 1000.0};
221  EXPECT_EQ(5.0, median(std::move(threeValues)));
222 
223  std::vector<WithWeight<double>> weightedValues{{0.0, 0.2},
224  {3.0, 0.1},
225  {5.0, 0.1},
226  {100.0, 0.3},
227  {1000.0, 0.3}};
228  EXPECT_EQ(100.0, weightedMedian(std::move(weightedValues)));
229 }
230 
231 TEST(TrackFindingCDCTest, covariance_simple_average)
232 {
233  // Simple average of two scalar values
234  const ParameterVector<1> par1{ -1};
235  const CovarianceMatrix<1> cov1{2};
236 
237  const ParameterVector<1> par2{1};
238  const CovarianceMatrix<1> cov2{2};
239 
240  ParameterVector<1> par{};
241  CovarianceMatrix<1> cov{};
242 
243  double chi2 = CovarianceMatrixUtil::average(par1, cov1, par2, cov2, par, cov);
244  EXPECT_NEAR(0, par[0], 1E-8);
245  EXPECT_NEAR(1, cov[0], 1E-8);
246  EXPECT_NEAR(1, chi2, 1E-8);
247 }
248 
249 TEST(TrackFindingCDCTest, covariance_half_projection_average)
250 {
251  // Same average of two scalar values but one is 'projected' stretched by a factor of 2
252  const JacobianMatrix<1, 1> ambi1{2};// Ambiguity is 2.
253  const ParameterVector<1> par1{ -2}; // Parameter scales with 2
254  const CovarianceMatrix<1> cov1{8}; // Covariance scales with 2 * 2
255 
256  const ParameterVector<1> par2{1};
257  const CovarianceMatrix<1> cov2{2};
258 
259  ParameterVector<1> par{};
260  CovarianceMatrix<1> cov{};
261 
262  double chi2 = CovarianceMatrixUtil::average(par1, cov1, ambi1, par2, cov2, par, cov);
263  EXPECT_NEAR(0, par[0], 1E-8);
264  EXPECT_NEAR(1, cov[0], 1E-8);
265  EXPECT_NEAR(1, chi2, 1E-8);
266 }
267 
268 TEST(TrackFindingCDCTest, covariance_kalman_update_average)
269 {
270  // Same average of two scalar values but one is 'projected' / stretched by a factor of 2
271  // Now as an inplace update using the Kalman formula.
272  const JacobianMatrix<1, 1> ambi1{2}; // Ambiguity is 2.
273  const ParameterVector<1> par1{ -2}; // Parameter scales with 2
274  const CovarianceMatrix<1> cov1{8}; // Covariance scales with 2 * 2
275 
276  ParameterVector<1> par{1};
277  CovarianceMatrix<1> cov{2};
278 
279  double chi2 = CovarianceMatrixUtil::kalmanUpdate(par1, cov1, ambi1, par, cov);
280  EXPECT_NEAR(0, par[0], 1E-8);
281  EXPECT_NEAR(1, cov[0], 1E-8);
282  EXPECT_NEAR(1, chi2, 1E-8);
283 }
Belle2::TrackFindingCDC::SinEqLine::getIPeriodFromIHalfPeriod
static int getIPeriodFromIHalfPeriod(int iHalfPeriod)
Helper function to translate the index of the half period to index of the containing period.
Definition: SinEqLine.h:161
Belle2::TrackFindingCDC::CovarianceMatrixUtil::average
static double average(const ParameterVector< N > &par1, const CovarianceMatrix< N > &cov1, const ParameterVector< N > &par2, const CovarianceMatrix< N > &cov2, ParameterVector< N > &par, CovarianceMatrix< N > &cov)
Averages two parameter vectors taking into account their respective covariances.
Definition: CovarianceMatrixUtil.h:156
Belle2::TrackFindingCDC::ESignUtil::ESign
ESign
Enumeration for the distinct sign values of floating point variables.
Definition: ESign.h:37
Belle2::TrackFindingCDC::ESignUtil::isValid
bool isValid(ESign s)
Returns true if sign is ESign::c_Plus, ESign::c_Minus or ESign::c_Zero.
Definition: ESign.h:56
Belle2::TrackFindingCDC::Line2D::fromSlopeIntercept
static Line2D fromSlopeIntercept(const double slope, const double intercept)
Constructs a line from its slope and intercept over the first coordinate (default forward orientation...
Definition: Line2D.h:80
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::ESignUtil::opposite
ESign opposite(ESign s)
Return the opposite sign. Leaves ESign::c_Invalid the same.
Definition: ESign.h:52
Belle2::Unit::deg
static const double deg
degree to radians
Definition: Unit.h:119
Belle2::TEST
TEST(TestgetDetectorRegion, TestgetDetectorRegion)
Test Constructors.
Definition: utilityFunctions.cc:18
Belle2::TrackFindingCDC::SinEqLine::getEIncDec
static EIncDec getEIncDec(const Vector2D &lower, const Vector2D &upper)
Determines if the function is increasing or decreasing in the intervall.
Definition: SinEqLine.h:143
Belle2::TrackFindingCDC::CovarianceMatrixUtil::kalmanUpdate
static double kalmanUpdate(const ParameterVector< N1 > &par1, const CovarianceMatrix< N1 > &cov1, const JacobianMatrix< N1, M > &ambiguity1, ParameterVector< M > &par2, CovarianceMatrix< M > &cov2)
Updates a parameter and its covariance by integrating information from parameter vector in a projecte...
Definition: CovarianceMatrixUtil.h:259