Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-13 03:24:05 

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 computeLocalTime,
0017                                  bool energyWeight) {
0018   LogDebug("TrackstersPCA_Eigen") << "------- Eigen -------" << std::endl;
0019 
0020   for (auto &trackster : tracksters) {
0021     Eigen::Vector3d point;
0022     point << 0., 0., 0.;
0023     Eigen::Vector3d barycenter;
0024     barycenter << 0., 0., 0.;
0025 
0026     auto fillPoint = [&](const reco::CaloCluster &c, const float weight = 1.f) {
0027       point[0] = weight * c.x();
0028       point[1] = weight * c.y();
0029       point[2] = weight * c.z();
0030     };
0031 
0032     // Initialize this trackster with default, dummy values
0033     trackster.setRawEnergy(0.f);
0034     trackster.setRawEmEnergy(0.f);
0035     trackster.setRawPt(0.f);
0036     trackster.setRawEmPt(0.f);
0037 
0038     size_t N = trackster.vertices().size();
0039     if (N == 0)
0040       continue;
0041     float weight = 1.f / N;
0042     float weights2_sum = 0.f;
0043 
0044     for (size_t i = 0; i < N; ++i) {
0045       auto fraction = 1.f / trackster.vertex_multiplicity(i);
0046       trackster.addToRawEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0047       if (std::abs(layerClusters[trackster.vertices(i)].z()) <= z_limit_em)
0048         trackster.addToRawEmEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0049 
0050       // Compute the weighted barycenter.
0051       if (energyWeight)
0052         weight = layerClusters[trackster.vertices(i)].energy() * fraction;
0053       fillPoint(layerClusters[trackster.vertices(i)], weight);
0054       for (size_t j = 0; j < 3; ++j)
0055         barycenter[j] += point[j];
0056     }
0057     float raw_energy = trackster.raw_energy();
0058     float inv_raw_energy = 1.f / raw_energy;
0059     if (energyWeight)
0060       barycenter *= inv_raw_energy;
0061     trackster.setBarycenter(ticl::Trackster::Vector(barycenter));
0062     std::pair<float, float> timeTrackster;
0063     if (computeLocalTime)
0064       timeTrackster = ticl::computeLocalTracksterTime(trackster, layerClusters, layerClustersTime, barycenter, N);
0065     else
0066       timeTrackster = ticl::computeTracksterTime(trackster, layerClustersTime, N);
0067 
0068     trackster.setTimeAndError(timeTrackster.first, timeTrackster.second);
0069     LogDebug("TrackstersPCA") << "Use energy weighting: " << energyWeight << std::endl;
0070     LogDebug("TrackstersPCA") << "\nTrackster characteristics: " << std::endl;
0071     LogDebug("TrackstersPCA") << "Size: " << N << std::endl;
0072     LogDebug("TrackstersPCA") << "Energy: " << trackster.raw_energy() << std::endl;
0073     LogDebug("TrackstersPCA") << "raw_pt: " << trackster.raw_pt() << std::endl;
0074     LogDebug("TrackstersPCA") << "Means:          " << barycenter[0] << ", " << barycenter[1] << ", " << barycenter[2]
0075                               << std::endl;
0076     LogDebug("TrackstersPCA") << "Time:          " << trackster.time() << " +/- " << trackster.timeError() << std::endl;
0077 
0078     if (N > 2) {
0079       Eigen::Vector3d sigmas;
0080       sigmas << 0., 0., 0.;
0081       Eigen::Vector3d sigmasEigen;
0082       sigmasEigen << 0., 0., 0.;
0083       Eigen::Matrix3d covM = Eigen::Matrix3d::Zero();
0084       // Compute the Covariance Matrix and the sum of the squared weights, used
0085       // to compute the correct normalization.
0086       // The barycenter has to be known.
0087       for (size_t i = 0; i < N; ++i) {
0088         fillPoint(layerClusters[trackster.vertices(i)]);
0089         if (energyWeight && trackster.raw_energy())
0090           weight = (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) * inv_raw_energy;
0091         weights2_sum += weight * weight;
0092         for (size_t x = 0; x < 3; ++x)
0093           for (size_t y = 0; y <= x; ++y) {
0094             covM(x, y) += weight * (point[x] - barycenter[x]) * (point[y] - barycenter[y]);
0095             covM(y, x) = covM(x, y);
0096           }
0097       }
0098       covM *= 1.f / (1.f - weights2_sum);
0099 
0100       // Perform the actual decomposition
0101       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::RealVectorType eigenvalues_fromEigen;
0102       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d>::EigenvectorsType eigenvectors_fromEigen;
0103       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(covM);
0104       if (eigensolver.info() != Eigen::Success) {
0105         eigenvalues_fromEigen = eigenvalues_fromEigen.Zero();
0106         eigenvectors_fromEigen = eigenvectors_fromEigen.Zero();
0107       } else {
0108         eigenvalues_fromEigen = eigensolver.eigenvalues();
0109         eigenvectors_fromEigen = eigensolver.eigenvectors();
0110       }
0111 
0112       // Compute the spread in the both spaces.
0113       for (size_t i = 0; i < N; ++i) {
0114         fillPoint(layerClusters[trackster.vertices(i)]);
0115         sigmas += weight * (point - barycenter).cwiseAbs2();
0116         Eigen::Vector3d point_transformed = eigenvectors_fromEigen * (point - barycenter);
0117         if (energyWeight && raw_energy)
0118           weight = (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) * inv_raw_energy;
0119         sigmasEigen += weight * (point_transformed.cwiseAbs2());
0120       }
0121       sigmas /= (1.f - weights2_sum);
0122       sigmasEigen /= (1.f - weights2_sum);
0123 
0124       trackster.fillPCAVariables(eigenvalues_fromEigen,
0125                                  eigenvectors_fromEigen,
0126                                  sigmas,
0127                                  sigmasEigen,
0128                                  3,
0129                                  ticl::Trackster::PCAOrdering::ascending);
0130 
0131       LogDebug("TrackstersPCA") << "EigenValues from Eigen/Tr(cov): " << eigenvalues_fromEigen[2] / covM.trace() << ", "
0132                                 << eigenvalues_fromEigen[1] / covM.trace() << ", "
0133                                 << eigenvalues_fromEigen[0] / covM.trace() << std::endl;
0134       LogDebug("TrackstersPCA") << "EigenValues from Eigen:         " << eigenvalues_fromEigen[2] << ", "
0135                                 << eigenvalues_fromEigen[1] << ", " << eigenvalues_fromEigen[0] << std::endl;
0136       LogDebug("TrackstersPCA") << "EigenVector 3 from Eigen: " << eigenvectors_fromEigen(0, 2) << ", "
0137                                 << eigenvectors_fromEigen(1, 2) << ", " << eigenvectors_fromEigen(2, 2) << std::endl;
0138       LogDebug("TrackstersPCA") << "EigenVector 2 from Eigen: " << eigenvectors_fromEigen(0, 1) << ", "
0139                                 << eigenvectors_fromEigen(1, 1) << ", " << eigenvectors_fromEigen(2, 1) << std::endl;
0140       LogDebug("TrackstersPCA") << "EigenVector 1 from Eigen: " << eigenvectors_fromEigen(0, 0) << ", "
0141                                 << eigenvectors_fromEigen(1, 0) << ", " << eigenvectors_fromEigen(2, 0) << std::endl;
0142       LogDebug("TrackstersPCA") << "Original sigmas:          " << sigmas[0] << ", " << sigmas[1] << ", " << sigmas[2]
0143                                 << std::endl;
0144       LogDebug("TrackstersPCA") << "SigmasEigen in PCA space: " << sigmasEigen[2] << ", " << sigmasEigen[1] << ", "
0145                                 << sigmasEigen[0] << std::endl;
0146       LogDebug("TrackstersPCA") << "covM:     \n" << covM << std::endl;
0147     }
0148   }
0149 }
0150 
0151 std::pair<float, float> ticl::computeLocalTracksterTime(const Trackster &trackster,
0152                                                         const std::vector<reco::CaloCluster> &layerClusters,
0153                                                         const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0154                                                         const Eigen::Vector3d &barycenter,
0155                                                         size_t N) {
0156   float tracksterTime = 0.;
0157   float tracksterTimeErr = 0.;
0158   std::set<uint32_t> usedLC;
0159 
0160   auto project_lc_to_pca = [](const std::vector<double> &point, const std::vector<double> &segment_end) {
0161     double dot_product = 0.0;
0162     double segment_dot = 0.0;
0163 
0164     for (int i = 0; i < 3; ++i) {
0165       dot_product += point[i] * segment_end[i];
0166       segment_dot += segment_end[i] * segment_end[i];
0167     }
0168 
0169     double projection = 0.0;
0170     if (segment_dot != 0.0) {
0171       projection = dot_product / segment_dot;
0172     }
0173 
0174     std::vector<double> closest_point(3);
0175     for (int i = 0; i < 3; ++i) {
0176       closest_point[i] = projection * segment_end[i];
0177     }
0178 
0179     double distance = 0.0;
0180     for (int i = 0; i < 3; ++i) {
0181       distance += std::pow(point[i] - closest_point[i], 2);
0182     }
0183 
0184     return std::sqrt(distance);
0185   };
0186 
0187   constexpr double c = 29.9792458;  // cm/ns
0188   for (size_t i = 0; i < N; ++i) {
0189     // Add timing from layerClusters not already used
0190     if ((usedLC.insert(trackster.vertices(i))).second) {
0191       float timeE = layerClustersTime.get(trackster.vertices(i)).second;
0192       if (timeE > 0.f) {
0193         float time = layerClustersTime.get(trackster.vertices(i)).first;
0194         timeE = 1.f / pow(timeE, 2);
0195         float x = layerClusters[trackster.vertices(i)].x();
0196         float y = layerClusters[trackster.vertices(i)].y();
0197         float z = layerClusters[trackster.vertices(i)].z();
0198 
0199         if (project_lc_to_pca({x, y, z}, {barycenter[0], barycenter[1], barycenter[2]}) < 3) {  // set MR to 3
0200           float deltaT = 1.f / c *
0201                          std::sqrt(((barycenter[2] / z - 1.f) * x) * ((barycenter[2] / z - 1.f) * x) +
0202                                    ((barycenter[2] / z - 1.f) * y) * ((barycenter[2] / z - 1.f) * y) +
0203                                    (barycenter[2] - z) * (barycenter[2] - z));
0204           time = std::abs(barycenter[2]) < std::abs(z) ? time - deltaT : time + deltaT;
0205 
0206           tracksterTime += time * timeE;
0207           tracksterTimeErr += timeE;
0208         }
0209       }
0210     }
0211   }
0212   if (tracksterTimeErr > 0.f)
0213     return {tracksterTime / tracksterTimeErr, 1.f / std::sqrt(tracksterTimeErr)};
0214   else
0215     return {-99.f, -1.f};
0216 }
0217 
0218 std::pair<float, float> ticl::computeTracksterTime(const Trackster &trackster,
0219                                                    const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0220                                                    size_t N) {
0221   std::vector<float> times;
0222   std::vector<float> timeErrors;
0223   std::set<uint32_t> usedLC;
0224 
0225   for (size_t i = 0; i < N; ++i) {
0226     // Add timing from layerClusters not already used
0227     if ((usedLC.insert(trackster.vertices(i))).second) {
0228       float timeE = layerClustersTime.get(trackster.vertices(i)).second;
0229       if (timeE > 0.f) {
0230         times.push_back(layerClustersTime.get(trackster.vertices(i)).first);
0231         timeErrors.push_back(1.f / pow(timeE, 2));
0232       }
0233     }
0234   }
0235 
0236   hgcalsimclustertime::ComputeClusterTime timeEstimator;
0237   return timeEstimator.fixSizeHighestDensity(times, timeErrors);
0238 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!