Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-03 05:27:16 

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 #include <vector>
0012 #include <functional>
0013 
0014 void ticl::assignPCAtoTracksters(std::vector<Trackster> &tracksters,
0015                                  const std::vector<reco::CaloCluster> &layerClusters,
0016                                  const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0017                                  double z_limit_em,
0018                                  const hgcal::RecHitTools &rhtools,
0019                                  bool computeLocalTime,
0020                                  bool energyWeight,
0021                                  bool clean,
0022                                  int minLayer,
0023                                  int maxLayer) {
0024   LogDebug("TrackstersPCA_Eigen") << "------- Eigen -------" << std::endl;
0025 
0026   for (auto &trackster : tracksters) {
0027     LogDebug("TrackstersPCA_Eigen") << "start testing teackster with size:" << trackster.vertices().size() << std::endl;
0028 
0029     Eigen::Vector3f point;
0030     point << 0., 0., 0.;
0031     Eigen::Vector3f barycenter;
0032     barycenter << 0., 0., 0.;
0033     Eigen::Vector3f filtered_barycenter;
0034     filtered_barycenter << 0., 0., 0.;
0035 
0036     auto fillPoint = [&](const reco::CaloCluster &c, const float weight = 1.f) {
0037       point[0] = weight * c.x();
0038       point[1] = weight * c.y();
0039       point[2] = weight * c.z();
0040     };
0041 
0042     // Initialize this trackster with default, dummy values
0043     trackster.setRawEnergy(0.f);
0044     trackster.setRawEmEnergy(0.f);
0045     trackster.setRawPt(0.f);
0046     trackster.setRawEmPt(0.f);
0047 
0048     size_t N = trackster.vertices().size();  ///< Number of layer clusters in trackster
0049     if (N == 0)
0050       continue;
0051     float weight = 1.f / N;
0052     float weights2_sum = 0.f;
0053 
0054     std::vector<float> layerClusterEnergies;
0055 
0056     for (size_t i = 0; i < N; ++i) {
0057       auto fraction = 1.f / trackster.vertex_multiplicity(i);
0058       trackster.addToRawEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0059       if (std::abs(layerClusters[trackster.vertices(i)].z()) <= z_limit_em)
0060         trackster.addToRawEmEnergy(layerClusters[trackster.vertices(i)].energy() * fraction);
0061 
0062       // Compute the weighted barycenter.
0063       if (energyWeight)
0064         weight = layerClusters[trackster.vertices(i)].energy() * fraction;
0065       fillPoint(layerClusters[trackster.vertices(i)], weight);
0066       for (size_t j = 0; j < 3; ++j)
0067         barycenter[j] += point[j];
0068 
0069       layerClusterEnergies.push_back(layerClusters[trackster.vertices(i)].energy());
0070     }
0071     float raw_energy = trackster.raw_energy();
0072     float inv_raw_energy = 1.f / raw_energy;
0073     if (energyWeight)
0074       barycenter *= inv_raw_energy;
0075     trackster.setBarycenter(ticl::Trackster::Vector(barycenter));
0076 
0077     LogDebug("TrackstersPCA_Eigen") << "cleaning is  :" << clean << std::endl;
0078 
0079     std::vector<unsigned> filtered_idx;  // indices of layer clusters to consider for cleaned PCA
0080     float filtered_energy = 0;
0081     float inv_filtered_energy = 0;
0082     if (clean) {
0083       // Filter layerclusters for the cleaned PCA
0084       auto maxE_vertex = std::distance(layerClusterEnergies.begin(),
0085                                        std::max_element(layerClusterEnergies.begin(), layerClusterEnergies.end()));
0086       auto maxE_layer = getLayerFromLC(layerClusters[trackster.vertices(maxE_vertex)], rhtools);
0087 
0088       auto vertices_by_layer = sortByLayer(trackster, layerClusters, rhtools);
0089 
0090       for (unsigned i = 1; i <= rhtools.lastLayer(); ++i) {
0091         auto vertices_in_layer = vertices_by_layer[i];
0092         if (vertices_in_layer.empty())
0093           continue;
0094 
0095         std::vector<float> energies_in_layer;
0096         for (auto vrt : vertices_in_layer)
0097           energies_in_layer.push_back(layerClusters[trackster.vertices(vrt)].energy());
0098 
0099         unsigned maxEid_inLayer = std::distance(energies_in_layer.begin(),
0100                                                 std::max_element(energies_in_layer.begin(), energies_in_layer.end()));
0101 
0102         // layer based filtering of what goes into the PCA
0103         if ((int)i >= (int)maxE_layer - minLayer && (int)i <= (int)maxE_layer + maxLayer) {
0104           auto filtered_vert = vertices_in_layer[maxEid_inLayer];
0105           filtered_idx.push_back(filtered_vert);
0106 
0107           const auto &maxE_LC = layerClusters[trackster.vertices(filtered_vert)];
0108           fillPoint(maxE_LC, maxE_LC.energy() * (1.f / trackster.vertex_multiplicity(filtered_vert)));
0109           for (size_t j = 0; j < 3; ++j)
0110             filtered_barycenter[j] += point[j];
0111           filtered_energy += maxE_LC.energy();
0112         }
0113       }
0114       inv_filtered_energy = 1. / filtered_energy;
0115       filtered_barycenter *= inv_filtered_energy;
0116     }
0117     LogDebug("TrackstersPCA_Eigen") << "min, max " << minLayer << "  " << maxLayer << std::endl;
0118 
0119     std::pair<float, float> timeTrackster;
0120     if (computeLocalTime)
0121       timeTrackster = ticl::computeLocalTracksterTime(trackster, layerClusters, layerClustersTime, barycenter, N);
0122     else
0123       timeTrackster = ticl::computeTracksterTime(trackster, layerClustersTime, N);
0124 
0125     trackster.setTimeAndError(timeTrackster.first, timeTrackster.second);
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 
0135     if (N > 2) {
0136       Eigen::Vector3f sigmas;
0137       sigmas << 0., 0., 0.;
0138       Eigen::Vector3f sigmasEigen;
0139       sigmasEigen << 0., 0., 0.;
0140       Eigen::Matrix3f covM = Eigen::Matrix3f::Zero();
0141       // Compute the Covariance Matrix and the sum of the squared weights, used
0142       // to compute the correct normalization.
0143       // The barycenter has to be known.
0144 
0145       auto calc_covM = [&](size_t i) {
0146         fillPoint(layerClusters[trackster.vertices(i)]);
0147         if (energyWeight && trackster.raw_energy()) {
0148           weight = (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) *
0149                    (clean ? inv_filtered_energy : inv_raw_energy);
0150           if (trackster.vertex_multiplicity(i) > 1)
0151             LogDebug("TrackstersPCA_Eigen")
0152                 << "trackster.vertex_multiplicity(i)   :" << trackster.vertex_multiplicity(i);
0153         }
0154         weights2_sum += weight * weight;
0155         for (size_t x = 0; x < 3; ++x) {
0156           for (size_t y = 0; y <= x; ++y) {
0157             covM(x, y) += weight * (point[x] - (clean ? filtered_barycenter[x] : barycenter[x])) *
0158                           (point[y] - (clean ? filtered_barycenter[y] : barycenter[y]));
0159             covM(y, x) = covM(x, y);
0160           }
0161         }
0162       };
0163       if (clean) {
0164         for (size_t i : filtered_idx) {
0165           calc_covM(i);
0166         }
0167       } else {
0168         for (size_t i = 0; i < N; ++i) {
0169           calc_covM(i);
0170         }
0171       }
0172 
0173       covM *= 1.f / (1.f - weights2_sum);
0174 
0175       // Perform the actual decomposition
0176       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f>::RealVectorType eigenvalues_fromEigen;
0177       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f>::EigenvectorsType eigenvectors_fromEigen;
0178       Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigensolver(covM);
0179       if (eigensolver.info() != Eigen::Success) {
0180         eigenvalues_fromEigen = eigenvalues_fromEigen.Zero();
0181         eigenvectors_fromEigen = eigenvectors_fromEigen.Zero();
0182       } else {
0183         eigenvalues_fromEigen = eigensolver.eigenvalues();
0184         eigenvectors_fromEigen = eigensolver.eigenvectors();
0185       }
0186 
0187       // Compute the spread in the both spaces.
0188       auto calc_spread = [&](size_t i) {
0189         fillPoint(layerClusters[trackster.vertices(i)]);
0190         sigmas += weight * (point - (clean ? filtered_barycenter : barycenter)).cwiseAbs2();
0191         Eigen::Vector3f point_transformed =
0192             eigenvectors_fromEigen * (point - (clean ? filtered_barycenter : barycenter));
0193         if (energyWeight && raw_energy)
0194           weight = (layerClusters[trackster.vertices(i)].energy() / trackster.vertex_multiplicity(i)) *
0195                    (clean ? inv_filtered_energy : inv_raw_energy);
0196         sigmasEigen += weight * (point_transformed.cwiseAbs2());
0197       };
0198 
0199       if (clean) {
0200         for (size_t i : filtered_idx) {
0201           calc_spread(i);
0202         }
0203       } else {
0204         for (size_t i = 0; i < N; ++i) {
0205           calc_spread(i);
0206         }
0207       }
0208 
0209       sigmas /= (1.f - weights2_sum);
0210       sigmasEigen /= (1.f - weights2_sum);
0211 
0212       trackster.fillPCAVariables(eigenvalues_fromEigen,
0213                                  eigenvectors_fromEigen,
0214                                  sigmas,
0215                                  sigmasEigen,
0216                                  3,
0217                                  ticl::Trackster::PCAOrdering::ascending);
0218 
0219       LogDebug("TrackstersPCA") << "EigenValues from Eigen/Tr(cov): " << eigenvalues_fromEigen[2] / covM.trace() << ", "
0220                                 << eigenvalues_fromEigen[1] / covM.trace() << ", "
0221                                 << eigenvalues_fromEigen[0] / covM.trace() << std::endl;
0222       LogDebug("TrackstersPCA") << "EigenValues from Eigen:         " << eigenvalues_fromEigen[2] << ", "
0223                                 << eigenvalues_fromEigen[1] << ", " << eigenvalues_fromEigen[0] << std::endl;
0224       LogDebug("TrackstersPCA") << "EigenVector 3 from Eigen: " << eigenvectors_fromEigen(0, 2) << ", "
0225                                 << eigenvectors_fromEigen(1, 2) << ", " << eigenvectors_fromEigen(2, 2) << std::endl;
0226       LogDebug("TrackstersPCA") << "EigenVector 2 from Eigen: " << eigenvectors_fromEigen(0, 1) << ", "
0227                                 << eigenvectors_fromEigen(1, 1) << ", " << eigenvectors_fromEigen(2, 1) << std::endl;
0228       LogDebug("TrackstersPCA") << "EigenVector 1 from Eigen: " << eigenvectors_fromEigen(0, 0) << ", "
0229                                 << eigenvectors_fromEigen(1, 0) << ", " << eigenvectors_fromEigen(2, 0) << std::endl;
0230       LogDebug("TrackstersPCA") << "Original sigmas:          " << sigmas[0] << ", " << sigmas[1] << ", " << sigmas[2]
0231                                 << std::endl;
0232       LogDebug("TrackstersPCA") << "SigmasEigen in PCA space: " << sigmasEigen[2] << ", " << sigmasEigen[1] << ", "
0233                                 << sigmasEigen[0] << std::endl;
0234       LogDebug("TrackstersPCA") << "covM:     \n" << covM << std::endl;
0235     }
0236   }
0237 }
0238 
0239 std::pair<float, float> ticl::computeLocalTracksterTime(const Trackster &trackster,
0240                                                         const std::vector<reco::CaloCluster> &layerClusters,
0241                                                         const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0242                                                         const Eigen::Vector3f &barycenter,
0243                                                         size_t N) {
0244   float tracksterTime = 0.;
0245   float tracksterTimeErr = 0.;
0246 
0247   auto project_lc_to_pca = [](const std::array<float, 3> &point, const std::array<float, 3> &segment_end) {
0248     float dot_product = 0.0;
0249     float segment_dot = 0.0;
0250 
0251     for (int i = 0; i < 3; ++i) {
0252       dot_product += point[i] * segment_end[i];
0253       segment_dot += segment_end[i] * segment_end[i];
0254     }
0255 
0256     float projection = 0.0;
0257     if (segment_dot != 0.0) {
0258       projection = dot_product / segment_dot;
0259     }
0260 
0261     std::array<float, 3> closest_point;
0262     for (int i = 0; i < 3; ++i) {
0263       closest_point[i] = projection * segment_end[i];
0264     }
0265 
0266     float distanceSquared = 0.f;
0267     for (int i = 0; i < 3; ++i) {
0268       distanceSquared += std::pow(point[i] - closest_point[i], 2);
0269     }
0270     return distanceSquared;
0271   };
0272 
0273   constexpr float c = 29.9792458;  // cm/ns
0274   for (size_t i = 0; i < N; ++i) {
0275     // Add timing from layerClusters not already used
0276     float timeE = layerClustersTime.get(trackster.vertices(i)).second;
0277     if (timeE > 0.f) {
0278       float time = layerClustersTime.get(trackster.vertices(i)).first;
0279       timeE = 1.f / pow(timeE, 2);
0280       float x = layerClusters[trackster.vertices(i)].x();
0281       float y = layerClusters[trackster.vertices(i)].y();
0282       float z = layerClusters[trackster.vertices(i)].z();
0283 
0284       if (project_lc_to_pca({{x, y, z}}, {{barycenter[0], barycenter[1], barycenter[2]}}) < 9.f) {  // set MR to 3
0285         float invz = 1.f / z;
0286         float deltaT = 1.f / c *
0287                        std::sqrt(((barycenter[2] * invz - 1.f) * x) * ((barycenter[2] * invz - 1.f) * x) +
0288                                  ((barycenter[2] * invz - 1.f) * y) * ((barycenter[2] * invz - 1.f) * y) +
0289                                  (barycenter[2] - z) * (barycenter[2] - z));
0290         time = std::abs(barycenter[2]) < std::abs(z) ? time - deltaT : time + deltaT;
0291 
0292         tracksterTime += time * timeE;
0293         tracksterTimeErr += timeE;
0294       }
0295     }
0296   }
0297   if (tracksterTimeErr > 0.f)
0298     return {tracksterTime / tracksterTimeErr, 1.f / std::sqrt(tracksterTimeErr)};
0299   else
0300     return {-99.f, -1.f};
0301 }
0302 
0303 std::pair<float, float> ticl::computeTracksterTime(const Trackster &trackster,
0304                                                    const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
0305                                                    size_t N) {
0306   std::vector<float> times;
0307   std::vector<float> timeErrors;
0308 
0309   for (size_t i = 0; i < N; ++i) {
0310     float timeE = layerClustersTime.get(trackster.vertices(i)).second;
0311     if (timeE > 0.f) {
0312       times.push_back(layerClustersTime.get(trackster.vertices(i)).first);
0313       timeErrors.push_back(1.f / pow(timeE, 2));
0314     }
0315   }
0316 
0317   hgcalsimclustertime::ComputeClusterTime timeEstimator;
0318   return timeEstimator.fixSizeHighestDensity(times, timeErrors);
0319 }
[ 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!