Belle II Software  release-05-02-19
SphericityEigenvalues.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Umberto Tamponi (tamponi@to.infn.it) *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 
12 #include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
13 #include <framework/logging/Logger.h>
14 #include <Eigen/Core>
15 #include <Eigen/Eigenvalues>
16 
17 #include <string>
18 
19 using namespace Belle2;
20 
22 {
23  Eigen::Matrix3f sphericityTensor;
24 
25  // elements of the matrix, in rows
26  // n = r*3+c, using all 0-based indexes:
27  // 00 = 0
28  // 01 = 1
29  // 10 = 3
30  // 12 = 5
31  // 22 = 8
32  // etc...
33  // diagonal = 0, 4, 8
34  double elements[9] = {0.};
35 
36  // normalization
37  double norm = 0;
38 
39 
40  if (m_momenta.size() < 2) {
41  B2WARNING("The particle list has less than 2 elements. The sphericity matrix will not be calculated");
42  return;
43  }
44 
45  for (const auto& p : m_momenta) {
46  elements[0] += p.X() * p.X(); // diag
47  elements[1] += p.X() * p.Y();
48  elements[2] += p.X() * p.Z();
49 
50  elements[3] += p.Y() * p.X();
51  elements[4] += p.Y() * p.Y(); // diag
52  elements[5] += p.Y() * p.Z();
53 
54  elements[6] += p.Z() * p.X();
55  elements[7] += p.Z() * p.Y();
56  elements[8] += p.Z() * p.Z(); // diag
57  norm += p.Mag2();
58  }
59 
60  for (short i = 0; i < 3; i++) {
61  for (short j = 0; j < 3; j++) {
62  sphericityTensor(i, j) = elements[i * 3 + j] / norm;
63  }
64  }
65 
66  auto eigenVals = sphericityTensor.eigenvalues();
67  Eigen::ComplexEigenSolver<Eigen::MatrixXcf> ces(sphericityTensor);
68 
69  // unfortunately Eigen does not provide the eigenvalues in
70  // any specific order, so we have to sort them and keep also the correct eigenvector-eigenvalue
71  // associations...
72 
73  short order[3] = {0, 1, 2};
74 
75  std::vector<float> tmpLambda;
76  for (short i = 0; i < 3; i++) {
77  tmpLambda.push_back(eigenVals[i].real());
78  }
79 
80  // position of the largest Eigenvalue
81  order[0] = std::distance(tmpLambda.begin(), std::max_element(tmpLambda.begin(), tmpLambda.end()));
82  // position of the smallest Eigenvalue
83  order[2] = std::distance(tmpLambda.begin(), std::min_element(tmpLambda.begin(), tmpLambda.end()));
84  // position of the middle eigenvalue
85  order[1] = (short)(3.1 - (order[0] + order[2]));
86 
87  for (short i = 0; i < 3; i++) {
88  short n = order[i];
89  m_lambda[i] = eigenVals[n].real();
90  auto eigenVector = ces.eigenvectors().col(n);
91  m_eVector[i].SetX(eigenVector[0].real());
92  m_eVector[i].SetY(eigenVector[1].real());
93  m_eVector[i].SetZ(eigenVector[2].real());
94  }
95 
96  return;
97 }
98 
Belle2::SphericityEigenvalues::m_eVector
TVector3 m_eVector[3]
The eigenvectors.
Definition: SphericityEigenvalues.h:89
Belle2::SphericityEigenvalues::calculateEigenvalues
void calculateEigenvalues()
Calculates eigenvalues and eigenvectors.
Definition: SphericityEigenvalues.cc:21
Belle2::SphericityEigenvalues::m_momenta
std::vector< TVector3 > m_momenta
The particles' momenta.
Definition: SphericityEigenvalues.h:90
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SphericityEigenvalues::m_lambda
double m_lambda[3]
The eigenvalues.
Definition: SphericityEigenvalues.h:88