Belle II Software light-2406-ragdoll
B2Vector3.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#include <framework/utilities/TestHelpers.h>
9#include <framework/geometry/B2Vector3.h>
10#include <TVector3.h>
11#include <TError.h> // to shut up stupid messages
12
13#include <gtest/gtest.h>
14
15using namespace std;
16using namespace Belle2;
17
18
19namespace {
20 struct PointRThetaPhi { double r; double theta; double phi; };
21 struct PointXYZ { double X; double Y; double Z; };
22 std::vector<PointRThetaPhi> test_vectors(const std::vector<double>& R = {1}, int nTheta = 3, int nPhi = 4, bool prune = true)
23 {
24 std::vector<PointRThetaPhi> result;
25 auto get_theta = [nTheta](int i) -> double {
26 if (nTheta < 2) return 0.;
27 return M_PI / (nTheta - 1) * i;
28 };
29 auto get_phi = [nPhi](int i) -> double {
30 if (nPhi < 2) return 0.;
31 return 2 * M_PI / nPhi * i;
32 };
33
34 for (double r : R) {
35 for (int i = 0; i < nTheta || (nTheta < 2 && i == 0); ++i) {
36 double theta = get_theta(i);
37 for (int j = 0; j < nPhi || (nPhi < 2 && j == 0); ++j) {
38 double phi = get_phi(j);
39 result.push_back(PointRThetaPhi{r, theta, phi});
40 // if theta is 0 or 180 degree phi does not have any effect except for zero signs, so disable those for now
41 if (prune && (theta == 0 || theta == M_PI)) break;
42 }
43 // if r is 0, theta and phi don't have any effect except for zero signs, so disable those for now
44 if (prune && r == 0) break;
45 }
46 }
47 return result;
48 }
49
50 TEST(B2Vector3, DoubleGetters)
51 {
52 gErrorIgnoreLevel = kError;
53 TVector3 tvec;
54 B2Vector3D bvec;
55 for (auto& rtp : test_vectors({0, 1, 1e20}, 64, 64)) {
56 tvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
57 bvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
58 EXPECT_EQ(bvec, tvec);
59
60 B2Vector3D bvec2(tvec);
61 EXPECT_EQ(bvec2, bvec);
62 TVector3 tvec2 = bvec;
63 EXPECT_EQ(tvec2, tvec);
64
65 EXPECT_DOUBLE_EQ(bvec.CosTheta(), tvec.CosTheta()) << bvec.PrintString();
66 EXPECT_DOUBLE_EQ(bvec.Mag(), tvec.Mag()) << bvec.PrintString();
67 EXPECT_DOUBLE_EQ(bvec.Mag2(), tvec.Mag2()) << bvec.PrintString();
68 EXPECT_DOUBLE_EQ(bvec.Perp(), tvec.Perp()) << bvec.PrintString();
69 EXPECT_DOUBLE_EQ(bvec.Perp2(), tvec.Perp2()) << bvec.PrintString();
70 EXPECT_DOUBLE_EQ(bvec.Phi(), tvec.Phi()) << bvec.PrintString();
71 EXPECT_DOUBLE_EQ(bvec.Pt(), tvec.Pt()) << bvec.PrintString();
72 EXPECT_DOUBLE_EQ(bvec.Px(), tvec.Px()) << bvec.PrintString();
73 EXPECT_DOUBLE_EQ(bvec.Py(), tvec.Py()) << bvec.PrintString();
74 EXPECT_DOUBLE_EQ(bvec.Pz(), tvec.Pz()) << bvec.PrintString();
75 EXPECT_DOUBLE_EQ(bvec.Theta(), tvec.Theta()) << bvec.PrintString();
76 EXPECT_DOUBLE_EQ(bvec.X(), tvec.X()) << bvec.PrintString();
77 EXPECT_DOUBLE_EQ(bvec.Y(), tvec.Y()) << bvec.PrintString();
78 EXPECT_DOUBLE_EQ(bvec.Z(), tvec.Z()) << bvec.PrintString();
79 EXPECT_DOUBLE_EQ(bvec.X(), tvec.X()) << bvec.PrintString();
80 EXPECT_DOUBLE_EQ(bvec.Y(), tvec.Y()) << bvec.PrintString();
81 EXPECT_DOUBLE_EQ(bvec.Z(), tvec.Z()) << bvec.PrintString();
82 //Since TVector3::Mag uses sqrt(x**2 + y**2 + z**2) instead of std::hypot
83 //there are some imprecisions in TVector3 so just compare this with float
84 //precision
85 EXPECT_FLOAT_EQ(bvec.Eta(), tvec.Eta()) << bvec.PrintString();
86 EXPECT_FLOAT_EQ(bvec.PseudoRapidity(), tvec.PseudoRapidity()) << bvec.PrintString();
87 }
88 }
89
90 TEST(B2Vector3, DoubleGettersVector)
91 {
92 gErrorIgnoreLevel = kError;
93 TVector3 tvec;
94 B2Vector3D bvec;
95 for (auto& rtp : test_vectors({0, 1, 1e20}, 16, 16)) {
96 tvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
97 bvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
98
99 TVector3 tvec2;
100 B2Vector3D bvec2;
101 for (auto& rtp2 : test_vectors({0, 1, 1e20}, 16, 16)) {
102 tvec2.SetMagThetaPhi(rtp2.r, rtp2.theta, rtp2.phi);
103 bvec2.SetMagThetaPhi(rtp2.r, rtp2.theta, rtp2.phi);
104 EXPECT_DOUBLE_EQ(tvec.Angle(tvec2), bvec.Angle(bvec2)) << bvec.PrintString() << ", " << bvec2.PrintString();
105 EXPECT_DOUBLE_EQ(tvec.DeltaPhi(tvec2), bvec.DeltaPhi(bvec2)) << bvec.PrintString() << ", " << bvec2.PrintString();
106 EXPECT_DOUBLE_EQ(tvec.Dot(tvec2), bvec.Dot(bvec2)) << bvec.PrintString() << ", " << bvec2.PrintString();
107 EXPECT_DOUBLE_EQ(tvec.Perp(tvec2), bvec.Perp(bvec2)) << bvec.PrintString() << ", " << bvec2.PrintString();
108 EXPECT_DOUBLE_EQ(tvec.Perp2(tvec2), bvec.Perp2(bvec2)) << bvec.PrintString() << ", " << bvec2.PrintString();
109 //ROOT uses sqrt(x**2+y**2+z**2) for the magnitude, we use
110 //std::hypot(std::hypot(x,y), z). This leads to imprecissions in ROOT
111 //so we have to adapt
112 EXPECT_NEAR(tvec.DeltaR(tvec2), bvec.DeltaR(bvec2), 1e-14) << bvec.PrintString() << ", " << bvec2.PrintString();
113 EXPECT_NEAR(tvec.DrEtaPhi(tvec2), bvec.DrEtaPhi(bvec2), 1e-14) << bvec.PrintString() << ", " << bvec2.PrintString();
114 }
115 }
116 }
117
118 TEST(B2Vector3, Rotation)
119 {
120 gErrorIgnoreLevel = kError;
121 TVector3 tvec;
122 B2Vector3D bvec;
123
124 // rotate \vec{e_x} by 90deg around \vec{e_z} to become \vec{e_y}
125 bvec.SetXYZ(1., 0., 0.);
126 tvec.SetXYZ(1., 0., 0.);
127 bvec.Rotate(M_PI_2, B2Vector3D(0., 0., 1.));
128 tvec.Rotate(M_PI_2, TVector3(0., 0., 1.));
129 EXPECT_NEAR(bvec.X(), 0., 1e-14) << bvec.PrintString();
130 EXPECT_NEAR(bvec.Y(), 1., 1e-14) << bvec.PrintString();
131 EXPECT_NEAR(bvec.Phi(), tvec.Phi(), 1e-14) << bvec.PrintString();
132 EXPECT_NEAR(bvec.Theta(), tvec.Theta(), 1e-14) << bvec.PrintString();
133 EXPECT_NEAR(bvec.X(), tvec.X(), 1e-14) << bvec.PrintString();
134 EXPECT_NEAR(bvec.Y(), tvec.Y(), 1e-14) << bvec.PrintString();
135 EXPECT_NEAR(bvec.Z(), tvec.Z(), 1e-14) << bvec.PrintString();
136
137 // rotate \vec{e_x} by 180deg around \vec{e_z} to become -\vec{e_x}
138 bvec.SetXYZ(1., 0., 0.);
139 tvec.SetXYZ(1., 0., 0.);
140 bvec.Rotate(M_PI, B2Vector3D(0., 0., 1.));
141 tvec.Rotate(M_PI, TVector3(0., 0., 1.));
142 EXPECT_NEAR(bvec.X(), -1., 1e-14) << bvec.PrintString();
143 EXPECT_NEAR(bvec.Y(), 0., 1e-14) << bvec.PrintString();
144 EXPECT_NEAR(bvec.Phi(), tvec.Phi(), 1e-14) << bvec.PrintString();
145 EXPECT_NEAR(bvec.Theta(), tvec.Theta(), 1e-14) << bvec.PrintString();
146 EXPECT_NEAR(bvec.X(), tvec.X(), 1e-14) << bvec.PrintString();
147 EXPECT_NEAR(bvec.Y(), tvec.Y(), 1e-14) << bvec.PrintString();
148 EXPECT_NEAR(bvec.Z(), tvec.Z(), 1e-14) << bvec.PrintString();
149
150 // just test a "random" case
151 bvec.SetXYZ(4., 3., 2.);
152 tvec.SetXYZ(4., 3., 2.);
153 bvec.Rotate(M_PI / 2.5, B2Vector3D(1., 2., 3.));
154 tvec.Rotate(M_PI / 2.5, TVector3(1., 2., 3.));
155 EXPECT_NEAR(bvec.Mag(), tvec.Mag(), 1e-14) << bvec.PrintString();
156 EXPECT_NEAR(bvec.Mag2(), tvec.Mag2(), 1e-14) << bvec.PrintString();
157 EXPECT_NEAR(bvec.Perp(), tvec.Perp(), 1e-14) << bvec.PrintString();
158 EXPECT_NEAR(bvec.Perp2(), tvec.Perp2(), 1e-14) << bvec.PrintString();
159 EXPECT_NEAR(bvec.Phi(), tvec.Phi(), 1e-14) << bvec.PrintString();
160 EXPECT_NEAR(bvec.Theta(), tvec.Theta(), 1e-14) << bvec.PrintString();
161 EXPECT_NEAR(bvec.X(), tvec.X(), 1e-14) << bvec.PrintString();
162 EXPECT_NEAR(bvec.Y(), tvec.Y(), 1e-14) << bvec.PrintString();
163 EXPECT_NEAR(bvec.Z(), tvec.Z(), 1e-14) << bvec.PrintString();
164
165 // rotate the test vectors all around the same axis by different angles
166 const B2Vector3D baxis(4., 3., -2.);
167 const TVector3 taxis(4., 3., -2.);
168
169 for (auto& rtp : test_vectors({0, 1, 1e20}, 64, 64)) {
170
171 double epsilon = (rtp.r < 1e10 ? 1e-10 : rtp.r * 1e-10);
172
173 // just test rotations around baxis / taxis by different angles in 0...2 M_PI
174 for (int i = 0; i < 28; i++) {
175 double angle = 2.*M_PI / 28 * i;
176 tvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
177 bvec.SetMagThetaPhi(rtp.r, rtp.theta, rtp.phi);
178 bvec.Rotate(angle, baxis);
179 tvec.Rotate(angle, taxis);
180
181 // Check for the single values to be +-PI and differing by 2 PI, which means both of them are
182 // basically equal, at least in terms of trigonometrical functions
183 if (fabs(fabs(bvec.Phi() - tvec.Phi()) - 2.*M_PI) < 1e-14) {
184 bvec.SetPhi(-bvec.Phi());
185 }
186
187 // Phi, Theta and CosTheta are not affected by large values of r, thus keep 1e-12
188 EXPECT_NEAR(bvec.CosTheta(), tvec.CosTheta(), 1e-12) << bvec.PrintString();
189 EXPECT_NEAR(bvec.Phi(), tvec.Phi(), 1e-12) << bvec.PrintString();
190 EXPECT_NEAR(bvec.Theta(), tvec.Theta(), 1e-12) << bvec.PrintString();
191// EXPECT_NEAR(bvec.Mag(), tvec.Mag(), (bvec.Mag() < 10 ? 1e-12 : bvec.Mag() * 1e-12)) << bvec.PrintString();
192// EXPECT_NEAR(bvec.Mag2(), tvec.Mag2(), (bvec.Mag2() < 10 ? 1e-12 : bvec.Mag2() * 1e-12)) << bvec.PrintString();
193// EXPECT_NEAR(bvec.Perp(), tvec.Perp(), (bvec.Perp() < 10 ? 1e-12 : bvec.Perp() * 1e-12)) << bvec.PrintString();
194// EXPECT_NEAR(bvec.Perp2(), tvec.Perp2(), (bvec.Perp2() < 10 ? 1e-12 : bvec.Perp2() * 1e-12)) << bvec.PrintString();
195 EXPECT_NEAR(bvec.X(), tvec.X(), epsilon) << bvec.PrintString();
196 EXPECT_NEAR(bvec.Y(), tvec.Y(), epsilon) << bvec.PrintString();
197 EXPECT_NEAR(bvec.Z(), tvec.Z(), epsilon) << bvec.PrintString();
198 }
199 }
200 }
201} // namespace
A fast and root compatible alternative to TVector3.
Definition: B2Vector3.h:42
DataType Phi() const
The azimuth angle.
Definition: B2Vector3.h:151
DataType Pz() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:441
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType Perp2() const
The transverse component squared (R^2 in cylindrical coordinate system).
Definition: B2Vector3.h:196
DataType Theta() const
The polar angle.
Definition: B2Vector3.h:153
DataType CosTheta() const
Cosine of the polar angle.
Definition: B2Vector3.h:155
void SetMagThetaPhi(DataType mag, DataType theta, DataType phi)
setter with mag, theta, phi
Definition: B2Vector3.h:259
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Eta() const
Returns the pseudo-rapidity.
Definition: B2Vector3.h:331
DataType DeltaPhi(const B2Vector3< DataType > &v) const
returns phi in the interval [-PI,PI)
Definition: B2Vector3.h:228
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
std::string PrintString(unsigned precision=4) const
create a string containing vector in cartesian and spherical coordinates
Definition: B2Vector3.h:481
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:159
DataType DeltaR(const B2Vector3< DataType > &v) const
return deltaR with respect to input-vector
Definition: B2Vector3.h:245
DataType Mag2() const
The magnitude squared (rho^2 in spherical coordinate system).
Definition: B2Vector3.h:157
DataType Dot(const B2Vector3< DataType > &p) const
Scalar product.
Definition: B2Vector3.h:290
DataType PseudoRapidity() const
Returns the pseudo-rapidity, i.e.
Definition: B2Vector3.h:319
DataType DrEtaPhi(const B2Vector3< DataType > &v) const
return DrEtaPhi with respect to input-vector
Definition: B2Vector3.h:253
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
Definition: B2Vector3.h:200
DataType Py() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:439
DataType Pt() const
The transverse component (R in cylindrical coordinate system).
Definition: B2Vector3.h:198
DataType Px() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:437
void SetXYZ(DataType x, DataType y, DataType z)
set all coordinates using data type
Definition: B2Vector3.h:464
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
void Rotate(DataType alpha, const B2Vector3< DataType > &v)
Rotation around an arbitrary axis v with angle alpha.
Definition: B2Vector3.h:399
void SetPhi(DataType phi)
Set phi keeping mag and theta constant.
Definition: B2Vector3.h:162
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.