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