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