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