12 #include <analysis/ContinuumSuppression/SphericityEigenvalues.h>
13 #include <framework/logging/Logger.h>
15 #include <Eigen/Eigenvalues>
23 Eigen::Matrix3f sphericityTensor;
34 double elements[9] = {0.};
41 B2WARNING(
"The particle list has less than 2 elements. The sphericity matrix will not be calculated");
46 elements[0] += p.X() * p.X();
47 elements[1] += p.X() * p.Y();
48 elements[2] += p.X() * p.Z();
50 elements[3] += p.Y() * p.X();
51 elements[4] += p.Y() * p.Y();
52 elements[5] += p.Y() * p.Z();
54 elements[6] += p.Z() * p.X();
55 elements[7] += p.Z() * p.Y();
56 elements[8] += p.Z() * p.Z();
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;
66 auto eigenVals = sphericityTensor.eigenvalues();
67 Eigen::ComplexEigenSolver<Eigen::MatrixXcf> ces(sphericityTensor);
73 short order[3] = {0, 1, 2};
75 std::vector<float> tmpLambda;
76 for (
short i = 0; i < 3; i++) {
77 tmpLambda.push_back(eigenVals[i].real());
81 order[0] = std::distance(tmpLambda.begin(), std::max_element(tmpLambda.begin(), tmpLambda.end()));
83 order[2] = std::distance(tmpLambda.begin(), std::min_element(tmpLambda.begin(), tmpLambda.end()));
85 order[1] = (short)(3.1 - (order[0] + order[2]));
87 for (
short i = 0; i < 3; i++) {
90 auto eigenVector = ces.eigenvectors().col(n);