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