Belle II Software development
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
38 TEST_F(eventShapeCoreAlgorithmTest, Thrust)
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
57 TEST_F(eventShapeCoreAlgorithmTest, CleoCones)
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
122 TEST_F(eventShapeCoreAlgorithmTest, FoxWolfram)
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);
138 FW.calculateBasicMoments();
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);
149 FW.calculateBasicMoments();
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);
165 FW.calculateBasicMoments();
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
178 TEST_F(eventShapeCoreAlgorithmTest, HarmonicMoments)
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);
206 MM.calculateAllMoments();
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
237 Sph.calculateEigenvalues();
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
Belle2::CleoCones
Class to calculate the Cleo clone variables.
Definition: CleoCones.h:22
Belle2::CleoCones::cleo_cone_with_roe
std::vector< float > cleo_cone_with_roe()
Returns calculated Cleo Cones constructed from only ROE tracks.
Definition: CleoCones.h:45
Belle2::CleoCones::cleo_cone_with_all
std::vector< float > cleo_cone_with_all()
Returns calculated Cleo Cones constructed from all tracks.
Definition: CleoCones.h:40
Belle2::FoxWolfram
Class to calculate the Fox-Wolfram moments up to order 8.
Definition: FoxWolfram.h:28
Belle2::FoxWolfram::getR
double getR(int i) const
Returns the i-th moment normalized to the 0th-order moment.
Definition: FoxWolfram.h:89
Belle2::FoxWolfram::calculateBasicMoments
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition: FoxWolfram.cc:14
Belle2::FoxWolfram::setMomenta
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
Belle2::HarmonicMoments
Class to calculate the Harmonic moments up to order 8 with respect to a given axis.
Definition: HarmonicMoments.h:28
Belle2::HarmonicMoments::calculateAllMoments
void calculateAllMoments()
Calculates the moments up to order 8.
Definition: HarmonicMoments.cc:42
Belle2::HarmonicMoments::getMoment
double getMoment(short i, double sqrts) const
Returns the moment of order i.
Definition: HarmonicMoments.h:94
Belle2::SphericityEigenvalues
Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-mom...
Definition: SphericityEigenvalues.h:24
Belle2::SphericityEigenvalues::calculateEigenvalues
void calculateEigenvalues()
Calculates eigenvalues and eigenvectors.
Definition: SphericityEigenvalues.cc:19
Belle2::SphericityEigenvalues::getEigenvalue
double getEigenvalue(short i) const
Returns the i-th Eigenvalue.
Definition: SphericityEigenvalues.h:62
Belle2::Thrust
Class to calculate the thrust axis.
Definition: Thrust.h:22
Belle2::Thrust::calculateThrust
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::XYZVector > &momenta)
calculates the thrust axis
Definition: Thrust.cc:71
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.