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 {
32
33 class eventShapeCoreAlgorithmTest : public ::testing::Test {
34 protected:
35 };
36
38 TEST_F(eventShapeCoreAlgorithmTest, Thrust)
39 {
40 std::vector<ROOT::Math::PxPyPzEVector> momenta;
41 // random generated numbers
42 momenta.emplace_back(0.5935352844151847, 0.28902324918117417, 0.9939000705771412, 2.);
43 momenta.emplace_back(0.7097025137911714, 0.5118418422879152, 0.44501044145648994, 2.);
44 momenta.emplace_back(0.6005771199332856, 0.12366608492454145, 0.7541373665256832, 2.);
45 momenta.emplace_back(0.8548902083897467, 0.6887268865943484, 0.34301136659215437, 2.);
46 momenta.emplace_back(0.26863521039211535, 0.011148638667487942, 0.96186334199901, 2.);
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::PxPyPzEVector> momenta;
62 std::vector<ROOT::Math::PxPyPzEVector> sig_side_momenta;
63 std::vector<ROOT::Math::PxPyPzEVector> roe_side_momenta;
64
65 // "Signal Side"
66 sig_side_momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432, 2.);
67 sig_side_momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105, 2.);
68 sig_side_momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818, 2.);
69 sig_side_momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237, 2.);
70 sig_side_momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986, 2.);
71
72 // "ROE Side"
73 roe_side_momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974, 2.);
74 roe_side_momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026, 2.);
75 roe_side_momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607, 2.);
76 roe_side_momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064, 2.);
77 roe_side_momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065, 2.);
78
79 // "All momenta"
80 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432, 2.);
81 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105, 2.);
82 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818, 2.);
83 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237, 2.);
84 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986, 2.);
85 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974, 2.);
86 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026, 2.);
87 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607, 2.);
88 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064, 2.);
89 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065, 2.);
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::PxPyPzEVector> momenta;
125
126 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432, 2.);
127 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105, 2.);
128 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818, 2.);
129 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237, 2.);
130 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986, 2.);
131 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974, 2.);
132 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026, 2.);
133 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607, 2.);
134 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064, 2.);
135 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065, 2.);
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., 1.);
146 momenta.emplace_back(-1., 0., 0., 1.);
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., 1.);
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<PxPyPzEVector> partMom;
183
184 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432, 2.);
185 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105, 2.);
186 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818, 2.);
187 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237, 2.);
188 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986, 2.);
189 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974, 2.);
190 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026, 2.);
191 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607, 2.);
192 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064, 2.);
193 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065, 2.);
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.Vect().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<PxPyPzEVector> partMom;
219
220 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432, 2.);
221 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105, 2.);
222 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818, 2.);
223 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237, 2.);
224 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986, 2.);
225 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974, 2.);
226 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026, 2.);
227 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607, 2.);
228 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064, 2.);
229 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065, 2.);
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:23
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:46
Belle2::CleoCones::cleo_cone_with_all
std::vector< float > cleo_cone_with_all()
Returns calculated Cleo Cones constructed from all tracks.
Definition CleoCones.h:41
Belle2::FoxWolfram
Class to calculate the Fox-Wolfram moments up to order 8.
Definition FoxWolfram.h:28
Belle2::FoxWolfram::setMomenta
void setMomenta(const std::vector< ROOT::Math::PxPyPzEVector > &momenta)
Sets the list of momenta used for the FW moment calculation, overwriting whatever list has been set b...
Definition FoxWolfram.h:73
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::HarmonicMoments
Class to calculate the Harmonic moments up to order 8 with respect to a given axis.
Definition HarmonicMoments.h:29
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:95
Belle2::SphericityEigenvalues
Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-mom...
Definition SphericityEigenvalues.h:25
Belle2::SphericityEigenvalues::calculateEigenvalues
void calculateEigenvalues()
Calculates eigenvalues and eigenvectors.
Definition SphericityEigenvalues.cc:17
Belle2::SphericityEigenvalues::getEigenvalue
double getEigenvalue(short i) const
Returns the i-th Eigenvalue.
Definition SphericityEigenvalues.h:63
Belle2::Thrust
Class to calculate the thrust axis.
Definition Thrust.h:23
Belle2::Thrust::calculateThrust
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::PxPyPzEVector > &momenta)
calculates the thrust axis
Definition Thrust.cc:71
Belle2::TEST_F
TEST_F(GlobalLabelTest, LargeNumberOfTimeDependentParameters)
Test large number of time-dep params for registration and retrieval.
Definition globalLabel.cc:69
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.