Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:40

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "DataFormats/Common/interface/ValueMap.h"
0003 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0004 #include "TrackstersPCA.h"
0005 
0006 #include <iostream>
0007 #include <set>
0008 
0009 #include <Eigen/Core>
0010 #include <Eigen/Dense>
0011 
0012 void ticl::assignPCAtoTracksters(std::vector<Trackster> &tracksters,
0013                                  const std::vector<reco::CaloCluster> &layerClusters,
0014                                  const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0015                                  double z_limit_em,
0016                                  bool energyWeight) {
0017   LogDebug("TrackstersPCA_Eigen") << "------- Eigen -------" << std::endl;
0018 
0019   for (auto &trackster : tracksters) {
0020     Eigen::Vector3d point;
0021     point << 0., 0., 0.;
0022     Eigen::Vector3d barycenter;
0023     barycenter << 0., 0., 0.;
0024 
0025     auto fillPoint = [&](const reco::CaloCluster &c, const float weight = 1.f) {
0026       point[0] = weight * c.x();
0027       point[1] = weight * c.y();
0028       point[2] = weight * c.z();
0029     };
0030 
0031     // Initialize this trackster with default, dummy values
0032     trackster.setRawEnergy(0.f);
0033     trackster.setRawEmEnergy(0.f);
0034     trackster.setRawPt(0.f);
0035     trackster.setRawEmPt(0.f);
0036 
0037     size_t N = trackster.vertices().size();
0038     float weight = 1.f / N;
0039     float weights2_sum = 0.f;
0040     Eigen::Vector3d sigmas;
0041     sigmas << 0., 0., 0.;
0042     Eigen::Vector3d sigmasEigen;
0043     sigmasEigen << 0., 0., 0.;
0044     Eigen::Matrix3d covM = Eigen::Matrix3d::Zero();
0045 
0046     std::vector<float> times;
0047     std::vector<float> timeErrors;
0048     std::set<uint32_t> usedLC;
0049 
0050     for (size_t i = 0; i < N; ++i) {
0051       auto fraction = 1.f / trackster.vertex_multiplicity(i);
0052       trackster.addToRawEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0053       if (std::abs(layerClusters[trackster.vertices(i)].z()) <= z_limit_em)
0054         trackster.addToRawEmEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0055 
0056       // Compute the weighted barycenter.
0057       if (energyWeight)
0058         weight = layerClusters[trackster.vertices(i)].energy() * fraction;
0059       fillPoint(layerClusters[trackster.vertices(i)], weight);
0060       for (size_t j = 0; j < 3; ++j)
0061         barycenter[j] += point[j];
0062 
0063       // Add timing from layerClusters not already used
0064       if ((usedLC.insert(trackster.vertices(i))).second) {
0065         float timeE = layerClustersTime.get(trackster.vertices(i)).second;
0066         if (timeE > 0.f) {
0067           times.push_back(layerClustersTime.get(trackster.vertices(i)).first);
0068           timeErrors.push_back(1. / pow(timeE, 2));
0069         }
0070       }
0071     }
0072     if (energyWeight && trackster.raw_energy())
0073       barycenter /= trackster.raw_energy();
0074 
0075     hgcalsimclustertime::ComputeClusterTime timeEstimator;
0076     std::pair<float, float> timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors);
0077 
0078     // Compute the Covariance Matrix and the sum of the squared weights, used
0079     // to compute the correct normalization.
0080     // The barycenter has to be known.
0081     for (size_t i = 0; i < N; ++i) {
0082       fillPoint(layerClusters[trackster.vertices(i)]);
0083       if (energyWeight && trackster.raw_energy())
0084         weight =
0085             (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) / trackster.raw_energy();
0086       weights2_sum += weight * weight;
0087       for (size_t x = 0; x < 3; ++x)
0088         for (size_t y = 0; y <= x; ++y) {
0089           covM(x, y) += weight * (point[x] - barycenter[x]) * (point[y] - barycenter[y]);
0090           covM(y, x) = covM(x, y);
0091         }
0092     }
0093     covM *= 1. / (1. - weights2_sum);
0094 
0095     // Perform the actual decomposition
0096     Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::RealVectorType eigenvalues_fromEigen;
0097     Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::EigenvectorsType eigenvectors_fromEigen;
0098     Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(covM);
0099     if (eigensolver.info() != Eigen::Success) {
0100       eigenvalues_fromEigen = eigenvalues_fromEigen.Zero();
0101       eigenvectors_fromEigen = eigenvectors_fromEigen.Zero();
0102     } else {
0103       eigenvalues_fromEigen = eigensolver.eigenvalues();
0104       eigenvectors_fromEigen = eigensolver.eigenvectors();
0105     }
0106 
0107     // Compute the spread in the both spaces.
0108     for (size_t i = 0; i < N; ++i) {
0109       fillPoint(layerClusters[trackster.vertices(i)]);
0110       sigmas += weight * (point - barycenter).cwiseAbs2();
0111       Eigen::Vector3d point_transformed = eigenvectors_fromEigen * (point - barycenter);
0112       if (energyWeight && trackster.raw_energy())
0113         weight =
0114             (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) / trackster.raw_energy();
0115       sigmasEigen += weight * (point_transformed.cwiseAbs2());
0116     }
0117     sigmas /= (1. - weights2_sum);
0118     sigmasEigen /= (1. - weights2_sum);
0119 
0120     // Add trackster attributes
0121     trackster.setBarycenter(ticl::Trackster::Vector(barycenter));
0122     trackster.setTimeAndError(timeTrackster.first, timeTrackster.second);
0123     trackster.fillPCAVariables(
0124         eigenvalues_fromEigen, eigenvectors_fromEigen, sigmas, sigmasEigen, 3, ticl::Trackster::PCAOrdering::ascending);
0125 
0126     LogDebug("TrackstersPCA") << "Use energy weighting: " << energyWeight << std::endl;
0127     LogDebug("TrackstersPCA") << "\nTrackster characteristics: " << std::endl;
0128     LogDebug("TrackstersPCA") << "Size: " << N << std::endl;
0129     LogDebug("TrackstersPCA") << "Energy: " << trackster.raw_energy() << std::endl;
0130     LogDebug("TrackstersPCA") << "raw_pt: " << trackster.raw_pt() << std::endl;
0131     LogDebug("TrackstersPCA") << "Means:          " << barycenter[0] << ", " << barycenter[1] << ", " << barycenter[2]
0132                               << std::endl;
0133     LogDebug("TrackstersPCA") << "Time:          " << trackster.time() << " +/- " << trackster.timeError() << std::endl;
0134     LogDebug("TrackstersPCA") << "EigenValues from Eigen/Tr(cov): " << eigenvalues_fromEigen[2] / covM.trace() << ", "
0135                               << eigenvalues_fromEigen[1] / covM.trace() << ", "
0136                               << eigenvalues_fromEigen[0] / covM.trace() << std::endl;
0137     LogDebug("TrackstersPCA") << "EigenValues from Eigen:         " << eigenvalues_fromEigen[2] << ", "
0138                               << eigenvalues_fromEigen[1] << ", " << eigenvalues_fromEigen[0] << std::endl;
0139     LogDebug("TrackstersPCA") << "EigenVector 3 from Eigen: " << eigenvectors_fromEigen(0, 2) << ", "
0140                               << eigenvectors_fromEigen(1, 2) << ", " << eigenvectors_fromEigen(2, 2) << std::endl;
0141     LogDebug("TrackstersPCA") << "EigenVector 2 from Eigen: " << eigenvectors_fromEigen(0, 1) << ", "
0142                               << eigenvectors_fromEigen(1, 1) << ", " << eigenvectors_fromEigen(2, 1) << std::endl;
0143     LogDebug("TrackstersPCA") << "EigenVector 1 from Eigen: " << eigenvectors_fromEigen(0, 0) << ", "
0144                               << eigenvectors_fromEigen(1, 0) << ", " << eigenvectors_fromEigen(2, 0) << std::endl;
0145     LogDebug("TrackstersPCA") << "Original sigmas:          " << sigmas[0] << ", " << sigmas[1] << ", " << sigmas[2]
0146                               << std::endl;
0147     LogDebug("TrackstersPCA") << "SigmasEigen in PCA space: " << sigmasEigen[2] << ", " << sigmasEigen[1] << ", "
0148                               << sigmasEigen[0] << std::endl;
0149     LogDebug("TrackstersPCA") << "covM:     \n" << covM << std::endl;
0150   }
0151 }