Belle II Software  release-06-02-00
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 <TVector3.h>
19 #include <TRandom3.h>
20 #include <TMath.h>
21 #include <boost/math/special_functions/legendre.hpp>
22 
23 #include <gtest/gtest.h>
24 
25 using namespace std;
26 
27 namespace Belle2 {
33  class eventShapeCoreAlgorithmTest : public ::testing::Test {
34  protected:
35  };
36 
38  TEST_F(eventShapeCoreAlgorithmTest, Thrust)
39  {
40  std::vector<TVector3> 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 TVector3 thrustB = Thrust::calculateThrust(momenta);
49 
50  EXPECT_FLOAT_EQ(0.925838, thrustB.Mag());
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<TVector3> momenta;
62  std::vector<TVector3> sig_side_momenta;
63  std::vector<TVector3> 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 TVector3 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<TVector3> 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<TVector3> 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  TVector3 axis(0., 0., 1.);
196  // repeats the calculation
197  double moment[9] = {0.};
198  for (auto& p : partMom) {
199  double pMag = p.Mag();
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<TVector3> 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_all
std::vector< float > cleo_cone_with_all()
Returns calculated Cleo Cones constructed from all tracks.
Definition: CleoCones.h:39
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:44
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< TVector3 > &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:13
Belle2::HarmonicMoments
Class to calculate the Harmonic moments up to order 8 with respect to a given axis.
Definition: HarmonicMoments.h:28
Belle2::SphericityEigenvalues
Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-mom...
Definition: SphericityEigenvalues.h:23
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::TEST_F
TEST_F(eventShapeCoreAlgorithmTest, Sphericity)
Test the calculation of the Sphericity eigenvalues and eigenvectors.
Definition: eventShapeCoreAlgorithm.cc:216
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17