TICL/src/SuperclusteringDNNInputs.cc

0001 /** Computation of input features for superclustering DNN. Used by plugins/TracksterLinkingBySuperClustering.cc and plugins/SuperclusteringSampleDumper.cc */
0002 // Author: Theo Cuisset - theo.cuisset@cern.ch
0003 // Date: 11/2023
0004 #include "RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h"
0005
0006 #include <algorithm>
0007 #include <numeric>
0008 #include <cmath>
0009
0010 #include <Math/Vector2D.h>
0011 #include <Math/Vector3D.h>
0012 #include <Math/Rotation3D.h>
0013 #include <Math/VectorUtil.h>
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "DataFormats/HGCalReco/interface/Trackster.h"
0017
0018 namespace ticl {
0019
0020   std::vector<float> SuperclusteringDNNInputV1::computeVector(Trackster const& ts_base, Trackster const& ts_toCluster) {
0021     /*  We use the barycenter for most of the variables below as that is what seems to have been used by Alessandro Tarabini,
0022         but using PCA might be better.
0023         (It would need retraining of the DNN)
0024     */
0025     return {{
0026         std::abs(ts_toCluster.barycenter().Eta()) - std::abs(ts_base.barycenter().Eta()),  //DeltaEtaBaryc
0027         ts_toCluster.barycenter().Phi() - ts_base.barycenter().phi(),                      //DeltaPhiBaryc
0028         ts_toCluster.raw_energy(),                                                         //multi_en
0029         ts_toCluster.barycenter().Eta(),                                                   //multi_eta
0030         (ts_toCluster.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())),         //multi_pt
0031         ts_base.barycenter().Eta(),                                                        //seedEta
0032         ts_base.barycenter().Phi(),                                                        //seedPhi
0033         ts_base.raw_energy(),                                                              //seedEn
0034         (ts_base.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())),              //seedPt
0035     }};
0036   }
0037
0038   // Helper functions for angles. Adapted from ROOT (3D vectors -> 2D vectors)
0039   template <class Vector1, class Vector2>
0040   float CosTheta2D(const Vector1& v1, const Vector2& v2) {
0041     float arg;
0042     float v1_r2 = v1.X() * v1.X() + v1.Y() * v1.Y();
0043     float v2_r2 = v2.X() * v2.X() + v2.Y() * v2.Y();
0044     float ptot2 = v1_r2 * v2_r2;
0045     if (ptot2 <= 0) {
0046       arg = 0.0;
0047     } else {
0048       float pdot = v1.X() * v2.X() + v1.Y() * v2.Y();
0049       using std::sqrt;
0050       arg = pdot / sqrt(ptot2);
0051       if (arg > 1.0)
0052         arg = 1.0;
0053       if (arg < -1.0)
0054         arg = -1.0;
0055     }
0056     return arg;
0057   }
0058   template <class Vector1, class Vector2>
0059   inline float Angle2D(const Vector1& v1, const Vector2& v2) {
0060     return static_cast<float>(std::acos(CosTheta2D(v1, v2)));
0061   }
0062
0063   std::vector<float> SuperclusteringDNNInputV2::computeVector(Trackster const& ts_base, Trackster const& ts_toCluster) {
0064     using ROOT::Math::XYVectorF;
0065     using ROOT::Math::XYZVectorF;
0066     using ROOT::Math::VectorUtil::Angle;
0067     XYZVectorF const& pca_seed_cmsFrame(ts_base.eigenvectors(0));
0068     XYZVectorF const& pca_cand_cmsFrame(ts_toCluster.eigenvectors(0));
0069     XYZVectorF xs(pca_seed_cmsFrame.Cross(XYZVectorF(0, 0, 1)).Unit());
0070     ROOT::Math::Rotation3D rot(xs, xs.Cross(pca_seed_cmsFrame).Unit(), pca_seed_cmsFrame);
0071
0072     XYZVectorF pca_cand_seedFrame = rot(pca_cand_cmsFrame);  // seed coordinates
0073
0074     float explVar_denominator = std::accumulate(
0075         std::begin(ts_toCluster.eigenvalues()), std::end(ts_toCluster.eigenvalues()), 0.f, std::plus<float>());
0076     float explVarRatio = 0.;
0077     if (explVar_denominator != 0.) {
0078       explVarRatio = ts_toCluster.eigenvalues()[0] / explVar_denominator;
0079     } else {
0080       edm::LogWarning("HGCalTICLSuperclustering")
0081           << "Sum of eigenvalues was zero for trackster. Could not compute explained variance ratio.";
0082     }
0083
0084     return {{
0085         std::abs(ts_toCluster.barycenter().Eta()) - std::abs(ts_base.barycenter().Eta()),  //DeltaEtaBaryc
0086         ts_toCluster.barycenter().Phi() - ts_base.barycenter().phi(),                      //DeltaPhiBaryc
0087         ts_toCluster.raw_energy(),                                                         //multi_en
0088         ts_toCluster.barycenter().Eta(),                                                   //multi_eta
0089         (ts_toCluster.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())),         //multi_pt
0090         ts_base.barycenter().Eta(),                                                        //seedEta
0091         ts_base.barycenter().Phi(),                                                        //seedPhi
0092         ts_base.raw_energy(),                                                              //seedEn
0093         (ts_base.raw_energy() * std::sin(ts_toCluster.barycenter().Theta())),              //seedPt
0094         static_cast<float>(Angle(pca_cand_cmsFrame, pca_seed_cmsFrame)),  // theta : angle between seed and candidate
0095         Angle2D(XYVectorF(pca_cand_seedFrame.x(), pca_cand_seedFrame.z()), XYVectorF(0, 1)),  // theta_xz_seedFrame
0096         Angle2D(XYVectorF(pca_cand_seedFrame.y(), pca_cand_seedFrame.z()), XYVectorF(0, 1)),  // theta_yz_seedFrame
0097         Angle2D(XYVectorF(pca_cand_cmsFrame.x(), pca_cand_cmsFrame.y()),
0098                 XYVectorF(pca_seed_cmsFrame.x(), pca_seed_cmsFrame.y())),  // theta_xy_cmsFrame
0099         Angle2D(XYVectorF(pca_cand_cmsFrame.y(), pca_cand_cmsFrame.z()),
0100                 XYVectorF(pca_seed_cmsFrame.y(), pca_seed_cmsFrame.z())),  // theta_yz_cmsFrame
0101         Angle2D(XYVectorF(pca_cand_cmsFrame.x(), pca_cand_cmsFrame.z()),
0102                 XYVectorF(pca_seed_cmsFrame.x(), pca_seed_cmsFrame.z())),  // theta_xz_cmsFrame
0103         ts_toCluster.eigenvalues()[0],                                     // explVar
0104         explVarRatio                                                       // explVarRatio
0105     }};
0106   }
0107
0108   std::unique_ptr<AbstractSuperclusteringDNNInput> makeSuperclusteringDNNInputFromString(std::string dnnInputVersion) {
0109     if (dnnInputVersion == "v1")
0110       return std::make_unique<SuperclusteringDNNInputV1>();
0111     else if (dnnInputVersion == "v2")
0112       return std::make_unique<SuperclusteringDNNInputV2>();
0113     assert(false);
0114   }
0115
0116 }  // namespace ticl