17 LogDebug(
"TrackstersPCA_Eigen") <<
"------- Eigen -------" << std::endl;
19 for (
auto &trackster : tracksters) {
23 barycenter << 0., 0., 0.;
32 trackster.setRawEnergy(0.
f);
33 trackster.setRawEmEnergy(0.
f);
34 trackster.setRawPt(0.
f);
35 trackster.setRawEmPt(0.
f);
37 size_t N = trackster.vertices().size();
39 float weights2_sum = 0.f;
43 sigmasEigen << 0., 0., 0.;
46 std::vector<float> times;
47 std::vector<float> timeErrors;
48 std::set<uint32_t> usedLC;
50 for (
size_t i = 0;
i <
N; ++
i) {
51 auto fraction = 1.f / trackster.vertex_multiplicity(
i);
60 for (
size_t j = 0;
j < 3; ++
j)
64 if ((usedLC.insert(trackster.vertices(
i))).second) {
65 float timeE = layerClustersTime.
get(trackster.vertices(
i)).
second;
67 times.push_back(layerClustersTime.
get(trackster.vertices(
i)).
first);
68 timeErrors.push_back(1. /
pow(timeE, 2));
72 if (energyWeight && trackster.raw_energy())
73 barycenter /= trackster.raw_energy();
81 for (
size_t i = 0;
i <
N; ++
i) {
83 if (energyWeight && trackster.raw_energy())
85 (
layerClusters[trackster.vertices(
i)].energy() / trackster.vertex_multiplicity(
i)) / trackster.raw_energy();
87 for (
size_t x = 0;
x < 3; ++
x)
88 for (
size_t y = 0;
y <=
x; ++
y) {
90 covM(
y,
x) = covM(
x,
y);
93 covM *= 1. / (1. - weights2_sum);
96 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::RealVectorType eigenvalues_fromEigen;
97 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::EigenvectorsType eigenvectors_fromEigen;
98 Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(covM);
99 if (eigensolver.info() != Eigen::Success) {
100 eigenvalues_fromEigen = eigenvalues_fromEigen.Zero();
101 eigenvectors_fromEigen = eigenvectors_fromEigen.Zero();
103 eigenvalues_fromEigen = eigensolver.eigenvalues();
104 eigenvectors_fromEigen = eigensolver.eigenvectors();
108 for (
size_t i = 0;
i <
N; ++
i) {
110 sigmas +=
weight * (
point - barycenter).cwiseAbs2();
112 if (energyWeight && trackster.raw_energy())
114 (
layerClusters[trackster.vertices(
i)].energy() / trackster.vertex_multiplicity(
i)) / trackster.raw_energy();
115 sigmasEigen +=
weight * (point_transformed.cwiseAbs2());
117 sigmas /= (1. - weights2_sum);
118 sigmasEigen /= (1. - weights2_sum);
122 trackster.setTimeAndError(timeTrackster.first, timeTrackster.second);
123 trackster.fillPCAVariables(
126 LogDebug(
"TrackstersPCA") <<
"Use energy weighting: " << energyWeight << std::endl;
127 LogDebug(
"TrackstersPCA") <<
"\nTrackster characteristics: " << std::endl;
128 LogDebug(
"TrackstersPCA") <<
"Size: " <<
N << std::endl;
129 LogDebug(
"TrackstersPCA") <<
"Energy: " << trackster.raw_energy() << std::endl;
130 LogDebug(
"TrackstersPCA") <<
"raw_pt: " << trackster.raw_pt() << std::endl;
131 LogDebug(
"TrackstersPCA") <<
"Means: " << barycenter[0] <<
", " << barycenter[1] <<
", " << barycenter[2]
133 LogDebug(
"TrackstersPCA") <<
"Time: " << trackster.time() <<
" +/- " << trackster.timeError() << std::endl;
134 LogDebug(
"TrackstersPCA") <<
"EigenValues from Eigen/Tr(cov): " << eigenvalues_fromEigen[2] / covM.trace() <<
", "
135 << eigenvalues_fromEigen[1] / covM.trace() <<
", "
136 << eigenvalues_fromEigen[0] / covM.trace() << std::endl;
137 LogDebug(
"TrackstersPCA") <<
"EigenValues from Eigen: " << eigenvalues_fromEigen[2] <<
", "
138 << eigenvalues_fromEigen[1] <<
", " << eigenvalues_fromEigen[0] << std::endl;
139 LogDebug(
"TrackstersPCA") <<
"EigenVector 3 from Eigen: " << eigenvectors_fromEigen(0, 2) <<
", "
140 << eigenvectors_fromEigen(1, 2) <<
", " << eigenvectors_fromEigen(2, 2) << std::endl;
141 LogDebug(
"TrackstersPCA") <<
"EigenVector 2 from Eigen: " << eigenvectors_fromEigen(0, 1) <<
", "
142 << eigenvectors_fromEigen(1, 1) <<
", " << eigenvectors_fromEigen(2, 1) << std::endl;
143 LogDebug(
"TrackstersPCA") <<
"EigenVector 1 from Eigen: " << eigenvectors_fromEigen(0, 0) <<
", "
144 << eigenvectors_fromEigen(1, 0) <<
", " << eigenvectors_fromEigen(2, 0) << std::endl;
145 LogDebug(
"TrackstersPCA") <<
"Original sigmas: " << sigmas[0] <<
", " << sigmas[1] <<
", " << sigmas[2]
147 LogDebug(
"TrackstersPCA") <<
"SigmasEigen in PCA space: " << sigmasEigen[2] <<
", " << sigmasEigen[1] <<
", "
148 << sigmasEigen[0] << std::endl;
149 LogDebug(
"TrackstersPCA") <<
"covM: \n" << covM << std::endl;