Belle II Software light-2406-ragdoll
eventShapeCoreAlgorithm.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 * This test files used to be the continuumSuppression test. *
10 ****************************************************************/
11
12#include <analysis/ContinuumSuppression/Thrust.h>
13#include <analysis/ContinuumSuppression/CleoCones.h>
14#include <analysis/ContinuumSuppression/FoxWolfram.h>
15#include <analysis/ContinuumSuppression/HarmonicMoments.h>
16#include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
17
18#include <TRandom3.h>
19#include <TMath.h>
20#include <boost/math/special_functions/legendre.hpp>
21
22#include <gtest/gtest.h>
23
24using namespace std;
25using namespace ROOT::Math;
26
27namespace Belle2 {
33 class eventShapeCoreAlgorithmTest : public ::testing::Test {
34 protected:
35 };
36
39 {
40 std::vector<ROOT::Math::XYZVector> momenta;
41 // random generated numbers
42 momenta.emplace_back(0.5935352844151847, 0.28902324918117417, 0.9939000705771412);
43 momenta.emplace_back(0.7097025137911714, 0.5118418422879152, 0.44501044145648994);
44 momenta.emplace_back(0.6005771199332856, 0.12366608492454145, 0.7541373665256832);
45 momenta.emplace_back(0.8548902083897467, 0.6887268865943484, 0.34301136659215437);
46 momenta.emplace_back(0.26863521039211535, 0.011148638667487942, 0.96186334199901);
47
48 const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(momenta);
49
50 EXPECT_FLOAT_EQ(0.925838, thrustB.R());
51 EXPECT_FLOAT_EQ(0.571661, thrustB.X());
52 EXPECT_FLOAT_EQ(0.306741, thrustB.Y());
53 EXPECT_FLOAT_EQ(0.660522, thrustB.Z());
54 }
55
58 {
59 const bool use_all = true;
60 const bool use_roe = true;
61 std::vector<ROOT::Math::XYZVector> momenta;
62 std::vector<ROOT::Math::XYZVector> sig_side_momenta;
63 std::vector<ROOT::Math::XYZVector> roe_side_momenta;
64
65 // "Signal Side"
66 sig_side_momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
67 sig_side_momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
68 sig_side_momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
69 sig_side_momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
70 sig_side_momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
71
72 // "ROE Side"
73 roe_side_momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
74 roe_side_momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
75 roe_side_momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
76 roe_side_momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
77 roe_side_momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
78
79 // "All momenta"
80 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
81 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
82 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
83 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
84 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
85 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
86 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
87 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
88 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
89 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
90
91
92 // Calculate thrust from "Signal Side"
93 const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(sig_side_momenta);
94
95 CleoCones myCleoCones(momenta, roe_side_momenta, thrustB, use_all, use_roe);
96
97 const auto& cleo_cone_with_all = myCleoCones.cleo_cone_with_all();
98 const auto& cleo_cone_with_roe = myCleoCones.cleo_cone_with_roe();
99
100 EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_all[0]);
101 EXPECT_FLOAT_EQ(4.7405558, cleo_cone_with_all[1]);
102 EXPECT_FLOAT_EQ(1.7517139, cleo_cone_with_all[2]);
103 EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_all[3]);
104 EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_all[4]);
105 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[5]);
106 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[6]);
107 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[7]);
108 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[8]);
109
110 EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_roe[0]);
111 EXPECT_FLOAT_EQ(1.9106253, cleo_cone_with_roe[1]);
112 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[2]);
113 EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_roe[3]);
114 EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_roe[4]);
115 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[5]);
116 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[6]);
117 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[7]);
118 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[8]);
119 }
120
123 {
124 std::vector<ROOT::Math::XYZVector> momenta;
125
126 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
127 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
128 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
129 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
130 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
131 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
132 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
133 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
134 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
135 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
136
137 FoxWolfram FW(momenta);
139 EXPECT_FLOAT_EQ(0.63011014, FW.getR(2));
140
141 // Tests using the analytic values provided by Fox & Wolfram
142 // for simple distributions
143 // first: Back-to-back versors
144 momenta.clear();
145 momenta.emplace_back(1., 0., 0.);
146 momenta.emplace_back(-1., 0., 0.);
147
148 FW.setMomenta(momenta);
150
151 EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1) < 0.001);
152 EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
153 EXPECT_TRUE(TMath::Abs(FW.getR(2) - 1) < 0.001);
154 EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
155 EXPECT_TRUE(TMath::Abs(FW.getR(4) - 1) < 0.001);
156
157 // second: equatorial line
158 momenta.clear();
159 TRandom3 rnd;
160 for (int i = 0; i < 10000; i++) {
161 double phi = rnd.Uniform(0., 2 * TMath::Pi());
162 momenta.emplace_back(TMath::Cos(phi), TMath::Sin(phi), 0.);
163 }
164 FW.setMomenta(momenta);
166
167 EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1.) < 0.001);
168 EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
169 EXPECT_TRUE(TMath::Abs(FW.getR(2) - 0.25) < 0.001);
170 EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
171 EXPECT_TRUE(TMath::Abs(FW.getR(4) - 0.14) < 0.001);
172
173 }
174
175
176
179 {
180
181 float dummySqrtS = 10.;
182 std::vector<XYZVector> partMom;
183
184 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
185 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
186 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
187 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
188 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
189 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
190 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
191 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
192 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
193 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
194
195 XYZVector axis(0., 0., 1.);
196 // repeats the calculation
197 double moment[9] = {0.};
198 for (auto& p : partMom) {
199 double pMag = p.R();
200 double cTheta = p.Dot(axis) / pMag;
201 for (short i = 0; i < 9; i++)
202 moment[i] += pMag * boost::math::legendre_p(i, cTheta) / dummySqrtS;
203 }
204
205 HarmonicMoments MM(partMom, axis);
207
208 for (short i = 0; i < 9; i++)
209 EXPECT_FLOAT_EQ(moment[i], MM.getMoment(i, dummySqrtS));
210 }
211
212
213
214
216 TEST_F(eventShapeCoreAlgorithmTest, Sphericity)
217 {
218 std::vector<XYZVector> partMom;
219
220 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
221 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
222 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
223 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
224 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
225 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
226 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
227 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
228 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
229 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
230
231
232 SphericityEigenvalues Sph(partMom);
233 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(0));
234 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(1));
235 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(2));
236
238
239 EXPECT_FLOAT_EQ(0.87272972, Sph.getEigenvalue(0));
240 EXPECT_FLOAT_EQ(0.095489956, Sph.getEigenvalue(1));
241 EXPECT_FLOAT_EQ(0.031780295, Sph.getEigenvalue(2));
242 EXPECT_FLOAT_EQ(1., Sph.getEigenvalue(0) + Sph.getEigenvalue(1) + Sph.getEigenvalue(2));
243 }
244
245
247} // namespace
Class to calculate the Cleo clone variables.
Definition: CleoCones.h:22
std::vector< float > cleo_cone_with_roe()
Returns calculated Cleo Cones constructed from only ROE tracks.
Definition: CleoCones.h:45
std::vector< float > cleo_cone_with_all()
Returns calculated Cleo Cones constructed from all tracks.
Definition: CleoCones.h:40
Class to calculate the Fox-Wolfram moments up to order 8.
Definition: FoxWolfram.h:28
double getR(int i) const
Returns the i-th moment normalized to the 0th-order moment.
Definition: FoxWolfram.h:89
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition: FoxWolfram.cc:14
void setMomenta(const std::vector< ROOT::Math::XYZVector > &momenta)
Sets the list of momenta used for the FW moment calculation, overwriting whatever list has been set b...
Definition: FoxWolfram.h:73
Class to calculate the Harmonic moments up to order 8 with respect to a given axis.
void calculateAllMoments()
Calculates the moments up to order 8.
double getMoment(short i, double sqrts) const
Returns the moment of order i.
Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-mom...
void calculateEigenvalues()
Calculates eigenvalues and eigenvectors.
double getEigenvalue(short i) const
Returns the i-th Eigenvalue.
Class to calculate the thrust axis.
Definition: Thrust.h:22
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::XYZVector > &momenta)
calculates the thrust axis
Definition: Thrust.cc:71
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.