Belle II Software light-2406-ragdoll
FoxWolfram.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#include <analysis/ContinuumSuppression/FoxWolfram.h>
10
11using namespace Belle2;
12using namespace ROOT::Math;
13
15{
16 // clear the momenta, in case someone calls this function twice
17 for (double& i : m_moment)
18 i = 0.;
19
20 const auto begin = m_momenta.begin();
21 const auto end = m_momenta.end();
22
23 // loops over the vector pairs
24 for (auto iter1 = begin; iter1 != end; iter1++) {
25 const XYZVector pVec1 = (*iter1);
26 double pMag1 = pVec1.R();
27
28 // avoids to iterate twice over the same pairs
29 for (auto iter2 = iter1; iter2 != end; iter2++) {
30 const XYZVector pVec2 = (*iter2);
31 double magProd = pMag1 * pVec2.R(); // product of the vector's magnitudes
32 double cTheta = pVec1.Dot(pVec2) / magProd; // costheta_ij
33
34 // Since the FW moment definition requires to double count all the
35 // pairs of different particles, but the smart loop implemented here doesn't,
36 // multiply each entry by 2.
37 if (iter1 != iter2) magProd *= 2;
38
39 // Fills the moments' list.
40 // This part is quite ugly, but hard-coding the Legendre polynomials makes the code
41 // much faster than using the boost libraries, which are implementing the recursive formulas.
42 // This implementation should also be faster than a switch...case one.
43 double cTheta2 = cTheta * cTheta;
44 double cTheta3 = cTheta2 * cTheta;
45 double cTheta4 = cTheta2 * cTheta2;
46 m_moment[0] += magProd;
47 m_moment[1] += magProd * cTheta;
48 m_moment[2] += magProd * 0.5 * (3.*cTheta2 - 1);
49 m_moment[3] += magProd * 0.5 * (5.*cTheta3 - 3.*cTheta);
50 m_moment[4] += magProd * 0.125 * (35.*cTheta4 - 30.*cTheta2 + 3.);
51 }
52 }
53 return;
54}
55
56
57
59{
60 // clear the momenta, in case someone calls this function twice
61 for (double& i : m_moment)
62 i = 0.;
63
64 const auto begin = m_momenta.begin();
65 const auto end = m_momenta.end();
66
67 // loops over the vector pairs
68 for (auto iter1 = begin; iter1 != end; iter1++) {
69 const XYZVector pVec1 = (*iter1);
70 double pMag1 = pVec1.R();
71
72 // avoids to iterate twice over the same pairs
73 for (auto iter2 = iter1; iter2 != end; iter2++) {
74 const XYZVector pVec2 = (*iter2);
75 double magProd = pMag1 * pVec2.R(); // product of the vector's magnitudes
76 double cTheta = pVec1.Dot(pVec2) / magProd; // costheta_ij
77
78 // Since the FW moment definition requires to double count all the
79 // pairs of different particles, but the smart loop implemented here doesn't,
80 // multiply each entry by 2.
81 if (iter1 != iter2) magProd *= 2;
82
83 // Fills the moments' list.
84 // This part is quite ugly, but hard-coding the Legendre polynomials makes the code
85 // much faster than using the boost libraries, which are implementing the recursive formulas.
86 // This implementation should also be faster than a switch...case one.
87 double cTheta2 = cTheta * cTheta;
88 double cTheta3 = cTheta2 * cTheta;
89 double cTheta4 = cTheta2 * cTheta2;
90 double cTheta5 = cTheta4 * cTheta;
91 double cTheta6 = cTheta3 * cTheta3;
92 double cTheta7 = cTheta6 * cTheta;
93 double cTheta8 = cTheta4 * cTheta4;
94
95 m_moment[0] += magProd;
96 m_moment[1] += magProd * cTheta;
97 m_moment[2] += magProd * 0.5 * (3.*cTheta2 - 1);
98 m_moment[3] += magProd * 0.5 * (5.*cTheta3 - 3.*cTheta);
99 m_moment[4] += magProd * 0.125 * (35.*cTheta4 - 30.*cTheta2 + 3.);
100 m_moment[5] += magProd * 0.125 * (63.*cTheta5 - 70 * cTheta3 + 15.*cTheta);
101 m_moment[6] += magProd * 0.0625 * (231.*cTheta6 - 315 * cTheta4 + 105 * cTheta2 - 5.);
102 m_moment[7] += magProd * 0.0625 * (429.*cTheta7 - 693.*cTheta5 + 315.*cTheta3 - 35.*cTheta);
103 m_moment[8] += magProd * 0.0078125 * (6435.*cTheta8 - 12012.*cTheta6 + 6930.*cTheta4 - 1260.*cTheta2 + 35.);
104 }
105 }
106 return;
107}
108
double m_moment[9]
The moments.
Definition: FoxWolfram.h:93
std::vector< ROOT::Math::XYZVector > m_momenta
The particle's momenta.
Definition: FoxWolfram.h:94
void calculateAllMoments()
Method to perform the calculation of the moments up to order 8.
Definition: FoxWolfram.cc:58
void calculateBasicMoments()
Method to perform the calculation of the moments up to order 4, which are the most relevant ones.
Definition: FoxWolfram.cc:14
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24